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

Pack only the essentials: Adaptive dictionary learning for kernel ridge regression
Apr 24, 2026Daniele Calandriello, Alessandro Lazaric, Michal ValkoKernel Ridge Regression Kernel Regression
SequeL team, INRIA Lille - Nord Europe, France
Analysis of Nystrom method with sequential ridge leverage scores
Apr 22, 2026Daniele Calandriello, Alessandro Lazaric, Michal ValkoKernel Ridge Regression Kernel Regression
SequeL team INRIA Lille - Nord Europe
Differentiable Kernel Ridge Regression for Deep Learning Pipelines
May 4, 2026Jean-Marc Mercier, Gabriele SantinKernel Ridge Regression Kernel Regression
1MPG-Partners, Paris, France · Department of Environmental Sciences, Informatics and Statistics, Ca’ Foscari University of Venice, Italy

Apr 24, 2026J/K move · Enter open · S save

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

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

Open paper AskPDF arXiv

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.

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