PaperlayerView exact revision on arXiv

Source-linked AI summary

Kernel methods through the roof: handling billions of points efficiently

Giacomo Meanti, Luigi Carratino, Lorenzo Rosasco, Alessandro Rudi

arXiv:2006.10350v2cs.LGstat.ML

TL;DR

Kernel methods are difficult to scale because standard implementations require costly kernel-matrix operations. The paper combines Nyström approximation and preconditioned gradient solvers with GPU and multi-GPU techniques; it reports state-of-the-art accuracy in minutes on datasets reaching billions of points.

  • Problem

    Kernel methods have strong theoretical foundations, but implementations that handle kernel matrices quadratic in sample size make large datasets difficult to address.

    Source: §1 · p. 2

  • Method

    The authors combine Nyström approximation and preconditioned conjugate-gradient solvers with GPU acceleration, multi-GPU parallelization, and out-of-core linear algebra.

    Source: §1 · p. 3, §3 · p. 6, §3.2 · p. 8

  • Results

    The approach achieves state-of-the-art accuracy in minutes on benchmarks with millions to billions of points, and the authors report favorable efficiency and accuracy against other solutions.

    Source: §1 · p. 3

  • Takeaways & Limitations

    The authors provide an easy-to-use library and conclude that kernel methods can be effectively deployed on large-scale problems.

    Source: §1 · p. 3, §5 · p. 12

Abstract

from arXiv · show

Kernel methods provide an elegant and principled approach to nonparametric learning, but so far could hardly be used in large scale problems, since na\"ive implementations scale poorly with data size. Recent advances have shown the benefits of a number of algorithmic ideas, for example combining optimization, numerical linear algebra and random projections. Here, we push these efforts further to develop and test a solver that takes full advantage of GPU hardware. Towards this end, we designed a preconditioned gradient solver for kernel methods exploiting both GPU acceleration and parallelization with multiple GPUs, implementing out-of-core variants of common linear algebra operations to guarantee optimal hardware utilization. Further, we optimize the numerical precision of different operations and maximize efficiency of matrix-vector multiplications. As a result we can experimentally show dramatic speedups on datasets with billions of points, while still guaranteeing state of the art performance. Additionally, we make our software available as an easy to use library.

1 Introduction

Kernel methods offer theoretically grounded nonlinear learning, but their kernel matrices make large datasets computationally difficult. This paper develops and evaluates GPU-based kernel solvers designed to handle billions of points efficiently.

1 Introduction

  • Kernel methods extend classical linear models through nonlinear feature maps and can support infinite-dimensional models when their feature-space inner products are computable.Despite their theoretical appeal, implementations commonly require handling a kernel matrix quadratic in sample size, making datasets of 10^4 to 10^5 challenging.Source: §1 · pp. 1–2
  • The approach achieves state-of-the-art accuracy in minutes on benchmarks spanning millions to billions of points.The authors report favorable efficiency and accuracy compared with other large-scale implementations.Source: Figure 1 · p. 2, §1 · p. 3
  • The solver uses Nyström approximation to reduce problem size and a preconditioned gradient method, with GPU acceleration and multi-GPU parallelization.The implementation uses out-of-core linear algebra, optimizes numerical precision and matrix-vector multiplication, and is tested on datasets ranging from millions to billions of points.Source: §1 · p. 3
  • The authors release the implementation as an easy-to-use library and report that kernel methods can be deployed effectively on large-scale problems.The library is developed on top of PyTorch.Source: §1 · p. 3

2 Background

The background frames kernel learning as regularized empirical risk minimization in a feature space, with kernels allowing finite coefficient problems even for infinite-dimensional features. Nyström approximation and preconditioned conjugate gradients address the resulting computational burden, with the paper also extending the approach to logistic and related losses.

