Joker: Joint Optimization Framework for Lightweight Kernel Machines
Authors: Junhong Zhang
Zhihui Lai
Venue: ICML 2025 poster (2025)
Keywords: unified optimization
large-scale kernel method
lightweight model
Primary area: general_machine_learning->kernel_methods
TL;DR: We propose a unified framework for kernel methods (not limited to KRR) with low computation requirements.

Abstract: Kernel methods are powerful tools for nonlinear learning with well-established theory.
  The scalability issue has been their long-standing challenge.
  Despite the existing success,
  there are two limitations in large-scale kernel methods:
  (i) The memory overhead is too high for users to afford;
  (ii) existing efforts mainly focus on kernel ridge regression (KRR), while other models lack study.
  In this paper, we propose **Joker**, a joint optimization framework for diverse kernel models, including KRR, logistic regression, and support vector machines.
  We design a dual block coordinate descent method with trust region (DBCD-TR) and adopt kernel approximation with randomized features,
  leading to low memory costs and high efficiency in large-scale learning.
  Experiments show that **Joker** saves up to 90% memory but achieves comparable training time and performance (or even better) than the state-of-the-art methods.

Joker: Joint Optimization Framework for Lightweight Kernel Machines
Junhong Zhang'
Zhihui Lai 12
Abstract
ability challenges with large datasets due the induced large-scale quadratic programming (Wen al., 2018).
Kernel methods are powerful tools for nonlinear learning with well-established theory. The scala- bility issue has been their long-standing challenge. Despite the existing success, there are two limita- tions in large-scale kernel methods: (i) The mem- ory overhead is too high for users to afford; (ii) existing efforts mainly focus on kernel ridge re- gression (KRR), while other models lack study. In this paper, we propose Joker, a joint optimization framework for diverse kernel models, including KRR, logistic regression, and support vector ma- chines. We design a dual block coordinate descent method with trust region (DBCD-TH) (T) and adopt kernel approximation with randomized features, leading to low memory costs and high efficiency in large-scale learning. Experiments show that Joker saves up to 90% memory but achieves com- parable training time and performance (or even better) than the state-of-the-ar methods.
Numerous researchers have dedicated their efforts to tack- ling the scalability issues. Nyström method (Williams & Seeger, 2000) and randomized features (Rahimi & Recht, 2007) are two representative techniques to reduce the com- putation of kernel methods. Nevertheless, researchers pre- ferred the former one (Yang et al., 2012). Based on Nyström method, Rudi et al. (2017) proposed KRR model, Falkon, for large-scale data. Its variants (Rudi et al., 2018; Mcanti et al., 2020) still hold the state-of-the-ar performance and efficiency for large-scale kernel methods. However, the outstanding performance is priced at the expense of high memory usage. Meanti et al. (2020) suggests that a good performance of Falkon- methods necessitates large number of Nyström centers. Yet, they require O(M)) mem- ory to store the preconditioner for acceleration, where M is the number of Nyström centers. This leads to dilemma be- tween memory consumption and performance. To perform Falkon with M = 1.2 x 10 on the HIGGS dataset as in (Meanti et al., 2020) using single precision, one needs at least 55GB of memory, which is unaffordable for most users. In other words, memory resources become a bottleneck for large kernel machines (Abedsoltan al., 2023).
1. Introduction
Kernel methods, formidable paradigm in nonlinear leamn ing, stand alongside deep learning as a prominent approach in the machine learning landscape. Recently, theoretical progress that connects the modern deep neural network with the kernel machines has been made (Jacot et al., 2018; Geif- man et al., 2020; Chen Xu, 2020), highlighting the great potential of kernel methods. In the big data era, kernel methods' scalability has become a concern. Kernel ridge re- gression (KRR) is the simplest kernel model. However, it is associated with linear system potentially requiring (nz)) operations, where n is the size of the training set. Another instance is the kernel support vector machine (SVM) (Cortes & Vapnik, 1995). Despite its prevalence, it also faces scal-
We also note that the existing large-scale kernel methods mainly focus on KRR while other models lack study. For classification models, to our knowledge, the recent refer- ences about large-scale kernel logistic regression (KLR) and SVM! ' which are two classifiers commonly used in prac- tice, are relatively sparse. (Marteau- Ferey et al., 2019) and (Wen al., 2018) the two latest influential works about KLR and SVM, respectively. However, they still meet the bottlenecks of either high memory or time consumption.
In short, this paper aims to address two issues of the large- scale kernel methods: expensive memory requirement and limited model diversity. We propose a joint optimization framework for kernel machines, Joker, which significantly reduces memory footprint while keeping comparable per- formance and efficiency to the state of the art, and covers a broad range of models with a unified optimization scheme. The contributions, some backgrounds, and a detailed review
'College of Computer Science and Software Engineering, Shenzhen University. Shenzhen, 518060, China 'Guangdong Provincial Key Laboratory of Intelligent Information Process- ing, Shenzhen, 518060, China. Correspondence to: Zhihui Lai <lai hui 163.com>.
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
'SVM includes SVC and SVR, where "C" and "R" stand for classification and regression, respectively. Here refer to SVC.

---

