BLOKS: Scaling Kernel Solvers with Sketched Block Preconditioning

Giordano Vitale; Giacomo Meanti; Marco Rando; Blanca Ayuso de Dios; Lorenzo Rosasco

cs.LGOct 9, 2026

BLOKS: Scaling Kernel Solvers with Sketched Block Preconditioning

Authors: Giordano Vitale, Giacomo Meanti, Marco Rando, Blanca Ayuso de Dios, Lorenzo Rosasco

Organizations: MaLGa Center, DIBRIS, Universit`a degli Studi di Genova, Genoa, Italy · Universit´e Cˆote d’Azur, INRIA, CNRS, LJAD, Nice, France · Dipartimento di Matematica e Applicazioni, Universit`a degli Studi Milano–Bicocca, Italy · Istituto di Matematica Applicata e Tecnologie Informatiche “E. Magenes”, CNR, Pavia, Italy ·

Abstract

Kernel Ridge Regression (KRR) is a theoretically sound algorithm for nonparametric learning, which is also often used to investigate properties of deep learning. However, the exact KRR formulation incurs time and memory costs preventing applicability to large datasets, motivating the development of approximate solvers. These algorithms rely on (i) iterative methods, to avoid the computational burden of direct solutions; (ii) sketching or batching, to reduce the problem size; and (iii) preconditioning, to accelerate convergence. While crucial for fast convergence, preconditioning incurs memory costs that can be managed by lowering the number of free parameters of the model, potentially at the expense of expressivity and hence accuracy. To address this challenge, we propose BLOKS, a Conjugate Gradient-based solver for Nyström-sketched KRR built around a novel preconditioner that combines a block Jacobi approximation with sketching. The BLOKS preconditioner has memory cost independent of both the dataset size and the number of parameters, allowing scaling to large models. We provide a memory-efficient, GPU-accelerated implementation and show that BLOKS achieves excellent accuracy while retaining favorable time/memory complexity. Indeed, on standard large-scale datasets used to benchmark kernel solvers, BLOKS achieves comparable or better performances than state-of-the-art methods.

Figures & tables

Panel of Figure 1 — Figure 1: First block of $H$ (above), and sketched approximation (below). Constants omitted for visual clarity.

	Taxi ( $n\approx 10^{9}$ )		Higgs ( $n\approx 10^{7}$ )		Airlines ( $n\approx 10^{6}$ )		CIFAR10-5M ( $n\approx 10^{6}$ )		COMET MC ( $n\approx 10^{7}$ )
Algorithm	time	RMSE	time	AUC	time	MSE	time	accuracy	time	AUC
Falkon	$2\text{\,}\mathrm{h}$	${307.9}$	$7\text{\,}\mathrm{min}$	${0.8194}$	$2\text{\,}\mathrm{min}$	${0.720}$	$14\text{\,}\mathrm{min}$	${87.7\%}$	$46\text{\,}\mathrm{s}$	${0.9911}$
Joker	$7\text{\,}\mathrm{h}$	${307.8}$	$47\text{\,}\mathrm{min}$	${0.8200}$	$1\text{\,}\mathrm{h}$	${0.718}$	$1\text{\,}\mathrm{h}$	${87.7\%}$	$56\text{\,}\mathrm{min}$	${0.9923}$
ASkotch	–	–	$4\text{\,}\mathrm{h}$	${0.8237}$	$30\text{\,}\mathrm{min}$	${0.717}$	–	–	$34\text{\,}\mathrm{min}$	${0.9941}$
EigenPro 4	–	–	$4\text{\,}\mathrm{h}$	${0.8104}$	$2\text{\,}\mathrm{h}$	${0.791}$	$1\text{\,}\mathrm{h}$	${86.9\%}$	$2\text{\,}\mathrm{h}$	${0.9759}$
Bloks	$5\text{\,}\mathrm{h}$	${307.0}$	$2\text{\,}\mathrm{h}$	${0.8254}$	$47\text{\,}\mathrm{min}$	${0.706}$	$2\text{\,}\mathrm{h}$	${88.0\%}$	$31\text{\,}\mathrm{min}$	${0.9943}$

Table 1: Large-scale dataset comparisons. For each method/dataset combination we report the training time and a relevant measure of error or accuracy on that dataset. Standard deviations (over 5 data splits) are reported in Appendix D for space reasons. Bold green cells indicate the best training time and the best accuracy on a dataset; light-green cells indicate accuracy within 1 standard deviation of the best.

	MillionSongs ( $n\approx 10^{5}$ )		Benzene ( $n\approx 10^{5}$ )		MiniBooNE ( $n\approx 10^{5}$ )		FashionMNIST ( $n\approx 10^{4}$ )		SVHN ( $n\approx 10^{4}$ )