2 Background

  • Kernel methods use feature maps to define nonlinear predictors, while a kernel inner product permits finite coefficient formulations even with infinite-dimensional features.Regularized empirical risk minimization adds a penalty on the RKHS norm; the representer theorem expresses the solution through kernel evaluations.Source: §2 · pp. 3–4
  • Handling the full kernel matrix becomes prohibitive near 10^5 samples, while Nyström approximation uses m sampled inducing points and can preserve the basic bound with m = O(√n).The passage notes that multi-GPU approaches had recently scaled to 10^6 points.Source: §2 · p. 4
  • For squared loss, kernel ridge regression reduces to a linear system; the presented algorithm combines a Nyström-based preconditioner with conjugate gradients.The pseudocode applies the preconditioner and matrix-vector operator within the conjugate-gradient solve.Source: §2 · pp. 4–5
  • The method reaches the stated statistical bound in O(log n) conjugate-gradient steps, with O(n√n log n) time and O(n) memory.The same ideas extend to logistic and other self-concordant losses, including softmax, using repeated preconditioned solves for approximate Newton steps.Source: §2 · p. 5

3 Reformulating kernel solvers for multi-core/multi-GPU architectures

The solver is reworked for GPU architectures with limited local memory and costly transfers, aiming to reduce memory use and increase operations per byte. It combines in-place storage, out-of-core computation, multi-GPU execution, and transfer/compute overlap.

3 Reformulating kernel solvers for multi-core/multi-GPU architectures

  • The design targets GPU constraints—small accelerator memory and slow transfers—by reducing memory use and increasing operations per byte.The implementation focuses on minimizing the solver’s memory footprint, fitting work into GPU memory, and improving accelerator utilization.Source: §3 · p. 6

3.1 Overcoming RAM memory bottleneck

  • One m × m matrix is allocated and overwritten through preconditioner construction, avoiding separate storage of Kmm during conjugate-gradient iterations.The preconditioner occupies approximately 90% of RAM; the matrix is no longer needed in the iterations because its occurrences cancel in the transformed system.Source: §3.1 · pp. 6–7, Figure 2 · p. 6, Figure 4 · p. 7

3.2 Fitting in GPU memory and dealing with multiple GPUs

  • Out-of-core chunks let kernel-vector products and Cholesky operations fit limited GPU memory, while block-cyclic row assignment distributes Cholesky work across GPUs.For the kernel product, data blocks move to the device while the kernel matrix is computed there; the matrix itself is not moved. Cholesky updates involve inter-GPU transfers.Source: §3.2 · pp. 8–9, Figure 5 · p. 8

3.3 Optimizing data transfers and other improvements.

  • Overlapping RAM loads, saves, and GPU computation improves utilization; 64-bit summation when computing Kmm blocks helps avoid loss of precision that can make kernels non-positive-definite.For equal-duration work, the paper describes t computations completing in t + 2 time units with overlap, rather than 3t serially; KeOps handles thin matrix blocks, and sparse routines support sparse inputs.Source: Figure 3 · p. 6, §3.3 · pp. 9–10

4 Large-scale experiments

Across the tested datasets, Falkon converges faster than the compared implementations, with accuracy generally stronger than GPflow and GPyTorch. The experiments also report a 20× improvement over the earlier base implementation and include comparisons across large-scale datasets.

Optimization tests

  • 20× improvement over the base implementation of [42] is reported for the HIGGS experiment with 1×10^5 centers and 10 conjugate-gradient iterations.The improvement is visible in equal parts for preconditioner computations and the iterative conjugate-gradient algorithm.Source: §4 · p. 10, Table 1 · p. 11

Results

  • Average speedup ranges from 6× versus EigenPro to > 10× versus GPyTorch; GPflow and GPyTorch both have worse accuracy than Falkon.LogFalkon gives a small but consistent accuracy boost on binary classification problems at the expense of higher running time.Source: §4 · p. 12
  • On HIGGS, Falkon reports slightly higher error than the original Falkon library [42], attributed to low-precision floating-point numbers; other datasets show no significant performance differences.GPyTorch could not be run on TIMIT because of difficulties setting up a multi-class benchmark, which the authors say is not a software limitation.Source: §4 · p. 12