Joker
Let I S [n]] be index sets, and K KI.J be the submatrix of K with rows indexed by I and columns indexed by J.
Table 1. Comparison of large-scale kernel methods on the HIGGS dataset. The "hybrid" means Joker supports both exact and inexact kernel models.
We use KI, == Kz,jn and _ ;= K for short. Let
((.,)) -) be a proper and convex loss function. A generic kernel machine can be written as:
Method
Type
Memory
Time
Models
Falkon
inexact
50GB
1 hour
KRR
LogFalkon
EigenPro3
inexact
inexact
Rs 7GB
50GB
> 15 hours
hour
KLR
KRR
min 1 + Σ Cy,, (0, (0, (((x))))) OE 2
A
(I)
LIBSVM
exact
< 2GB
1 week
SVC, SVR
ThunderSVM
exact
8GB
week
SVC, SVR
where ((x)) is a nonlinear map satisfying K(x, = (y(()), (((x')xx The Representer Theorem (Schölkopf & Smola, 1998) states that the optimal solution 0 satisfies 0- -i= pfz() for some ax E RR. Exact models seeks such a" without kernel approximation. For exam- ple, KRR uses e(y, 9) (y 8)² /2 and corresponds to closed-form solution a' = (K + XI) y, however, which is computationally expensive for large-scale data. Inexact models approximate the kernel functions with less computa- tional cost. For example, the Nyström based method utilizes approximation K K K KI.: with a specific index set I [n]] satisfying I M and M < It equivalent to restricting 0 in subspace: 0- = Eiez 0x;4(x:)) making dimension of parameters largely reduced.
Joker (ours)
hybrid
x 2GB
x 1 hour
see Table 2
of related work are presented in the rest of this section.
1.1. Contributions
Our contributions are summarized three breakthroughs:
We develop Joker, a joint optimiza- tion framework for diverse kernel models beyond KRR, presenting a general scheme support a wide range of large-scale kernel machines.
Unified scheme:
Low consumption: Joker is lightweight. We propose novel solver, dual block coordinate descent with trust region (DBCD-TR). Owing to its low space and time complexity, the hardware requirement of large-scale learming is significantly reduced.
1.3. Related work
1.3.1. LARGE-SCALE KERNEL MACHINES
We start from KRR. Ma & Belkin (2017) proposed Eigen- Pro to solve the exact KRR using preconditioned gradi- ent descent. Despite its success, its computational cost is still prohibitive. The Nyström method has become a dom- inant technique for inexact kernel models. Among them, Falkon (Rudi et al., 2017) leverages the conjugate gradient descent with Cholesky- -based preconditioner, which made significant progress large-scale KRR with remarkable performance. The falkon-based method (Rudi et al., 2018; Carteau-Ferey et al., 2019) usually requires ((M) mem- ory to store the Cholesky factor, where M is the number of Nyström centers, resulting in memory limitation in practice. To this end, EigenPro3 avoids high space complexity by projected gradient descent (Abedsoltan et al., 2023), how- ever, compromising with high time complexity per iteration. Hence, Abedsoltan et al. (2024) proposed delayed projec- tion technique to improve its efficiency. Finally, despite the relevance of KRR, we do not elaborate on works about Gaus- sian Process (de G. Matthews et al., 2017; Gardner et al., 2018) in this paper as the techniques are quite different.
We implemented KRR, KLR, SVM, etc. based on Joker, and conducted experiments with a single RTX 3080 (I0GB). It shows that Joker achieves state-of-the-art performance with a low mem- ory budget and moderate training time.
Superior performance:
1.2. Preliminaries
This paper focuses on supervised learning problems. Let X' be a compact sample space, and y' be a label set. The training dataset {(x,, y:) includes n samples, where 2; E X is the feature vector of the i-th sample and yA e y is its label. For a function f, its domain is dom f = {ax E co), and its gradient and Hessian æ are f(())) and V² respectively. The Fenchel conjugate of f is fz) := suPaeddom f (y. x) f(z). The infimal convolution of f and g is (u) == inf f(u) + g(x - u). Let [m] == {1, 2, ......... n}. For two vectors a, E IR", a < b means a, S by for all i E [M]. Notation (a, b) represents the inner product between a and b.
X f(x)
Another focus of this paper lies in kernel models for clas- sification. ThunderSVM implements the kernel SVM ac- celerated with GPU for large-scale data (Wen et al., 201 8), emerging as a competent alternative to LIBSVM (Chang & Lin, 2011). However, their solver, Sequential Minimization Optimization (SMO) (Platt, 1998), becomes out of date for coping with big data in modern applications, causing high
A Mercer's kernel K((,.): : X x X R uniquely induces reproducing kernel Hilbert space (RKHS) H with the endowed inner product (5-) such that K(x,.) :) E H for all x E X, and (f, (T,(K())))))) = f(z) for all f E 7/ and x E X (Schölkopf & Smola, 1998). We denote K E Ruxn the kernel matrix with K,y- = K( 2;)), where j E [n].
2

---

Joker
Table Fenchel conjugate associated with losses. bÆnt(a)) == xr log + xx) log(1 x) denotes binary entropy, and ( log 0 = 0. In the regressss,,,, have p,g > 1 and p= + +
1.
Task
Model
Loss functions
y
Ey(u)
(7(w)
dom Eg
Regression
KRR
Square loss
IR
(v w)2/2
(u) + ey
R
Regression
Lp-Reg Reg.
Lp loss
18
lv
u/p/p
|u/9/g + ay
Regression
Li-Reg
Absolute loss
I
ey
-1 S V -6 S
Regression
Huber Reg.
Huber loss
IR
((-) /2L16 y -1) (u)
(u) + ey
Regression
SVR
e-insensitive
R
()- (i-e (u)
c/ul + vy
-1 US
Classification
SVC
Hinge
{-1.1]
max(0, yu}
ay
-1 S "y S 0
Classification
L2-SVC
Squared hinge
{ {-1.1}
max(0,1 max(0,1-yull 101-yyu)2 //
(u) + bEnt(-ey)
by S 0
Classification
KLR
Logistic
1.1
log(1 + exp(-yu)
-1 S vy <
time consumption. Based on the Newton method, (Marteau- Ferey et al., 2019; Meanti al., 2020) implements Log- Falkon, fast large-scale KLR, yet has the same memory issue with Falkon. In this paper, delve into joint op- timization scheme for lightweight kernel models, reaching a remarkable balance between performance, memory us- age, and training time. We present an intuitive comparison between the prevalent kernel methods and Joker in Table
We begin with its dual problem (Section 2.1) and present an optimization roadmap (Section 2.2). Inexact Joker is then proposed in Section 2.3.
2.1. Joint Optimization Problem by Duality
We first present direct result for the dual problem of (2). Theorem 2.1. Let E-(-) R + R+ defined as Ev( := l(y, u). Then the optimal solution of (2) is given by
1.3.2. DUAL COORDINATE DESCENT ALGORITHMS
0
a (p)))
(3)
Coordinate descent methods (CD) iteratively select one vari- able for optimization while keeping all other variables fixed, aiming to decrease the objective function incrementally. In machine learning, CD has succeeded in the fast training of linear SVM (Hsieh et al., 2008; Dai & Qiu, 2023). Moreover, Shalev-Shwartz & Zhang (2013) proposed dual optimiza- tion framework with coordinate ascent methods for linear models. A Newten-based dual CD method is investigated for unconstrained optimization (Qu et al., 2016). In Joker, we design a unified optimization scheme also using duality. Nevertheless, our work has crucial differences from theirs in two aspects. (i) We focus on kernel models, which are usually ill ill-conditioned and more challenging in optimiza- tion than the linear models. (ii) We employ block coordinate descent (BCD) for more efficient optimization. Note that Tu et al. (2016) and Rathore et al. (2024) also proposed to address KRR with the BCD algorithm. Notably, (Rathore et al., 2024) introduces the Nesterov acceleration technique for BCD and proves its convergence under proper condi- tions. Beyond their scope, we investigate a general class of kernel models, including but not limited KRR, KLR, and SVM, and tackle the intricate constrained optimization.
a= arg min = Ka 1 I
E
n
E (-dos);
(4)
aEl? 2
where St = far: -da; E dom Eyg i E [n]]] is the feasible
region, (.) is the Fenchel conjugate of (.))), and K
kernel matrix with Ki₃ ((((x,), ((xx)))
The proof is shown in Appendix A. We first present some important properties. (i) Problem (4) is convex due to the convexity of Eg (-). (ii) The strong duality holds according to Slater's condition (Boyd et al., 2004), meaning that if a is optimal for (4), then 0* is optimal for (2). (iii) Based on the closeness and the convexity of the conjugate function (Beck, 2017, Theorem 4.3), the domain of Ey. (.) should be closed interval. Therefore, the constraints in (4) simply box constraints, i.e., ag E
with -oo S Tf < T!" S cc. (iv) Problem (4) can be better conditioned than the primal form (2), as the dual Hessian is linearly dependent on K and the primal Hessian is a quadratic form of K. Since K is usually ill-conditioned, dual optimization should converge faster than the primal one. This benefit is also noted by Tu et al. (2016) and Rathore et al. (2024).
2. Joker
Theorem 2.1 covers a wide range of kemne models using different loss functions, We list representative ones in Table 2. Notably, problem (4) can easily adapt to some sophisti- (2) cated loss functions, whereas its primal problem (2) may be tricky. To state this, we introduce the following proposition:
Joker focus on convex problem (1) reformulated as:
min 0 2
1
Σ Elv. (0, ((xx)))).
3

---

Joker
Algorithm 1: Trust region dwwiceddifeerettillllee f) Input :Initial point æg, block kernel matrix KB,B,
2.2. Dual Block Coordinate Descent with Trust Region
Coordinate descent (CD) is particularly suitable for the op- timization of (4) for its separable structure (Nutini et al., 2022). The key idea of CD is optimize one variable while fixing others in each iteration, exhibiting inexpensive com- putation and storage. Despite its simplicity, it suffers from slow convergence for large-scale kernel machines since it only updates one variable at a time. Thus, we propose up- dating multiple variables (i.e., a block) simultaneously to leverage the merit of parallel computing, leading to our core solver, Dual Block Coordinate Descent with Trust Region (DBCD-TR). We show the main idea and leave some details in Appendix B.
kernel gradient g = Ks. a, function f()).
Max region size Amax. threshold 17 E [0, 1/4),
tolerance > 0, and max iteration TTR-
0g,0 f agg, Ao+ f Amax/ /4;
T [(ag).upper,
TB 4 (agg).ooor:::
for k = 1 to TTR do
gs < g V((BB))), Qx + K + 42 f(aB,))
8k_- TCG-Stel haug (hiaaugiuia0 gh, B, TB rb),
ie. Algorithm 2;
P& + (J(aBB.)) J (a + sk))/ ()) pu (s:));
if Pk < 0.5 then
Ax+1 + Ax/4:
In each iteration of block coordinate descent, we pick an index set B = (ii, ... 11 } C [n] and Bg == |m]]] [m]]/B is the index set associated with the fixed variables. We define f(ag) == ee fy (-Xax)/ / for simplicity, where og is the subvector of a indexed by B. Minimizing (4) w.r.t ons while fixing au we have:
else if pk > 0.75 and - - ||s|| < e then
Ak+1 t (24 Amax): ):
else
AK+1 f AR:
end
if Pk: > then
0g < 0g + 8;
1
T min J(as) = ag Kg + agg ag
+ f(as),
g - 9 + KBB8;
2
end
8.1. T S ag < T
(6)
end
There are two challenges in solving subproblem (6): (i) The smoothness of f(as) is unknown; (ii) even with smooth ((ss) the box constraints are still tricky. For (ii), the projected Newton method (Gafni & Bertsekas, 1984) is a good choice to cope with the box constraints. However, its difficulty is the step size tuning. An improper step size could degrade the projected Newton step (Schmidt et al., 2011; Nutini et al., 2022).
Proposition 2.2. If E.) can be written as (u) 2
2... ... (g,,)) then the dual problem of (2) is
1 T Ka+ + min a aEl?
Su.r (-Xos)-
(5)
with Q (a: -ao; E N. dom Sg (1, E [m] [s)).
To solve (6), we employ the trust region method (Sorensen, 1982), which is an iterative algorithm for efficient optimiza- tion. Denote Org, the k-th iterate of the trust region proce- dure. The next iterate is given by 08, == OB,, + 8 with a proper step 8% restricted in the trust region {s : || < Dx]. where A is the radius. Define a quadratic model function:
Proof. The result is simply by induction and the property of infimal convolution: (f) = Ei + E (Beck, 2017, Theorem 4.16).
L
For example, Huber loss is used for robust regression, which is defined as
1 TQs..
(=())==((<BB)) + g s +
(7)
2
1/2 (u 2)2 8 u - yl 2 Otherwise.
if Lu wa < 8,
such that (x((s) =2 J(OBB + +) for || S As. Then we find the step 8% by minimizing an easier quadratic function (s). The trust region method can overcome two chal- lenges in our problem (6) naturally. Firstly, it is compatible with non-smooth optimization (Baraldi & Kouri, 2025), and secondly, it implicitly tunes the step size by adjusting the radius A in each iteration, giving a crucial safeguard for convergence (Baraldi & Kouri, 2024).
((y,w) w)
may be uneasy to optimize its primal problem directly. However, as it can be written as the infimal convolution form: (x)) = (- -1) (u), the Fenchel conjugate ()) is easily derived via Proposition 2.2. Interestingly, Huber loss and square loss (KRR) have almost the same dual problem, with the only difference in the feasible region. The commonalities of the two models imply that they can be solved similarly.
We first study the most common scenario where f is twice differentiable, which is applicable to most losses listed in Table 2. The procedure described below is summarized in
4

---

Joker
Algorithm 2: Truncated CG-Steihaug
heuristic truncated CG-Steihaug method, shown in Al- gorithm 2. The key is to terminate the CG procedure if violates the box constraints or goes beyond the trust region boundary, and finally project s back to the feasible region. Compared with the projected Newton method suggested in (Gafni & Bertsekas, 1984; Nutini et al., 2022), Algorithm computes a truncated CG step instead of the exact inverse Q 9. and thus is more efficient than their projected New- ton step. Moreover, Algorithm 2 elegantly embeds the step size tuning into the projected Newton by the nature of the trust region, which also eases the implementation.
Input Quadratic model Q, 9, initial guess ag,
region size A, bounds TV Tg. tolerance e. Output: Truncated CG step
O <- -9, <- r; r2old c rT;; while not converged do w +- (=2old)((d (d) nnxt t s + wd:
if || |next] > A then Determines w' > 0 such that ||s + w'dl A:
8 - 8 w'd and break;
else if Smext violates box constraints then break;
Now we consider the complexity of DBCD-TR per iteration. Its space complexity is only 2 lying in the storage of Kg,8. 8. each iteration, Algorithm is repeated TTR times, resulting in a time complexity of O(TTRTcG [B])), where TeG is the number of CG iterations. In our implementation, TTR S 50 and TcG is generally tiny (S 10) due to CG's fast convergence and the truncations. The most expensive com- putation lies in K5,: with time complexity of ((nd (B)) supposing a single kernel evaluation costs ((d) time. This brings us to the next major issue be overcome.
end
8 f Suext; r wQd, r2new t- rr:
Trnnew S then
| break; end
v + r2new/:2old, t r vd;
end
return max (min =g- as), as):
2.3. Inexact Joker via Randomized Features
Algorithm 1. In this case,
construct the quadratic model (+) using Taylor expansion. That is, let Qs K KB, + V (OxBx)) and 9s = Kg.0 a f((BB))) in (7).
A blueprint to solve the pivotal problem (4) is presented in Section 2.2. Kg,8 and ,,, are computed exactly in DBCD-TR, we obtain exact Joker. However, its bot- tleneck occurs in computing K,..a for the (Bd) time complexity, which becomes heavy computational burden when 2 10 To alleviate it, we propose inexact Joker. The goal is to approximate the exact kernel evaluations with a umite-dimensional mapping (-)) : X + R M with M < n such that ((xx, x') s <b(x) *b('). The following discussion and our implementation are based on the Ran- dom Fourier feature (RFF) (Rahimi & Recht, 2007), a well- studied kernel approximation approach. Using Bochner's theorem, one can write a shift-invariant kernel (i.e, the value of KC x²) only depends on x x') as
Then the "next step" 85 is given by:
8k= argmin Sag S
2 Q&s + I 8, s.t.
| SA
However, 8y may not be the good enough step when (s) does not approximate J(OBB, + 8) well. Considering this, we should evaluate the quality of sæ and only apply 0B,6411 = @B,. + 8& for the qualified step, and keep un- moved otherwise, i.e., 0B = 0g, A generic trust region procedure evaluates 85 by the ratio:
Pk: = x(BB,)) J(OB, + 8k) |r (0) (*k)
(9)
K(xx, x') = / c c (e- dpp = Ew
A large Pk suggests that the objective J(-) is decreased sufficiently, and we tend to accept 8t. Specifically, 8 is qualified if >> 17 where E (0. 1/4] is the acceptance threshold. On the other hand, we can know that (s) cannot approximate J(ag, + s) when Pk is small. In this case, the radius of the trust region should be reduced, e.g., A&+1 == A+/4. Oppositely, we can enlarge the radius in the next iteration to allow a larger step when Pk: is large.
where = V-1 denotes imaginary unit, Pie is a proper probability distribution associated with the kernel K., and Sw (x) = exp(-juu T x)) Considering IC is real-valued, we can further derive
E KC(,x') w~Pp Up.
[2 coss x+b) cos(w" x' +b)].
where Uj0,2x) denotes the uniform distribution on [0, 2m).
The subsequent issue is to find an effective solver for (8). The vanilla trust region problem can be solved efficiently with the conjugate gradient method (CG) proposed by Stei- haug (1983). However, for (8), extra consideration should be taken on the box constraints. To this end, we propose a
Based on the Monte Carlo method, (xx)) can be defined as:
(<z) = v 2 M
cos (os(Wx + b), by Up.227) wi PK.
(10)
5