Algorithm	time	rel. error $\times 10^{3}$	time	MSE $\times 10^{2}$	time	AUC	time	accuracy	time	accuracy
Falkon	$40\text{\,}\mathrm{s}$	${4.492}$	$27\text{\,}\mathrm{s}$	${1.60}$	$18\text{\,}\mathrm{s}$	${0.9733}$	$30\text{\,}\mathrm{s}$	${91.2\%}$	$1\text{\,}\mathrm{min}$	${80.4\%}$
Joker	$9\text{\,}\mathrm{min}$	${4.510}$	$4\text{\,}\mathrm{min}$	${1.72}$	$1\text{\,}\mathrm{min}$	${0.9717}$	$32\text{\,}\mathrm{s}$	${90.4\%}$	$52\text{\,}\mathrm{s}$	${76.7\%}$
ASkotch	$13\text{\,}\mathrm{min}$	${4.475}$	$11\text{\,}\mathrm{min}$	${0.16}$	$5\text{\,}\mathrm{min}$	${0.9741}$	$1\text{\,}\mathrm{h}$	${91.0\%}$	$1\text{\,}\mathrm{h}$	${79.9\%}$
EigenPro 4	$19\text{\,}\mathrm{min}$	${4.502}$	$7\text{\,}\mathrm{min}$	${1.99}$	$5\text{\,}\mathrm{min}$	${0.9689}$	$1\text{\,}\mathrm{min}$	${90.8\%}$	$2\text{\,}\mathrm{min}$	${77.8\%}$
Bloks	$4\text{\,}\mathrm{min}$	${4.471}$	$4\text{\,}\mathrm{min}$	${0.15}$	$40\text{\,}\mathrm{s}$	${0.9718}$	$14\text{\,}\mathrm{s}$	${91.2\%}$	$1\text{\,}\mathrm{min}$	${80.2\%}$

Table 3: Medium/small-scale dataset comparisons. Cell colors as in Table 1 ; standard deviations in Appendix D .

Appendix figures & tables11 assets

Supplementary material from the paper’s appendix.

Appendix

Algorithm	$n$	$m$	$m_{0}$	Memory
Falkon	data size	model size (Nyström)	N/A	$\mathcal{O}(m^{2})$
Krill	data size	model size (Nyström)	embedding dimension	$\mathcal{O}(m^{2})$
EigenPro 4	data size	model size (Nyström)	preconditioner rank	$\Omega(m_{0}^{2})$
Joker	data size	model size (RFF)	block size	$\mathcal{O}(m_{0}^{2})$
ASkotch	data and model size	N/A (Full KRR)	rank of Nys. approx.	$\mathcal{O}(m_{0}^{2})$
Bloks	data size	model size (Nyström)	block size	$\mathcal{O}(m_{0}^{2})$

Appendix

Table 4: Parameters and memory requirements of state of the art algorithms.

Panel of Figure 5 — Figure 5: First block of $\widetilde{H}^{\prime}$ (above), and sketched approximation (below). Constants omitted for visual clarity.

Preconditioner	$j$ -th block	FLOPS	Interactions captured by the $j$ -th block
$M_{\text{BJ}}$ ( 28 )	$\left(K_{nm_{0}}^{(j)\top}K_{nm_{0}}^{(j)}+\lambda n\widetilde{K}_{j}\right)^{-1}$	$\mathcal{O}(nm_{0}^{2}+m_{0}^{3})$	$\{x_{1},\dots,x_{n}\}$ and $\{\widetilde{x}_{(j-1)m_{0}+1},\dots,\widetilde{x}_{m_{jm_{0}}}\}$
$M_{\text{BJF}}$ ( 30 )	$\left(\frac{n}{m}K_{mm_{0}}^{(j)\top}K_{mm_{0}}^{(j)}+\lambda n\widetilde{K}_{j}\right)^{-1}$	$\mathcal{O}(mm_{0}^{2}+m_{0}^{3})$	$\{\widetilde{x}_{1},\dots,\widetilde{x}_{m}\}$ and $\{\widetilde{x}_{(j-1)m_{0}+1},\dots,\widetilde{x}_{m_{jm_{0}}}\}$
$M$ ( 13 )	$\left(\frac{nq}{m_{0}}\widetilde{K}_{j}^{2}+\lambda nq\widetilde{K}_{j}\right)^{-1}$	$\mathcal{O}(m_{0}^{3})$	$\{\widetilde{x}_{1}^{(j)},\dots,\widetilde{x}_{m_{j}}^{(j)}\}$ and $\{\widetilde{x}_{(j-1)m_{0}+1},\dots,\widetilde{x}_{m_{jm_{0}}}\}$

Appendix

Table 5: Difference in information and time cost between block preconditioners. FLOPS refers to the number of operations required to compute and invert a single block. The total cost of the preconditioner is given by multiplying this cost by the number of blocks

