Daniele Calandriello
email: daniele.calandriello@inria.fr

Alessandro Lazaric
email: alessandro.lazaric@inria.fr

Michal Valko
email: michal.valko@inria.fr

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

Introduction

One of the major limits of kernel ridge regression (KRR) is that for n samples storing and manipulating the kernel matrix K n requires O(n 2 ) 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 K n , 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 K n and the construction of low-rank approximations. Both of these choices significantly affect the accuracy of the resulting approximation [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF]. Bach [START_REF] Bach | Sharp analysis of low-rank kernel matrix approximations[END_REF] 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 [START_REF] Gittens | Revisiting the Nyström method for improved large-scale machine learning[END_REF] (i.e., a kernel matrix with weakly correlated columns). Alternatively, Alaoui and Mahoney [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] 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 [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] 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. [START_REF] Calandriello | Analysis of Nyström method with sequential ridge leverage scores[END_REF] 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-ESTIMATE) does not depend on the minimal eigenvalue anymore, it introduces a dependency on the largest eigenvalue of K n , 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-ESTIMATE, 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.

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 e n,i ∈ R n 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 = {(x t , y t )} n t=1 , with input x t ∈ X ⊆ R d and output y t = f (x t ) + η 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 K t ∈ R t×t is obtained as [K t ] ij = K(x i , x j ) for any i, j ∈ [t] and we denote by y t , f t ∈ R t the vectors with components y i and f (x i ), i ∈ [t]. Whenever a new point x t+1 arrives, the kernel matrix K t+1 ∈ R t+1×t+1 is obtained by bordering K t as

K t+1 = K t k t+1 k T t+1 k t+1 (1)
where

k t+1 ∈ R t is such that [k t+1 ] i = K(x t+1 , x i ) for any i ∈ [t] and k t+1 = K(x t+1 , x t+1 ).
At any time t, the objective of kernel regression is to find the vector w t ∈ R t that minimizes the regularized quadratic loss

w t = arg min w y t -K t w 2 + µ w 2 = (K t + µI) -1 y t , (2)
where µ ∈ R is a regularization parameter. If µ is properly tuned, then w t achieves a near-optimal risk R(

w t ) = E η ||f t -K t w t || 2 2 .
Nonetheless, the computation of the final w n requires O(n 3 ) time and O(n 2 ) space, which is infeasible for large datasets. Nyström approximation. A common approach to reduce the complexity is to (randomly) select m columns of K t according to some distribution p t = {p t,i } t i=1 and construct the dictionary

I t = {(i j , k t,ij , p t,ij )} m
j=1 , which contains the set of indices i j ∈ [t], the corresponding columns and their weights. Given a dictionary I t , the regularized Nyström approximation of K t is obtained as

K t = K t S t (S T t K t S t + γI m ) -1 S T t K t , (3)
where the selection matrix S t ∈ R t×m is defined as S t = [(q p t,i1 ) -1/2 e t,i1 , . . . , (q p t,im ) -1/2 e t,im ], q is a constant, and γ is a regularization term (possibly different from µ). At this point, K t can be used to compute w t = ( K t + µI t ) -1 y t efficiently using block inversion, reducing the complexity from O(n 3 ) to O(nm 2 + m 3 ) time and from O(n 2 ) to O(nm) space.

Ridge leverage scores. The accuracy of K t is strictly related to the distribution p t used to construct the dictionary I t . In particular, Alaoui and Mahoney [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] showed that sampling according to the γ-ridge leverage scores (RLS) of K t leads to an accurate Nyström approximation. Definition 1.

Given K t = U t Λ t U T t , the γ-ridge leverage score (RLS) of column i ∈ [t] is τ t,i = k T t,i (K t + γI t ) -1 e t,i = e T t,i K t (K t + γI t ) -1 e t,i , (4)
Furthermore, the effective dimension of the kernel is defined as

d eff (γ) t = t i=1 τ t,i .
Similar to standard leverage scores (i.e., j [U ] 2 i,j ), RLSs measure the importance of each point x i for the kernel regression. Furthermore, the sum of the RLSs is the effective dimension d eff (γ) t , which measures the intrinsic capacity of the kernel K t 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 [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF]). Let ε ∈ [0, 1] and I n be the dictionary built with m columns selected proportionally to RLSs {τ n,i }.

