  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.