---

Joker
Algorithm 3: DBCD-TR for problem (4)
Table 3. Complexity comparison, |8| < M < n. (Log)Falkon
Input- :Kernel K,, feasible region Ω, function C-),, parameter X, block size b, max iteration T.
has extra setup time of ((M²) and post-process time of (M²))
Mapa is the memory cost of EigenPro2 (Ma & Belkin, 2019).
Output: multiplier a (predictor 0)
Methods
Space
Operations per epoch
Initialize a E ft partition [x] into blocks B,,."" BB if using inexact model then
(Log)Falkon
M + Mdd
n.dd
Sample W E Rxx PK. b RM,, U/0,2a) Define map (b(xr) == cos War + b):
EigenPro3
Maga + Md
nMd + [ड] 0(M2)
Exact Joker
[B
n²d [8) (1+
Initialize f such that 0 = Ei= (za)(
Inexact Joker
[B/ + Md
nAdd + (ड| (M] + Md
end
for / = 1,2, .... do
Randomly pick a block B E [BI,... .. BmIi
Inexact Joker needs extra memory of ((Md) to store W. Notably, increasing M generally produces better approxima- tion and performance, but also costs more time and storage. One can find theoretical guide for setting M in (Lanthaler & Nelsen, 2023). Up to this point, we can present the com- plete DBCD-TR procedure in Algorithm 3.
Let f(oss == Eies Sy, (-Xa;)/X: X;
using inexact model then
K8 4- (xxB) T 4b (XB):
K + (xBB) 0:
end
ag now + TrustRegionnass K8,8, Ks. f) to B solve problem (6), i.e., Algorithm! 1;
Finally. we justify why we do not consider the Nyström method in the proposed Joker, although it is the preferred approximation method in many works such as (Yang et al., 2012: Rudi et al.. 2017). In short, Nyström method is un- suitable in our scenario due its heavy computation and storage in each iteration. Assume Z E X are the Nyström centers and M is the number of centers. Then the non- linear map becomes xb(x) == L L ((,,x), where L is the Cholesky factor of K(Z, Z). Therefore, to utilize the Nyström method, one should first compute L with ((3) time complexity and store with O(M) space complexity. When updating a block B during training, the computation of eb( (xx) costs O(M Md) time, where ( is from the inverse of the triangular matrix L. In other words, the time complexity per block update is at least O(M + Md), so it is too expensive when M is large. Unfortunately, the existing works (Rudi et al., 2017: Abedsoltan et al., 2023) have shown that a large M is necessary to yield a satisfying performance. Therefore, the Nyström method does not fit Joker. In contrast, RFF is more efficient and scalable. We compare the complexity of different methods in Table 3.
using inexact model then
6 + 4b (Xx)) new ars): -
end
ag 8 end
return (and 0 if using inexact model);
where cos(-) is applied element-wise, W E Mxd b E IR M and each row w; b, are sampled from the probability distributions px((w) and Up. (.23) respectively. Gaussian ker- nel K(x, x²) = exp(-Ile -| - x' (212) is a widely used kernel in RFF, which corresponds to Px (w) NN I)),. i.e., the Gaussian with zero mean and covariance I²I. Note that although RFF is initially proposed for shift-invariant kernels, it has been sufficiently developed to diverse kernels, such as dot- product kernels and additive kernels. One can find a comprehensive summary of RFF for various kernels in the latest report of (Dai et al., 2014).
Now we use new kernel (x, x²) == (1x) (ax))) to replace the exact one IC. Then the kernel matrix becomes Ki (xx;) =b(z)), Let (xxx)) == vx(x) in Theorem 2.1,
3. Practical Instances
This section presents some example models of Joker and the associated implementation issues. Additionally, SVR is also discussed in Appendix B.2.
we obtain 0= = Ei= (;,b(ar)) leading to
KB. = Σ 1=1
xb (Xx) T (a)))) = (XB)) To.
(11)
3.1. Simple cases: KRR, Huber regression and Ly-SVC
where 2b (XB))== : (x:)))) Therefore, time complexity of evaluating K a is reduced to ( To imple- ment this, we must maintain the weight vector 0 during the optimization. Once ag is updated, 0 is updated with only BB) time complexity:
Based on Table 2 and Theorem 2.1, the dual problems of KRR, Huber, and L2C share the same form:
min - a [K + AI) y aa
(13)
aen 2
gnew == gold Σ (a - a. (2))),
(12)
where Ω = R for KRR, Ω = {a /|allo S 6/2} for
Huber, and = {a azy. 2 0) for L2-SVC That is, the
iEB
6

---