q=m/m_{0}

	HIGGS			AIRLINES
	$1-\mathrm{AUC}$	time (s)	time w.r.t. Bloks	MSE	time (s)	time w.r.t. Bloks
Block Jacobi	$0.2358$	$103.17$	$\times 12.98$	$0.9627$	$87.36$	$\times 15.29$
Block Jacobi - Falkon	$0.2349$	$17.56$	$\times 2.21$	$0.9629$	$6.35$	$\times 1.11$
Bloks	$0.2335$	$7.95$	$\times 1.00$	$0.9629$	$5.72$	$\times 1.00$

Appendix

Table 6: Comparison of different block preconditioning strategies on subsamples of the HIGGS and AIRLINES datasets. Means over

5

-fold cross validation, standard deviations omitted for space.

Figure 6: Subsample of the HIGGS dataset,

5

-fold cross-validation. (a) relative error

\relerr_{H}

(solid) and relative residual

\relres_{M}

(dashed); (b) test

1-\mathrm{AUC}

. Curves are averages over the folds, shaded bands are the standard deviations.

Figure 7: Subsample of the AIRLINES dataset,

5

-fold cross-validation. (a) relative error

\mathrm{relerr}_{H}

(solid) and relative residual

\mathrm{relres}_{M}

(dashed); (b) test MSE. All curves are averages over the folds. The shaded bands in (a) are the standard deviations across the folds; in (b) we do not report them, as they are of order

4\times 10^{-2}

, much larger than the gaps between the curves, and would overlap into a single transparent block, collapsing the curves in a single flat line.

Table 12

Table 13

dataset	algorithm	kernel	$\sigma$	$\lambda$	$m$	$m_{0}$	epochs	other
Higgs	Falkon	Gaussian	3.5	$3\text{\times}{10}^{-8}$	$100\,000$	N/A	20
	Joker	Laplacian	median ( $10.7$ )	$\frac{1}{n}$	$200\,000$	512	2	SVM
	ASkotch	Gaussian	3.5	$3\text{\times}{10}^{-8}$	Full	$100$	20
	EigenPro 4	Laplacian †	15	$0$	$500\,000$	$10\,000$ / $1000$	10
	Bloks	Gaussian	3.5	$3\text{\times}{10}^{-8}$	$400\,000$	$100\,000$	$20$	recompute blocks
COMET MC	Falkon	Gaussian	0.25	$1\text{\times}{10}^{-8}$	$50\,000$	N/A	20

Appendix

Table 7: Parameters used for each dataset and algorithm. Here

\sigma

denotes the kernel bandwidth and

\lambda

the regularization parameters. For each algorithm, the meaning of the parameters

m

and

m_{0}

is explained in Table 4 .

		train	test	test			train	test	test
dataset	algorithm	time (s)	time (s)	error	dataset	algorithm	time (s)	time (s)	error
Taxi	Falkon	$8970.89$	$61.20$	$307.9$\pm$0.1$	MillionSongs	Falkon	$39.69$	$0.12$	$0.004\,492$
	Joker	$25\,131.92$	$541.45$	$307.8$\pm$0.5$		Joker	$537.77$	$2.04$	$0.004\,510$
	ASkotch	Exceeds time limits				ASkotch	$809.03$	$1.14$	$0.004\,475$
	EigenPro 4	Exceeds time limits				EigenPro 4	$1151.67$	$1.78$	$0.004\,502$
	Bloks	$18\,330.89$	$121.54$	$307.0$\pm$0.1$		Bloks	$254.47$	$0.46$	$0.004\,471$

Appendix

Table 8: Training time, test time and test error for each dataset and algorithm. For the test error we report the mean and standard deviation over

5

data splits, except for MillionSongs. The test error metric depends on the dataset: RMSE (Taxi), MSE (Airlines, Benzene),

1-\mathrm{AUC}

(Higgs, COMET MC, MiniBooNE), classification error (CIFAR10-5M, FashionMNIST, SVHN) and relative error (MillionSongs).

$m$	$200\mathrm{k}$	$400\mathrm{k}$	$800\mathrm{k}$	$1.6\mathrm{M}$	$3.2\mathrm{M}$
train time (s)	$3140.71$	$6268.30$	$12\,574.60$	$25\,940.69$	$52\,901.50$
test time (s)	$3.46$	$6.89$	$13.76$	$27.48$	$54.92$
AUC	$0.8232$\pm$0.0003$	$0.8254$\pm$0.0003$	$0.8266$\pm$0.0003$	$0.8272$\pm$0.0002$	$0.8277$\pm$0.0003$

Appendix

Table 9: Training time, test time and AUC for each model size

m

of Bloks on the HIGGS dataset. For the AUC we report the mean and standard deviation over

5

data splits.

Nyström Method