5 Conclusions

The paper presents an easy-to-use library for large-scale kernel methods, reporting strong accuracy and computational costs. It identifies other loss functions, optimization approaches, and more structured kernels as possible further developments.

5 Conclusions

  • The GPU-adapted library for large-scale kernel methods achieves what the authors describe as excellent accuracy and computational costs.The paper frames the software as a contribution toward flexible, easy-to-use machine-learning libraries.Source: §5 · p. 12
  • Further development could consider other loss functions, optimization approaches, and more structured kernels that may improve efficiency.Source: §5 · p. 12
  • Table 2 presents accuracy and running-time comparisons on large-scale datasets.Source: Table 2 · p. 13

A.1 Relative impact of performance optimizations

The ablation shows that combining precision, preconditioner, multi-GPU, and matrix-vector optimizations yields a nearly 20× improvement over the baseline. Individual steps contribute speedups across CPU, GPU, and preconditioner computations.

A.1 Relative impact of performance optimizations

  • Nearly 20× improvement over the baseline came from the combined performance optimizations.The baseline resembles the original Falkon implementation, with CPU preconditioning, float64 computations, and GPU-accelerated matrix-vector products.Source: §A.1 · p. 20
  • Using float32 for all computations gave 2× speedup on the CPU part and 3× on the GPU part.The Cholesky decomposition required care to avoid errors.Source: §A.1 · p. 20
  • The GPU preconditioner reduced its runtime from more than 20 minutes to just under 3 minutes.Adding a second GPU gave a perfect 2× speedup for CG iterations and a 1.5× speedup for the preconditioner.Source: §A.1 · p. 20
  • For HIGGS's thin data matrix, using KeOps for matrix-vector multiplication further improved its speed.The experiment used m = 1×10^5 and 10 epochs.Source: §A.1 · p. 20

A.2 Multi-GPU scalability

Falkon scales across multiple GPUs, with scaling depending on preconditioner size and operation. Matrix-vector performance also depends on data dimensionality, and the implementation switches between its routine and KeOps accordingly.

A.2 Multi-GPU scalability

  • At preconditioner size 140,000, two GPUs reached a peak speedup of around 1.8×; at low sizes, the gain was negligible.For n = 40,000, the matrix fits in GPU memory, allowing an in-place decomposition on one GPU.Source: §A.2 · p. 20, Figure 7 · p. 22
  • Falkon scales remarkably well on TAXI even with four GPUs.Figure 6 reports full-algorithm multi-GPU scalability on TAXI; the result depends on scaling both the preconditioner and conjugate-gradient iterations.Source: §A.2 · p. 20, Figure 6 · p. 21
  • LAUUM's CPU-side operation makes both GPU threads share CPU resources, while Cholesky scaling is limited by intrinsic data dependencies.The LAUUM implementation was out-of-place; porting it to the GPU was left as future work.Source: §A.2 · p. 21
  • As n increases, KeOps is approximately 10 times faster than the specialized routine, while that routine scales linearly with data dimensionality d.KeOps scales polynomially with d, so the final algorithm switches implementation based on a dimensionality threshold; the matrix-vector operation scales almost perfectly across multiple GPUs.Source: §A.2 · pp. 21–22

A.3 Additional information on the datasets

The experiments span datasets with varied sizes, feature dimensions, target types, and learning tasks. Most use standardized data and random train/test splits, with dataset-specific preprocessing and evaluation procedures.