If m = O( 1 ε 2 d eff (γ) n log( n δ )), the Nyström approxi- mation K n is a γ-approximation of K t , that is 0 K t -K t γ 1-ε K t (K t + γI) -1 γ 1-ε I and the risk of w t is R( w t ) ≤ (1 + γ µ 1 1-ε )R( w t ).
Unfortunately, computing exact RLS requires storing K n , and has the same O(n 2 ) 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-ESTIMATE algorithm [START_REF] Calandriello | Analysis of Nyström method with sequential ridge leverage scores[END_REF] with the major algorithmic difference that the sampling probabilities are computed directly on estimates τ t,i without renormalizing them by an estimate of d eff (γ) 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 I t .

1) Estimation of RLS. We introduce an RLS estimator that improves on [START_REF] Calandriello | Analysis of Nyström method with sequential ridge leverage scores[END_REF], showing that it can be efficiently computed. At any time t, let Q t = i Q t,i be the number of columns |I t | contained in the dictionary at time t, and S t ∈ R t×Qt the selection matrix constructed so far. Let S t+1 ∈ R (t+1)×(Qt+q) be constructed as [S t , (q) -1/2 e t+1,t+1 , . . . , (q) -1/2 e t+1,t+1 ] by adding q copies of e t+1,t+1 to the selection matrix. Denoting α = (1 + ε)/(1 -ε), we define the RLS estimator as

τ t+1,i = 1 + ε αγ k i,i -k t+1,i S S T K t+1 S + γI -1 S T k t+1,i . (5)
Algorithm 1 The SQUEAK algorithm Input: Dataset D, regularization γ, µ, q Output: K n , w n 1: Initialize I 0 as empty, p 1,0 = 1 2: for t = 0, . . . , n -1 do

3: Receive new column [k t+1 , k t+1 ] 4:
Compute α-approximate RLS { τ t+1,i : i ∈ I t ∪ {t + 1}}, using I t , [k t+1 , k t+1 ], and Eq. 5 5:

Set p t+1,i = max {min { τ t+1,i , p t,i } , p t,i /2}

6:

Initialize

I t+1 = ∅ 7:
for all j ∈ {1, . . . , t} do 8:

Q t,j = |{i = j : i ∈ I t }| 9: if Q t,j = 0 then 10: Q t+1,j ∼ B( p t+1,j / p t,j , Q t,j ) 11:
Add Q t+1,j copies of (j, k t+1,j , p t+1,j ) to I t+1 .

12:

end if 13:
end for 14:

Q t+1,t+1 ∼ B( p t+1,t+1 , q)

15:

Add Q t+1,t+1 copies of (t + 1, k t+1,t+1 , p t+1,t+1 ) to I t+1 16: end for 17: Compute K n using I n and Eq. 3 18: Compute w n using K n , y n

DICT-UPDATE SHRINK

EXPAND If Q t ≥ q, then τ t+1,i can be computed in O(Q 3 t ) time (O(Q t ) to compute k t+1,i S and O(Q t ) 3
to invert the inner matrix) and O(Q 2 t ) space. If Q t < q the same applies with q replacing Q t . Furthermore, we have the following guarantee. Lemma 1. Assume that the dictionary I t induces a γ-approximate kernel K t . Then for all i such that i ∈ {I t ∪ {t + 1}}, τ t+1,i computed using Eq.5 is an α-approximation of the RLS τ t,i , that is τ t+1,i (γ)/α ≤ τ t+1,i ≤ τ t+1,i (γ).

