Pack only the essentials: Adaptive dictionary learning for kernel ridge regression Daniele Calandriello Alessandro Lazaric Michal Valko SequeL team, INRIA Lille - Nord Europe, France
arXiv:2604.22386v1 [stat.ML] 24 Apr 2026
{daniele.calandriello, alessandro.lazaric, michal.valko}@inria.fr
1
Introduction
One of the major limits of kernel ridge regression (KRR) is that for n samples storing and manipulating the kernel matrix Kn requires O(n2 ) space, which becomes rapidly unfeasible for large n. Many solutions focus on how to scale KRR by reducing its space (and time) complexity without compromising the prediction accuracy. A popular approach is to construct low-rank approximations of the kernel matrix by randomly selecting a subset of m columns from Kn , thus reducing the space complexity to O(nm). These methods, often referred to as Nyström approximations, mostly differ in the distribution used to sample the columns of Kn and the construction of low-rank approximations. Both of these choices significantly affect the accuracy of the resulting approximation [5]. Bach [2] showed that uniform sampling preserves the prediction accuracy of KRR (up to ε) only when the number of columns m is proportional to the maximum degree of freedom of the kernel matrix. This may require sampling O(n) columns in datasets with high coherence [4] (i.e., a kernel matrix with weakly correlated columns). Alternatively, Alaoui and Mahoney [1] showed that sampling columns according to their ridge leverage scores (RLS) (i.e., a measure of the influence of a point on the regression) produces an accurate Nyström approximation with only a number of columns m proportional to the average degrees of freedom of the matrix, called effective dimension. Unfortunately, the complexity of computing RLS is comparable to solving KRR itself, making this approach unfeasible. However, Alaoui and Mahoney [1] proposed a fast method to compute a constant-factor approximation of the RLS and showed that accuracy and space complexity are close to the case of sampling with exact RLS at the cost of an extra dependency on the inverse of the minimal eigenvalue of the kernel matrix. Unfortunately, the minimal eigenvalue can be arbitrarily small in many problems. Calandriello et al. [3] addressed this issue by processing the dataset incrementally and updating estimates of the ridge leverage scores, effective dimension, and Nyström approximations on-the-fly. Although the space complexity of the resulting algorithm (INK-E STIMATE) does not depend on the minimal eigenvalue anymore, it introduces a dependency on the largest eigenvalue of Kn , which in the worst case can be as big as n. This can potentially reduce the advantage of the method. In this paper we introduce SQUEAK, a new algorithm that builds on INK-E STIMATE, but uses unnormalized RLS and an improved RLS estimator. As a consequence, the algorithm is simpler, does not need to compute an estimate of the effective dimension for normalization, and it achieves a space complexity that is only a constant factor worse than sampling according to the exact RLS.
2
Background
Notation. We use curly capital letters A for collections and |A| for the number of entries in A, upper-case bold letters A for matrices and lower-case bold letters a for vectors. We denote by [A]ij and [a]i the (i, j) element of a matrix and i-th element of a vector respectively. We use en,i ∈ Rn for the i-th indicator vector of dimension n. Finally, the set of the first n integers is [n] := {1, . . . , n}. Kernel regression. We consider a regression dataset D = {(xt , yt )}nt=1 , with input xt ∈ X ⊆ Rd and output yt = f ⋆ (xt ) + ηt , where f ⋆ is an unknown target function and ηt is a zero-mean i.i.d. noise. We denote by K : X × X → R a positive definite kernel function. Given the first t samples in D, the kernel matrix Kt ∈ Rt×t is obtained as [Kt ]ij = K(xi , xj ) for any i, j ∈ [t] and we denote by yt , ft⋆ ∈ Rt the vectors with components yi and f ⋆ (xi ), i ∈ [t]. Whenever a new point xt+1 arrives, the kernel matrix Kt+1 ∈ Rt+1×t+1 is obtained by bordering Kt as " # Kt kt+1 Kt+1 = (1) T kt+1 kt+1
where kt+1 ∈ Rt is such that [kt+1 ]i = K(xt+1 , xi ) for any i ∈ [t] and kt+1 = K(xt+1 , xt+1 ). b t ∈ Rt that minimizes the At any time t, the objective of kernel regression is to find the vector w regularized quadratic loss b t = arg min ∥yt − Kt w∥2 + µ∥w∥2 = (Kt + µI)−1 yt , w
(2)
w
b t achieves a near-optimal where µ ∈ R is a regularization parameter. If µ is properly tuned, then w b t ) = Eη ||ft⋆ − Kt w b t ||22 . Nonetheless, the computation of the final w b n requires O(n3 ) risk R(w 2 time and O(n ) space, which is infeasible for large datasets. Nyström approximation. A common approach to reduce the complexity is to (randomly) select m columns of Kt according to some distribution pt = {pt,i }ti=1 and construct the dictionary It = {(ij , kt,ij , pet,ij )}m j=1 , which contains the set of indices ij ∈ [t], the corresponding columns and their weights. Given a dictionary It , the regularized Nyström approximation of Kt is obtained as e t = Kt St (ST Kt St + γIm )−1 ST Kt , K t t
(3)
where the selection matrix St ∈ Rt×m is defined as St = [(qe pt,i1 )−1/2 et,i1 , . . . , (qe pt,im )−1/2 et,im ], e t can be q is a constant, and γ is a regularization term (possibly different from µ). At this point, K −1 e t + µIt ) yt efficiently using block inversion, reducing the complexity e t = (K used to compute w from O(n3 ) to O(nm2 + m3 ) time and from O(n2 ) to O(nm) space. e t is strictly related to the distribution pt used to construct Ridge leverage scores. The accuracy of K the dictionary It . In particular, Alaoui and Mahoney [1] showed that sampling according to the γ-ridge leverage scores (RLS) of Kt leads to an accurate Nyström approximation. Definition 1. Given Kt = Ut Λt UT t , the γ-ridge leverage score (RLS) of column i ∈ [t] is −1 −1 τt,i = kT et,i = eT et,i , (4) t,i (Kt + γIt ) t,i Kt (Kt + γIt ) Pt Furthermore, the effective dimension of the kernel is defined as deff (γ)t = i=1 τt,i . P Similar to standard leverage scores (i.e., j [U ]2i,j ), RLSs measure the importance of each point xi for the kernel regression. Furthermore, the sum of the RLSs is the effective dimension deff (γ)t , which measures the intrinsic capacity of the kernel Kt when its spectrum is soft-thresholded by a regularization γ. Using RLS in constructing a Nyström approximation leads to the following result. Proposition 1 (Alaoui and Mahoney [1]). Let ε ∈ [0, 1] and In be the dictionary built with m columns selected proportionally to RLSs {τn,i }. If m = O( ε12 deff (γ)n log( nδ )), the Nyström approximation e n is a γ-approximation of Kt , that is 0 ⪯ Kt − K e t ⪯ γ Kt (Kt + γI)−1 ⪯ γ I and the risk K 1−ε 1−ε γ 1 e t is R(w e t ) ≤ (1 + µ 1−ε )R(w b t ). of w
Unfortunately, computing exact RLS requires storing Kn , and has the same O(n2 ) space requirement as solving Eq. 2. In the next section, we introduce SQUEAK, an RLS-based incremental algorithm able to preserve the same accuracy of Prop. 1 without requiring to know the RLS in advance, and that generates a dictionary only a constant factor larger than exact RLS sampling.
3
Incremental Nyström approximation with ridge leverage scores
SQUEAK (Alg. 1) builds on the INK-E STIMATE algorithm [3] with the major algorithmic difference that the sampling probabilities are computed directly on estimates τt,i without renormalizing them by an estimate of deff (γ)t . SQUEAK introduces two key elements: 1) an improved, accurate estimator of the RLS and 2) an incremental sampling scheme for the construction of the dictionary It . 1) Estimation of RLS. We introduce an RLS estimatorPthat improves on [3], showing that it can be efficiently computed. At any time t, let Qt = i Qt,i be the number of columns |It | contained in the dictionary at time t, and St ∈ Rt×Qt the selection matrix constructed so far. Let St+1 ∈ R(t+1)×(Qt +q) be constructed as [St , (q)−1/2 et+1,t+1 , . . . , (q)−1/2 et+1,t+1 ] by adding q copies of et+1,t+1 to the selection matrix. Denoting α = (1 + ε)/(1 − ε), we define the RLS estimator as T −1 T 1+ε τet+1,i = ki,i − kt+1,i S S Kt+1 S + γI S kt+1,i . (5) αγ 2
Algorithm 1 The SQUEAK algorithm Input: Dataset D, regularization γ, µ, q e n, w en Output: K 1: Initialize I0 as empty, pe1,0 = 1 2: for t = 0, . . . , n − 1 do 3: Receive new column [kt+1 , kt+1 ] 4: Compute α-approximate RLS {e τt+1,i : i ∈ It ∪ {t + 1}}, using It , [kt+1 , kt+1 ], and Eq. 5 5: Set pet+1,i = max {min {e τt+1,i , pet,i } , pet,i /2} 6: Initialize It+1 = ∅ 7: for all j ∈ {1, . . . , t} do 8: Qt,j = |{i = j : i ∈ It }| 9: if Qt,j ̸= 0 then 10: Qt+1,j ∼ B(e pt+1,j /e pt,j , Qt,j ) S HRINK D ICT-U PDATE 11: Add Qt+1,j copies of (j, kt+1,j , pet+1,j ) to It+1 . 12: end if 13: end for 14: Qt+1,t+1 ∼ B(e pt+1,t+1 , q) E XPAND 15: Add Qt+1,t+1 copies of (t + 1, kt+1,t+1 , pet+1,t+1 ) to It+1 16: end for e n using In and Eq. 3 17: Compute K e n , yn e n using K 18: Compute w
If Qt ≥ q, then τet+1,i can be computed in O(Q3t ) time (O(Qt ) to compute kt+1,i S and O(Q3t ) to invert the inner matrix) and O(Q2t ) space. If Qt < q the same applies with q replacing Qt . Furthermore, we have the following guarantee. e t . Then for all i such Lemma 1. Assume that the dictionary It induces a γ-approximate kernel K that i ∈ {It ∪ {t + 1}}, τet+1,i computed using Eq. 5 is an α-approximation of the RLS τt,i , that is τt+1,i (γ)/α ≤ τet+1,i ≤ τt+1,i (γ). 2) Sequential sampling. At each time step t, SQUEAK receives a new column [kt+1 , kt+1 ]. This can be implemented either by having a separate algorithm that constructs each column sequentially and streams it to SQUEAK, or by storing just the samples (with an additional O(td) space complexity) and computing the column once. Adding a new column to the matrix can either decrease the importance of columns already observed (i.e., if they are correlated to the new column) or leave it unchanged (i.e., if they are orthogonal) and thus the RLS evolves as τt+1,i ≤ τt,i [3, App. A, Lem. 4]. In the D ICT-U PDATE loop, the dictionary is updated to reflect the change in importance of old columns (e.g., pt,i = τt,i may decrease) and to add the new column proportionally to its RLS τt+1,t+1 . The dictionary It , and the new column are used to compute new approximate RLS τet+1,i as in Eq. 5, which in turn define the new sampling probabilities pet+1,i . The D ICT-U PDATE phase is composed of two steps. For each index i ∈ [t], the S HRINK step counts the number of copies Qt,i present in It , and then draws a sample from the binomial B(e pt+1,i /e pt,i , Qt,i ), where taking pet+1,i = min {e τt+1,i , pet,i } ensures that the binomial probability at L10 is well defined. The more pet+1,i is lower than pet,i , the more Qt+1,i will be lower than Qt,i . If the probability pet+1,i continues to decrease over time, it is also possible that Qt+1,i is decreased to zero, and column i is completely dropped from the dictionary. Intuitively, the S HRINK step stochastically reduces the size of the dictionary to reflect the reductions of the RLSs. Conversely, the E XPAND step adds the new column to the dictionary with a number of copies (from 0 to q) which depends on its estimated relevance pet+1,t+1 . Unlike in [3], the approximate probabilities pet,i are not obtained by normalizing the approximate τet,i by an estimate of the effective dimension and thus they do not necessarily sum to one. Yet, we guarantee that pet,i ≤ pt,i ≤ 1 by construction. Note that SQUEAK never estimates again the RLS of a columns dropped from It . Moreover, computing Eq. (5) requires only to construct the kernel sub-matrix for samples whose indices are in It . Therefore, if we are only interested in estimating the approximate e t , SQUEAK is the first RLS sampling algorithm that can RLS τet,i and not the regression weights w operate in a single pass over the dataset (store and access only the samples in It instead of the whole Dt ), without ever constructing the whole matrix. Thm. 1 guarantees that SQUEAK succeeds in e n with high probability. returning a γ-approximate matrix K 3
|In | (Total space = O(n|In |))
Acc. loss
Increm.
3
n
n
1
N/A
3 ndmax 2 d n n + max ε ε n(|In |)2
dmax,n ε
(1 + 4ε)
No
2
(1 + 2ε)
No
2
(1 + 2ε)
Yes
(1 + 2ε)2
Yes
2
N/A
Time E XACT Bach [2] Alaoui and Mahoney [1] Calandriello et al. [3] SQUEAK RLS- SAMPLING
λmin +nµε λmin −nµε
2 2 λ2 max n deff (γ)n γ2 ε2 n2 deff (γ)2 n ε2 ndeff (γ)2 n ε2
n) deff (γ)n + Tr(K µε
λmax deff (γ)n γ ε2 deff (γ)n ε2 deff (γ)n ε2
(1 + 2ε)
Table 1: Comparison of Nyström methods. λmax and λmin refer to largest and smallest eigenvalues of Kn . 1+ε Theorem 1. Let α = 1−ε and γ > 1. For any 0 ≤ ε ≤ 1, and 0 ≤ δ ≤ 1, if we run Alg. 1 with parameter q = O( εα2 log( nδ )) to compute a sequence of random dictionaries It each with a random number of entries |It |, then with probability 1 − δ, for all iterations t ∈ [n] e t (Eq. 3) associated with It is a γ-approximation of Kt . (1) The Nyström approximation K P (2) The number of stored columns is |It | = i Qt,i ≤ O(qdeff (γ)t ) ≤ O( εα2 deff (γ)n log( nδ )). 1 e t satisfies R(w e t ) ≤ (1 + µγ 1−ε b t ). (3) The solution w )R(w
As the previous theorem holds for any t ∈ [n], SQUEAK has any-time guarantees on its space complexity, approximation, and risk performance. In fact, (1) combined with Lem. 1 shows that, at all steps, τet,i are α-approximate RLSs estimates. Since adding a column to Kt can only increase the effective dimension (i.e., deff (γ)t ≤ deff (γ)t+1 ) [3, App. A, Lem. 5], from (2) we see that the number of columns stored by SQUEAK over iterations never exceeds the budget O(deff (γ)n log(n)) required by sampling columns according to the exact RLS computed over the whole dataset. Notice that this is obtained by automatically increasing the dictionary size (and space occupation) over time to adapt to the growth in effective dimension of the data, which does not need to be known in advance. Furthermore, if the size of the dictionary grows too large w.r.t. the memory available, we can still terminate the algorithm knowing that the intermediate dictionary returned is a good approximation of the part of dataset processed. We can also restart the process with a larger γ, since deff (γ)n is inversely proportional to γ. The tradeoffs of this approach are quantified by (3), which shows that e t incur a risk only a factor roughly (1 + γ/µ) away from the corresponding exact all solutions w b t . This means that choosing a small γ < µ allows to achieve a risk close to the exact solution w solution for a large range of µ, at the cost of increasing the space, while larger γ require less space but it may prevent from tuning µ optimally. Finally, it is important to notice that even in the worst case deff (γ)n = n, SQUEAK requires only log(n) more space than storing the whole matrix.
4
Discussion
Table 1 compares several Nyström approximation methods w.r.t. their space complexity and risk. For all methods, we omit O(log(n)) factors. The space complexity of uniform sampling [2] scales P with the maximal degree of freedom dmax . Since dmax = n maxi τn,i ≥ i τn,i = deff (γ)n , uniform sampling is often outperformed by RLS sampling. While Alaoui and Mahoney [1] also sample according to RLS, their two-pass estimator is not very accurate. In particular, the first pass requires to sample O (nµε/(λmin − nµε)) columns, which quickly grows above n2 when λmin becomes small. Finally, [3] require that the maximum dictionary size is fixed in advance, which implies some knowledge of the effective dimensions deff (γ)n , and requires estimating both τet,i and deeff (γ)t . In particular, this extra estimation effort causes an additional λmax /γ factor to appear in the space complexity. This factor cannot be easily estimated, and causes a space complexity of n3 in the worst case. We also include RLS- SAMPLING, a fictitious algorithm that receives the exact RLS in input, as an ideal baseline for all RLS sampling algorithms. From the table, we can therefore see that SQUEAK achieves the same space complexity (up to constant factors) as knowing the RLS in advance. Moreover, although in this paper we only considered fixed design KRR, γ-approximation e n are commonly used in similar problems such as random design KRR, or Kernel guarantees for K PCA. Finally, with a more careful analysis, we can generalize SQUEAK and its guarantees to the distributed setting, where multiple machines construct dictionaries in parallel on separate datasets, and then recursively merge them to construct a dictionary for the union of the datasets. 4
References [1] Ahmed El Alaoui and Michael W. Mahoney. Fast randomized kernel methods with statistical guarantees. In Neural Information Processing Systems, 2015. [2] Francis Bach. Sharp analysis of low-rank kernel matrix approximations. In Conference on Learning Theory, 2013. [3] Daniele Calandriello, Alessandro Lazaric, and Michal Valko. Analysis of Nyström method with sequential ridge leverage scores. In Uncertainty in Artificial Intelligence, 2016. [4] Alex Gittens and Michael W Mahoney. Revisiting the Nyström method for improved large-scale machine learning. In International Conference on Machine Learning, 2013. [5] Alessandro Rudi, Raffaello Camoriano, and Lorenzo Rosasco. Less is more: Nyström computational regularization. In Neural Information Processing Systems, 2015.
5