A.3 Additional information on the datasets

  • The dataset suite was chosen to represent varied data sizes, types, and learning tasks; the default split was 80% training and 20% testing.Preprocessing generally involved cleaning and standardizing to zero mean and unit standard deviation, except where specified otherwise.Source: §A.3 · p. 22
  • HIGGS has 1.1 × 10^7 rows, 28 features, and a binary target; TIMIT has 1.2 × 10^6 rows, 440 features, and 144 classes.HIGGS reports 1 minus AUC and binary classification error; TIMIT reports classification error on a standardized subset of 57,242 samples.Source: §A.3 · p. 22
  • YELP contains 1.5 × 10^6 text reviews encoded as sparse 3-gram count vectors with 6.52 × 10^7 features.The data were not normalized, and RMSE was measured on a random 20% sample.Source: §A.3 · p. 23
  • TAXI has 1.1 × 10^9 rows, nine features, and a continuous target, with RMSE measured on a random 20% subsample.The data were standardized; taxi trips longer than five hours and trips with pickup or drop-off locations outside NYC were excluded.Source: §A.3 · p. 23
  • AIRLINE has 5.93 × 10^6 rows and eight features; its continuous-target task reports MSE on random test sets containing 33% of the data.AIRLINE-CLS thresholds the target at 0 and reports classification error and 1 minus AUC on 100,000 random test points.Source: §A.3 · p. 23
  • MSD has 5.1 × 10^5 rows and 90 features, with relative error reported on a standard test set of 51,630 samples.Its target is continuous, and the data were standardized.Source: §A.3 · p. 24
  • SUSY has 5 × 10^6 rows and 18 features, with a binary target and classification error reported on 20% of the data.The data were standardized to zero mean and unit standard deviation.Source: §A.3 · p. 24

A.4 Additional information on the experimental settings

The experimental settings describe tuning and implementation choices for EigenPro2, GPFlow, GPyTorch, Falkon, and Logistic Falkon. The GP baselines use different optimization and covariance choices, while Falkon variants tune kernel and regularization parameters.

A.4 Additional information on the experimental settings

  • EigenPro2 tunes the preconditioner's complexity but is limited to datasets fitting entirely in GPU memory.Some experiments therefore used subsampled data; some also required manually lowering the automatically inferred learning rate by a fixed integer.Source: §A.4 · p. 24
  • GPFlow used SVGP with likelihoods matched to regression, binary, and multiclass tasks, and all final experiments used diagonal covariance.Full covariance was rarely beneficial and increased training times slightly; using float32 made natural-gradient optimization unstable, leaving the double-precision/natural-gradient tradeoff unresolved.Source: §A.4 · p. 24
  • GPFlow whitening cut per-epoch runtime by about 2× while slowing convergence by around the same amount, so it was used only for HIGGS.The overall running-time difference was not strongly affected.Source: §A.4 · p. 24
  • GPyTorch's SVGP could not be run on TIMIT because of issues handling multiple outputs.The experiments used single-GPU execution; unwhitened training was around 3× faster without hampering convergence, so it was selected throughout.Source: §A.4 · pp. 24–25
  • GPyTorch used natural gradients for variational parameters and Adam for other hyperparameters, with distribution parameterizations chosen by task.The natural-gradient variational distribution was used for regression, while classification used a lower-triangular parameterization.Source: §A.4 · pp. 24–25
  • Falkon tuned kernel length-scale, inducing-point count, and regularization amount using a coarse-to-fine length-scale search.Source: §A.4 · p. 25
  • Logistic Falkon tuned length-scale, inducing-point count, and regularization path, but was not very sensitive to the exact path.Setting the final λ was sufficient; many paths leading to that value worked similarly.Source: §A.4 · p. 25

A.5 Additional benchmarks

The benchmarks report that Falkon scales better across GPUs on larger datasets, while comparisons with ThunderSVM show speedups at comparable accuracy.