Joker
three models only have differences in the feasible region. Their similarity eases the practical implementation. Due to the simplicity of the quadratic functions, their convergence is usually fast. The empirical results (Table 4) show that the elapsed time of these three models is close, suggesting that they can achieve comparable speed. That is, Joker fills the potential efficiency gap between different kernel models.
implementation of Joker is based PyTorch without extra acceleration libraries.
The used datasets cover both regression (MSD, HEPC) and classification tasks (SUSY, HIGGS, CIFAR-5M) with the sample size ranging from 10 to 10". They are frequently used in the literature of large-scale kernel methods, and the details are summarized in Appendix C. The largest dataset, CIFAR-5M, has 10 classes and 5 million samples, each with 3072 dimensions, barely fitting within the machine's ' RAM. All datasets are normalized using the z-score trick.
3.2. A Complicated case: KLR
The case of KLR is more complicated and has some practical problems. From now on, we define ã :- y O a, where C denotes the element-wise product. We obtain the dual problem for KLR::
Using Joker's framework, we implement three regressors: KRR, Huber, and SVR, and two classifiers: (L2))SVC and KI... The compared methods are their state-of-the-art counterparts, including Falkon (KRR). EigenPro3 (KRR), LogFalkon (KLR), ThunderSVM (SVC, SVR). We uti- lize one-versus-rest (OVR) in Joker to support multi-class problems. OVR is also available in ThunderSVM, while LogFalkon only supports binary classification (Meanti et al., 2020) and thus is not applicable on CIFAR-5M, a 10-class dataset. The kernels used in the experiments include the Gaussian and Laplacian kernels, which are two prevalent choices in practice. The regularization parameter À in Joker is tuned from {2 i = -6, .7} via grid search. According to (Nutini et al., 2022), increasing the block size reduces the needed iterations for the convergence of BCD, In most cases, we employ the inexact Joker models with the block size [ = 512 considering the tradc-off between total training time and memory consumption. Except for MSD, we can employ the exact Joker with block size [6] = 2048 due to its relatively small sample size. Exact Joker needs not to store RFF random matrix W and can afford larger block size than the inexact Joker. We also increase the block size to 1024 for Joker- to accelerate convergence. Falkon and LogFalkon can be run with M = 2.5 x 104 within the limited memory, Details of further parameter settings are shown in Appendix C.
1 T Ka+ + Σ re log σ3 + -A-1- -ās) log(A- -2.))))) a.). min x ael? 2
The feasible region is Ω = fa: : a. < 1/A}. Due to lim x log x = 0, it can be extended to 0 S ãg S 1/A by defining 0 log 0 = 0. However, a problem arises when āg near the boundary. Consider the gradient and Hessian:
vf(āB), = log((.)) - Log(A a.). f(ās)) [Q. Aa.))
We can that they unbounded 0 and 1 1/A. This poses two issues. First, the quadratic model in (8) becomes ill-conditioned near the boundary, causing a rapid shrinkage of the trust region (i.e., A + 0) and slow convergence (Baraldi & Kouri, 2025). The second issue is the potential catastrophic cancellation ( (Yu et al., 2011). When ã A O, the result of 1/A ā; may be inaccurate because of the limited precision of the computer, further leading to incorrect logarithms in the gradient computation.
To mitigate these issues, we first redefine the feasible region as & < ã; < 1/\- - e, where 6 is the distance from 1/A to its next smaller floating-poit number. In this way, 1/A - â will always be precise for all feasible a; Moreover, we uti- lize a modified Hessian H., = min ( 1/2) to avoid the ill-conditioned model. This also alleviates possible catastrophic cancellation when computing the frequent op- eration Qs. Finally, we increase the block size as suggested in (Nutini et al., 2022) to overcome the slow convergence issue. We found that these strategies significantly improve the numerical stability and convergence speed.
4.1. Results and Analysis
Table 4 shows the result of the performance comparison. Note that ThunderSVM and Joker-SVM both include SVR and SVC cases, where SVR is applied to MSD and HEPC datasets, and SVC is applied to others. We observe that Joker- methods reach the lowest memory usage on all tested datasets. To obtain equivalent performance in (Meanti et al., 2020) on HIGGS (z 74.22% with M = 100s, unaf- fordable for the used machine). Joker only needs 1.9GB GPU memory, saving at least 95% storage costs compar- ing (Meanti et al., 2020). Despite low memory, Joker still outperforms the state-of-the-an methods in most cases, pre- senting almost no accuracy sacrifice. Joker's time is signifi-
4. Experiments
To highlight that Joker can obtain promising performance under limited computational budget, we conduct experi- ments on a machine with a single consumer GPU (NVIDIA RTX 3080, 10GB) and 64GB RAM. The experiments al- ways use single precision unless otherwise specified. The
*Code available at GitHub: Apple- peplee-hang/ae 2hang/3oker-paper
7

---

Joker
Table The performance comparison on regression (MSD, HEPC) and classification (SUSY, HIGGS. CIFAR-5M datasets. "!": lower is better, and vice versa for "f" "NA": meaning not applicable. "Timeout": the running time exceeds the limit of I week. "f": using the double floating-point precision. Data below each term indicates the time and peak GPU memory consumption. ThunderSVM and JoKer-SVM apply the SVR model on MSD and HEPC datasets and the SVC model on others. The top-2 results are highlighted in bold.
MSD (n 0.5M)
HEPC (n == 2M)
SUSY (n Rz 5M)
HIGGS (n 11M)
CIFAR-5M
Methods
rel. error (x (x10-9,0)
RMSE (x10-9,4)
AUC (%; f)
ACC (%,1)
AUC (%; 1)
ACC (%,
ACC (%, 1)
0.4984=0000113 (6min, 9.8GB)
5.4642+40.0188/
87.61=0.00
80.38:+0.00
80.90=0.09
73.42++0.07
68.24:10.33
Falkon
(19min, 9.9GB)
(9min, 6.0GB)
(31min, 9.9GB)
(1.9h, 9.9GB)
87.77 10.05
80.4910.000
80.43
73.0410.00
LogFalkon
NA
NA
NA
(13min, 6.5GB)
(45min, 9.9GB)
EigenPro3
4.5512-10.0047
5.0417:10.0014
86.99_0.00
80.08:10.01
79.7410.13
72.46110.06
72.94.10:00
(1.0h, 1.6GB)
(2.2h, 188B)
(2.2h, 1.7GB)
(18h, 7.0GB)
(80h, 6.9GB)
4.6431=0.0257
6.083420.0847
79.32=0.00
80.22=00001
ThunderSVM
Timeout
Timeout
(3.2h, 5.0GB)
(3.2h, 8.0GB)
(15h, 7.8GB)
4.4868 10.0011
1.7000.0.0007
87.63 0.00
80.41.10.01
81.9410.17
74.03
73.32.10.01
Joker-KRR
(35min, 0.7GB)
(31min, 1,5GB)
(25min, 1.2GB)
(1.0h, 1.9GB)
(2. 5.3GB)
4.5058=00.0100
4.7160=10.004
87.64=0.00
80.41=0.011
81.83+0.27
74.01+0.21
73.66=00.11
Joker-Huber
(36min, 0.7GB)
(36min, L5GB)
(23min, 1.2GB)
(57min, 1.9GB)
(h. 5,3GB)
4.600410.000
4.837611000342
0.01
80.44.0.02
82.40 12.00
74.41.10.05
74.47
Joker-SVM
(35min, 0.7GB)
(27min, 1,5GB)
(25min, L.2GB)
(56min, 1.9BB)
(2.0h, 5.3GB)
87.73
80.42:10.01
82.11 10000
74.17:10001
74.88
Joker-KI.R
NA
NA
(1.1h, 1.7GB)
(1.6h, 2.6GB)
(3.2h, 5.9GB)
cantly lower than EigenPro3 and ThunderSVM. Falkon- based methods are the fastest and have substantial gap com- pared to EigenPro3 and ThunderSVM. However, Joker alleviates this gap. Specifically, EigenPro3 and Thunder- SVM use at least 10x training time compared to Falkon on MSD, and Joker reduces it to 5x time. On the other hand, EigenPro3 needs 36x time (18 hours) of Falkon (0.5 hour), and Joker reduces such gap to 2x time. Thus, Joker obtains comprehensively better results under the same hardware conditions, and achieves a good trade-off between mem- ory, time, and accuracy, demonstrating that Joker is highly scalable for large-scale learning.
10
Joker-Svi SVM (ime)
Testaccuracy 196
EigenPro3 (tme)
Elapsed time
Loge a kon (tme)
loke-SVV (ACC)
EigenPro3 (ACC)
Logi e lonn (ACC)
65
60
8
Model
x10*
Figure 1. Perforrance versus the model size on HIGGS.
The performance of various models in Joker is worth noting. The exact Joker-KRR obtains the best performance on MSD using less than 1GB of memory, showing the efficacy of the exact Joker on relatively small datasets. Joker-Huber surpasses other regression models on the HEPC dataset, which may benefit from its robustness. JoK shows outstanding performance on the classification tasks. Thun- derSVM is the most time-consuming method. Particularly, converges slowly when its penalty parameter (equiv- alent to 1/A in (2)) is large. Such inefficiency may stem from the SMO solver, which updates only two variables at each iteration, significantly slower than DBCD-TR. In com- parison, Joker-SVM converges fast when utilizing large penalty e = 1/A and always outperforms ThunderSVM with less time and memory. We use a larger block size, [B] = 1024, in Joker-KLR rather than B = 512 in other
inexact Joker models to overcome the slow convergence caused by the ill- ll-connitioned Hessian. Thus, Joker-KLR takes more time and memory than other Joker models, De- spite the ill-conditioning issue, Joker-KLR never meets numerical error degrading the performance in single preci- sion, and still obtains promising accuracy.
Figure/ illustrates the time and performance versus the model size. In general, a larger model size yields better performance. Joker obtains the best accuracy even with the smallest model. The time of EigenPro3 increases rapidly with the model size, while Joker's 's time is almost unaffected. In spite of the lowest elapsed time, Falkon fails to scale up the model size due to its expensive memory cost.
A comparison of the training progress of different models is shown in Figure 2. The curves of Joker are not mono-
8

---

Joker
75
KRR
3.25
0.2
Falkon
1ggnnro33
Accursry (%)
1st epoch ends
15
primal dej
vol. loss
-0442
Oqjectvv velue
dual ozj. primal cbj
Valddeion los value
0.4
RMSE
0.15
60
Objectire
15
Valisatoon
35
SW Logf akkn
0.38
0.05
50
0.5
0.34
EigenPro3
0.32
0
500
Time (sec)
1000
1500
1000
Time (soc)
2000
3000
number teratoons
The number of ttattions =10
(a) HEPC
(b) HIGGS
(a) MSD, exact KRR
(b) HIGGS, inexact KLR
Figure 3. The primal, dual objectives, and validation loss of Joker versus the iteration steps.
Figure2 2. Test performance versus time.
tonical since Joker optimizes the dual problem, which does not guarantee the monotonic decrease of the primal objec- tive. Nonetheless, Joker eventually converges with the best performance. HEPC is a challenging regression dataset, where Falkon obtains suboptimal results using the single precision (RMSE s 17x x 10- It is mitigated by using double precision (RMSE Rs 5.4 x 10-3). The reason may be the loss of precision caused by matrix decomposition, while Joker and EigenPro3 are free of this issue and obtain better results than Falkon. Nonetheless, the slow conver- gence of EigenPro3 makes it less competitive than Joker. EigenPro3 just finished the first epoch on HIGGS when Joker and LogFalkon are nearly convergent.
& Singer, 2001) and softmax regression, whose equality- constrained dual problems are the major issue. In addition, the convergence speed of DBCD- TR is still unclear. Exist- ing theoretical results, e.g.. (Richtárik & Takáč, 2014; Nutini et al., 2015; 2022), suggest that DBCD-TR has at least a linear convergence rate. A sharper rate worth exploring. We feel it is non-trivial and leave it as future work.
A Acknowledgment
This work was supported in part by the National Natural Science Foundation of China under Grant 62476175 and Grant 62272319, and in part by the Natural Science Foun- dation of Guangdong Province (Grant 1023A151 2024A 15101117711677,, 2023B1212060078 and Technology Plamning Project of Shenzhen Mu- nicipality under Grant JCYJ2022081 (CY12024011114220009
We illustrate the convergence pattern of Joker in Figure 3. The proposed method optimizes (4), i.e., minimizing the negative of the dual objective. Therefore, the dual objective keeps increasing, which fits the pattern shown the figures. The primal objective is always larger than or equal to the dual objective due to weak duality. As mentioned before, the primal objective and loss may not decrease monotonically. Nonetheless, Figure 3 shows their major trend of decreasing. With sufficient iterations, dual and primal objectives tend to get close and merge, indicating the convergence of Joker.
and Science
110000800099 and
Impact Statement
This work improves the optimization framework of kernel machines. It significantly reduces the computational require- ment and expense of large-scale kernel methods. Compar- ing the existing methods that may require top-line hardware (c.g., A100, V100 GPU), the proposed method is more friendly for weaker hardware and could be widely applied. Moreover, the proposed techniques support various kernel models beyond KRR, presenting a more comprehensive solution for kernel learning and related applications. To conclude, this work lowers the cost of large-scale learning, widens the application scope of kernel methods, and makes it more accessible to the public.
5. Conclusion and Future Directions
Scalability is crucial issue for kernel methods. In this paper, we propose a novel optimization framework for large- scale kernel methods named Joker, breaking the memory bottleneck and pushing the development of models beyond KRR. The proposed solver, DBCD-TR, provides a modern and efficient solution to dual optimization in kernel ma- chines. We show the effectiveness of Joker on variety of kernel methods, including KRR, SVM, KLR, etc. Even with consumer hardware and limited memory, Joker obtains state-of-the-an performance within acceptable training time, making the low-cost kernel methods possible in practice.
References
Abedsoltan, A., Belkin, M., and Pandit, P. Toward large kernel models. In the 40th International Conference on Machine Learning (ICML), volume 202, pp. 61-78. PRMI., 2023.
Regarding future work, generalizations of model (2) can be explored, e.g., multi-class SVM proposed by (Crammer

---

Joker
Abedsoltan, A., Ma, S., Pandit, P., and Belkin, M. Fast training of large kernel models with delayed projections, November 2024, URL https:://// / arxiv. org/ abs/ 2411.16658
Dai, B., Xie, B., He, N., Liang, Y., Raj, A., Balcan, M.- F. F., and Song. L. Scalable Kernel Methods via Dou- bly Stochastic Gradients. In Proceedings of the 28th International Conference on Neural Information Process- ing Systems (NeurIPS), volume 27. Curran Associates, Inc., 2014. URL https: /ttp///axivvv/..///iv/.../...avm www.oorg///ao... absp//bs//// 11007/10.. 5599v4.
Aravkin, A. Y., Baraldi, R., and Orban, D. A proximal quasi- Newton trust-region method for nonsmooth regularized optimization. SIAM Journal on Optimization, 32(2):900- 929, 2022. ISSN 1052-6234. doi: 10.1 10.1137721MM4409536
de G. Matthews, A. G., van der Wilk, M., Nickson, T., Fujii, K., Boukouvalas, A., León-Villagrá, P., Ghahramani, Z., and Hensman, J. GPflow: A Gaussian process library us- ing TensorFlow. Journal of Machine Learning Research, 18(40): 2017.
Baldi, P., Sadowski, P., and Whiteson, D. Searching for exotic particles in high-energy physics with deep learning. Nature communications,, 5(1):4308, 2014.
Baraldi, R. J. and Kouri, D. P. Local convergence analysis of an inexact trust-region method for nonsmooth optimiza- tion. Optimization Letters, (63))663-688) April 2024. ISSN 1862-4480. doi:
Gafni, E. M. and Bertsekas, D. P. Two-metric projection methods for constrained optimization. SIAM Journal on Control and Optimization, 22(6))33-964,, 1984. doi: 10.1 137/0322061.
Baraldi, R. and Kouri, D. Efficient proximal sub- problem solvers for a nonsmooth trust-region method. Computational Optimization and Applications, 90(1):193- 226, January 2025. ISSN 1573-2894. $10589-024-006688x.
Gardner, J., Pleiss, G., Weinberger, K. Q.. Bindel, D., and Wilson, A. G. GPyTorch: Blackbox matrix-matrix gaus- sian process inference with GPU acceleration. In Proceed- ings of the 32nd Conference on Neural Information Pro- cessing Systems (NeurIPS), volume 31, pp. 7576-7586. Curran Associates, Inc., 2018
doi: 10.1007/
Beck, A. First-order methods in optimization. SIAM, 2017.
Garreau, D., Jitkrittum,, W., and Kanagawa, M. Large sam- ple analysis of the median heuristic, October 2018. URL https: : aarxiv.vvvv/////a///..///aaa org/ 0ttp/17770007666
Bertin-Mahieux, T., Ellis, D. P., Whitman, B., and Lamere, P. The million song dataset. In Proceedings of the 12th International Conference on Music Information Retrieval (ISMIR), 2011.
Geifman, A., Yadav, A., Kasten, Y., Galun, M., Jacobs, D., and Ronen, B. On the similarity between the Laplace and neural tangent kernels. In Proceedings of the 34th Conference Neural Information Processing Systems (NeurlPS), volume 33, pp. 1451-1461, Curran Associates, Inc., 2020.
Boyd, S., Boyd, S. P., and Vandenberghe, L. Convex opti- mization. Cambridge University press, 2004.
Chang, C.-C. and Lin, C.-J. LIBSVM: A library for sup- port vector machines. ACM Transactions on Intelligett Systems and Technology, 2:27:1 -2::2,7 2011.
Hsieh, C.-J,, Chang, K.-W., Lin, C.-J,, Keerthi, S. S., and Sundararajan, S. A dual coordinate descent method for large-scale linear SVM. In the 25th International Confer- ence Machine Learning (ICML), pp. 408-415. ACM Press, 2008
Chen, L. and Xu, S. Deep neural tangent kernel and laplace kernel have the same RKHS. In the 8th International Conference on Learning Representations (ICLR), October 2020.
Jacot, A., Gabriel, F., and Hongler, C. Neural tangent ker- nel: Convergence and generalization in neural networks. In Proceedings of the Conference on the 32nd Neural Information Processing Systems (NeurIPS), volume 31, pp. 8571 15880. Curran Associates, Inc., 2018.
Cortes, C. and Vapnik, V. Support- Networks. Ma- chine Leaming, 20(3):273-22) 297, 1995. ISSN 1573-0565.
Crammer, K. and Singer, Y. On the algorithmic implementa- tion of multiclass kernel-based vector machines. Journal of Machine Learning Research, 2:265-292, 2001.
Lanthaler, S. and Nelsen, N. H. Error Bounds for Learning with Vector-Valued Random Features. In Proceedings of the 37th Conference on Neural Information Processing Systems (NeurIPS), volume 36, pp. 71834-71861, 2023.
Dai, B. and Qiu, Y. ReHLine: Regularized composite ReLU- ReHU loss minimization with linear computation and lin- ear convergence. In Proceedings of the 37th Conference on Neural Information Processing Systems (NeurIPS), volume 36, pp. 14366-1434 Curran Associates, Inc., December 2023.
Lin, J. A., Padhy, S., Antoran, J., Tripp, A., Terenin, A., Szepesvari, C.,, Hemmandez-Lobato, M., and Janz, D. Stochastic gradient descent for gaussian processes done
10

---

Joker
right. In the 12th International Conference on Learning Representations (ICIR), 2024.
Rathore, P., Frangella, Z.,, and Udell, M. Have ASkotch:
A neat solution for large-scale kernel ridge regression. arXiv preprint arXiv:2407.100700 2024. URL https:
Ma, S. and Belkin, M. Diving into the shallows: A com- putational perspective on large-scale shallow learning. In Proceedings of the 31st Conference on Neural Infor- mation Processing Systems (NeurIPS), volume 30, pp- 3778-3787, 2017.
/ wwtp://..wgaab////4//////
Richtárik, P. and Takáë, M. Iteration complexity of random- ized block-coordinate descent methods for minimizing a composite function. Mathematical Programming. 144(1): 1-38, 2014.
Ma, S. and Belkin, M. Kernel machines that adapt to gpus for effective large batch training. In Proceedings of Ma- chine Learning and Systems (MISys), volume 1, pp. 360- 373, 2019.
Rudi, A., Carratino, L.,, and Rosasco, L. FALKON: An optimal large scale kernel method. In Proceedings of the 31st Conference on Neural Information Processing Systems (NeurIPS), volume 30, pp. 3888-38989 Curran A Associates, Inc., 2017.
Marteau- Ferey, U., Bach, F., and Rudi, A. Globally conver- gent newton methods for ill-conditioned generalized self- concordant losses. In Proceedings of the 33rd Conference on Neural Information Processing Systems (NeurIPS), volume 32, pp. 7636-7646. Curran Associates, Inc., 2019.
Rudi, A., Calandriello, D., Carratino, L., and Rosasco, L. On fast leverage score sampling and optimal learning. In Proceedings of the 32nd Conference on Neural Infor- mation Processing Systems (NeurIPS), volume 31, pp. 5672-5682. Curran Associates, Inc., 2018.
Meanti, G., Carratino, L., Rosasco, L., and Rudi, A. Kernel methods through the roof: Handling billions of points effi- ciently. In Proceedings of the 34th Conference on Neural Information Processing Systems (NeurIPS), volume 33, pp. 14410-14422. Curran Associates, Inc., 2020.
Schmidt, M. Graphical model structure learning with 11- regularization. PhD Thesis, The University of British Columbia, pp. 26, 2010.
Nakkiran, P., Neyshabur, B., and Sedghi, H. The deep bootstrap framework: Good online learners are good of- fline generalizers. the 9th International Conference on Learning Representations (ICLR), 2021.
Schmidt, M., Kim, D., and Sra, S. Projected newton-type methods in machine learning. In Optimization for Ma- chine Learning, pp. 305-329, The MIT Press, September 2011.
Nutini, J., Schmidt, M., Laradji, I., Friedlander, M., and Koepke, H. Coordinate descent converges faster with the Sounhwll rule than random selection. In the 32nd International Conference on Machine Learning (ICML), pp. 1632-1641. PMLR, June 2015.
Schölkopf, B. and Smola, A. J. Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT press, 1998.
Shalev-Shwards S. and Zhang, T. Stochastic dual coordi- nate ascent methods for regularized loss minimization. Journal of Machine Learning Research, 14(16))567-59, 2013.
Nutini, J., Laradji, I., and Schmidt, M. Let's make block coordinate descent converge faster: Faster greedy rules, message-passingg, active-set complexity, and superlinear convergence. Journal of Machine Learning Research, 23 (131):1-74, 2022.
Sorensen, D. C. Newton's Method with a Model Trust Region Modification. SIAM Journal on Numerical Anal- ysis, 19(2):40- 426, April 1982. ISSN 0036-1429 doi: 10. 10.1137/771902.
Platt, J. C. Sequential minimal optimization: A fast algo- rithm for training support vector machines. Technical report, Advances in Kernel Methods Support Vector Learning, 1998.
Steihaug, T.
The conjugate gradient method and trust re- gions in large scale optimization. SIAM Journal on Nu- merical Analysis, 20(3):626-637, 1983. ISSN 0036-1429.
Qu, Z., Richtarik, P., Takac, M., and Fercoq, O. SDNA: Stochastic dual newton ascent for empirical risk mini- mization. In the 33rd International Conference on Ma- chine Learning, pp. 1823-1832. PMLR, June 2016.
Tu, S., Roelofs, R., Venkataraman, S., and Recht, B. Large scale kernel learning using block coordinate de- scent. arXiv preprint arXiv:1602.05311, 2016. URL https: wwwp:/v..wwr//////..eee//e// 053100....0..............11....8
Rahimi, A. and Recht, B. Random features for large-scale kernel machines. In Proceedings of the 2 Conference on Neural Information Processing Systems (NeurIPS), volume 14 F, pp. 467-470, 2007.
Wen, Z., Shi, Li, Q., He, B., and Chen, J. ThunderSVM: A fast SVM library on GPUs and CPUs. Journal of Machine Learning Research, 19(21):1-5, 2018.

---

Joker
Williams, C. and Seeger, M. Using the Nyström method to speed up kernel machines. In Proceedings of the 14th Conference on Neural Information Processing Systems (NeurlPS), volume 13, pp. 682-688. MIT Press, 2000.
Yang, T., Li, Y--., Mahdavi, M., Jin, R., and Zhou, Z.-H. Nyström method random Fourier features: A theoreti- cal and empirical comparison. In Proceedings of the 26th Conference on Neural Information Processing Systems (NeurIPS), volume 25. Curran Associates, Inc., 2012.
Yu, H.-F., Huang, F.-L,, and Lin, C.-J. Dual coordi- nate descent methods for logistie regression and max- imum entropy models. Machine Learning, 85(1):41- 75, October 2011. ISSN 1573-0565. doi: 10.1007/ 110994-010-5221-8
12

---

Joker
Appendix
A. Proof of Theorem 2.1 Proof. Note that problem (2) is equivalent to:
min 00 0,u
Σ En (u:), s.t. (0, (p(x))) u;-
(A.1)
Introducing the Lagrange multipliers a E RR, we can obtain the Lagrangian function:
L(O, u a = 10,00 0) + Σ (us) + Σ a0 sp(x))) u:).
(A.2)
The first-order Karush-Kuhn-Tucck (KKT) condition gives
0" arg min L(O, u, a)= Σ a:p(()))
(A.3)
1
u arg min L(0, u, 03) = arg min
Sn, (1:) + 0:4;
arg max
Aayu Sy. (u:).
(A.4)
1=1
Recall the definition of the Fenchel conjugate, the maximum in (A.4) is Eist Ey. (-da;). Therefore, the dual problem is
max min L(0, u. a.) max ΣΣ m (,au(p(x)).
n
Σ S (-Xas), s.t. Aas E dom E
(A.5)
a u,6
Or equivalently, by applying negative sign to the objective function, have the minimization problem:
1 TI Á
well arg min a
T Ka
Σ eg (-Aas), where fa e R =>aa E dom E.
(A.6)
Equations (A.3) and (A.6) give the primal-dual relationship, which completes the proof.
B. Details of Dual Block Coordinate Descent with Trust Region (DBCD-TR)
B.1. Strategies of Block Coordinate Descent
As pointed out by Nutini al. (2022), there are several aspects consider when designing the BCD algorithm:
(i) How to obtain the block candidates?
(ii) How select the block size?
(ii) How choose the block to update?
(iv) How to update the block?
In fact, (iv) has been discussed in the main paper. Here we elaborate on the other three aspects.
(i): We adopt the fixed block strategy, i.e., the potential selected block is static during the optimization. Specifically, the block candidates
obtained by partitioning the index set (1.... n... Although Nutini et al. (2022) suggests "variable blocks" (the block is dynamically constructed during the optimization) for its faster convergence, its computational expense higher (e.g., the gradient of all samples should be maintained during optimization). So we choose the fixed block strategy for its simplicity, which also shows satisfactory efficiency in practice.
(ii): Increasing block size generally decreases the iterations to reach convergence of BCD (Nutini et al., 2022), but it may increase the computation cost and memory usage of each iteration. So we should find a trade-off between the storage
13

---

Joker
S | and the time inner- where TCC is the number of BCD iterations and /nnnr the time complexity of each inner iteration. We have shown that linner ( in the main paper. For TECD, Nutini al. (2022) points out that larger block sizes generally lead to faster convergence. (Richtárik & Takáč, 2014; Nutini et al., 2015) show a linear convergence rate of proximal gradient descent for o-strongly convex functions:
( Ln O
log )
TECD S
(A.7)
where L is the sum of Lipschitz constants of all blocks. However, this result is not tight. Borrowing the analysis from (Nutini et al., 2022), DBCD-TR adopts the "optimal updates" paradigm, which should be faster than (A.7), even possibly superlinear, but the precise rate is still unclear. Despite this, we still use (A.7) as a tool to analyze the block size setting. In our experiments, we find that [B] 512 is balanced choice for KRR, SVC, SVR and Huber. However, because of the ill-conditioned problem, i..,, extremely large L in (A.7), the convergence of Loker-KLR is slower than other Joker models. To overcome it, except for the numerical techniques mentioned in the main paper, increasing the block size is also feasible (by reducing the upper bound of TBCD). Regarding this, we set [ = 1024 for Joker-Kl. in our experiments.
(iii): Based on the fixed block strategy, we randomly select the block from the candidates with equal probability. Its advantage over the cylic approach is that it can avoid the potential bias of the fixed selection order. The alternative and probably better strategy is the greedy approach, e.g. Gauss-Southwell rules presented (Nutini al., 2022). However, they need evaluate the gradient (and possibly the Hessian) of all blocks, which is expensive in the large-scale data scenario.
B.2. Discussion of Nonsmoothness: SVR as an Example
Nonsmoothness is an intrinsic issue for trust region methods. The universal solutions include proximal gradient descent (Aravkin et al., 2022) and proximal Newton projected step (sometimes intractable) (Nutini al., 2022). é (-normm and piecewise linear functions are typical examples of nonsmooth functions. The dual problem of SVR is representative one of the former:
2 L T Ka - y a + =||ael1,,
min a
s.t. ||]||00 S X
(A.8)
commonly used trick was used (Chang & Lin, 2011) and (Schmidt, 2010). That is, take at, a 2 such that at a and obtain the following problem:
min I at.a- 2 aaat a- K(at -a ----aaaa-- -aaaa= +1 +a=) s.t. 0 at aj S X
1
(A.9)
It can be proven that the optimal solution to (A.8) and (A.9) satisfies at = max(0, oj max(0, - aa,; }. Then the SVR problem becomes a quadratic programming problem that can be solved with DBCD-TR.
On the other hand, we provide a simpler strategy. Minimizing J J(x + 8) is equivalent to:
1 T K + (Kg.a - ys) + e as + allı, s.t. |||||2 - A,
min $E 2
(A.10)
where 52 {s> |lag + 8||= I X-1 is the feasible set. Now we assume ||arg + a||2 2 ||ar/]] + sign(cg)) T The equality holds when all elements ors and og + s have the same sign, and the errors occur on the different signed elements. Therefore, the trust region size should be sufficiently small to keep the sign consistency as much possible. Then the trust region subproblem becomes the following:
min = K + s sT (KB..a a - ys + sign(ars)) s.t. ||s||2 S A, $E 2"
(A.11)
which falls back into the framework discussed in the main paper with Q K5,8,9 KB ys + sign(ag).
C. Details of Experiments
C.1. Datasets
We evaluate the models with the datasets that are commonly used in the literature on kernel methods. In these datasets, MSD, SUSY, and HIGGS were used in (Rudi et al., 2017; Meanti et al., 2020), etc. The HEPC dataset was used in (Lin
14

---

Joker
al., 2024), and the CIFAR-5M dataset was benchmarked in (Abedsoltan et al., 2023). Like most of the previous works, the z-score normalization is always applied make data with zero mean and unit variance. The information the datasets is summarized in Table A.I.
Table A.I. Summary of datasets used in experiments
Dataset
Task
d
Data Split
Metrics
MSD
Regression
5.1 10
90
90% train, 10% test
Relative Error (L)
НEPС
Regression
2.05 x 106
11
90% train, 10% test
RMSE (L)
SUSY
Binary Classification
5.0 10'
18
80% train, 20% test
AUC (t), Accuracy (t)
HIGGS
Binary Classification
1.1 107
28
80% train, 20% test
AUC (t), Accuracy (t)
CIFAR-5M
10-class Classification
5.0 x 106
3072
80% train, 20% test
Accuracy (1)
Million-song dataset (MSD) (Bertin-Mahieux et al., 2011). This dataset contains audio features for year prediction. Available at attps://archive ics httpa///atttttttte////arrrreeeee
Household Electric Power Consumption (HEPC) dataset: is the same as the "HouseElec" dataset used (Lin et al., 2024), available using the python package https://githee www/to//////oo/o twweeoreeaanssccccattttaa.a..
Supersymmetric particle classification (SUSY) dataset (Baldi et al., 2014): This a binary classification task to distinguish between supersymmetric particles and background process. Available at https://archive. uci wcc.eedu/dtatasse/27///
HIGGS dataset (Baldi et al., 2014): This a binary classification task to distinguish between the Higgs boson and the background process. Available https: /arrciiieeeeeeeeee ics uci. edu/ (attaset/280//it
CIFAR-55M dataset (Nakkiran et al., 2021): This is a generated dataset based on CIFAR-10. Available at https: //github.bb//////// wwt.://ottpp.cceeeeeeeee ci far5m.
C.2. Implementation details
We implement RFF of two widely-used kernels, the Gaussian and Laplacian, in inexact Joker. RFF is constructed by the following formulas:
/ 2 (xx) = V M
cos Wx++ b),
(A.12)
where W E Mxd and E RM are random matrices. Each element of b is always sampled from the uniform distribution U [0,20] independently. The distribution of W depends on the kernel:
Gaussian kernel: Each element Wiy is sampled independently from Gaussian distribution with zero mean and variance 2:
(-1ee 20
| x'112)
K(z, = cxp
p(w) = N(0,
(A.13)
Laplacian kernel: Each element Wij is sampled independently from the Cauchy distribution with scale 1/o:
x x (-1 p(w)=
K(x x') = exp
(A.14)
Tw + +1)
In our experiments, we implement five Joker models: KRR, Huber, SVC, SVR, and KLR. the main paper, we merge SVC and SVR into JoM. The summary is shown in Table A.2.
15

---

Joker
Table A.2. Lookup table Joker models experiments. bEnt(x) = log x + - x)) log(1 - z)) the binary entropy function
and + -).
Model
Primal problem
Dual problem
Constraint
Joker-KRR
min _1011*
1
(v (0, (ax))))
min 2a (K +1ja ay'a vJa
-00 as $00 0c
2
2A
==1
Joker-Hude
min
|1011
1
(yt (0, (ax))))
min 2
a T (K +11a- yyaa
A az
2x
==1
Joker-SVC
|1011
1
(0 (mm)))
min a T (K A1ja -vaa 2
S ag 00
2
2X
Joker-SVR
=110 min 2
Σ
[1:: (0, ((xx)),)))) | - =|-
min 2a
Ka =||aa| 1 -
- S a $ x
1
1=1
Joker-KI LL.
min
|||
Σ
log(1 + - (....)))))))))
min 50 T Ka + Σ bEnt(-Aaup
1
0 0. S A
C.3. Parameter settings
The major hyperparameters of the models are shown in the experiments in Table A.3. We tend to use the best parameters reported in literature. Median heuristic is used the bandwidth parameter of kernels, denoted "median" in Table A.3. Our empirical results show that it always produces satisfying performance, and its rationality is also supported by (Garreau et al., 2018).
16

---

Joker
Table A.3. The major hyperparameters of the models in the experiments. "NA" means applicable.
Dataset
Methods
Parameters
MSD
HЕРС
SUSY
HIGGS
CIFAR-5M
Ası
10-6
10-9
10-6
10-a
10-3
M
25000
25000
10000
25000
20000
Falkon
precision
float32
float64
float32
float32
float32
kernel,
Gaussiun, 6
Gaussian,
Gaussiun, 4
Gaussian,
Gaussian, median
epochs
20
50
20
10
50
A
10-9
10-9
M
25000
LogFalkon
#Newton step
NA
NA
NA
kernel, epochs
float32
float32
15
M
80000
80000
50000
10
EigenPro3
kernel, epochs
Laplacian, median
Laplacian, median
Laplacian, median
Laplacian, median
Laplacian, median
50
30
16
32
32 Gaussian, median
ThunderSVM kernel.
Gaussian, median
Gaussian, median
Gaussian, median
Gaussian, median
SVR
0.25
0.25
NA
NA
NA
A
1
2-7
2
2
2-7
M
Exact model
50000
10
10
2 x 10
Joker-KRR
kernel, block size
Gaussian, median
Laplacian, median
Laplacian, median
Laplacian, median
Gaussian, median
2048
512
512
512
512 40000
literations
10000
25000
10000
50000
2-7
2-5
2-7
2-7
M
Exact model
50000
10a
10a
2 x 10
5 kernel, block size
2
Joker-Huber
Gaussian, median
Laplacian, median
Laplacian. median
Laplacian, median
Gaussian, median
2048
512
512
512
512 40000
#iterations
10000
25000
10000
50000
2-7
2-5
2-7
2-7
M kernel, Gaussian, median
Exact model
10000
10a
10a
2 x 10
Laplacian. median
Laplacian. median
Laplacian. median
Gaussian, median
Joker-SVM
in SVR
0.25
0.25
NA
NA
NA
block size
2048
512
512
512
512
#iterations
10000
25000
10000
50000
40000
2
2-3
2-7
M kernel, block size
105
10a
2 x 10a
KLR
NA
NA
Laplacian. mediam
Laplacian. median
Gaussian, median
1024
1024
1024
#iterations
40000
50000
50000
17


Reviews:
[{"review_id": "fPDQ4Tdpxm", "reviewer": "Reviewer_1dCw", "summary": "The paper proposes Joker, a novel optimization scheme that aims at scaling kernel methods beyond current computational limitations. It is versatile and can handle several objective functions in a similar manner. The core idea of Joker is to solve the dual problem with a block coordinate descent with trust region. An approximate version based on random fourier features is also developed; its time and space complexity are low and Joker exhibits good performance on large scale (for kernel methods) datasets.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weak. 1: precise writing.\nWe appreciate your suggestions to make the expression more precise.\nThe following are our responses and the plan of revision.\n\n(line 94): Your expression is rigorous. We should add that $K$ satisfying Mercer condition is positive definite (Mercer condition is equivalent to positive-definiteness).\n\n(line 105, 107): Thank you for pointing out the typos. We have fixed them.\n\n(line 132): Yes it should be convexity. We have corrected it.\n\n(line 171): Your understanding on robust regression and KRR are precise.\nHowever, what we want to emphasize is the efficiency, i.e.,\nwe can solve kernel Huber regression and KRR within a similar time.\nWe rewrite this sentence to make it clearer.\n\n(line 189): Thank you. We noticed that it is imprecise expression.\nWe move this sentence after we finish the discussion on the twice-differentiable $f$.\n\n(line 230):\nYou may misunderstand $T_{TR}$.\nIn fact, $T_{TR}$ is the iteration number of trust region procedure in solving eq.(6)(*maxIter* in Algorithm 3),\nwhich is small because the trust region method converges fast (typically superlinear).\n\nThe number of block update iterations ($T$ in Algorithm 2) is large, typically many times of $n/|B|$.\nIn your example with 512000 samples and a block size of 512, needing 1000 times to go through all data (one pass),\nwe may set $T=10000$ and exactly go through 10 passes.\nAnother example in our experiments, we use $T=50000$ with block size 512 for Joker-SVM on HIGGS (see Table A.3), roughly complete $40000*512/5M\\approx 4$ passes.\n\nWe clarify three iteration procedures:\n- Algorithm 2: DBCD-TR outer-loop, iteration times $T\\sim O(n/|B|)$.\n- Algorithm 3 (in Appendix B): called by Algorithm 2, iteration times $T_{TR}\\leq50$.\n- Algorithm 1: called by Algorithm 3, iteration times $T_{CG}\\leq10$.\n\n(line 254): We regard \"directly searched for on the feature space\" as optimizing eq.(2) directly using the approximation feature $\\varphi(x):=\\psi(x)$.\nThis is related to the primal-based methods involving $M$ variables and can be related to [1] and Falkon (using Nystrom feature).\nTo obtain a promising convergence rate these methods usually utilize second-order algorithms.\nThe challenges occurs when handling the Hessian or precondition matrix that may cause $O(M^2)$ computation.\nIn contrast, the proposed Joker optimizes the dual variables but not eq.(2) itself.\nwe only maintains the KKT condition $\\theta=\\sum_{i=1}^n\\alpha_i\\psi(x_i)$ once the dual variables updated using eq.(11).\nThis procedure aims to reduce the computation complexity of $K_{B,:}\\alpha_B$, as stated in Section 2.3.\nCompared to the direct method,\nthe proposed dual optimization leverages the separable structure of eq.(4),\nallowing the block coordinate descent to process a small working set in each iteration.\nIn this way, the computation of processing Hessian is reduced to $O(|B|^2)$.\n\n(line 272): We sincerely thank your suggestion.\nThe provided reference proves that a generalization error of $O(n^{-1/2})$ can be obtained using $O(n^{1/2})$-dimensional RFF.\nThis can be more useful in practice to guide the selection of $M$.\nWe will update the our citation.\n## Weak. 2: limitation of RFF?\nRFF is proposed for shift-invariant kernels initially.\nHowever, it has been developed to diverse kernels like dot-product kernels and additive kernels.\nTable 1 of [2] gives a good summary.\nFastfood [3] also presents a principled way to construct random features beyond shift-invariant kernels.\nRegarding NTKs, [4] also provides a fast algorithm to obtain their explicit features.\nThese results allow for obtaining RFF of a broad range of kernels and it is no longer a limitation of Joker. \n\n# Other\n## Access to Code.\nResponse to \"However no code is provided...\":\nWe actually provided the code.\nThe anonymous GitHub link is at the bottom of page 6 of the PDF.\nWe welcome you to review the code and give suggestions.\n## Related work.\nThank you for your approval and suggestions on related work.\nIndeed, we noticed that ASkotch [5] was updated after our submission.\nWe found many new results, and we will update our review of this paper.\n## Theorem 2.1\nIndeed, many results similar to Theorem 2.1 was mentioned in the literature.\nWe did not list Theorem 2.1 as our contribution.\nThis theorem aims to highlight the dual formulation of the kernel methods,\nwhich is the key problem of this paper.\n# Ref.\n[1] Hsia et al. Preconditioned conjugate gradient methods in truncated Newton frameworks for large-scale linear classification, ACML2018.\n\n[2] Dai et al. Scalable Kernel Methods via Doubly Stochastic Gradient, Neurips2014, latest report: arxiv.org/pdf/1407.5599.\n\n[3] Le et al. Fastfood: Approximate Kernel Expansions in Loglinear Time, ICML2013.\n\n[4] Han et al. Fast Neural Kernel Embeddings for General Activations, Neurips2022.\n\n[5] Rathore et al. Have ASkotch: A Neat Solution for Large-scale Kernel Ridge Regression, arxiv.", "claims_and_evidence": "Claims of better or even performance under smaller memory budget are convincing.", "methods_and_evaluation_criteria": "Yes, the evaluation criteria makes sense.", "theoretical_claims": "Theoretical claims from section 2 are correct.", "experimental_designs_or_analyses": "The experimental designs are correct. However no code is provided, and given the practical nature of the paper, I would have liked to audit the code.", "relation_to_broader_scientific_literature": "- The Eigenpro series of papers could be better discussed in the introduction.\n\n- The most recent addition to large scale KRR might be the arxiv paper [1], from July 2024, with a revised version from February 2025. However, given that it is recent and likely not through with the review process itself, I understand that it might be touchy to ask for a comparison. Still, it deserves more than the current hard-to-notice citation.\n\nOtherwise the related work section is well organized. The idea of using trust regions on the dual problem is novel and promising, especially the part about the quadratic extension at (8).\n\nThe contribution from theorem 2.1 is minimal. The dualization of the regularized empirical risk minimization problem in a RKHS has been extensively studied over two decades. Similarly, the contribution of proposition 2.2 is again not original, it is a straightforward extension of the case of the infimal convolution of two loss functions.\n\n[1] Have ASkotch: A Neat Solution for Large-scale Kernel Ridge Regression, Rathore et al.", "essential_references_not_discussed": "None.", "other_strengths_and_weaknesses": "+: The paper is very relevant to ICML and the practical improvement over existing methods is valuable.\n\n-: The writing is often not precise, see comments or suggestions box.\n\n-: Inexact Joker relies on random fourier features, thus is limited to shift-invariant kernel (as proposed at least).", "other_comments_or_suggestions": "- [Line 94] A Mercer kernel is typically a kernel that satisfies the assumptions needed to apply Mercer's theorem: continuous on a compact domain. Any positive definite kernel uniquely induces a RKHS.\n\n- [Line 105] Typo $\\phi(x_i)_{\\mathcal{H}}$\n\n- [Line 107] where $\\phi(x)$ is the linear map -> a linear map, as there can be several.\n\n- [Line 132] the closeness and the convexity\n\n- [Line 171] \"It implies that if the constraints are properly handled,\nkernel Huber regression can be as efficient as KRR\" -> what this implies is that if the constraints are satisfied by the dual solution to the KRR, then both solutions coincide. There exist cases of contaminated data where the Huber loss estimator outperforms the KRR precisely because the dual coefficients must pertain to a smaller ball and cannot be too influenced by outliers.\n\n- [Line 189] \"in some rare cases\" -> either $f$ is assumed twice differentiable, or it is not. Pick one but do not claim full generality while proposing ad-hoc modifications.\n\n- [Line 230] \"In our implementation, $T_{TR} \\leq 50$\" -> ok but that does not mean that $n$ does not appear in the complexity. For example if you use blocks of size $512$ and have $512k$ samples you would approximately need to go through 1k times to make a pass on all the dual variables. And this would not even guarantee convergence. Or am I wrong here ?\n\n- [Line 254] \"we obtain $\\theta = \\sum_{i=1}^n \\alpha_i \\psi(x_i)$\" -> usually when using random features, one big advantage is that the parameter space is reduced, so that $\\theta$ can be directly searched for on the feature space. Could you comment on the difference with your approach here ?\n\n- [Line 272] The results from [Rahimi and Recht, 2007] are known to be suboptimal and certainly do not represent a theoretical guide for setting $M$. See e.g. [Error Bounds for Learning with Vector-Valued Random Features, Lanthaler and Nelsen, Neurips 2023].", "questions_for_authors": "My main criticism about the paper is that I find that it lacks rigor. Improving on the points raised in \"other comments or suggestions\" is critical to me.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "b0SDiUpjxY", "reviewer": "Reviewer_MsAQ", "summary": "This paper explores a joint optimization framework for diverse kernel models, including KRR, logistic regression, and support vector machines.\nThe authors employed a dual block coordinate descent method with trust region (DBCD-TR) and kernel approximation with randomized features to solve the proposed model, which makes the algorithms have low memory costs and high efficiency in large-scale learning.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Q1. bias term.\nIt is for simplicity and keeping consistency with the recent literature of kernel methods, where they also do not consider it.\nA simple way to include the bias is to append a constant after $\\varphi(x)$.\nHowever, adding this term generally has no impact on the final performance,\nso our paper as well as the related work tends to disregard it.\nThis is good feedback, and we will briefly explain it in our paper.\n\n## Q2. benefit of trust region.\nThe main reason is that projected gradient descent (PGD) converges slowly.\nThe best rate of PGD is linear, specifically, $O(\\kappa\\log(1/\\epsilon))$, where $\\kappa$ is the condition number.\nBut due to the kernel matrix, eq.(6) usually has large $\\kappa$ and the linear rate may be insufficient.\nSo many latest kernel methods, such as [1] use second-order information for acceleration.\nThe proposed trust region (TR) has the same purpose and can reach a superlinear rate [2],\nsignificantly faster than PGD.\n\n## Q3. approximation of eq.(7).\nThis is a misunderstanding of TR method.\nIt aims to find multiple steps toward the optimal point,\nand each point is restricted in a small region (i.e., TR) wherethe  approximation error of eq.(7) is tiny.\nThe procedure of TR method guarantees that the steps always decrease the objective function and finally reach the optimality (especially, with superlinear convergence rate [2]).\nIf the approximation error of eq.(7) is large and yields a suboptimal step, TR's procedure will reject such a step automatically and decrease the region size.\nIt may be clearer to see Algorithm 3 in Appendix B.\n\n## Q4. Nonconvex loss.\nOur paper focuses on the large-scale kernel method with the convex losses, but not the nonconvex ones.\nWe have stated that \"Joker focuses on convex problem...\" (line 143),\nand nonconvexity is not within the scope of our paper.\nDespite this, Joker still significantly develops kernel models with diverse loss functions,\nand improves the speed and memory footprint compared to the existing kernel methods.\n\nKernel methods with nonconvex losses should be a future work.\nThe approaches may be significantly different from the proposed Joker.\nThe feasible approaches may include Difference of Convex functions (DC) programming, convex relaxation, and approximation\n(e.g., piecewise linear approximation method [1] mentioned by reviewer Zf7v).\n\n## Q5. Convergence.\nWe decided not to include convergence for the following reasons:\n\n1. Our work is more practical, as stated by reviewer 1dCw, it has \"the practical nature\".\nThe proposed method has remarkable contributions:\nlargely reducing the memory footprint of large-scale kernel methods, and unifying the optimization scheme of different kernel models.\nWe believe the two advances are significant enough in practice.\n\n2. We feel that the rigorous proof of the convergence is non-trivial and deserves in-depth study in future work.\nThe proposed DBCD-TR algorithm is a relatively complicated algorithm.\nBlock coordinate descent and the truncated CG-Steihaug method (Algorithm 1) may become two major challenges in investigating the convergence rate.\nIndeed, a simple linear convergence rate $O(n\\log(1/\\epsilon)/|B|)$ can be obtained easily using the Polyak-Lojasiewicz condition.\nHowever, it is vacuous since the first-order block coordinate descent methods also have this rate.\nSo the linear rate cannot highlight the improvement of DBCD-TR.\nConsidering the usage of second-order information, it is reasonable to guess that DBCD-TR can have a superlinear convergence rate, as hinted in Section 6 of [2].\n\nNonetheless, we still can give an outline of the convergence proof:\n- Sufficient descent. We first prove that truncated CG-Steihaug satisfies the sufficient descent condition in [3].\n- Global Convergence to a stationary point. The limit point of the proposed trust region is a stationary point.\n- Local Superlinear convergence: The trust region process (Algorithm 2) converges with a superlinear rate same as the projected Newton.\n- Giving the iteration complexity using the above results, and investigate the influence of the Hessian in different cases.\n\nThe above illustrates the convergence of one block eq.(6).\nFor the outer loop (the block coordinate),\nwe may establish the iteration complexity result following the analysis in [4] and combining the decrease made by the trust region.\n\n# Ref.\n[1] Teo et al. Bundle Methods for Regularized Risk Minimization. JMLR, 2009.\n\n[2] Nutini et al. Let's Make Block Coordinate Descent Converge Faster: Faster Greedy Rules, Message-Passing, Active-Set Complexity, and Superlinear Convergence. JMLR, 2022.\n\n[3] Baraldi et al. Efficient proximal subproblem solvers for a nonsmooth trust-region method. Computational Optimization and Applications, 2025.\n\n[3] Li et al. On Faster Convergence of Cyclic Block Coordinate Descent-type Methods for Strongly Convex Minimization. JMLR, 2018.", "claims_and_evidence": "The proposed approach shows a good performance on some tasks.", "methods_and_evaluation_criteria": "The evaluation metric used in this paper is reasonable.", "theoretical_claims": "I checked the proofs of  some theorems, and they sound correct.", "experimental_designs_or_analyses": "Experimental results demonstrate a good performance compared to baselines.", "relation_to_broader_scientific_literature": "The proposed techniques can be employed to deal with various kernel models beyond KRR, thereby giving a more comprehensive solution for kernel learning and related applications.", "essential_references_not_discussed": "The paper contains sufficient discussed references", "other_strengths_and_weaknesses": "This paper explores a joint optimization framework for diverse kernel models, including KRR, logistic regression, and support vector machines. For weaknesses, please refer to the problem below.", "other_comments_or_suggestions": "No suggestions.", "questions_for_authors": "1)Classical kernel-based methods usually contain biases. Why are the biases not considered in the proposed model?\n\n2)For Eq.(6), one may employ projected gradient methods to solve them. What is the benefit of the rust-region method adopted in this article?\n\n3)The authors employ Taylor expansions in Eq.(7). This approximation will introduce errors such that the solution deviates from the optimal solution.\n\n4)When nonconvex loss functions are employed, there exists the dual gap between the primal problem and the dual problem. How to deal with this problem?\n\n5)It would be much better to prove the convergence of the algorithm under proper assumptions.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Qj6J1GdLrm", "reviewer": "Reviewer_Zf7v", "summary": "The authors propose a general optimization scheme for kernel machines. Similar to Teo et al. (2009), conjugate loss functions allow for unified representations that are solved with block coordinate descent in dual space (Sorensen, 1982). The authors report on empirical results showing that their optimizer is more efficient than its competitors.\n\n[Slightly updated score after rebuttal]", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Questions\n## Q1. Convergence\nThe linear convergence rate of the proposed DBCD-TR can be proven using Polyak-Lojasiewicz condition.\nHowever, it is a loose bound and does not highlight the advantage of DBCD-TR.\nThe tighter bound (a reasonable guess is a superlinear rate) is a challenging work, considering that DBCD-TR itself is a quite complicated algorithm.\nTherefore, we tend not to present the convergence result in this paper, and focus on the practical design of the algorithm.\nWe have explained this in response to Q5 of the reviewer MsAQ.\n\n## Q2. Comparison on OCAS [1]\nWe argue that OCAS and the BMRM are unsuitable to be benchmarks in our experiment for the following reasons.\n\n1. They are designed for linear models but not kernel methods.\nThis distinction is important because kernel methods are usually associated with an ill-conditioned problem,\nwhich means that *a fast solver for a linear model could be slow in the kernel regime*.\nThe recent kernel methods, including EigenPro series, Falkon, LogFalkon aim to overcome the ill-conditioning via preconditioning or Newton methods.\nOur solution, DBCD-TR, uses **trust region** method to incorporate truncated Newton step and optimizes **multiple variables (a block)** at a time to leverage the merit of parallel computing.\n\n2. They are outdated since LIBLINEAR is faster than them, as shown in the experiments of [3].\nTheoretically, LIBLINEAR is based on dual coordinate descent [2] and has a linear convergence rate $O(\\log(1/\\epsilon))$ superior to $O(1/\\epsilon)$ of OCAS.\nSo if we must make a comparison, we should choose LIBLINEAR but not [2] and [3].\n\n3. The proposed DBCD-TR is much faster than LIBLINEAR.\nAs mentioned before, OCAS, BMRM and LIBLINEAR are for linear models.\nBut it does not mean that they cannot handle the kernel models.\nAs mentioned by [3], we can first obtain the random Fourier feature (RFF) of data and apply the linear models,\nthus obtaining an inexact kernel machine.\nUsing this approach, we run LIBLINEAR (GPU implementation for fairness) on HIGGS ($M=10^5$, same as Joker-SVM in our experiment) during rebuttal.\nHowever, it hardly converges and only obtains an accuracy of 64.7% after 3 days. \nFor OCAS and BMRM, they will cost more time.\nTherefore, we think they are unsuitable to be benchmarks in our experiment.\n\n## Q3. Relation to BMRM [3]\nBMRM is also for linear SVM as OCAS and LIBLINEAR.\nBMRM minimizes the piecewise lower bound (PLB) of the primal objective via the Fenchel conjugate.\nThen it minimizes PLB by optimizing a series of dual problems.\nIn contrast, Joker does not approximate the primal objective function,\nbut optimizes the dual objective directly.\nWe have presented more discussion in Q2.\n\n# Other Comments\n## 1. Typos in Appendix\nWe have fixed the typos you mentioned.\nThank you for pointing them out.\n\n## 2. Notation of infimal convolution\nYes, we note the infimal convolution as \"$\\square$\" following the notation of most textbooks.\n\n## 3. Interpretation of Table 4\nThank you for your comments.\nOur words of \"substantial\" may be imprecise.\nWhat we want to state is that the time gap between other methods like EigenPro and Falkon is large,\nand Joker alleviate this gap significantly.\nWe realize that the analysis combined with data may be convincing:\nEigenPro3 and ThunderSVM use at least 10x training time than Falkon on MSD,\nand Joker reduces it to 5~6x time.\nEspecially on HIGGS, EigenPro3 needs 36x time (18 hours) of Falkon (0.5 hour),\nand Joker reduces the gap to 2x time (1 hour). \n\n# Ref.\n[1] Franc et al. Optimized cutting plane algorithm for support vector machines. ICML 2008.\n\n[2] Hsieh et al. A dual coordinate descent method for large-scale linear SVM. ICML 2008.\n\n[3] Teo et al. Bundle Methods for Regularized Risk Minimization. JMLR, 2009.", "claims_and_evidence": "The authors claim generality of their approach due to the use of conjugate loss functions. The idea is not novel and always lead to unified representations of quadratic objectives (e.g., regularized empirical risks). The empirical results show that the present approach is actually saving a great deal of memory when compared with other optimizers. I am not up to date in this domain anymore but back than the fastest solver was OCAS and a similar versatility were contained in Teo et al. (see below for references). Unfortunately, both are not cited or compared in the paper.\n\nI checked the proof of Theorem 2.1. The proof seems to be correct, but may suffer from some smaller inaccuracies. \nFrom my understanding, the last sum in Equation A.2. (line 616/617) should be :\n$$ \\sum_{i=1}^{n}\\alpha_i(u_i-\\theta^T\\varphi(x_i))$$\nchanging the sign of the expression to be coherent with Equations A.3 and A.4, and replacing the $\\omega$ with  $\\theta$.", "methods_and_evaluation_criteria": "The approach is evaluated against four baselines (two are missing in my view, see additional literature below) on five large data sets. The proposed method appears to beat the baselines in terms of memory footprint and often also predictive accuracies.", "theoretical_claims": "The authors show hardly any theoretical claims. It would have been nice to learn more about convergence of the inexact variant with randomization of features.", "experimental_designs_or_analyses": "See above.", "relation_to_broader_scientific_literature": "The paper seems a bit outdated although faster computation is always appreciated, even for kernel machines that are kind of outruled by networks at the moment. There has been quite a great deal of optimization approaches for kernel machines end of the nineties until perhaps 2010. I am citing two papers that came to my mind but I have to admit that I don't remember much from what I know once...", "essential_references_not_discussed": "V. Franc and S. Sonnenburg. Optimized cutting plane algorithm for support vector machines. In A. McCallum and S. Roweis, editors, Proceedings of the International Conference on Machine Learning, pages 320–327. Omnipress, 2008. (OCAS used to be the fastest SVM solver back then)\n\nC. H. Teo, S.V.N. Vishwanathan, A. Smola, Q. Le: Bundle Methods for Regularized Risk Minimization. Journal of Machine Learning Research 11 (2010) 311-365, 2009.", "other_strengths_and_weaknesses": "The proposed method is actually very memory efficient and performs well compared with the baselines. There is a lack of theory (eg convergence) for such a paper in my view.", "other_comments_or_suggestions": "In my PDF, the infimal convolution is denoted by a square that seems identical to the end of a proof environment. Is that they symbol you wanted us to see? The authors claim that \"Falkon based methods are the fastest, but the gap between Joker and them is not substantial\" (line 370-373). This seems to be a rather subjective interpretation of the results in Table 4, where Falkon performs comparatively at only about 1/6 of the time requirement on the MSD dataset and LogFalkon achieving the best results at about half the time of the best Joker based approach. There seems to be a typo in Equation A.5, where it should be  $\\alpha_i\\alpha_j$ instead of $\\alpha_i\\alpha_k$.", "questions_for_authors": "What about convergence proofs or bounds? How does performance and efficiency compare to OCAS? How does the approach relate to Teo et al?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "BstCstPlUX", "reviewer": "Reviewer_aGkm", "summary": "This paper proposes a novel algorithm for extremely large-scale kernel machines. Its main contribution lies in Theorem 2.1, which reformulates the objective function in the dual problem of kernel methods into a form based on decoupled conjugate functions. This ensures convexity and strong duality, making it possible to use the block coordinate descent method for optimization. In each subproblem, the conjugate function is further approximated by performing a Taylor expansion and optimizing the approximation. Experimental results show that this method outperforms some previous algorithms for very large-scale kernel machines.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weak. 1: \"Exact\" and \"Inexact\"\nIn the context of this paper, the \"inexact\" and \"exact\" are not about the approximation of eq.(7),\nwhich has been stated in Section 1.2.\n\"Exact\" refers to solving the problem eq.(1) without approximating the kernel function $K(\\cdot,\\cdot)$, and these methods usually involves $n$-dimensional variable $\\alpha$.\n\"Inexact\" approximates the kernel function with a map $\\psi(\\cdot)$ (e.g., Nystrom and RFF) such that $K(x,x')\\approx \\psi(x)^\\top\\psi(x')$ and reduces computation burden.\nThese terminologies are consistent with [1].\n\nBesides, the trust region method in solving eq.(6) is not similar to Nystrom.\nAs a well-studied optimization technique,\nthe trust region method has a theoretical guarantee of convergence to the optimal point.\nWe have explained it in response to Q3 of reviewer MsAQ.\n\n## Weak. 2: impact of other hyperparameters\nThank you for the suggestion.\nWe have listed the optimal hyperparameters in Table A.3.\nThe block size is emphatically discussed because it directly influences the memory footprint and the convergence,\nwhile the max region size $\\Delta_{max}$ has a lower impact. \nThis is because the trust region procedure will adaptively tune the region size to ensure a sufficient decrease of the objective function.\nSo the final performance is usually insensitive to $\\Delta_{max}$ as long as it is not extremely tiny or huge. \nEmpirically, $\\Delta_{max}\\in[4,64]$ is moderate.\n\n\n# Other comments\nThank you for your careful review.\nWe have corrected the typos you mentioned.\n\n# Ref.\n[1] Rahimi et al. Random Features for Large-Scale Kernel Machines, Neurips2007.", "claims_and_evidence": "Strength:\n\n1. The formulation provided in Theorem 2.1 is ingenious. The subsequent ideas, including block decomposition optimization, approximate optimization, and the use of random features for approximation, are sound.\n\nWeakness:\n\n1. I believe the author's definition of \"exact\" is somewhat biased. The core of the proposed method is to focus only on one subset during each optimization step (as in Equation (6)). For this subproblem, approximate optimization is required (as in Equation (7)). Essentially, this is similar to methods like Nystrom, which select a sub-dataset and are inherently inexact. If the authors wish to claim that their method is \"exact\", they need to theoretically guarantee that the iterative optimization of the subproblem with box constraints can reach the optimal solution of Problem (4) or provide some error bounds. The authors could refer to related proofs in block coordinate descent methods to explain this.\n\n2. In the experiments, the authors discussed the impact of the choice of $|\\mathcal{B}|$ on the results, but I did not notice some discussion regarding the setting of other parameters, such as the maximum region size $\\Delta$.", "methods_and_evaluation_criteria": "The compared methods and criteria are reasonable.", "theoretical_claims": "There are some typos in the proof of Theorem 2.1. Please refer to \"Other Comments or Suggestions\".", "experimental_designs_or_analyses": "Good.", "relation_to_broader_scientific_literature": "None.", "essential_references_not_discussed": "None.", "other_strengths_and_weaknesses": "See Claims And Evidence.", "other_comments_or_suggestions": "1. In Equation (1), it should be written as $\\langle \\theta, \\phi(x_i) \\rangle_{\\mathcal{H}}$.\n\n2. In Equation (A.2), $w$ should be replaced with $\\theta$. Additionally, to align with Equations (A.3) and (A.4), it would be better to use $\\alpha_i(u_i - \\theta^\\top \\phi(x_i))$.\n\n3. Due to the notation $\\xi_y(u)$ in the main context, it is preferable to use something like $(f \\square g)(u) \\coloneqq \\inf_p f(p) + g(u - p)$.", "questions_for_authors": "None.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper addresses the scalability of kernel methods by proposing a unified optimization framework for various kernel models—including KRR, logistic regression, and support vector machines.  Joker uses a dual block coordinate descent method with trust region and randomized feature-based kernel approximation. It significantly reduces memory usage (up to 90%) while maintaining or surpassing the training speed and performance of state-of-the-art methods.

The paper is well written, and the empirical performance is impressive. The convergence proof is lacking, but the authors gave a reasonable outline of the proof in the rebuttal. Overall, the paper is a good addition to the conference.
Decision: Accept (poster)