See all 15 papers

Apr 24, 2026

Pack only the essentials: Adaptive dictionary learning for kernel ridge regression

One of the major limits of kernel ridge regression (KRR) is that storing and manipulating the kernel matrix K_n for n samples requires O(n^2) space, which rapidly becomes unfeasible for large n. Nystrom approximations reduce the space complexity to O(nm) by sampling m columns from K_n. Uniform sampling preserves KRR accuracy (up to epsilon) only when m is proportional to the maximum degree of freedom of K_n, which may require O(n) columns for datasets with high coherence. Sampling columns according to their ridge leverage scores (RLS) gives accurate Nystrom approximations with m proportional to the effective dimension, but computing exact RLS also requires O(n^2) space. (Calandriello et al. 2016) propose INK-Estimate, an algorithm that processes the dataset incrementally and updates RLS, effective dimension, and Nystrom approximations on-the-fly. Its space complexity scales with the effective dimension but introduces a dependency on the largest eigenvalue of K_n, which in the worst case is O(n). In this paper we introduce SQUEAK, a new algorithm that builds on INK-Estimate but uses unnormalized RLS. As a consequence, the algorithm is simpler, does not need to estimate the effective dimension for normalization, and achieves a space complexity that is only a constant factor worse than exact RLS sampling.

Daniele Calandriello, Alessandro Lazaric, Michal Valko

SequeL team, INRIA Lille - Nord Europe, France

Kernel Ridge Regression Kernel Regression Importance Sampling Nyström Method Low-Rank Approximation Ridge Regression

Ask PDF

Apr 22, 2026

Analysis of Nystrom method with sequential ridge leverage scores

Large-scale kernel ridge regression (KRR) is limited by the need to store a large kernel matrix K_t. To avoid storing the entire matrix K_t, Nystrom methods subsample a subset of columns of the kernel matrix, and efficiently find an approximate KRR solution on the reconstructed matrix. The chosen subsampling distribution in turn affects the statistical and computational tradeoffs. For KRR problems, recent works show that a sampling distribution proportional to the ridge leverage scores (RLSs) provides strong reconstruction guarantees for the approximation. While exact RLSs are as difficult to compute as a KRR solution, we may be able to approximate them well enough. In this paper, we study KRR problems in a sequential setting and introduce the INK-ESTIMATE algorithm, that incrementally computes the RLSs estimates. INK-ESTIMATE maintains a small sketch of K_t, that at each step is used to compute an intermediate estimate of the RLSs. First, our sketch update does not require access to previously seen columns, and therefore a single pass over the kernel matrix is sufficient. Second, the algorithm requires a fixed, small space budget to run dependent only on the effective dimension of the kernel matrix. Finally, our sketch provides strong approximation guarantees on the distance between the true kernel matrix and its approximation, and on the statistical risk of the approximate KRR solution at any time, because all our guarantees hold at any intermediate step.

Daniele Calandriello, Alessandro Lazaric, Michal Valko

SequeL team INRIA Lille - Nord Europe

Kernel Ridge Regression Kernel Regression Nyström Method Randomized Sketching Spectral Methods Ridge Regression+1

Ask PDF

May 4, 2026

Differentiable Kernel Ridge Regression for Deep Learning Pipelines

Deep neural networks dominate modern machine learning, while alternative function approximators remain comparatively underexplored at scale. In this work, we revisit kernel methods as drop-in components for standard deep learning pipelines. We introduce Sparse Kernels (SKs), a differentiable, localized, and lazy variant of kernel ridge regression (KRR) that defers training to inference time and reduces to the solution of small local systems. We integrate SKs into PyTorch as modular layers that preserve end-to-end trainability, and we show that they expose three distinct sets of parameters -- feature representations, target values, and evaluation points -- each of which can be fixed or learned. This decomposition broadens the design space available to practitioners, enabling, in particular, training-free transfer, nonlinear probing, and hybrid kernel-neural models. Across convolutional networks, vision transformers, and reinforcement learning, SK-based modules serve two complementary roles: in some settings, they match the performance of trained neural readouts with substantially less training; in others, they augment existing models and improve their performance when used as additional components. Our results suggest that kernel methods, once made scalable and differentiable, can be readily integrated with deep learning rather than treated as a separate paradigm.

Jean-Marc Mercier, Gabriele Santin

1MPG-Partners, Paris, France · Department of Environmental Sciences, Informatics and Statistics, Ca’ Foscari University of Venice, Italy

Kernel Ridge Regression Kernel Regression Kernel Methods Ridge Regression

Ask PDF

Abstract

Nyström Method

Explore similar work

Pack only the essentials: Adaptive dictionary learning for kernel ridge regression

Analysis of Nystrom method with sequential ridge leverage scores

Differentiable Kernel Ridge Regression for Deep Learning Pipelines