A.5 Additional benchmarks

  • Larger datasets scale better from one to two GPUs as startup costs and the preconditioner’s lower scaling ratio are amortized.Table 4 presents Falkon performance for one and two GPUs across the considered datasets.Source: §A.5 · p. 25, Table 4 · p. 26
  • Falkon runs 2–10 times faster than ThunderSVM, depending on problem size, with comparable accuracy on three image datasets.The benchmarks use smaller datasets because ThunderSVM solves the full SVM problem and has cubic time scaling.Source: §A.5 · p. 25, Table 5 · p. 27
  • InCoreFalkon averages a further 2× speedup over standard Falkon, but requires a dataset that fits in GPU memory with room for other computations.Table 5 compares Falkon, InCoreFalkon, and ThunderSVM; hyperparameters were tuned for approximately equal accuracy.Source: §A.5 · p. 25, Table 5 · p. 27
  • Table 3 summarizes learning-rate, training-subsampling, and Newton-step settings for the algorithm–dataset combinations.Newton steps count separate runs of the main Falkon algorithm for Logistic Falkon.Source: Table 3 · p. 26

A.6 Performance comparisons in a literature review

The literature review finds the benchmark results in line with previously reported kernel-method results, while distinguishing deep Gaussian-process results from the paper’s intended comparisons.

A.6 Performance comparisons in a literature review

  • Previously reported results with both accuracy and running time are in line with the paper’s benchmarks, as summarized in Table 6.The review omits literature results that do not report running time.Source: §A.6 · p. 27
  • Some Table 6 results exceed Falkon’s accuracy because deep Gaussian processes can learn better data representations through many parameters.The paper reports these models for completeness, not as direct kernel-method comparisons.Source: §A.6 · p. 27

B Logistic Falkon Algorithm

Logistic Falkon extends Falkon to generalized self-concordant losses using an approximate Newton method. The out-of-core algorithms tile linear algebra operations to process matrices beyond GPU memory while coordinating work across GPUs.

B Logistic Falkon Algorithm

  • The generalized self-concordant loss class includes logistic regression, robust regression, multiclass losses, and generalized linear models with bounded features.The cited examples distinguish these losses from self-concordant functions.Source: B Logistic Falkon Algorithm · pp. 27, 29
  • Table 6 surveys results on the considered datasets, listing Falkon beside other implementations with runtime and hardware where available.Source: Table 6 · p. 28
  • Generalized self-concordant losses achieve learning-rate bounds similar to kernel ridge regression under standard regularity conditions and refined regularity conditions.The result concerns the regularized empirical-risk minimizer.Source: B Logistic Falkon Algorithm · p. 29
  • The generalized self-concordant-loss solver uses Nyström reduction and an approximate Newton step whose weighted preconditioner approximates the Hessian.The approximate Hessian is not built explicitly; its factorization is used to precondition conjugate gradient.Source: B Logistic Falkon Algorithm · p. 29
  • GSC-Falkon needs only a small number of WeightedFalkon calls with few conjugate-gradient iterations to achieve optimal learning rates.The paper’s goal is to make this theoretically fast algorithm practical on available computational architectures.Source: B Logistic Falkon Algorithm · p. 30

C Out-Of-Core Algorithms

  • The out-of-core Cholesky and triangular-multiplication algorithms divide matrices into tiles that fit in GPU memory and manage tile dependencies.The triangular multiplication exploits structure with specialized operations such as LAUUM, SYRK, and TRMM.Source: C Out-Of-Core Algorithms · pp. 30–31
  • Multi-GPU work is assigned block-cyclically to balance triangular-matrix workloads, while smaller tiles improve workload balance but can reduce utilization.Source: C Out-Of-Core Algorithms · p. 30
  • Cholesky updates coordinate tile state to avoid coarse barriers, but dependencies and the serial first-tile decomposition prevent perfect parallelism.The implementation tracks tile update counts and can overlap computation with transfers, at the cost of more GPU memory and smaller possible tiles.Source: C Out-Of-Core Algorithms · pp. 31–32

Related papers

Loading 2006.10350v2…