2) Sequential sampling. At each time step t, SQUEAK receives a new column [k t+1 , k t+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 DICT-UPDATE loop, the dictionary is updated to reflect the change in importance of old columns (e.g., p t,i = τ t,i may decrease) and to add the new column proportionally to its RLS τ t+1,t+1 . The dictionary I t , and the new column are used to compute new approximate RLS τ t+1,i as in Eq. 5, which in turn define the new sampling probabilities p t+1,i . The DICT-UPDATE phase is composed of two steps. For each index i ∈ [t], the SHRINK step counts the number of copies Q t,i present in I t , and then draws a sample from the binomial B( p t+1,i / p t,i , Q t,i ), where taking p t+1,i = min { τ t+1,i , p t,i } ensures that the binomial probability at L10 is well defined. The more p t+1,i is lower than p t,i , the more Q t+1,i will be lower than Q t,i . If the probability p t+1,i continues to decrease over time, it is also possible that Q t+1,i is decreased to zero, and column i is completely dropped from the dictionary. Intuitively, the SHRINK step stochastically reduces the size of the dictionary to reflect the reductions of the RLSs. Conversely, the EXPAND step add the new column to the dictionary with a number of copies (from 0 to q) which depends on its estimated relevance p t+1,t+1 . Unlike in [START_REF] Calandriello | Analysis of Nyström method with sequential ridge leverage scores[END_REF], the approximate probabilities p t,i are not obtained by normalizing the approximate τ t,i by an estimate of the effective dimension and thus they do not necessarily sum to one. Yet, we guarantee that p t,i ≤ p t,i ≤ 1 by construction. Note that SQUEAK never estimates again the RLS of a columns dropped from I t . Moreover, computing Eq. ( 5) requires only to construct the kernel sub-matrix for samples whose indices are in I t . Therefore, if we are only interested in estimating the approximate RLS τ t,i and not the regression weights w t , SQUEAK is the first RLS sampling algorithm that can operate in a single pass over the dataset (store and access only the samples in I t instead of the whole D t ), without ever constructing the whole matrix. Thm. 1 guarantees that SQUEAK succeeds in returning a γ-approximate matrix K n with high probability. Acc. loss Increm.

EXACT (1) The Nyström approximation K t (Eq. 3) associated with I t is a γ-approximation of K t .

n 3 n 1 N/A Bach [2] ndmax 2 n ε + dmax 3 n ε dmax,n ε (1 + 4ε) No Alaoui and Mahoney [1] n(|In|) 2 λ min +nµε λ min -nµε deff(γ)n + Tr(Kn) µε (1 + 2ε) 2 No Calandriello et al. [3] λ 2 max γ 2 n 2 d eff (γ) 2 n ε 2 λmax γ d eff (γ)n ε 2 (1 + 2ε) 2 Yes SQUEAK n 2 d eff (γ) 2 n ε 2 d eff (γ)n ε 2 (1 + 2ε) 2 Yes RLS-SAMPLING nd eff (γ) 2 n ε 2 d eff (γ)n ε 2 (1 + 2ε) 2 N/A
(2) The number of stored columns is

|I t | = i Q t,i ≤ O(qd eff (γ) t ) ≤ O( α ε 2 d eff (γ) n log( n δ )). (3) The solution w t satisfies R( w t ) ≤ (1 + γ µ 1 1-ε )R( w t ).
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, τ t,i are α-approximate RLSs estimates. Since adding a column to K t can only increase the effective dimension (i.e., d eff (γ) t ≤ d eff (γ) t+1 ) [3, App. A, Lem. 5], from [START_REF] Bach | Sharp analysis of low-rank kernel matrix approximations[END_REF] we see that the number of columns stored by SQUEAK over iterations never exceeds the budget O(d eff (γ) 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 d eff (γ) n is inversely proportional to γ. The tradeoffs of this approach are quantified by [START_REF] Calandriello | Analysis of Nyström method with sequential ridge leverage scores[END_REF], which shows that all solutions w t incur a risk only a factor roughly (1 + γ/µ) away from the corresponding exact solution w t . This means that choosing a small γ < µ allows to achieve a risk close to the exact 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 d eff (γ) n = n, SQUEAK requires only log(n) more space than storing the whole matrix.

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 [START_REF] Bach | Sharp analysis of low-rank kernel matrix approximations[END_REF] scales with the maximal degree of freedom d max . Since d max = n max i τ n,i ≥ i τ n,i = d eff (γ) n , uniform sampling is often outperformed by RLS sampling. While Alaoui and Mahoney [START_REF] Alaoui | Fast randomized kernel methods with statistical guarantees[END_REF] 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 n 2 when λ min becomes small. Finally, [START_REF] Calandriello | Analysis of Nyström method with sequential ridge leverage scores[END_REF] require that the maximum dictionary size is fixed in advance, which implies some knowledge of the effective dimensions d eff (γ) n , and requires estimating both τ t,i and d eff (γ) 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 n 3 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 guarantees for K n are commonly used in similar problems such as random design KRR, or Kernel 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.

  Time |In| (Total space = O(n|In|))

Table 1 :

 1 Comparison of Nyström methods. λmax and λmin refer to largest and smallest eigenvalues of Kn. Theorem 1. Let α = 1+ε 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 I t each with a random number of entries |I t |, then with probability 1 -δ, for all iterations t ∈ [n]