Fast Kernel Methods for Generic Lipschitz Losses via $p$-Sparsified Sketches
Authors: Tamim El Ahmad
Pierre Laforgue
Florence d'Alché-Buc
Venue: Accepted by TMLR ()

Abstract: Kernel methods are learning algorithms that enjoy solid theoretical foundations while suffering from important computational limitations. Sketching, which consists in looking for solutions among a subspace of reduced dimension, is a well-studied approach to alleviate these computational burdens. However, statistically-accurate sketches, such as the Gaussian one, usually contain few null entries, such that their application to kernel methods and their non-sparse Gram matrices remains slow in practice. In this paper, we show that sparsified Gaussian (and Rademacher) sketches still produce theoretically-valid approximations while allowing for important time and space savings thanks to an efficient \emph{decomposition trick}. To support our method, we derive excess risk bounds for both single and multiple output kernel problems, with generic Lipschitz losses, hereby providing new guarantees for a wide range of applications, from robust regression to multiple quantile regression. Our theoretical results are complemented with experiments showing the empirical superiority of our approach over state-of-the-art sketching methods.

Published in Transactions on Machine Learning Research (09/2023)
Fast Kernel Methods for Generic Lipschitz Losses via p-Sparsified Sketches
Tamim El Ahmad LTCI, Télécom Paris, Paris, Fronce
kamim.edabmad@@elecommmoro
Pierre Laforgue
picrre.laforuueehmmii..
Department of Computer Science, University of Milan, Itoly
Florence d'Alché-Buc
LTCI, Télécom Paris, IP Paris, France
Reviewed on OpenReview: https: openreview. net/ net/foree forum? id- ddrggggRRTTOu
Abstract
Kernel methods are learning algorithms that enjoy solid theoretical foundations while suffering from important computational limitations. Sketching, which consists in looking for solutions among a subspace of reduced dimension, is a well-studied approach to alleviate these computational burdens. However, statistically-
sketches, such as the Gaussian one, usually contain few null entries, such that their application kernel methods and their non-sparse Gram matrices remains slow in practice. In this paper, we show that sparsified Gaussian (and Rademacher) sketches still produce theoretically-valie approximations while allowing for imporlant time and space savings thanks to elficient decomposilion trick. To support our method, we derive excess risk bounds for both single and multiple output kernel problems, with generic Lipschitz losses, hereby providing new guarantees for a wide range of applications, from robust regression to multiple quantile regression. Our theoretical results are complemented with experiments showing the empirical superiority of our approach over state-of-the-ate sketching methods.
1 Introduction
Kernel methods hold a privileged position in machine learning, as they allow to tackle large variety of learning tasks in a unique and generic framework, that of Reproducing Kernel Hilbert Spaces (RKHSs), while enjoying solid theoretical foundations (Steinwart & Christmann, 2008b; Scholkopf & Smola, 2018). From scluer-valued to multiple output regression (Micchelli & Pontil, 2005; Carmeli et al., 2006; 2010), these approaches play a central role in nonparametric learning, showing a great flexibility. However, when implemented naively, kernel methods raise major issues in terms of time and memory complexity, and are often thought of as limited to "fat data", i.., datasets of reduced size but with large number of input features. One way to scale up kernel methods are the Random Fourier Features (Rahimi & Recht, 2007; Rudi & Rosasco, 2017; Sriperumbudur & Szabó, 2015; Li et al., 2021), but they mainly apply to shift-invariant kernels. Another popular approach is to use sketching methods, first exemplified with Nyström approximations (Williams & Seeger, 2001; Drineas et al., 2005; Bach, 2013; Rudi et al., 2015). Indeed, sketching has recently gained lot of interest in the kernel community due to its wide applicability (Yang al., 2017; Lacotte et al., 2019; Kpotufe& Sriperumbudur, 2020; Lacotte Pilanci, 2022; Gazagnadou et al., 2022) and its spectacular successes when combined to preconditioners and GPUs (Meanti et al., 2020).

---

Published in Transactions on Machine Learning Research (09/2023)
Sketching as random projection method (Mahoney et al., 2011; Woodruff, 2014) is rooted in the Johnson- Lindenstrauss lemma (Johnson & Lindenstrauss, 1984), and consists in working in reduced dimension subspaces while benefiting from theoretical guarantees. Learning with sketched kernels has mostly been studied in the case of scalar-valued regression, in particular in the emblematic case of Kernel Ridge Regression (Alaoui & Mahoney, 2015; Avron el al., 2017; Yang el al., 2017; Chen & Yang, 2021a). For several identified sketching types (e.g., Gaussian, Randomized Orthogonal Systems, adaptive sub-sampling), the resulting estimators come with theoretical guarantees under the form of the minimax optimality of the empirical approximation error. However, an important blind spot of the above works is their limitation to the square loss. Few papers go beyond Ridge Regression, and usually exclusively with sub-sampling schemes (Zhang et al., 2012; Li et al., 2016; Della Vecchia et al., 2021). In this work, we derive excess risk bounds for sketched kernel machines with generic Lipschitz-continuours losses, under standard assumption on the sketch matrix, solving an open problem from Yang et al. (2017). Doing so, we provide theoretical guarantees for a wide range of applications, from robust regression, based either on the Huber loss (Huber, 1964) or c-inscnsitive losses (Steinwart& & Christmann, 2008a), to quantile regression, tackled through the minimization of the pinball loss (Koenker, 2005). Further, we address this question in the general context of single and multiple output regression. Learning vector-valued functions using matrix-valued kernels (Micchelli & Pontil, 2005) has been primarily motivated by multi-task learning. Although equivalent in functional terms to scalar-valued kernels on pairs of input and tasks (Hein & Bousquet, 2004, Proposition 5), matrix-valued kernels (Álvarez et al., 2012) provide a way to define a larger variety of statistical learning problems by distinguishing the role of the inputs from that of the tasks. The computational and memory burden is naturally heavier in multi-task/ taaltiitaskkmulttouutppii multi-output regression, as the dimension of the output space plays an inevitable role, making approximation methods for matrix-valued kernel machines crucial issue. To our knowledge, this work is the first to address this problem under the angle of sketching. It is however worth mentioning Baldassarre et al. (2012), who explored spectral filtering approaches for multiple output regression, and the generalization of Random Fourier Features to operator-valued kernels by Brault et al. (2016).
An important challenge when sketching kernel machines is that the sketched items, e.g., the Gram matrix, are usually dense. Plain sketching matrices, such as the Gaussian one, then induce significantly more calculations than sub-sampling methods, which can be computed by applying mask over the Gram matrix. Sparse sketching matrices (Clarkson & Woodruff, 2017; Nelson & Nguyên, 2013; Cohen, 2016; Derezinski et al., 2021) constituue an important line of research to reduce complexity while keeping good statistical properties when applied to sparse matrices (e.g., matrices induced by graphs), which is not the case of a Gram matrix. Motivated by these considerations, we analyze a family of sketches, unified under the name of p-sparsified sketches, that achieve interesting trade-offs between statistical accuracy (Gaussian sketches can be recovered as a particular case of p-sparsified sketches) and computational efficiency. The p-sparsified sketches are also memory-efficient, as they do not require computing and storing the full Gram matrix upfront. Besides theoretical analysis, we provide extensive experiments showing the superiority of p-sparsified sketches over SOTA approaches such as accumulation sketches (Chen & Yang, 2021a).
Contributions. Our goal is to provide framework to speed-up both scalar and matrix-valuch kernel methods which is as general as possible while maintaining good theoretical guarantees. For that purpose, we present three contributions, which may be of independent interest.
We derive excess risk bounds for sketched kernel machines with generic Lipschitz continuots losses, both in the scalar and multiple output cases. We hereby solve an open problem from Yang et al. (2017), and provide first analysis to the sketching of vector-valued kernel methods.
We show that sparsified Gaussian and Rademacher sketches provide valid approximations when applied to kernel methods. They maintain theoretical guarantees while inducing important space and computation savings, as opposed to plain sketches.
We discuss how to learn these new skotched kernel machines, by mcans of an approximated feature map. We finally present experiments using Lipschitz-contimious losses, such as robust and quantile regression, on both synthetic and real-world datasets, supporting the relevance of our approach.
2

---

Published in Transactions on Machine Learning Research (09/2023)
Notation. For any matrix A E Rmxp, At is its pseudo-inverse, ||Allop its operator norm, A, R' its i-th row, and Ag E IR" its j-th column. The identity matrix of dimension d is la. For couple of random variables (X, Y) E X x y with distribution P. Px is the marginal distribution of X. For f : X - V, we use E [f] = Epx [f(x),, B = Ep (f()), and Enleyl = - ((f(xs), ye) for any function
y: x IR.
2 Sketching Kernels Machines with Lipschitz-Continuunn Losses
In this section, we derive excess risk bounds for sketched kernel machines with generic Lipschitz losses, for both scalar and multiple output regression.
2.1 Scalar Kernel Machines
We consider a general regression framework, from an input space X' to some scalar output space ) C R. Given a loss function & y y + IR such that = + ((z,y) y) is proper, lower semi-continuous and convex for every y, our goal is to estimate f = arg inf SER E(xy)~ [e (f(X), Y)], where H C )d a hypothosis set, and is a joint distribution over X' x V. Since P usually unknown, we assumne that we have access to a training dataset {(1, y))==1 composed of i..d. realisations drawn from P. We recall the definitions of a scalar-valued kernel and its RKHS (Aronszajn, 1950).
Definition (Scalar-valued kernel). A scalar-valued kernel a symmetric function : X - such that for all n E N, and any (x:))=1 X", (a:) E ER, we have Eij=t 0g k 24) aj 2 0.
Theorem 1 (RKHS). Let k be a kernel on X. Then, there exists a unique Hilbert space of functions He c. Rx such that k (,,2) Hs for all X, and such that we have h (x) (h, k (1,2)) for any (h, H4 X.
A kernel machine computes proxy for fe by solving
An 2 |1 ll
Σ (((a)))))) +
min JEHL
where Ay > 0 is a regularization parameter. By the representer theorem (Kimeldorf & Wahba, 1971; Schölkopf et al., 2001), the solution to Problem (1) is given by Jn = Σi= & k(;,21),, with & E R the solution to
n
on a Kaa
min RER
Σ
{([a)).))))
(2)
where E RaxR is the kernel Gram matrix such that Kiy= = x)))
Definition 2 (Regularized Kernel-based Sketched Estimator). Given a matria S E Raxn, with 8 & n, sketching consists in imposing the substitution a = ST in the empirical risk minimization problem stated in Equation (2). We then obtain an optimisation problem of reduced size on y, that yields the sketched estimator f. Σl (ST (TT)): k(), where 7 E Rs is solution to
min 1 yER
n Σ
A TSKST 2
(3)
In practice, usually obtains the matrix S by sampling it from random distribution. The literature is rich in examples of distributions that can be used to generate the sketching matrix S. For instance, the sub-sampling matrices, where each line of S is sampled from In, have been widely studied in the context. of kernel methods. They are computationally efficient from both time and space perspectives, and yield the so-called Nyström approach (Williams & Seeger, 2001; Rudi al., 2015). More complex distributions, such as Randomized Orthogonal System (ROS) sketching or Gaussian sketch matrices, have also been considered (Yang et. al., 2017). In this work, we first give general theoretical analysis of regularized kernel-based
3

---

Published in Transactions on Machine Learning Research (09/2023)
sketched estimators for any K-satisffable sketch matrix (Definition 3). Then, we introduce the p-sparsified sketches and prove their K-satisfiablity, as well as their relevance for kernel methods in terms of stat and computational trade-off.
Works about sketched kernel machines usually assess the performance of fe by upper bounding its squared 1,2 (Pxx) error, i.e., (1/n) Σ ((.:))) (ax))), where SH. is the minimizer of the true risk over He, supposed to be attained (Yang et al., 2017, Equation 2), or through its (relative) recovery error II.e - fu l /llfull see Lacotte& & Pilanci (2022, Theorem 3). In contrast, we focus on the excess risk of fs, the original quantity of interest. As revealed by the proof of Theorem 2, the approximation error of the excess risk can be controlled in terms of the (P((N) error, and we actually recover the results from Yang et al. (2017) when we particularize to the square loss with bounded outputs (second bound in Theorem 2). Furthermore, studying the excess risk allows to better position the performances of fs among the known off- off-the-sheil kernel- based estimators available for the targeted problem. To achieve this study, we rely the key notion of K K-satisfiability for a sketch matrix (Yang et al., 2017; Liu et al., 2019; Chen & Yang, 2021a).
Let K/n = U T be the cigendecomposstion of the Gram matrix, where D= diag (Hi,...,,p)) stores the eigenvalues of K/n in decreasing order. Let 6% be the critical radius of K/n, ie., the lowest value such that (6,,) ( E"= min(6 (i))1/ I 82. The existence and uniqueness of 82 is guaranteed for any RKHS associated with a positive definite kernel (Bartlett et al., 2006; Yang et al., 2017). Note that 52 is similar to the parameter E2 used in Yang et al. (2012) to analyze Nystrõm approximation for kernel methods. We define the statistical dimension of K as dn = min {i E (1, (I.......: Pj S { with dn = n if no such index j exists.
Definition 3 (K-satisfiabiilit, Yang et al. 2017). Let c > O be independent of n, U1 E Raxd, and U2 E (=x(n=d,) be the left and right blocks of the matriz U previously defined, and D2 = diag (b).)......,,)). A matric S is said to be K K-satisfiable for e if we have
(SU:)] T SU - Idn 1.pp op 1/2, and | op
S c5p
(4)
Roughly speaking, a matrix is K-satisfiable if it defines an isometry on the largest eigenvectors of K, and has a small operator norm on the smallest eigenvectors. For random sketching matrices, it is common to show K-satisfiability with high probability under some condition on the sketch size s, see c.g., Yang et al. (2017, Lemma 5) for Gaussian sketches, Chen & Yang (2021a, Theorem 8) for Accumulation sketches. In Section 3, we show similar results for p-sparsified sketches.
To derive our excess risk bounds, we place ourselves in the framework of Li et al. (2021), see Sections 2.1 and 3 therein. Namely, we assume that the true risk is minimized over H4 at JHH = arg min seH, E [ (X), Y)]- The existence of fu, is standard in the literature (Caponnetto & De Vito, 2007: Rudi & Rosasco, 2017; Yang et al., 2017), and implies that fH,, has bounded norm, see e.g., Rudi & Rosasco (2017, Remark 2). Similarly to Li et al. (2021), we also assume that estimators returned by Empirical Risk Minimization have bounded norms. Hence, all estimators considered in the present paper belong to some ball of finite radius R. However, we highlight that our results do not require prior knowledge on R. and hold uniformly for all finite R. As a consequence, we consider without loss of generality as hypothesis set the unit ball B (He) in HA, up to an posteriori rescaling of the bounds by R to recover the general case.
Assumption 1. The true risk is minimized at fHu- Assumption 2. The hypothesis set considered is B (Hz). Assumption 3. For all y E y, 2 + f(z.g) is L-Lippchitz. Assumption 4. For all x, x' E X, we have k(I,2') S K.
Assumption 5. The sketch S is K --satisfiable with constant c > 0.
Note that we discuss some directions to relax Assumption 2 in Appendix B. Many loss functions satisfy Assumption 3, such as the hinge loss (L= = 1), used in SVMs (Cortes & Vapnik, 1995), the e-insensitive & (Drucker et al., 1997), the k-Huber loss, known for robust regression (Huber, 1964), the pinball loss, used in quantile regression (Steinwart & Christmann, 2011), or the square loss with bounded outputs. Assumption 4 is standard (e.g., for the Gaussian kernel). Under Assumptions 1 to 5 we have the following result.
4

---

Published in Transactions on Machine Learning Research (09/2023)
Theorem 2. Let fo in Definition 2, suppose that Assumptions 1 to 5 hold, and let C = + v6c, with the constant from Assumplion 5. Then, for any å E (0,1) with probabilily at least 1 - 8 we have
An + 88
8 log (4/6)
B [l s L [e + LC + + +
+2V
(5)
2
n
Furthermore, ( v) = (2 - -)²² 12 and y c 10, with probability at least 1 - we have
I [tg,]] s [m + (02 + 2) An + C
(c²+2)
K + VR +
8 log (4/5)
Vn + 2/2
(6)
Proof sketch. The proof relies on the decomposition of the excess risk into two genoralization error terms and an approximation error term, i.e.,
E.l B El, - - B, +.,1,,11 -..07,-En..... In, 1 E Enlesne sn, BIL sm
(7)
The two generalization errors (of f. and (HH)) can be bounded using Bartlett & Mendelson (2003, Theorem 8) together with Assumptions 1 to 4. For the last term, we can use Jensen's inequality and the Lipschitz contimuity of the loss to upper bound this approximation crror by the square root of the sum of the square residuals of the Kernel Ridge Regression with largets the (x.).. The latter can in turn be upper bounded using Assumptions and 5 and Lemma 2 from Yang et al. (2017). When considering the square loss, Jensen's 's inequality is not necessary anymore, leading to the improved second term in the right-hand side of the last inequality in Theorem 2.
Recall that the rates in Theorem 2 are incomparable is to that of Yang et al. (2017, Theorem 2), since we focus on the excess risk while the authors study the squared (IP error. Preciscly, we recover their results as particular case with the square loss and bounded outputs, up to the generalization errors. Instead, note that we do recover the rates of Li et al. (2021, Theorem 1), based on a similar framework. Our bounds feature two different terms: a quantity related to the generalization errors, and a quantity governed by ån, deriving from the K-satisfiabilie bility analysis. The behaviour of the critical radius 8, crucially depends on the choice of the kernel. In Yang et al. (2017), the authors compute its decay rate for different kernels. For instance, we have 5% /lov (n)/n) for the Gaussian kernel, 52 O (1/n) for polynomial kernels, or 52 (n-3//3 for first-order Sobolev kernels. Note finally that by setting An 1/vn we attain rate of O (n), that is minimax for the kernel ridge regression, see Caponnetto & De Vito (2007).
Remark 1. Note that standard additional assumption the second order moments of the functions in Hi (Bartlett et al., 2005) allows to derive refined leorning rates for the generalization errors. These refined rales are expressed in lerms of Fite: the fixed point of new sub-rool function V.. In order to make the approcimation error of the same order, it is then necessary to prove the K K-satisfiability of S with respect to Pie instead of 82,- Whether it is possible to prove such a K K-satisfiability for standard sketches is however a nontrivial question, left as fulure work.
2.2 Matrix-valued Kernel Machines
In this section, we extend our results to multiple output regression, tackled in vector-valued RKHSs. Note that the output space y now a subset of Rdd with d 2 2. We start by recalling important notions about Matrix-Vahun Kernels (MVKs) and vector-valued RKHSs (vv-RKHSs).
Definition (Matrix-valued kernel). A MVK is an application K. X x X - L(RR), where L(Rd) is the set of bounded linear operators on Rd, such that K. (a, x') = KK (x',xx) for all (H, E X2, and such that for all n E N and any (x,, vi) E (X x y)) we have Ei U, IC (I), 2) D3 > O.
Thcorem 3 (Vector-valued RKHS). Let K. be a MVK. There is unique Hilbert space Hx. C F((,R") Rd), the vo- RK-S of K, such that for all E X, y Bd and f He we have x' + K. (xxxx') E Hx. and (J, (.,z) u) = f(x) yy-
5

---

Published in Transactions on Machine Learning Research (09/2023)
Note that we focus in this paper the finite-dimensionaa case, i.e., y C R such that for all x, E X, we have K(xx E Rdxdd For training sample { we define the Gram matrix as K= (K(,,
E
Brdxad, A common assumplion cousists in considering decomposable kernels: we assume that there exist a scalar kernel & and positive semidefinite matrix M Rxdd such that for all x, x, X we have x(x,
The Gram matrix can then be written K K &,,, where K R is the scalar Gram matrix, and @ denotes the Kronecker product. Decomposable kernels are widely spread in the literature they provide a good compromise between computational simplicity and expressivity -note that in particular they encapsulate independent learning, achieved with M = La- We now discuss two examples of relevant output matrices.
Example In joint quantile regression, one is interested in predicting different conditional quantiles of an outpul y given the impul x. If (Ti)ssa (0,1) denole the d different quantile levels, has been shown in - -Ty)²) favors close predictions for close quantile levels,
Sangnier et al. (2016) that choosing My = while limiting crossing effects.
In multiple output regression, it is possible leverage prior knowledge on the task relationships to design a relevant oulput matrix M. For instance, let P be the d d adjacency matrix of a graph in which the vertices are the tasks and edge erists between two tasks if and only if they (thought to be) related. Denoting by T.p the graph Laplacian associaled to P, Engenion et al. (2005) and Sheldon (2008) have proposed to use M (μLp + (1 _ with μ E [o, When μ== independent. When μ 1, we only rely on the prior knowledge encoded in. P.
Example 2.
have M Ia and all tasks considered
Given a sample (x;,::) E (x) and a decomposable kernel K. = kM (its associated vv-RKHS is Hx:), the penalized empirical risk minimisation problem is
min Σ ((f(x;) 1=1
+ An
An 1f
(8)
FEHx
where é Rd Fd + R is a loss such that 2 + ((2,8) y) is proper, lower semi-continuout and convex for all y Rdd By the vector-valued representer theorem (Micchelli & Pontil, 2005), we have that the solution to Problem (8) writes fn = Σラ= (((,x))a = Ej=1 Mā,, where A = an,) Raxd is the solution the problem
An Σ (KAM)) M. + 2
Tr (KAMAT).
min AERuxd
In this conlext, sketching consists in making the substitution = ST where S RRxx is a sketch matrix and I e Raxd is the parameter of reduced dimension to be learned. The solution to the sketched problem is then f. = - M(. ( 1 with I't E Rxd minimizing
n Σ
([KKSTMM v.) 2
An Tr (SKSTTMTT).
Thcorem Suppose that Assumptions to hold, that = KM decomposable kernel with M invertible, and let C in Theorem 2. Then for any 8 E (0,1) with probability at least we have
An + SLV
Tr (M) + 2/
8 Slog (4/8)
E [L4. < [e + LC LLC An + 1|M 52
(9)
n
Purthermore, if ( v) 11=-1112/2 and y B (RR), with probability at least 1 - 8 we have that
E B(G,1 S [ + (o² (c²+b) An 1 of +8 8 TI (M2 ||M/1/2 + 1/12
8 8log (4/8) (10) + +2
n
Proof sketch. The proof follows that of Theorem The main challenge to adapt Yang et al. (2017,
Lemma 2) to the multiple output setting. To do so, we leverage that I. is decomposable, such that the K K-satisfiability of S is sufficient, where K the scalar Gram matrix.

---

Published in Transactions on Machine Learning Research (09/2023)
Note that for M = Ia (independent prior), the third term of the right-hand side of both inequalities becomes of order Vd/n, that is typical of multiple output problems. If moreover we instantiate the bound for d = 1, we recover exactly Theorem 2. Finally, similarly to the scalar case in Theorem 2, looking at the least square case (Equation (10)). by setting A x 1/vn, we attain the minimax rate of D(1/vn), as stated in Caponnetto & De Vito (2007) and Ciliberto et al. (2020, Theorem 5). To the best of our knowledge, Theorem the first thcoretical result about sketched vector-valued kernel machines. We highlight that applies to generic Lipschitz losses and provides a bound directly on the excess risk.
2.3 Algorithmic details
We now discuss how to solve single and multiple output optimization problems. Let ((u v.), E [s]] be the eigempairs of SKST in descending order, Ü = [Uiglaxe (Vi,..., V.), rank(SKS"), and K, = U, and U, = (V1,..., v,).
where D = diag(E1.....
Proposition 1. Solving Problem (3) is equivalent to solving
1
Σ e 2w Es (2x)) m.) + An 1wll
min WERT TL
(11)
2
where zs (x) =
, (z, 2n)) E R..
Problem (11) thus writes as a linear problem with respect to the feature maps induced by the sketch, generalizing the results established Yang et al. (2012) for sub-sampling sketches. When considering multiple outputs, it is also possible to derive a linear feature map version when the kernel decomposable. These feature maps are of the form 2g 1/2 yielding matrices of size nd x rd that are prohibitive in terms of space, Appendix E. Note that alternative way is to see sketching R" (Chatalic el al., 2021). Instead, we directly learn r. For both single and multiple output problems, we consider losses not differentiable everywhere in Section 4 and apply ADAM Stochastic Subgradient Descent (Kingma& & Ba, 2015) for its ability to handle large datasets.
projection of the k(-,x,) into
Remark 2. In the previous sections, sketching is always leveraged in primal problems. However, for some of the loss funclions we consider, dual problems are usually more albractive (Cortes & Vopnik; 1995; Laforgue et al,, 2020). This naturally raises the question of investigating the interplay between sketching and duality on the algorithmic level. More details can be found in Appendia F.
3 p-Sparsified Sketches
We now introduce the p-sparsified sketches, and establish their K-satisfiability. The p-sparsified sketching matrices are composed of i.i.d. Rademacher or centered Gaussian entries, multiplied by independent Bernoulli variables of parameter p (the non-zero entrics are scaled to ensure that S defines an isometry in expectation). The sketch sparsity is controlled by p, and when the latter becomes small enough, contains many columns full of zeros. It is then possible to rewrite S as the product of mbb-Gaussian and sub- samplingg sketch of reduced size, which greatly accelerates the computations.
Definition 5. Let n, and p (0, 1]. A p-Sparsified Rademacher (p-SR) sketching matrix is random malria S E Rx whose entries Sij are independent and identically distribuled (f.i.d) as follows
with probability
y
Vap
Siz
0
with probability
1 p
(12)
with probability a
A p-Sparsified Gaussian (p-SG) sketching matric is a random matrix S E Raxn whose entries Sig i.d.d. follows Sig =
Gas
with probability p
Vpp 0
(13)
with probability
7

---

Published in Transactions on Machine Learning Research (09/2023)
where the Gij are i..dd standard normal random variables. Note that standard Gaussian sketches special case of p-SG sketches, corresponding to p=1.
Several works partially addressed p-SR sketches in the past litcrature. For instance, Baraniuk et al. (2008) establish that p-SR sketches satisfy the Restricted Isometry Property (based on concentration results from Achlioptas (2001)), but only for p 1 and 1/3. In Li et al. (2006), the authors consider generic p-SR sketches, but do not provide any theoretical result outside of a moment amalysis. The i.í.d. sparse embedding matrices from Cohen (2016) are basically m/s-SR sketches, where m 2 1, leading cach column to have exactly m nonzero elements in expectation. However, we were not able to reproduce the proof of the fohnson-Linddrstrauss property proposed by the author for his sketch (Theorem 4.2 in the paper, equivalent to the first claim of K-satisfiability, left-hand side of (4)). More precisely, we think that the assumptions considering "each entry is independently nonzero with probability m/s" and "each column has fixed mumber of nonzero entries" (m here) are conflicting. As far as we know, this is the first time p-SG sketches are introduced in the literature. Note that both (12) and (13) can be rewritten as Sty = (1/vsp) B., Rıj. where the Bij are i..dd Bernouilli random variables of parameter p, and the Rag are i.i.d. random variables, independent from the Bi,, such that [[[,s] = O and Ryy. = if = 1' and j = j', and 0 otherwise. Namely, for p-SG sketches R,y = Gy is a standard Gaussian variable while for p-SR sketches it a Rademacher random variable. It is then easy to check that p-SR and p-SG sketches define isometries in expectation. In the next. theorem, we show that p-sparsified sketches are K-satisfiable with high probability.
Theorem 5. Let S be a p-sparsified sketching matria. Then, there are some universal constants Co, C1 > 0 and a constant c(p), increasing with p, such that for s 2 max (Cod,/p² 83n) and with probability least 2 1 Cie Ce(ps) the sketch is K K-satisfiable for c = vF (1 t /lg (5) +1.
Proof skelch. To prove the left-hand side of (4), we use Boucheron et al. (2013, Theorem 2.13), which shows that any i..d. sub-Gaussian sketch matrix satisfies the To prove the right-hand side of (4), we work conditionally on a realization of the Bij, and use concentration results of Lipschitz functions of Rademacher or Gaussian random variables (Tao, 2012). We highlight that such concentration results do not hold for sub-Gaussian random variables in general, preventing from showing K-satisfiabilitt of generic sparsified sub-Gaussian sketches. Note that having Sij x BiyBi is key, and that sub-sampling uniformly at random non-zero entries instead of using i.i.d. Bernoulli variables would make the proof significantly more complex. We highlight that Theorem 5 strictly generalizes Yang et al. (2017, Lemma 5), recovered for p = 1, and extends the results to Rademacher sketches
lemma with high probability.
J
Computational property of p-sparsified sketches. In addition to be statistically accurate, p-sparsified sketches are computationally efficient. Indeed, recall that the main quantity one has to compute when sketching a kernel machine is the matrix SKST. With standard Gaussian sketches, that are known to be theoretically accurate, this computation takes Osn²)) operations. Sub-sampling sketches are notoriously less precise, but since they act as masks over the Gram matrix K, computing SKST can be done in (ss2) operations only, without having to store the entire Gram matrix upfront. Now, let S E R be a p-sparsifiod sketch, and 8' = Σラ= I(S, * 0.} be the number of columns of S with at least one nonzero element. The crucial observation that makes S computationally efficient is that we have
S = SsG Sss.
(14)
where SsG E IR* is obtained by deleting the null columns from S, and Sas E R xm is a sub-Sampling sketch whose sampling indices correspond to the indices of the columns in S with at least one non-zero entry". We refer to (14) as the decomposition trick. This decomposition is key, as we can apply first a fast sub-sampling sketch, and then sub- aub-Gaussian sketch on the sub-sampled Gram matrix of reduced size. Note that s' is a random variable. By independence of the entries, each column is null with probability (1 - p1)))" Then, by the independence of the columns we have that s' follows a Binomial distribution with parameters n and 1 - (1 - p)"; such that E [s/] = n(1 - (1-p)*)
Preeisely, Sas is the identity matrix Tz,, augmented with n s' mull columns inserted at the indices of the rull columns of S.
8

---

Published in Transactions on Machine Learning Research (09/2023)
Hence, the sparsity of the p-sparsified sketches, controlled by parameter p. is interesting degree of freedom to add: it preserves statistical guarantees (Theorem 5) while speeding-up calculations (14). Of course, there is no free lunch and one looses on one side what is gained on the other: when p decreases (sparser sketches), the lower bound to get guarantecs Z dn/p² increases, but the expected number of non-null columns s' decreases, thus accelerating computations (nole that for p = 1 we exactly recover the lower bound and number of non-null columns for Gaussian sketches).
Corollary 1. Let E Rx be a p-sparsified sketching matrix, and Co, Ct and c(p) as in Theorem 5. Then, setting p 0.7 and = (o.z/(0.79), S is K-satisfiable for c 9, with a probobility at least C2e ()000) These values of p and s minimize compulations while maintoining the guarantees.
Proof. By substituting Codn/p² into E[s'), one can show that it is optimal to set p z 0.7, independently from Co and dn.
Corollary gives the optimal values of and that ensure K-satisfiability of a p-sparsified sketching matrix while having some complexity reduction. However, the lower bound in Theorem 5 is sufficient condition, that might be conservative. Looking at the problem of setting 8 and p from the practitioner point of view, we also provide more aggressive empirical guidelines. Indeed, although this regime is not covered by Theorem 5, experiments show that setting 8 as for the Gaussian sketch and p smaller than 1/8 yield very interesting results, see Figure I(c). Overall, p-sparsified sketches (i) generalize Gaussian sketches by introducing sparsity as a new degree of freedom, (ii) enjoy a regime which theoretical guarantees are preserved and computation (slightly) accelerated, (iii) empirically yield competitive results also in aggressive regimes not covered by theory, thus achieving wide range of intesting accuracy cetnracy/coopppe /computations tradcoffs.
Related works. Sparse sketches have been widely studied in the literature, Clarkson & Woodruff (2017); Nelson & Nguyên (2013); Derczinski et al. (2021). However these sketches are well-suited when applied to sparse matrices (e.g., matrices induced by graphs). In fact, given a matrix A, computing SA with these types of sketching has a time complexity of the order nnz (A), the number of nonzero elements of A. Besides, these sketches usually are constructed such that each column has at least one nonzero element (e.g. CountSketch, OSNAP), hence no decomposition trick is possible. Regarding kernel methods, since a Gram matrix is typically dense (e-g., with the Gaussian kernel, mnz (K) n²), and since no decomposition trick can be applied, one has to compute the whole matrix K and store it, such that time and space complexity implied by such sketches are of the order of n² In practice, we show that we can set p small enough to computationally outperform classical sparse sketches and still obtain similar statistical performance. Note that an important line of research is devoted to improve the statistical performance of Nyström's 's approximation, either by adaptive sampling (Kumar et al., 2012; Wang & Zhang, 2013; Gittens & Mahoney, 2013), or leverage scores (Alaoui & Mahoney, 2015; Musco & Musco, 2017; Rudi et al., 2018; Chen & Yang, 2021b). We took the opposite route, as p-SG sketches accelerated but statistically degraded versions of the Gaussian sketch.
4 Experiments
We now empirically compare the performance of p-sparsified sketches against state-of-he-art approaches, namely Nyström approximation (Williams & Seeger, 2001), Gaussian sketch (Yang et al., 2017), Accumulation sketch (Chen & Yang, 2021a), CountSketch (Clarkson & Woodruff, 2017) and Random Fourier Features (Rahimi & Recht, 2007). We chose not to benchmark ROS sketches as CountSketch has equivalent statistical accuracy while being faster to compute. Results reported are averaged over 30 replicates.
4.1 Scalar regression
Robust regression. We generate dataset composed of 10, 000 training datapoints: 9, 900 input points drawn i.i.d. from 24 ([010, Lill) and 100 other drawn i.idd. from W (1.51 0.25/1o). The outputs are generated as y = fx) + e, where e N~N (N,00) 1) and
4
(*(x) = 421
33r + 204 + 36.
1 + e 0.5)
9

---

Published in Transactions on Machine Learning Research (09/2023)
Nystrom Gaussian
(p=5mm)
(p-30/))
Training time (in seconds)
0.2
8.0
-Nystram
0.1
(p-S/n)
10/)
(p-30m) p-5R (p=50m)
0.09
6.0
0.08
0.04
00.0
00.0
80.0
100
120
140
(a)
Testt relative MSE w.r.b. with K- K-Huber.
(b) Training time (seconds) w.r.b. 8 with R- K-Huber.
0.19
Nystrom
Gaussian
0.18
RFF
p-SR (p=30/n)
Accumulation (m=30)
CountSketch
0.13
E
MSE
0.12
0.1
0.08
0.06
0.01
5.0
Training time (in seconds) 16.0
11.0
(c) Test relative MSE w.r.t. training times with k-Huber.
Figure 1: Trade-off between Accuracy and Efficiency for p-SR sketches with x-Huber loss on synthetic dataset.
as introduced Friedman (1991). We generate a test set of Ate = 10, 000 points in the same way. We use the Gaussian kernel and select its bandwidth -as well as parameters X, and k. (and for eSVR) via. 5-folds cross-validation. We solve this 1D regression problem using the s-Huubrr loss, described in Appendix G. We learn the sketched kernel machines for different values of (from 40 to 140) and several values of p, the probability of being non-null in P-S sketch. Figure ((a) presents the test error a function of the sketch size s. Figure 1(b) shows the corresponding computational training time. All methods reduce their test error, measured in terms of the relative Mean Squared Error (MSE) when s increases. Note that increasing p increases both the precision and the training time, as expected. This behaviour recalls the Accumulation sketches, since we observe a form of interpolation between the Nyström and Gaussian approximations. The behaviour of all the different sketched kernel machines shown in Figure 1(c), where each of them appears point (training time, test MSE). We observe that p-SR sketches attain the smallest possible error (MSE < 0.05) at the lowest training time budget (mostly around 5.6 < túme < 6.6). Moreover, p-SR sketches obtain a similar precision range as the Accumulation sketches, but for smaller training times (both approaches improve upon CountSketch and Gaussian sketch in that respect). Nyström sketching, which similarly to our approach does not need computing the entire Gram matrix, is fast to compute. The method is however known to be sensitive the non-homogeneity the marginal distribution of the input data (Yang al., mixing matrix Ssc in (14) makes p-sparsified sketches more
2017, Section 3.3). In contrast, the robust, empirically shown in Figure 1(c). See Appendix H.1 for results on p-S sketches.
10

---

Published in Transactions on Machine Learning Research (09/2023)
4.2
Vector-valued regression
Table 1: Test pinball and crossing loss and training times (in seconds) with and without sketching (s = 50).
Dataset
Metrics
w/o Sketch
20/m,r-SR
20/m,, SG
Acc. m = 20
CountSketch
Pinball loss
51.28 = 0.67
54.75 + 0.74
54.78 = 0.72
54.73 + 0.75
54.60 0.72
Boston
Crossing loss
0.31 + 0.13
0.26 + 0.08
0.11 + 0.07
0.15 + 0.07
0.10 + 0.05
Training time
6.97 + 0.25
1.43 0.07
1.38 - 0.08
1.48
1.23 + 0.07
Pinball loss
2.78
2.66 1 0.02
2.64 0.02
2.67
2.65 + 0.02
otoliths
Crossing loss
5.18
5.46 0.06
5.43 0.05
5.46 0.066
5.44 1 0.05
Training time
606.8
20.4 + 0.5
20.0 0.3
22.1
20.9 + 0.3
Joint quantile regression. We choose the quantile levels as follows T = (0.1, 0.3, 0.5, 0.7, We apply subgradient algorithm to minimize the pinball loss described in Appendix G with ridge regularization and a kernel K== kM with M discussed in Example 1, and & a Gaussian kernel. We select regularisation parameter An and bandwidth of kernel o² via a 5-fold cross-validation We showcase the behaviour of the proposed algorithm for Joint Sketched Quantile Regression on two datasets: the Boston Housing dataset (Harrison Jr & Rubinfeld, 1978), composed of 506 data points devoted to house price prediction, and the Fish Otoliths dataset (Moen et al., 2018; Ordoñez et al., 2020), dedicated to fish age prediction from images of otoliths (calcium carbonate structures), composed of train and lest sets of size 3780 and 165 respectively. The results are averages over 10 random 70% 30% train-test splits for Boston dataset. For the Otoliths dataset we kept the initial given train-test split. The results are reported in Table 1. Skotching allows for a massive reduction of the training times while preserving the statistical performances. As a comparison, according lo the results of Sangnier et al. (2016), the best benchmark result for the Boston dataset in terms of test pinball loss is 47.4, while best test crossing loss is 0.48, which shows that our implementation does not compete in terms of quantile prediction but preserves the non-crossing property.
Table 2: ARRMSE and training times (in sec) with square loss and 8 = 100 when using Sketching,
Dataset
Metrics
w/o Sketch
20/
S0/n--SSG
Acc. m = 20
CountSketch
ARRMSE
0.575
0.584 = 0.003
0.583 = 0.003
0.592 + 0.001
0.575 = 0.0005
rfl
Training time
1.73
0.22 1 0.025
0.25 1 0.005
0.60 1 0.0004
0.66 1 0.013
ARRMSE
0.578
0.671 = 0.009
0.656 + 0.006
0.796 0.006
0.715 0.011
rf2
Training time
1.77
0.28 = 0.003
0.27 : 0.003
0.82 + 0.003
0.62 + 0.001
ARRMSE
0.418
0.422 + 0.002
0.123 + 0.001
0.423 + 0.001
0.420 0.001
scmld
Training time
9.36
0.45 + 0.022
0.45 - 0.019
0.86 - 0.006
2.49 = 0.035
ARRMSE
0.755
0.754 + 0.003
0.754 + 0.003
0.753 + 0.001
0.754 + 0.002
scm20d
Training time
6.16
0.38 + 0.016
0.38 0.017
0.70 0.032
1.91 + 0.047
Multi-output regression. We finally conducted experiments on multi-output kernel ridge regression. We used decomposable kernels, and took the largest datasets introduced in Spyromitros- Xioufis et al. (2016). They consist in four datasets, divided in two groups: River Flow (rfl and rf2) both composed of 4108 training data, and Supply Chain Management (semld and scm20d) composed of 8145 and 7463 training data respectively (more details and additional results can be found in Appendix H.2). We compare our non-sketched decomposable matrix-valued kernel machine with the sketched version. For the sake of conciseness, we only report here the Average Relative Root Mcan Squared Error (ARRMSE), see Table 2 and Appendix H.2. For all datasets, sketching shows strong computational improvements while maintaining the accuracy of non-sketched approaches.
Note that for both joint quantile regression and multi-output regression the results obtained after sketching (no matter the sketch chosen) are almost the same as that attained without sketching. It might be explained by two factors. First, the datasets studied have relatively small training sizes (from 354 training data for Boston to 8145 for semld). Second, prodicting jointly multiple outputs is a complex task, so that it appears more natural to obtain less differences and variance using varions types of sketches (or no sketch). However, in all cases sketching induces huge time saver.
11

---

Published in Transactions on Machine Learning Research (09/2023)
5 Conclusion
We proposed excess-risk bounds for sketched kernel machinos in the context of Lipschitz-contimoous losses, with results valid for both scalar and matrix-valued kernels. We introduced a novel sketching scheme that leverages the good empirical statistical guarantees of the Gaussian Sketching while combining them with the low cost of Nyström sketching. Numerical experiments show that this novel scheme opens the door to many applications beyond the squared loss. Improvements on multi-output regression can certainly be obtained by applying low-rank considerations in the output space as well.
Acknowledgments
The authors thank Olivier Forcoq for insightful discussions. This work was supported by the Télécom Paris research chair on Data Science and Artificial Intelligence for Digitalized Industry and Services (DSAIDIS) and by the French National Research Agency (ANR) through the LNR-18-CE23-06 APi project.
References
Dimitris Achlioptas. Database-friendly random projections. In Proceedings of the twentieth ACM SIGMOD- SIGACT-SIGARY symposium on Principles of database systems, pp. 274-281, 2001.
Ahmet Alacaoglu, Olivier Fercoq, and Volkan Cevher. Random extrapolation for primal-dual coordinate descent. In Proc of the International Conference on Machine Learning (ICML), pp. 191-201. PMLR, 2020.
Ahmed Alaoui and Michael W Mahoney. Fast randomized kernel ridge regression with statistical guarantecs. In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett (eds.), Advances in Neural Information Processing Systems (NeurIPS), volume 28, 2015.
Mauricio A. Álvarez, Lorenzo Rosasco, and Neil D. Lawrence, Kernels for vector-valued functions: a review. Foundations and Trends in Machine Learning, 4(3):195-266, 2012.
Nachman Aronszajn. Theory reproducing kernels. Transactions of the American Mathematical Society, pp- 337-404, 1950.
Haim Avron, Kenneth L Clarkson, and David P Woodruff. Faster kernel ridge regression using sketching and preconditioning. SIAM Journal on Matriz Analysis and Applications, 38(4):1116 1138, 2017.
Francis Bach. Sharp analysis of low-rank kernel matrix approximations. In Proc. of the 26th annual Conference Learning Theory, pp. 185 259- PMLR, 2013.
Luca Baldassarre, Lorenzo Rosasco, Aunalisa Barla, and Alessandro Verri. Multi-oulput learning via spectral filtering. Machine Learning, 87(3):259 301, 2012.
Richard Baraniuk, Mark Davenport, Ronald DeVore, and Michael Wakin. A simple proof of the restricted isometry property for random matrices. Constructive Approcimation, 283))2533-260) 2008
Peter L. Bartlett and Shahar Mendelson. Rademacher and gaussian complexities: Risk bounds and structural results. J. Mach. Learn. Res., 3:463-482 2003.
Peter L. Bartlett, Olivier Bousquet, and Shahar Mendelson. Local rademacher complexities. Ann. Statist., 33 (4):1497 1537, 2005.
Poter I. Bartlett, Michael Jordan, and Jon D McAuliffe. Convexity, classification, and risk bounds. Journal of the American Statistical Association, 101(473):135 138-1156, 2006.
Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration inequalities: A nonasymptotic theory of independence. Oxford university press, 2013.
Romain Brault, Markus Heinonen, and Florence Buc. Random fourier features for operator-valued kernels. In Asian Conference Machine Learning, pp. 110-125. PMT.R, 2016.
12

---

Published in Transactions on Machine Learning Research (09/2023)
Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squares algorithm. Founda- lions of Compulational Muthematics, 7(3):331-308 2007.
Claudio Carmeli, Ernesto De Vito, and Alessandro Toigo. Vector valued reproducing kernel hilbert spaces of integrable functions and mercer theorem. Analysis and Applications, 4(04))377-408, 2006.
Claudio Carmeli, Ernesto De Vito, Alessandro Toigo, and Veronica Umanitá. Vector valued reproducing kernel hilbert spaces and universality. Analysis and Applications, 8(01):19-61, 2010.
Antonin Chambolle, Matthias J Ehrhardt, Peter Richtárik, and Carola-Bibiane Schonlieb. Stochastic primal- dual hybrid gradient algorithm with arbitrary sampling and imaging applications. SIAM Journal on Oplimization, 28(4):2783 2808, 2018.
Antoine Chatalic, Luigi Carratino, Ernesto De Vito, and Lorenzo Rosasco. Mean nyström embeddings for adaplive compressive learning, 2021.
Yifan Chen and Yun Yang. Accumulations of projections unified framework for random sketches in kernel ridge regression. In International Conference on Artificial Intelligence and Stalislics, pp. 2953 2961. PMLR, 2021a.
Yifan Chen and Yun Yang. Fast statistical leverage score approximation in kernel ridge regression. Tn International Conference on Artificial Intelligence and Statistics, pp. 2935 2943. PMLR, 2021b.
Carlo Ciliberto, Lorenzo Rosasco, and Alessandro Rudi. A general framework for consistent structured prediction with implicit loss embeddings. J. Mach. Learn. Res., 21(98):1 67, 2020.
Kenneth L. Clarkson and David P. Woodruff. Low-rank approximation and regression in input sparsity time. J. ACM, 63),,, jan 2017. ISSN 0004-5411. doi: 10.1145/3019134. URL https: http:://dooo.ooo// 1145/3019144
Michael B. Cohen. Nearly tight oblivious subspace embeddings by trace incqualitics. In Proceedings of the 2016 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pp.
Laurent Condat. A primal-dusl splitting method for convex optimization involving lipschitzian, proximable and linear composite terms. Journal of optimization theory and applications, 158(2):460 479, 2013.
Corinna Cortes and Vladimir Vapnik. Support-vector networks. Machine learning, 30(3))273-22) 1995.
Andrea Della Vecchia, Jaouad Mourtada, Ernesto De Vito, and Lorenzo Rosasco. Regularized erm on random subspaces. In International Conference on Artificial Intelligence and Statistics, pp. 4006-40144 PMLR, 2021.
Michal Derezinski, Zhenyu Liao, E. Dobriban, and Michael W. Mahoney. Sparse sketches with small inversion bias. In COLT, 2021.
Petros Drineas, Michael W Mahoney, and Nello Cristianini. On the nyström method for approximating gram matrix for improved kernel-based learning. JMLR, 6(12), 2005.
Harris Drucker, Christopher JC Burges, Linda Kaufman, Alex Smola, and Vladimir Vapnik. Support vector regression machines. In Advances in neural information processing systems, pp. 155 161, 1997.
Thcodoros Evgeniou, Charles A. Micchelli, and Massimiliano Pontil. Learning multiple tasks with kernel methods. Journal of Machine Learning Research, 6(21):615 637, 2005.
Olivier Fercoq and Pascal Bianchi. A coordinate-descent primal-dual algorithm with large step size and possibly nonseparable functions. SIAM Journal on Oplimization, 29(1):100 134, Jan 2019.
Jerome H Friedman. Multivariate adaptive regression splincs. The annals of statistics, pp. 1-67, 1991.
Nidham Gazagnadou, Mark Tbrahim, and Robert M. Gower. Ridgeskotch: A fast sketching basod solver for large scale ridge regression. SIAM Journal on Matrir Analysis and Applications, 43(3):1440-1468 2022. doi: (..1377/11M1422 URL https://doi. htt/ddi.oorgloe org/10.0.//
13

---

Published in Transactions on Machine Learning Research (09/2023)
Alex Gittens and Michael Mahoney. Revisiting the nystrom method for improved large-scale machine learning. Proceedings of the 30th International Conference Machine Learning, 28(3):567-55) 575, 17-19 Jun 2013.
William Groves and Maria Gini. On optimizing airline ticket purchase timing. ACM Transactions on
Intelligent Systems and Technology (TIST), 7(1):1-28, 2015.
David Harrison Jr and Daniel L. Rubinfeld. Hodonic housing prices and the demand for clean air. Journal of environmental economics and management, 5(1):81 102, 1978.
Matthias Hein and Olivier Bousquct. Kernels, associated structures and gencralizations. Technical Report 127, Max Planck Institute for Biological Cybernetics, Tübingen, Germany, July 2004.
Peter J Huber. Robust estimation of location parameter. The Annals of Mathemalicnl Stalistics, pp. 73 101, 1964.
William B Johnson and Joram Lindenstrauss. Extensions of lipschitz mappings into a hilbert space 26. Contemporary mathematics, 26:28, 1984.
George Kimeldorf and Grace Wahba. Some results on tchebycheffian spline functions. Journal of mathematioal analysis and applications, 33(1):82-99) 1971.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In Yoshua Bengio and Yann LeCun (eds.), 3rd International Conference on Learning Representations, ICL. 2015, Son Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
Roger Koenker. Quantile regression. Cambridge university press, 2005.
Samory Kpotufe and Bharath K. Sriperumbudur. Gaussian sketching yiclds J-L lemma in RKHS. In Silvia Chiappa and Roberto Calandra (eds.), AISTATS 2020, volune 108 of Proceedings of Machine Learning Research, pp. 3928 3937. PMLR, 2020.
Sanjiv Kumar, Mehryar Mohri, and Ameet Talwalkar. Sampling methods for the nyström method. Mach. Learn. Res., 13:981-1006, 2012.
Jonathan Lacotte and Mert Pilanci. Adaptive and oblivious randomized subspace methods for high high-ddmensiionl optimization: Sharp analysis and lower bounds. IEEE Transactions on Information Theory, (8(5))3281-3303 2022. doi: 10.1109/TII ALL109/TTT202221 202..3166666
Jonathan Lacotte, Mert Pilanci, and Marco Pavone. High-dimensional optimization in adaptive random subspaces. In Proc. of the 33rd International Conference Neural Information Processing Systems, pp. 10847-10857, 2019.
Pierre Laforgue, Alex Lambert, Luc Brogat-Motte, and Florence d' Alché Buc. Duality in rkhss with infinite dimensional outpuls: Application to robust. losses. International Conference on Machine Learning, pp- 5598 5607. PMLR, 2020.
Ping Li, Trevor J Hastie, and Kenneth W Church. Very sparse random projections. In Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 287-296, 2006.
Zhe Li, Tianbao Yang, Lijun Zhang, and Rong Jin. Fast and accurate refined nyström- based kernel svm. Thirtieth AAAI Conference on Artificial Intelligence, 2016.
Zhu Li, CaanFFaaccois Ton, Dino Oglic, and Dino Sejdinovic. Towards unified analysis of random fourier features. Journal of Machine Learning Research, 22(108):1-51, 2021.
Meimei Liu, Zuofeng Shang, and Guang Cheng. Sharp theoretical analysis for nonparametric testing under
random projection. In Alina Beygelzimer and Daniel Hsu (ods.), Proceedings of the Thirty-Second Conference on Learning Theory, volume 99 of Proceedings of Machine Learning Research, pp. 2175-2209. PMLR, 25-28 Jun 2019.
14

---

Published in Transactions on Machine Learning Research (09/2023)
Michael W Mahoney et al. Randomized algorithms for matrices and data. Foundations and Trends& in Machine Learning, 3(2):123-22) 2011.
Jifí Matoušck. Lectures on Discrete Geometry. Graduate Texts in Mathematics, Springer, 2013.
Andreas Maurer. A vector-contraction incquality for rademacher complexitics. In International Conference on Algorithmic Learning Theory, pp. 3-17. Springer, 2016.
Giacomo Mcanti, Luigi Carratino, Lorenzo Rosasco, and Alessandro Rudi. Kernel methods through the roof: Handling billions of points efficiently. Advances in Neural Information Processing Systems (NeurIPS), 33, 2020.
Charles A Micchelli and Massimiliano Pontil. On learning vector-valued functions. Neural computation, 17 (1):177-204, 2005.
Endre Moen, Nils Olav Handegard, Vaneeda Allken, Ole Thomas Albert, Alf Harbitz, and Ketil Malde. Automatic interpretation of otoliths using deep learning. PLoS One, 13(12):c0204713, 2018.
Cameron Musco and Christopher Musco. Recursive sampling for the nyström method. Advances in Neural Information Processing Systems, 2017:3834-388 3846, 2017.
Jelani Nelson and Huy L. Nguyên. Osnap: Faster numerical linear algebra algorithms via sparser subspace embeddings. In 2013 IEEE 54th Annual Symposium on Foundations of Computer Science, pp. 117-126, 2013. doi: 10.1109/ N0.1109/FOCS.D0 FOCS.2033.1
Alba Ordoñez, Line Eikvil, ment-Bereeeeee Salberg, Alf Harbitz, Sean Meling Murray, and Michael C Kampffmeyer. Explaining decisions of deep neural networks used for fish age prediction. PloS one, (5(6):00233013, 2020.
Ali Rahimi and B. Recht. Random features for large scale kernel machines. NIPS, 20:1177 1184, 01 2007.
Alessandro Rudi and Lorenzo Rosasco. Generalization properties of learning with random features. In Advances on Neural Information Processing Systems (NeurIPS), pp. 3215 3225, 2017.
Alessandro Rudi, Raffaello Camoriano, and Torenzo Rosasco. Less is more: Nyström computational regularization. Advances in Neural Information Processing Systems, 28, 2015.
Alessandro Rudi, Daniele Calandriello, Luigi Carratino, and Lorenzo Rosasco. On fast leverage score sampling and optimal learning. In NewrIPS, 2018.
Maxime Sangnier, Olivier Fercoq, and Florence d'Alché Buc. Joint quantile regression in vector-valued RKHSs. In Advances in Newral Information Processing Systems (NeurIPS), Barcelona, France, December 2016.
Bernhard Scholkopf and Alexander J Smola. Learning with kernels: support vector machines, regularization, optimization, and beyond. Adaptive Computation and Machine Learning Series, 2018.
Bernhard SchölkopI, Ralf Herbrich, and Alex Smola. A generalized representer theorem. In Compulational Learning Theory: 14th Annual Conference on Computational Learning Theory, COLT 2001 and 5th European Conference on Computational Learning Theory, EuroCOLT 2001 Amsterdam, The Netherlands, July 16-19, 2001 Proceedings 14, pp. 416-426. Springer, 2001.
Daniel Sheldon. Graphical multi-task learning, 2008
Eleftherios cinoomitros-iioufin Grigorios Tsoumakas, William Groves, and Ioannis Vlahavas. Multi-target regression via input space expansion: treating targets as inputs. Machine Learning, 104(1):55 98, 2016.
Bharath K Sriperumbudur and Zoltán Szabó. Optimal rates for random fourier features. In NIPS, 2015.
Ingo Steinwart and Andreas Christmann. Sparsity of svms that use the cpsilon-insensitive loss. In Daphne Koller, Dale Schuurmans, Yoshua Bengio, and Léon Bottou (eds.), Advances in Neural Information Processing Systems 21 (NeurIPS), pp. 1569-1576. Curran Associates, Inc., 2008a.
15

---

Published in Transactions on Machine Learning Research (09/2023)
Ingo Steinwart and Andreas Christmann. Support vector machines. Springer Science & Business Media, 2008b.
Ingo Steinwart and Andreas Christmann. Estimating conditional quantiles with the help of the pinball loss. Bernoulli, 17(1):211-223 2011.
Michel Talagrand. Concentration of measure and isoperimetric inequalities in product spaces. Publications Mathématiques de l'Institut des Hautes Etudes Scientifiques, 81(1):73 205, 1995.
Terence Tao. Topics in random matrir theory, volume 132. American Mathematical Soc., 2012.
Bang Cong Vu. A splitting algorithm for dual monotone inclusions involving cocoercive operators, 2011.
Shusen Wang and Zhihua Zhang, Improving cur matrix decomposition and the nyström approximation via adaptive sampling. J. Mach. Learn. Res., 14(1):2729-25) 2769927679 jan 2013.
Christopher Williams and Matthias Seeger. Using the nyström method to speed up kernel machines. In T. Leen, T. Dietterich, and V. Tresp (eds.), Advances in Neural Information Processing Systems, volume 13, pp. 682-688. MIT Press, 2001.
David P. Woodruff. Sketching as tool for numerical linear algebra. Found. Trends Theor. Comput. Sci., 10 (1-2):1-157, 2014. doi: 10. 10.1561 00.166000400000000 URL https://doi.org/ 15t1//.0..1100000000000000
Tianbao Yang, Yu-feng Li, Mchrdad Mahdavi, Rong Jin, and Zhi-Hua Zhou. Nyström method vs random fourier features: A theoretical and empirical comparison. In F. Pereira, C. C. Burges, L. Bottou, and K. Q. Weinberger (eds.), Advonces in Neurol Information Processing Systems, volume 25, Curran Associates, Inc., 2012.
Yun Yang, Mert Pilanci, Martin J Wainwright, et al. Randomized skotches for kernels: Fast and optimal nonparametrie regression. The Anmals of Slalistics, 45(3):991- 1023, 2017.
Kai Zhang, Liang Lan, Zhuang Wang, and Fabian Moerchen. Scaling kernel svm on limited resources: A low-rank linearization approach. In Neil D. Lawrence and Mark Girolami (eds.). Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics, volume 22, pp. 1425-1433 2012.
16

---

Published in Transactions on Machine Learning Research (09/2023)
A Technical Proofs
In this section are gathered all the technical proofs of the results stated in the article.
Notation. We recall that we assume that training data (#6 We are i.i.d. realisations sampled from a joint probability density P(x,8) y). We define
En[ey] 1 Σ {(f(z)).u.
= [r()))(((((()))))
For class of functions F, the empirical Rademacher complexity (Bartlett, & Mendelson, 2003) is defined
2
R = E
sup fEF
0. (x:)
[1
En
where 01....,0, Un are independent Rademacher random variables such that Pfa; = = Po; Poo, -11- The corresponding Rademacher complexity is then defined as the expectation of the empirical Rademacher complexity
(F. = B R. (F))
A.1 Proof of Theorem 2
We first prove the first inequality in Theorem 2 for generic Lipschitz losses.
Proof. The proof follows that of Li al. (2021, Theorem 3). We decompose the expected learning risk
B B[e E[e,,] 1,, + E,le,,1 EnIL Enlesnn -E(L 1. We then use Bartlett & Mendelson (2003, Theorem 8) to bound Ele,,] E, [, and E, le (ee@ Ele
(7)
Lemma 1. (Bartlett & Mendelson, 2003, Theorem 8) Let (Ii, yil be i.d.d samples from P and let H. be the space of functions mopping from X to R. Denote loss function with 1 3 R [0, 1] and recall the learning risk function for all f E H is B [1,, together with the corresponding empirical risk function En [ = (1/n) Ei= I (Vi, (x:)). Then, for sample of size n, for all f E H and 8 E (0, 1), with probability 5/2, we have that
8log
E[4y < En [1y] + Rn
n
where 1 H ((.,1))) y) Ml. ((x)) - 1yy, 0) H).
Thus, since T. lies in the unit ball B (H²) of He by Assumption 2, we obtain thanks to the above lemma, with a probability at least 1 8
8 llog(2/))
E
[ En 1 S R, ( +
n
Then, by the Lipschitz continuity of (Assumption 3) and point of Theorem 12 from Bartlett & Mendelson (2003), we have that
R2 (f B (Hx)) - 2LR, (B (H+)).
Finally, Assumption 4 combined with Lemma 22 from Bartlett & Mendelson (2003) then yiclds
2 R₁ (5 (H+))= <
Σ k (x,, S /// K
17

---

Published in Transactions on Machine Learning Research (09/2023)
As consequence, we obtain
E (r.]
En [et,] S
4LVR vn
8 Blog(4/0)
(15)
m
and the exact same result applies [e Ellsnell by Assumption 2 and the opposite side of Lemma
We now focus on the last quantity to bound. Let Hs = {f = Ei= [S *) I E R²). By Assump- tions and and Jensen's inequality we have
En [1.] IEm
(fA(t.)) - 1 Σ e ((iH)
+ An ||
-
1
(JHe (x)).
n
2
n
inf JEHs 1rlbg <1
((f(x)) 1 Σ
(Sn. ((:)
An
115
L (x4) - SHe (2x) + An
<
inf SENt
2
i=1
1 n
fuu (x:))+ An
S L inf JEHs 11lb,
Σ
1f (x,)
n
_ 11fx -fp-- fiell2 112e+
An 2
=L
inf JeHs
1i 51
where, for any f E Hd. fx = (f(z)).... f(In)) E R". Let f.t = Ei=1 [ST k (-,z,,, where & is a solution to
inf yER
IKST - + Xn SKST
It casy to check that fā also a solution
1 1
- 1155-
- 11
inf JEHs 1/ slir? lbx
Since we have | 1x < by Assumption 2, it holds
inf 1 FEHg n lillng 21
||
inf
1 ||5x fit/ 1|2
JEHs 11 Ellffll
n
inf 1 7ER* n
|KKS 7 - 1xx 112+ >nn SKSTT
As a consequence,
E, [ En [enma] s L
inf 1 11x + An A SKST y+ An yER
2
Finally, since K-satisfiable sketch matrix, using Lemma 2 from Yang et al. (2017),
E, [ E, [emm.] < SLCVAn+ +
Am
(16)
where C 1 V6c and is a universal constant coming from K-satisfiable property. The desired bound is obtained by combining Equations (7), (15) and (16).
18

---

Published in Transactions on Machine Learning Research (09/2023)
We now give the proof of the second claim ie., the excess risk bound for kernel ridge regression.
Proof. We now assume that the outputs are bounded, hence, without loss of generality, y C [0, 1]. First, we prove Lipschitz-contimuity of the square loss under Assumptions 2 and 4. Let (x = (f (z) f E H&, E }} and g E Ha (2) + 1 (z y)², for y E [0. Hence, g'(z)
y and by
Assumptions 2 and 4
lo((:) S ||| |l S [f(x)) S 10 ()))) + S |1 pe, x2/2 1= 1,
for some f E Hu. and I E X since 2 E Hu (2). As a consequence, we obtain that g is (x112 + 1) -ippschhi, i.e.
(e 9) - 1 (f (z') v) - (2 + ( (x)
We can then obtain the same gencralisation bounds as above. Finally, looking at the approximation term, 112 2n
En [ - En [
1
|1
|| 112
2n
1
|1
11
2n
inf
||
X 112
+ An 2
SEHs
2n
1
1
An TKST
inf yERt n
4+
2
$ C (2+2))= An + c
We oblain the same bound as above without the square root since we do not need to invoke Jensen inequality as the studied loss is the square loss. Gathering all arguments, we obtain last inequality in Theorem 7
A.2 Refined analysis in the scalar case
As said in Remark 1, and similarly to Li et al. (2021), we can conduct a refined analysis, leading to faster convergence rates for the generalization errors, with the following additional assumption.
Assumption 6. There is a constant B such that, for all f E H4. we have
Blf - fnal 12 < BE ey e
(17)
It has be shown that many loss functions satisfy this assumption such as Hinge loss (Steinwart & Christmann, 2008b; Bartlett et al., 2006), truncated quadratic or sigmoid loss (Bartlett et al., 2006). Under Assumptions to and 6, the following result holds:
Theorem 6. We define, for 5 (0,1), the following sub-root function
Snn
1/2
Bn(r) 21C1
2
min (b2r, (orr,r)) mi)
C2 n
log 8
and let fits the fixed point of Sn, i.e, On (Phan) = fi Then, we have for all > 1 and ô E (0, 1) with probability greater than 1 1-8,
8,
[
< E [fmm
D-1 [I.CVin+ + + An
D
+
12D
H. +
2Ca log 5'
(18)
B
where C is in Theorem 2 and Ci, Cz, Co and b2 are some constants and Pie can be upper bounded by
h bo
H4
<
min oSh
+
Σm
Mi
n
where B and by are some constonts.
19

---

Published in Transactions on Machine Learning Research (09/2023)
Hence, that, in order to obtain faster learning rates than Theorem as Li et al. (2021), need to replace og by f However, according the expression of 2h and its dependencies to non-explicit constants,
1/2 it appears very difficult to prove that (1 En= min(fri. m4) S AH which is a necessary condition prove that sketch matrix S is K K-satisfiable. We still prove the above result following the proof of Theorem 4 in Li et al. (2021), and leave as an open problem to find faster rates than 6n.
Proof. We use the decomposition of the expected learning risk from Li et al. (2021)
D En f..1
E Blesnn]
E
D-1
D (Bulkg.) Enlesny
D
D
En Enlesny
Bleinel)
D-1
for D 1. The generalization errors can be bounded as in Li al. (2021) by Theorem 6. The approximation error is bounded using (16).
A.3 Proof of Theorem 4
Here, the proof uses the same decomposition of the excess risk (Equation (7)) as in single output settings. Since some works (Maurer, 2016) exist to easily extend generalisation bounds of functions in scalar-valued RKHS functions in vector-valued RKHS, the main challenge here is derive an approximation error for the multiple output settings. Hence, let us first state needed intermediate results that will prove later.
Lemma 2. For all f E Hx, such that || & 1, we have z (K-1 & M-1) 2 < 1, where = ((,,),...,,,(EE)) E RRdd
We are now equipped to state the main result that generalises Lemma 2 from Yang et al. (2017).
Lemma 3. Let Z = (I'(x),, (*m))) E Raxd for any J* E Hx. such that || S 1, where K. = kM, and S E paxn K- satisfiable matrix. Then we have
inf 1 LKKSTMM- - 1I + An Ir (KS PMT 9) s C² (11M +
(KSTAM
(19)
FERoxd n
where C C=1+ Vbc and e is the universal consland from Defivition 3.
Proof. We adapt the proof of Lemma 2 from Yang et al. (2017) to the multidimensional case. If we are able to find a I E Rxd such that
IKSTM 1 + X, Tr (KS 9) < C² (ll *m),
then in particular it also holds true for the minimizer. We recall the eigendecompositions 1K Knorm T and M = VAVT. Then the above problem rewrites
IDS - - + An Tr (ŠDSTTMT) S C² (l1Ml Am)
where $= SU and 0
1 T Z.V. We can rewrite H* (0) ....,... 0t) 2// (UT 8 (VT)=² 2", hence
n1/9
(D-1/2 1/2 & A 11/2) 0 11 = x" (K 2 M-1) x", with x" = (Zi..... 2 = =s(((z)) (..., ((.,) By
Lemma. 2, we have that 1/2 & A 112 0*11a < 1, and using the notation 7 (Ii....., Iw.) E gdd we can rewrile the above problem as finding 7 such that
1 (DST & AVT) r12 + A (SDST so M) 5 C2 (llM + Aaa)
20

---

Published in Transactions on Machine Learning Research (09/2023)
As in (Yang et al., 2017), we partition vector H E Rnd into two sub-vectors, namely 8 E Rddd and 02 E R(inddadd, the diagonal matrix D into two blocks Dr E RRdx and D2 E (((-))))(-- and finally, under the condition > dn. we let S E Raxd, and S2 E (n-da) denote the left and right block of $ respectively. By the K K-satisfiablility of S we have
S Td,,lop S and 2|
1lop S 2cd
(20)
By the first inequality, we have that S. S1 is invertible. In fact, assuming that there exists x E IR such that ||al|2 = and ST S X = O, then ||(5] $1 - I.. )all = > タ Then, we can define
3 = (5) (ST DI @VA-\ 09.
(21)
Hence,
110- - (DST & AV T)ll avT)) = |0 - (D, @ AVTT + (D2S] @ A
and we have
||0 SD @ AV M 2 |10 (DIST @
DJ' & VA 001112
|00
(DIST S S)) AV 08:11
|10 - 0111 0,
and
10: - (D2SI @ AVT) 1/12
8/ (D2S] S, (ST 5.) -1 DJ @ 12) 0:
S |0:] 122 +
() 1, -,- 1/2 & 1//2 0i 2
= 105112 +
( DJ 8 A412) (D 8 1011280-111 0-112) 01 2
< |10 22 +
(D2SI 3, (5T5,) 1/2 8 122
(D 1/2 8 A A-12
1/2
01
0p
105112 +
D25} S, (S) -1 DJ 1/2 1
A
[ 1/2
& A 1/2
1/2)
0
< |0511|
|D/L
/.DD//1
|| 1lop
())
D1
1/2
1/2
[1D-
112
83 11/2) 07
op
op
op
op
(22)
We now bound all lerms involved in (22). Since (D- 6 -1/3 0*1|2 S 1, then || and,
(Dr 1/2 3 A -1/3) 0ī S
10112
02
1=1 j=d +1
(3 2
$ 2 op
Asi
H4j
1=1
j=d +1
05
S M op
1=1 1=1
M Ilon
(D-1/2 @ A A1122) 0-
MMMM
21

---

Published in Transactions on Machine Learning Research (09/2023)
since My 8 for all j 2 dn + and An S |oM for all i d. Moreover, since ||ST S1 - Ian llop -
|| Sillop S 3 then 2; || lopp S 2 Besides, for a E Rd such that || 1, we have
||ST S1 - 11 - |IST S 12 - | - (1( S1 -Id.) x 2=1222 2.2
1 2
Then, we obtain that ||| S - and then |||] S 2 1,2 taking the eigenvector of S corresponding its smallest eigenvalue, obtain that || (ST Si) Moreover we have
Ilap' 2 1 and finally II (i) -1 llop s 2.
|0 1/2 1lop 1 pn¹
1
|D! 1lop S 6n:
|| 1/ 1lop I cda-
Thus,
1105 (D25g @ AV ( - (8.1 1/2 + 6n,
3
1/2
2-
|
2
8n
(1)))))) 1/2 (1 (1+c1) (v66
Finally,
0 (DS & AV a11 s MM (1+ (1+oV))²
(23)
Furthermore, looking into the second term,
ST (ŠDST @ M) 4- = D1/2ST & A
8)1 112
(D! & A
+ (D!/35 8 T )11/
(D, 1/2 & 1/2
2
(Dy ST3, S, (ST -1 DJ 8 11/2 A112 01
01
+
2
2
S + |3DD] 1/2 /212
|
1
2
(ST 5.) -1 -1
D D 1/2 2
(D. 1/2 - 11/2 01
2
op
2
21+ e 822
2 0%
= =1 6c
(1 + v6c)
2/
< (1 voc)
Finally, we obtain that
|10- (DST @ AVV 1112+ NNAT (ŠDST &M) s (1+ voe) [looo + +2n),
(24)
and as a conclusion
inf 1 IKSTTM-- - 2 A, Tr (KSTTMTT 8) s (11 An), P n
(25)
where C = Vbc.
Now, for the proof of Theorem 2, let us prove equation first inequality in Theorem 4.
22

---

Published in Transactions on Machine Learning Research (09/2023)
Proof. For any function in B (Hc)) {f E Ni llf $ 1), Lemma still holds, then
(lgg(2/))
E [ez) S En [ey]+ (1 B (Hb)) +
(26)
Then, using Corollary from Maurer (2016), we have that:
lf, B (Hxc) I V2LR, (Bxc))
(27)
where
d
R₂, = E
sup
sj (ax)))
I1
JEF
1= j=1
= E
sup SEP
Σ 0 / (x:)))
11 I0,,.., In
where 01...... 011....,0p are nd independent Rademacher variables, and for all 1 S sis o; = Os Hence
2 Σ i=1
(n(B(()))) E
sup
(o,, f(x)))
, En
|If |lf $1
E
sup
Krn, oz:
x1,-.. In
||
S1
1=1
Hx
1/2
2 E
K
F1,
Fn
t=1
Hx
1/2
2 E
(a,, (C(1,, Z1.
In
n
iリ=1
1/2
2
E
A(I,,, (ar,, Mar) (.....,, Tm
n
1/2
2
k(),, 2.)
E My (n,,..
n
f='==1
1/2
2
Σ k( 2.))
Σ Mye
n
2
[Tr (K @ M/)
2
(B S
1/2 Tr (M2
1/2
Finally, for any function f E B (Hx)), for all 5 (0, 1), we have for probability least - 8,
2x 4L
8 log (2/0)
IE [Cy] E., [Cr]
Tr (M) +2
(28)
n
n
23

---

Published in Transactions on Machine Learning Research (09/2023)
Now, for the approximation error term, proceed as in the proof of Theorem 2. Let Hs E ; 1 E Rexd]. By Assumptions 2 and 3 and Jensen's 's inequality,
n 7
n
En (%/
Bn
(f.())))))
Σ (fri (21),
1
e (f..)))
+ An ||
n E (friee
1
1 ((f(a)).B1) - Σ
1
e (frne + An
inf fEHs 1/llg s1
llfl
==1
Σ 1f (x:) SHC (2:)l12 + An
$
inf JEHs
n
2
1/llg $1
iml
Σ (x,)) fre (2t) An
L
inf JEHs 1ill $1
=||x 1 112 +
An
=
inf JEHs
2
1 21
where, for any f E Hk, fx (f(zi),..... f(xn)) E Ruud Let f. solution to
=
x.), [STi"], where [ is
inf TERoxd
1
JKSTMM. TM f l|r An Tr TrS
It is easy to check that f is also a solution to
1 ||| - lite l1
inf JEHs Bflg 2llfRll
Since we have [1 |1xx S 1 by Assumption 2, it holds
inf _| 1 JEHs n HH
inf
|| - fxx l1
JEHo
11lg S1
1 11 N
inf FERex n
1
IKSTTM - I + Xn (KS TMTT
As consequence, we have
An
E. - E. En.lrn S V
inf Pe
|KSTTM - SAA ll An, TKSTTTMMT S) +
2
Finally, by Lemma 3 and Equation (28), we obtain the result stated.
Furthermore, we give the proof of the second claim, i.c. the excess risk bound for kernel ridge multi-output regression.
Proof. We now assume that the outputs are bounded, hence, without loss of generality, y C B (Rd). First, we prove Lipschito-contimiity of the square loss under Assumptions and 4. Let g 1222 Hx (X) z ||= v11
24

---

Published in Transactions on Machine Learning Research (09/2023)
[=f(2)]2 1, for some f He and X. By We have that Vg(z) 2 y, and hence [VV( Assumptions 2 and 4 and Cauchy-Schward inequality, it is easy to check that
(1) < << [1l ( 1/2
which gives that (x))l /// o and then [|Vo(=)]|2 - 4/2 M11pp + 1. We finally obtain that
()) v) (f (z') 3)1 (x2/22 || || + f (2')l2
We can then obtain the same generalisation bounds as above. Finally, looking at the approximation term,
En [/.]
En
[mmm]]
fx
112
|| - Y112
2n
2n
- 2n
fx
fir.
12
1
1|2 An
<
inf
||
l 2
/EMs
2m
1/|N,51 $1
$ inf 7E82 n
IKST 7- --itt 11 + An T
TSKT
An 7+
2
C² + 1 An C
Here again, as in second claim of Theorem 2, we can directly use bound (19) and then, in combination with (28), we obtain the stated second claim in Theorem
Finally, we here prove Lemma
Proof. Let f Hx such that llf S and (f(x1), ....,, f(nn))) Rnd We define the lincar
operator 5x Hx + and such that (x) ((( f(r)) for all f Hx. Then for all f E Hx
and z (al 2,)) E Rdd we have
(8x(f), (f(z))
(1,Kn
Kin f,
= -ffSS(
i=1
1=1
Hx
Hence
2
=- (
- -(5x (K ( .s) HK
Rnd
We recall the eigendecompossitions K and M
n
K U (nD) UT
M VAVT Σ p
25

---

Published in Transactions on Machine Learning Research (09/2023)
Then,
n K M
( Emogur
Tilii (u.u,) @ (1
thll (u,) 2 (v²) (u, @ (v))
1444 (u,) 2 (0,) ((u) 8 (vy)))
and for all 1 S i, n and 1 - j.J 2 d, if (i,t)) (i.5), then (u,) & (ry) ((u)) @ (uy.) = 0 and otherwise (u+) & (uj)² ((ux) @ (vy)) Then, this allows to show that the operator norm of Kronecker product is the product of the operator norms, and that
(K & M) = Σ
(mpane) (u.) 2 (1z) ((()) @ (v)))²
(29)
i= j=
We define, for all i n and 1 S
1=1 Pi
ui Kmn Vz-
(30)
By definition, span (ssn)ssss.))sss) C span mal,
(KCis U/ ) ")isiiiniis and we show that the iPys are orthonor-
(pig. Pvy))
Ui Kaz Uj*
mup Hy Σ un Key U2
i=1
Hx
Σ
Us miy (Kur, M8g Kan 7exc
Mmpstti Aapepy
1.1'
1
1
Us uy, (0z, Kas )es
Vapals Vapety
12
Us, ue, k( my) (uy, Muy
Smpiiy Vapepy Hy
Un uz, ((01) (155 2y)
Mupllty Papppy -0 if jt1
1.2"
Otherwise, if 9',
(Pi5, /npli /
1
Σ
ua ue K( me
1.1'
1
1
(Kus, u. gm
Vnnu /nips
nap (us ust yg
0 if i
26

---

Published in Transactions on Machine Learning Research (09/2023)
Hence, (Pi). 0 (i,2) # (d.0") and if (i, (i,J),
(pig, =.
Finally, span (sss)isssaisssa) c span (KKs, 12) and dim ((s)))sssm nd dim (Kaa (isimnssss)) hence (pu)sssmassid yields orthonormal basis of span (Kx, V)ssssnasss) As consequence, all Hy can be decomposed as fi f2. with fi span (K2, (isisnss))))) and f2 E span (Ku: ()sisniss))))))) Thus, for all IRd, y can be written as y Σ 2==1 U3 and
(5),
(f(z)), 2)
Bi (f(x:). uy)Rs
B4y (f. Kz, Hy)
Big (fi, K2, 25)1 +
(fe. Ka, U5) He 23
5=1
(fi, Kr, Koi
1=1
(8x((1)) (5),
Hence, let f He such that ||5 S 1, written as f= Eis E fospis f with fiz E R for all ssisn n and 1 j S d and such that Σ Σj= Pij S and f E span (Kox, (Ka,)))ssssassss)) such that | I 1 (since Ilfl| = Ei= E/ f2 + |lf ll $ 1, we have that
d 5x(f) = Σ Σ
Si (1)1)),
=1
and, for all 151 S n,
(i)(x))
Jmpats Σ
Wip klar, xz) My
P=1
v
Má KT. U;
ni,
and then
Sx (1p1)) = v
Pj (Ku;) & = Vmpaniu & uj. nju
(31)
Finally,
12
8x(f) - ΣV Σ
fis (mpin us @ U/
(32)
27

---

Published in Transactions on Machine Learning Research (09/2023)
Besides,
d
(K& Sx(1)- =
(mphpt))) (u;) @ (1z) ((u,) @ (0;))
fey (mpuuy)) us Buy
fü (muint))) 112 us 8 Uz.
2= j=
Then,
Sx (K () M) -)-1sx(11)) = Σ
fo (mpp)1s) 11/2 5xx (us 1g)
i=1 1=1
fog (mppsts)) 1/2 Σ Knn U5
and finally,
(5x Sx(1) 1) =
Σ
fyfey (mpmpps) 1/2 us (Kar uy sprr)) me
Σ
faler
(mmm)) Σ
ui Kni 9u8
W-=1j if'=1
1=1
Σ
Σ
Jufer (eigs
W-=1ド==
1=1
Sxx (KK (x()))
Sx(1)) .))
1.
Thus, we do have the ellipse constraint
(K-1/2 -12 2 M 1112)
z112 5 1.
(33)
A.4 Proof of Theorem 5
A.4.1 First claim K-satisfiability
Let us now prove the first claim (lh.s. of equation (4)) of the K-satisfiability for p-SR and p-SG sketches. It. is articulated around the following two lemmas.
Lemma 4. Let M Idxd be a symmetric matrix, 5 E (0, 1), and C. be an E-cover of Bd. Then we have
1
|| 5
sup vec.
(0, Mo)|
1 1 25 x²
Proof. Let M, and C. in Lemma Let a Bd. By definition, there exist C. and w Bd such that v + We thus have
(u, Mu) Mu) + 2e(u, Mw) (w0, Mw),
(34)
28

---

Published in Transactions on Machine Learning Research (09/2023)
Taking the supremum both sides of (34) we obtain
sup B
(u, Muyl
sup
Muyi + 2e| Mwyl + + kw, Mw)l)
veC,, WEB
I sup (vv, Mo)l + 25 vec,
sup
(1, Mw)) [)] + E²
sup wEB
(ux,
Mw)l
wEB
S sup vEc,
((v, Mb)1 + 25
sup
1(v', Mw> M
d'EB",WEE
sup ((v, Mo)) + (2e + 22) Mllop: vEC,
or again
|M - 1 22 e²
sup
Me)].
wEC,
Lemma 5. Let S E Raxu be p-SR or p-SG sketch. Let v E B", then for every t >> 0, we have
p 1 - 1|012
41 < 2
Proof. The proof of Lemma 5 is largely adapted from the proof of Theorem 2.13 in Boucheron et al. (2013).
Let S E Rexn be a p-SR or a p-SG sketch, and v E B It is easy to check that for all i I 8 we have
1 |oll2 such that E [Su]]] a
||Su||3 - |101 = Σ [So]]
||
The proof Uhen consists in applying Bernstein's inequality (Boucheron et al., 2013, Theorem 2.10) to the random variables [Su]]- We now have to find some constants v and c such that Σi= E [[So]] I v and
[Su
20
d! 2 I
vucr2 for all 2 3.
From (12) and (13), it is casy to chock that the Sij are independent and 1/(sp) sub-Gaussian. Then, for all R, we have
n
E exp (A [Su,,)]
E
exp A Σ Sijty)
E
exp
(ASt)ug)]
J=1
X²
< exp
|001|
2sp
X²
A exp
2sp
The random variable [Sv], is therefore 1/(sp) sub- nab-Gaassian, and Theorem 2.1 from Boucheron et al. (2013) yields that for every integer 2 it holds
[So]
2q
<
g!
4
q!
4
9
2
sp
sp
29

---

Published in Transactions on Machine Learning Research (09/2023)
Choosing 9 2, we obtain
Σ IE [SSU]]]
Σ
4
16
sp
sp2
such that we can choose / = 1sp("pp)) and c == 4/(s)) Applying Theorem 2.10 from Boucheron et al. (2013) to the random variables [St]] finally gives that for any holds
4
2t
4t
P ||
p
I
2c_
sp
Proof of the first claim of the K-satisfiability. Let K R be Gram matrix, and R be a p-SR or p-SG sketch. Recall that we want to prove that there exists co such that
P |Io S SU2 - lopp 2 2e 2c
where K/n = UDUT is the SVD of K; and U1 E Bxdd contains the left part of U. Let & E (0, and C. = (u'....,,"" Vn be an E-cover of Bdd, We know that such a covering exists with cardinality Ns. I (1 + 2) 2)dr,
see c.g., Matoušek (2013). Let Q = U.S S SU1 Id., applying Lemma 4, we have
P 111lllop
1 2
1 28 g²
PP
sup iSA.
(10 Qu4)
2
1 1 28 - E2
Σ
Ku's Qu") >
2
S.
{
1 2e -2
||S
-010-1 ||
>
(35)
where w' = -U E EB. Now, by Lemma 5, for any w B",
we have
4
21
4/
P
||S 1
+
2e
p
sp
Let 2 32t/(a²p²), for some a $ 1. Then, we have
21
4t
<
a
g²p
S 2a, and therefore
8
8
PP [|]|[ - ||
[
20
$ 25 -
If we take a = -1-22- -2)/4 we obtain
P ||Su|| - || 1- 52
S 2e
as long as 2
p²(1 2 2y Now, let t= (F(1-2--67) + log(Ne), and 8 > 1024 log(1+2/s)
512
p²(1 (1112c-c2jaa dz. We do have
512t
512
2 +
log(NV)
+ =
p1 2= - 22)2
p² 2 222
2
such that
2
2
2e -0
P [|S -||w||3
< 26 -t
N
where co
22(1-25-e2) 1024
Plugging this result into (35), we get that as soon as s 2 1024 log! 1+2/s)
221- d, it holds
1 1 s 22
P op
Finally, we can tune optimize the lower bound on s. If we take 0.1, we obtain 2 5520dd.pp' and Co 2 /2/2500.
30

---

Published in Transactions on Machine Learning Research (09/2023)
A.4.2 Second claim of K-satisfiability
We now turn to the proof of the second claim (r.h.s. of equation (4)) of the K-satisfiability for p-SR and p-SG sketches. It builds upon the following two intermediate results, about the concentration of Lipschitz functions of Rademacher or Gaussian random variables.
Lemma 6. Let > 0, and let XI...., Xn be independent real random variables with |Xx & K for all isiSi n. Let F R + R be L-Lipschitz convex function. Then, there exist C, > 0 such that for any X one has
(MMXX)EFFC -XF(X))FF(X)22 KX} S C" exp (-dX2/13)-
Lemma 7. Let XI,..., X, be i.i.d. standard Gaussion random variables. Let F R R be L-TTipschtzz function. Then, there erist C, > 0 such that for any X one has
P(F(() - BF(X) L} S C' exp exp (-dX2/13).
The above two lemmas are taken from Tao (2012), see Theorems 2.1.12 and 2.1.13 therein, but are actually well known results in the literature. In particular, Lemma 6 is adaptated from Talagrand's inequality (Talagrand, 1995), while Lemma 7 is stated as Theorem 5.6 in Boucheron et al. (2013), with explicit constants. We however choose the writing by Tao (2012) in order be consistent with the Rademacher case.
Remark 3. Note that thanks to Lemma 6, we are even able to prove K K-satisfiability for any sketch matrix S whose entries are i.i.d. centered and reduced bounded random variables.
Proof of the second claim of the K-satisfiability. Let K E R"*" be Gram matrix, and S E Rax be p-SR or p-SG sketch. Recall that we want to prove that there exist positive constants c, C,, > 0 such that
|SDDD
2/1 llop > cōn S cie c28
where K/n UDUT is the SVD of K, U2 Rux(n-da) is the right part of U, and D2 (R(--d))((ndd))) is the right bottom part of D. Note that we have SU2D, = SU 11,2, where 5 = diag (0d.. D2) E RR" Following Yang et al. (2017), we have
SUD 1lop = sup
(u,, NEB*.
where & = E E : Ew E 8m-1, v = U.) Now, let u!.... N be 1/2-cover of B We know that such a covering exists with cardinality N < 5". We then have
SUD
sup w WER
SuYl
op
max sup (u', Su) vE
+
sup
Ku,
SS
2 wEB vEE vEE
max
sup S 3E
Sv)
SU D1/2
op
and rearranging implies that
SUD!/2
< 2 max sup iS 33R
(ut
Hence, for every > we have
P
(|SSnDD!
> cdn
SP
max sup SSN oEt
c (ut,
<
sup JEE
(u',Sii
_.
2
(36)
31

---

Published in Transactions on Machine Learning Research (09/2023)
Now, recall that
1 B BoR,
S
Vsp
where B E R is filled with i.id. Bernoulli random variables with parameter R E R is filled with i..dd Rademacher or Gaussian random variables for y-SR and p-SG sketches respectively, and o denotes the Hadamard (termwise) matrix product. The next step of the proof consists in controlling the right-hand side of (36) by showing that, conditionally on B, we have Lipschitz functions of Rademacher or Gaussian random variables, whose deviations can be bounded using Lemmas 6 and Therefore, from now on we assume B be fixed, and only consider the randomness with respect. to R. Let B*, and define F: Raxe R
F(R)-
sup D3E
(u. (B o (B+R)o)
sp
It is direct to check that F is a convex function. Morcover, we have
F(H = sup (u, (B o (B)))
wEē
Ku, (B R)UD)?3))
sup
vESn-1
( (B o R) T u)
sup
vESA-1
D ( R) u
Thus, for any R, R' we have
sp (R) - F (R') =
11/20 (B R) RTTU|
1/2
D1/2UT (B (R-R))
D1/2 L opp
||B (R RR) lop ||u||2
< 6n || B o (R 1R-HR RF)
6n R R'I
(37)
such that is (Sppsittt Moreover, we have
pp [ E
2 (Bo R) R)Tn| m1
S v E Iu (B o R2U (B R
uguu [[ o RUUDDDUJ (B R)T] kx
A,
uuu (UDDU Jur [[B Rl o R[wy]
A
11==1
A=1
Tr((2))
which implies
E [R() S π
Σ n
Σ d M
1
min(uy, 02) S
82
(38)
n
+1
32

---

Published in Transactions on Machine Learning Research (09/2023)
where we have used the definition of 62 and the assumption 2 8g Coming back to (36), obtain
p P|IISUDD?? olop cdn S E
P
sup NEE
(u, Sv) 8n B
5m }
55' E
F(R)
c
S 5
P
F(R) -E(F(R))
8n
(39)
2
vp
CC5 cxp
-c'
c
1
2
8p
(40)
2
vD
1
I Cexp
2 2
1
2
log(5) 8
where (39) comes from the upper bound on [R)) we derived in (38), and (40) derives from Lemmas and 7 applicd to the function whose Lipschitz constant has been established (37). Thercfore, taking
vp 1 + (lg (5) 1, we have
P 21|SU2D2/2110 olop cdn S GiE cie
with = O and C2 = d' (Vplog (5) + +
A.5 Proof of Proposition
We prove Proposition 1 thanks to duality properties.
Proof. Since problems (3) and (11) are convex problems under Slater's constraints, strong duality holds and will show that they admit the same dual problem
min CER
[ + 2A, 2A S- KST (SKS
(53)
First, we compute dual problem of (3), that can be rewritten
min Σ (w:) + Ann 7 TKS
n
(u:)
PER",UER 1=1
8.t. u = KSTT. 7.
Therefore the Lagrangian writes
Σ (w:) Ann SKST 7+ Σ Ciuu [KST
C(7, u),
Ann T SKST +
Σ l (m:)
Σ Cru S KS 7.
1=1
Differentiating with respect to and using the definition of the Cenchel-Legende transform, one gets
s(5) = inf "ER"
u(),u))
inf {ls (w:) + +Gui++ WER
+ inf TER
Ann 7
T - -
-lg (-Gi) 5 (SKSTTSSSKKK
2An
1=
33

---

Published in Transactions on Machine Learning Research (09/2023)
together with the equality SKST xn SKC, implying 2. (SKST))KK? where Ç E R² is the solution of the following dual problem
1 B KS KS (SKS
min BER
E (-
(53)
2
Then, compute dual problem of (11), that can be rewritten
An 1
E (Us, Uk) +
min WER" MER"
i=1
s.t. u = KST K.w.
Therefore the Lagrangian writes
Σ 4 (ui) + An
C(w, 5) =
1
2 +
G(u [KST Kiwli)
2
1=1
Σ 4 (u.) +
Am
2
Σ
c u, - w R. SKC.
i=1
Differentiating with respect to and using the definition of the Fenchel-Legendee transform, one gets
s(S)
inf "EER"
L(w, 5)
Σ inf {ce (u;) + Gui3 inf wERr
Ay 2
| w TT-1122 SKC}.
HER
We have that
(1ww113) = 2w
aw
0
(w K. SKC) = K,T SKC,
Ow
Then, setting the gradient to zero, we obtain
R,T SK?.
(41)
Ann
Hence, putting it into the Lagrangian,
KST (SKST)* SKG.
Ann
Ann
and
2A, S KS K,K,J SKÇ = 2A,
KST (SKST)² SKG.
Hence, ? E R is the solution the following dual problem
min
(-3++ + 2A,
KST KST (
(53)
BER
1=1
Finally, since both problems are convex and strong duality holds, we obtain through KKT conditions
(SKST) 1/2 -1/2" (SKST) (B./20 VTg
and
1 7ER* π 1=1
min
Σ A(KST 7/. v1) + An, T TSSS > min
Σ l (w zg (I.).
.1) + An
2 1wmm
2
DER!
1=l
34

---

Published in Transactions on Machine Learning Research (09/2023)
B On relaxing Assumption 2
In this section, we detail the discussion about relaxing Assumption 2, i.c. the restriction of the hypothesis set to the unit ball of the RKHS. Assumption is a classical assumption in kernel literature to apply generalisation bounds based on Rademacher complexity of a bounded ball of a RKHS. Moreover, it also useful in our case derive an approximation error bound, describing how K-satisfiability of sketch matrix allows to obtain good approximation of the minimiser of the risk. However, let discuss the consequences of relaxing this assumption. Indeed, all need is a bound on the norm of the estimators fs minimiser of the regularised ERM sketched problem and fi minimiser of the regularised ERM sketched denoised KRR problem. By dcfinition, noting Hs = (f = - [s (,, I 1 E RR we have that
f. arg min
Σ ((f(z)), Mi) + 2 1l n
JEHs
Hence,
An 2
1
(f.())))))) + An 1
1 Σ n
((0.))) $
n
if we assume that gaxigien 8(0, y:) $ to simplify the derivations. As consequence, we obtain that
|| 1 v Am
(42)
Similarly, we have that
1
|| - - 11 + An 1lf
= = arg min fEMs
that gives
1/2
Me |
S
11 llzee
Ann
By Assumptions 2 and 4,
1
[l 11
1
Σ
n
(Jren do (1 (1.2)) Hhx
7
n
1=1
I1 li & (x))
n
K,
and finally
11 He S
V
K Xm
(43)
Remark 4. Note that in the multiple output settings, we obtain
K TT(M) An
1 S
(44)
We are now equipped to derive the generalisation error bound IE [Cr.] E, [e..] and the approximation error bound E,, [ E,, [emma] We first focus the generalisation bound, and following the proof given in Appendix A.1 and given the new norm upper bound V Am. for any 8 (0. 1), with probability 1 1- 18// we have that
4L 2K
(4/6)
E f En
[eg]]
(45)
35

---

Published in Transactions on Machine Learning Research (09/2023)
This dependence in 1/VA, shows that, as expected by regularisation penalty, with fixed n, when An increases, the generalisation bound decreases and then we obtain a better generalisation performance. However, this behaviour does not reflect completely the role of An, since there exists a tradeoff between overfitting and underfitting, and then should not be set too large. We now focus on the approximation error bound. As in Appendix A.1, we oblain that
a
(f.(t.). 1
n n
En
(3.
En
(46)
n
inf JEHe I1 xy Sllfall
Σ (f((z))
==1
2=1
L n
inf JEHg 11l sllf-llng
If (ax) fuu (++)
S L
inf SEHs 11ll, sllfall
Σ 1f (xt) fHee (x.))²
=L
inf JEHs
1
|| fta 11
(47)
n
where, for any ʃ E Hu. f = (f(z))..., [(xn)) E ER. Let jth = Ei= [sT k (-,2;,,, where 5 is a solution to
inf yER
+ SKS T
It is easy to check that fR is also a solution to
inf JEHs 1/ <|14 1l
| f 1
(48)
Now, comparing (47) and (48), done in Appendix A.1, essentially boils down to comparing || =1, and || Il which is highly nontrivial question. In particular, the upper bounds (42) and (43) are not informative enough. Another solution could consist in adding -/ || to (46). However, the upper bound (42) then transforms this term into a constant bias. This can be explained (42) is very crude. Instead, having || bounded by A for a 2 -1/2 would be enough to exhibit bias term that vanishes An goes to O. Note that it would still degrade the tradeoff with the genealisation term. Hence, if generalization errors can be dealt with when removing Assumption 2, it is much more complex to control (46).
C Some background on statistical properties when using a p-sparsified sketch
In this section, we focus on p-sparsified sketches, and give, according to different standard choices of kernels (Gaussian, polynomial and first-order Sobolev), the different learning rates obtained for the excess risk as well as the condition on
We first derive the following corollaries for the excess risk of p-sparsified sketched estimator in both single and multiple output settings.
Corollary 2. For p-sparsified sketch matrix S with 2 max (Cudn /p2 (2,n), we have with probability greater than Cie-scly in the single oulput setting for generic Lipschitz loss,
E [1p.]] S R [esnn] + ILC An +
Xn
K + 8T
+
(V
(49)
2
n
36

---

Published in Transactions on Machine Learning Research (09/2023)
and if (2, (z y)² 12 and V [0, 11, with probability at least - ô we have that
B [ B [esra] + (c²+ A. C +
+ vi
Vn + +0 +0 (VV))
(50)
Corollary 3. For a p-sparsified sketch matrix S with max (Codn/p² 8gn), we have with probobility greater than Cie sc(pp) in the multiple output setting for a generic Lipschitz loss,
E [ S E [6 LCVA An + |M 82 + An
K Ir (M)
(V=))
3
+II
(51)
n
and (2, ||= v12 12 and Vc c 1 (R*), with probability at least
6 we have that
E [ S [eee + (
+2 (n 8n + 8TD (M) & *| 1/2
up + 1/2
(V
(52)
Vn
We summarize in Table 3 the different behaviours of 8 and d, in the different spectrum regimes considered, in order to explicit the exact condition on 8 each case, More specifically, for a Dth-order polynomial kernel, d,, for any n is at most D 1, leading s of order D 1 to be sufficient. Finally, we can derive the learning rate obtained as well as the exact condition on s for each scenario, see Table 3. Compared with Random Fourier Features (Li et al., 2021), we that obtain slightly degraded learning rates for Gaussian and first-order Sobolev kernels, in comparison with the O (1/vTi) rale the authors obtain. Our rates remain however very close.
Table 3: Statistical dimension, lower bound obtained on 8, and learning rate obtained for excess risk with p-sparsified sketches for different kernels.
Kernel
dn
Learning rate
Gaussian
Ving(n)
α v Iog(n)
S (log(nn//²)
(log(n)) 1/4 1//2
Polynomial
x
Ω (
Vn
Sobolev
xc 1/3
Ω (n!13/p²)
n//3
1/3
D Detailed algorithm of the generation and the decomposition of p-sparsified sketch
In this section, we detail the process of generating a p-sparsified sketch and decomposing it as a product of sub-Gaussian sketch Ssc and a sub-Sampling sketch Sgs.
Algorithm Generation of a p-sparsified sketch
input:s, n and p
Generate a 8 x matrix B whose entries are i.i.d. Bernouilli random variables of parameter p-
indices indices of non-null columns of B.
B' B where all mull columns have been deleted.
Generate matrix MsG of the same size as B' whose entries are either i.i.d. Gaussian or Rademacher random variables.
SsG Msg B', where denotes the component-wir wise Hadamard matrix product.
return Ssc and indices
37

---

Published in Transactions on Machine Learning Research (09/2023)
E Some background on complexity for single and multiple output regression
In this section, we detail the complexity in time and space of various matrix operations and iterations of stochastic subgradient descent in both single and multiple output settings for various sketching types.
We first recall the time and space complexities for elementary matrix products. The main advantage of using Sub-Sampling matrices is that computing SK is equivalent to sampling training inputs and construct a Gram matrix, hence gain huge time complexity since we do not compute matrix multiplication, as well as space complexity since we do not compute n x n Gram matrix. As consequence, the main advantage of our p-sparsified sketches is their ability to be decomposed as S = Ssa where Ssc E Raxs' is a sparse sub- Gaussian sketch and Sss E IR' x is a sub-Sampling sketch, as explained in Section 3 and Appendix D. This decomposition trick is particularly interesting when p is small, and since s' follows a Binomial distribution of parameters n and 1 (1 - p)* and we assume our settings that n large, hence we have that s' E [s'] -= nsp. In the following, we take s' = nsp. We recall that Accumulation matrices
from Chen & Yang (2021a) writes S = Ex= Sty), where the Siys are sub-sampling matrices whose each row is multiplied by independent Rademacher variables. Hence, both p-sparsified and Accumulation sketches are interesting since it completely benefits from computational efficiency of sub-sampling matrices. See table for complexity analysis of matrix multiplications.
Going into the complexity of the learning algorithms,, the main difference between single and multiple output settings are the computation of feature maps, relying on the construction of SKST and the computation of the square root of its pseudo-inverse for the single output setting which is not present in the multiple output settings. We assume in our framework that d and even d² are typically very small in comparison with n. Hence, we have that the complexity in the single output case is dominated by the complexity of the operation SKST, whereas in the multiple output case it is dominated by the complexity of the operation SK. We see that from a time complexity perspective, p-sparsified sketches outperform Accumulation sketches in single output scttings as long as p S m/n//s, and in multiple output scttings as long as p S m/ns. From space complexity perspective, Accunulation is always belter nsp is lypically greater than s; otherwise it shows poor performance. However, p is usually chosen such that nsp is not very large compared with
Table 4: Complexity of matrix operations for each skelching lype.
Sketching type
Complexity type
SK
SKST
time
0 (n2s)
O (n²s)
Gaussian
O (n²)
O (n²)
space
time
O (n²s²p
0 (n2
p-sparsified
(n2sp)
O pp))
space
time
O (nsm)
ms²
Accumulation
space
(ns)
0 (s
time
0 (n²)
0 (n²)
Skent
O (n²)
0 (n²)
space
time
O (ns)
0 (s²)
Sub-Sampling
space
O (ns)
0 (s²
F Discussion with dual implementation
In this section we detail the discussion about duality of kernel machines when using sketching. A first idea consists in computing the dual problem to the sketched problem (3). It writes
min CER
Σ c (-G++
2A, ST KST (SKST)
(53)
38

---

Published in Transactions on Machine Learning Research (09/2023)
where & E(;, (k), and f denotes the Fenchel-Legendre transform of such that sup f(z). First note that sketching with a subsampling matrix in the primal is thus equivalent to using Nyström approximation in the dual. This remark generalizes for any loss function the observation made in Yang et al. (2017) for the kernel Ridge regression. However, although the & might be casier to optimize, solving (53) seems not a meaningful option, as duality brought us back to optimizing over IR", what we initially intended to avoid. The natural alternative thus appears to use duality first, and then sketching. The resulting problem writes
min BER
Σ E (-(sT
0 SKS 0.
(54)
2A,
It is interesting to note that (54) is also the sketched version of Problem (53), which we recall is itself the dual to the sketched primal problem. Hence, sketching in the dual can be seen as a double approximation. As consequence, the objective value reached by minimizing (54) is always larger than that achieved by minimizing (3), and theoretical guarantees for such an approach are likely to be harder to obtain. Another limitation of (54) regards the (- Indeod, these terms generally contain the non-differentiable part of the objective function (for the e-insensitive Ridge regression instance), and are usually minimized by proximal gradient descent. However, using similar approach for (54) is impossible, since the proximal operator of (ST.) .) is only computable if S S = In, which is never the case. Instead, one may use a primal-dual algorithm (Chambolle et al., 2018; Vu, 2011: Condat, 2013), which solves the saddle-point optimization problem of the Lagrangian, but maintain dual variable in R Coordinate descent versions of such algorithms (Fercoq & Bianchi, 2019; Alacaoglu et al., 2020) may also be considered, as they leverage the possible spansity of S to reduce the per-iteration cost. In order to converge, these algorithms however require number of iteration that is of the order of n, making them hardly relevant in the large scale sctting we consider.
have (0.) 1 0|01 + (0, v) + + <||0[]= for
For all the reasons listed above, we thus believe that minimizing (53) or (54) is not theoretically relevant nor computationally attractive, and that running stochastic (sub-)gradien) descent on the primal problem, as detailed at the beginning of the section, the bost way to proceed algorithmically despite the possibly more elegant dual formulations. Finally, we highlight that although the condition ST S L, is almost surely not verified (we have S *x* with < n), still have B[ST $) I, for most sketching matrices. An interesting research direction could thus consist in running a proximal gradient descent assuming that S S = In, and controlling the error incurred by such an approximation.
G Examples of Lipschitz-contimuous losses
Robust losses for multiple output regression: For p 2 1, 1 R", and for all v.v E V, Ely,v')
g y'), where g is:
For Huber: For k > 0:
{
1 Hulb k (llvll> - 5)
if Iall S * otherwise
Vy E V,g(v))
-
For e-SVR (i.e. einsensitive & loss): For >>00 O:
= ||vlly - if |l 2
Vy E V.g(y)
otherwise
The pinball loss (Koenker, 2005) for joint quantile regression: For quantile levels, 71 < T2 < ... < Ta with 5 (0, 1), we define:
()(()),)) = Lf((xx) ylla),
39

---

Published in Transactions on Machine Learning Research (09/2023)
with the following definition for L, the extension of pinball loss to R (Sangnier et al., 2016):
For E :
TJTJ (T)
if 2 0, 1)r, ry < 0.
(r(r)=
H Additional experiments
In this section, bring some additional experiments and details.
H.1 Simulated dataset for single output regression
First, we report the plots obtained with k-Huber for p-SG sketches (see Figure 2) and note that we observe a behaviour similar lo p-SR sketches when varying p and in comparison to other lypes of sketching and RFFs. However, that the MSE obtained is slightly worse than p-SR sketches. An explanation might be that, in very sparse regime, i.c. very low p, p-SG sketch is too different than Gaussian sketch, making it lose some good statistical properties of Gaussian sketches. We however observe that the larger is, the smaller is the statistical performance between y-SR and p-SG sketches.
We then report in the following the corresponding plots obtained with e-VRR, that witnesses the same phenomenon observed earlier with H-Huber about the interpolation between Nyström method and Gaussian sketching while varying the probability of being different than 0 p, with p-SR sketches (see Figure 3) and p-SG sketches (see Figure 4).
H.2 Multi-Output Regression on real datasets
We here first brief presentation on River Flow and Supply Chain Management:
River Flow datasets aim predicting the river network flows for 48 hours in the future specific locations. These locations are 8 sites in the Mississippi River network in the United States and were obtained from the US National Weather Service. Dstaset rf2 extends rf1 since it contains additional precipitation forecast information for each of the 8 sites.
2.
The datasets semld and sm20d come from the Trading Agent Competition in Supply Chain Manage- ment (TAC SCM) tournament from 2010. More details about data preprocessing can be found in Groves & Gini (2015). The dataset contains prices of products at specific days, and the task is to predict to the next day mean price (semld) or mean price for 20-days in the future (scm20d) for each product in the simulation.
To conduct these experiments, the train-test splits used are the ones available at http://mulan, sourceforge. net/ datasets-mttr. html. Besides, we used multi-output Kernel Ridge Regression framework and an input Gaussian kernel and an operator M Ig. We selected regularisation parameter Xn and bandwidth of kernel a² via 5-fold cross-validation. Results are averages over 30 replicates for sketched models.
Table 5: Numbers of training samples (ng)), test samples (me)), features (4) and targets (d).
Dataset
Ner
Ric
d
Boston
354
152
13
5
otoliths
3780
165
4096
5
rfl
4108
5017
64
8
rf2
4108
5017
576
8
semld
8145
1658
280
16
scm20d
7463
1503
16
40

---

Published in Transactions on Machine Learning Research (09/2023)
Nystrom
Nystrom
0.13
0.12
P-5G (p-10/m)
P-SG (p=30/m)
P-SG (p-5/n
Gaussian
(p-50m)
Training time (in seconds)
8.5
P-5G (p-10m)
P-5G (0= 30m)
P-5G (p-5/n)
Gaussian
PG (p=50m)
8.0
0.1
7.0
0.08
6.5
0.00
5.5
0.04
40.0
60.0
80.0
100.
120.
140.
(a) Test relative MSE w.r.t. 8 with x- x-Huber.
(b) Training time (seconds) w.r.t. 8 with A- H-Huuer.
Nystrom Gaussian
0.2
RFF
0.17
p-SG (p=30/n)
Accumulation (m=30)
0.15
CountSketch
0.12 0.12
10
0.07
0.05
5.0
7.0
8.0
9.0
Training time (in seconds)
(c) Test relative MSE w.r.t. training times with k-Huber.
Figure 2: Trade-off between Accuracy and Efficiency for p-SG sketches with k-Huber loss on synthetic dataset.
We compare our non-sketched framework with the sketched one, and we furthermore compare our p-sparsified sketches with Accumulation sketch from Chen Yang (2021a) and CountSketch from Clarkson & Woodruff (2017). Moreover, we report the range of results obtained by SOTA methods available at Spyromitros- sirromitrossXooufis et al. (2016). All results in terms of Test RRMSE are reported Table we see that our p-sparsífied sketches allow ally statistical and computational performance, since we maintain an accuracy of the same order without sketching, and these sketches outperform Accumulation in terms of training times (see Table 2). Tn comparison to SOTA, our framework does not compete with the best results obtained in Spyromitros-Xionn Xioufis et al. (2016), but almost always remains within the range of results obtained with SOTA methods.
41

---

Published in Transactions on Machine Learning Research (09/2023)
Table 6: Test RRMSE and ARRMSE for different methods on the MTR datasels. For decomposable kernel-based models, loss here is square loss and s = 100 when performing Sketching.
Dalaset
Targets
w/o Sketch
20/w.-SSR
SS/ 0.583 + 0.003
Acc. m 20
CountSketch
SOTA
Mean
0.575
0.584 + 0.003
0.592 + 0.001
0.575 + 0.0005
[0.091, 0.983
chsi2
0.351
0.356 + 0.005
0.357 + 0.004
0.361 + 0.002
0.350 + 0.002
[0.033, 0.797
nasi2
1.085
1.124 = 0.003
1.124 + 0.003
1.082 = 0.0004
1.110 = 0.0003
[0.376, 1.951
eadm7
0.388
0.397 + 0.004
0.398 + 0.003
0.387 +: 0.001
0.387 + 0.001
[0.039. 1.019
rfl
sclm?
0.660
0.659 + 0.008
0.661 0.005
0.681 1: 0.002
0.648 0.002
[0.047, 1.503]
clkm?
0.283
0.281 + 0.001
0.282 + 0.001
0.293 + 0.0005
0.281 = 0.0004
[0.031, 0.587
vali2
0.614
0.633 + 0.008
0.635 0.010
0.656 + 0.003
0.611 0.003
[0.037, 0.571
napm?
0.593
0.628 + 0.020
0.614 + 0.016
0,627 + 0.003
0.601 + 0.003
[0.038, 0.909)
dldi4
0.629
0.597 0.004
0.597 0.003
0.646 + 0.001
0.614 0.002
[0.029, 0.534
Mean
0.578
0.671 = 0.009
0.656 0.006
0.796 = 0.006
0.715 0.011
[0.095, 1.103
chsi2
0.318
0.382 + 0.016
0.358 - 0.010
0.478 + 0.006
0.426 + 0.013
00.03,, 0.737
nasi2
1.099
1.084 = 0.005
1.092 + 0.006
1.018 = 0.003
1.036 = 0.002
00.384, 3.143
eadm7
0.342
0.390 1 0.013
0.369 + 0.007
0.456 1: 0.004
0.417 1 0.010
[0.040, 0.737
rf2
sclm?
0.610
0.719 + 0.030
0.672 0.021
0.948 1 0.014
0.852 0.030
[0.049, 0.970
clkm7
0.311
0.328 + 0.009
0.330 + 0.009
0.614 + 0.005
0.436 + 0.006
0.044, 0.891
vali2
0.712
0.960 + 0.044
0.894 0.043
0.890 + 0.017
0.939 0.028
[0.047, 0.956
napm?
0.589
0,812 + 0.014
0.831 + 0.017
1,100 + 0.023
0.856 = 0.032
[0.039, 0.617]
dldi4
0.646
0.696 0.010
0.701 0.011
0.855 = 0.004
0.761 0.007
[0.032, 0.770
Mean
0.418
0.422 = 0.002
0.423 0.001
0.423 = 0.001
0.420 0.001
00.330, 0.457
lbl
0.358
0.365 + 0.003
0.364 = 0.002
0.367 = 0.001
0.363 = 0.001
0..294, 0.409
mtlp2
0.352
0.360 = 0.003
0.362 + 0.003
0.362 = 0.001
0.358 = 0.001
00.008, 0.436
mttp3
0.409
0.419 + 0.003
0.416 + 0.002
0.417 1 0.001
0.416 1 0.002
[0.315, 0.442
mulp4
0.417
0.427 + 0.002
0.426 0.003
0.426 1 0.001
0.423 0.002
[0.325, 0.461
mtlp5
0.495
0.491 + 0.006
0.492 + 0.006
0.502 + 0.002
0.492 + 0.003
00.349, 0.530
mtlp6
0.534
0.524 + 0.008
0.527 0.006
0.537 + 0.002
0.527 0.002
[0.347, 0.540
mtlp7
0,531
0.519 + 0.008
0.523 + 0.006
0,534 + 0.002
0.523 + 0.003
00.338, 0.526
scmld
mtlp8
0.542
0.536 + 0.010
0.540 0.008
0.547 = 0.002
0.537 = 0.003
00.345, 0.504
mtlp9
0.385
0.395 + 0.003
0.395 + 0.002
0.390 +: 0.001
0.390 + 0.002
[0.323, 0.6
mtlp10
0.389
0.398 1 0.003
0.397 0.003
0.394 + 0.002
0.394 1 0.001
[0.339, 0.456
mtlp11
0.424
0.430 1 0.003
0.429 + 0.003
0.426 1 0.001
0.426 1 0.001
[0.327, 0.445
mtlp12
0.420
0.422 + 0.003
0.421 + 0.004
0.423 + 0.001
0.421 = 0.002
0..50, 0.466
mtlp13
0.349
0.358 + 0.004
0.354 + 0.004
0.351 + 0.001
0.351 + 0.001
0..32,, 0.419
mtlp14
0.347
0.364 + 0.004
0.363 + 0.003
0.350 1 0.001
0.355 + 0.002
00.356, 0.472
mtlp15
0.361
0.371 + 0.004
0.370 + 0.003
0.363 = 0.001
0.364 + 0.001
00.314, 0.406
mtlp16
0.376
0.382 + 0.003
0.384 + 0.003
0.376 + 0.001
0.378 + 0.001
00.322, 007
Mean
0.753
0.754 = 0.003
0.754 0.003
0.753 = 0.001
0.754 0.002
0..99,, 0.763
lbl
0.613
0.618 +: 0.002
0.618 1 0.002
0.614 1: 0.001
0.613 F. 0.001
[0.356, 0.678
mtlp2a
0.628
0.635 1 0.002
0.634 0.003
0.632 1. 0.001
0.631 1 0.002
[0.352, 0.688
mtlp3a
0.603
0.608 + 0.002
0.608 + 0.003
0.607 + 0.001
0.605 + 0.002
00.663, 0.683
mtlp4a
0.635
0.645 + 0.002
0.645 = 0.003
0.644 + 0.001
0.638 = 0.002
00.374, 0.730
mtlp5a
0.974
0.977 + 0.008
0.977 + 0.007
0.978 + 0.003
0.975 + 0.006
00.413, 0.846
mtlp6a
0.981
0.986 = 0.009
0.992 = 0.008
1.002 = 0.004
0.989 = 0.008
(0.424, 0.843
mtlp7a
0.996
1.001 = 0.008
1.004 0.007
L.005 # 0.006
1.000 = 0.009
00.404, 0.833
scm20d
mtlp8a
0.995
0.997 + 0.010
0.997 + 0.011
1.008 + 0.005
0.994 + 0.005
[0.107, 0.851
mtlp9a
0.708
0.704 + 0.003
0.702 0.003
0.698 1: 0.001
0.705 0.002
[0.382, 0.737
mtlp10a
0.718
0.722 + 0.004
0.722 + 0.004
0.716 + 0.001
0.723 = 0.003
00.413, 0.753
mtlp1la
0.729
0.730 + 0.003
0.729 0.003
0.725 + 0.001
0.728 0.002
00.402, 0.769
mllp12a
0.720
0.718 + 0.004
0.717 + 0.004
0.712 1 0.002
0.716 + 0.003
[0.429, 0.787
mtlpp3a
0.711
0.703 1 0.005
0.699 0.004
0.697 1 0.001
0.705 1 0.003
[0.400, 0.751
mtlp14a
0.683
0.673 + 0.004
0.670 + 0.003
0.668 1 0.001
0.675 + 0.002
[0.411, 0.779
mtlp15a
0.684
0.674 + 0.004
0.671 0.004
0.666 + 0.001
0.678 + 0.002
[0.384, 0.727
mt
0.689
0.677 + 0.005
0.676 + 0.005
0.672 + 0.001
0.682 + 0.003
[0.386, 0.754
42

---

Published in Transactions on Machine Learning Research (09/202))
Nystrom
Gaussian
Nystrom
P-SR (p=5/))
Gaussian
P-SR [p=10/n]
PRSR (p=5/))
p-SR (p-30/n)
P.SR [p-10/n]
P.SR (p-50mm)
Training time (in seconds)
P-5R (0-50m)
P-5A (0- 30ml
8.0
0.12
MSE
0.08
6.0
0.06
5.0
40.0
60.0
80.0
100.
120.
140
40.0
60.0
80.0
100.
120.
140.
(a) Test relative MSE w.r.t. with eVR..
(b) Training time (seconds) w.r.t. 8 with eSSV..
Nystrom
Gaussian
0.18
RFF
p-SR (p=30/n)
0.16
Accumulation (m=30)
CountSketch
0.14
0.14
MSE
0.11
0.1
0.08
0.06
6.0
7.0
Training time (in seconds)
(c) Test relative MSE w.r.t. training times with e-SVR.
Figure 3: Trade-off between Accuracy and Efficiency for p-SR sketches with c-SVR on synthetic dataset.
43

---

Published in Transactions on Machine Learning Research (09/202))
Nystrom Gaussian
P-SG (p=10/n)
p-SG (p=5/n)
12.0
P-SG (p-30/n)
0.14
P-5G (p-50m)
Training seconds
0.12
Nystrom
Gaussian
MSE
p-SG (p=S/n)
P.SG (p=10mm)
p.SG (p-30m)
p-5G (p-50mm)
8.0
0.08
7.0
0.06
40.0
60.0
00.0
120
140
(a) Test relative MSE w.r.t. with eV....
(b) Training time (seconds) w.r.t. with e-SVR.
0.22
Nystrom
Gaussian
RFF
p-SG (p=30/n)
0.17
Accumulation (m=30)
CountSketch
MSC
0.12
0.07
0.05
6.0
7.0
8.0
10.0 Training time (in seconds)
9.0
11.0
12.0
(c) Test relative MSE w.r.t. training times with e-SVR.
Figure 4: Trade-off between Accuracy and Efficiency for p-SG sketches with e-SVR on synthetic dataset,
44


Reviews:
[{"review_id": "M1ghZkhUjC", "reviewer": "Reviewer_sWZq", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n1. The paper is well-organized and easy to read. \n2. The risk bound is derived in a theoretical way.\n3. The performance is evaluated with both synthetic and real datasets, which provides empirical support for Contribution 2. \n\nWeaknesses\n1. The tightness of the derived bound is not discussed. \n2. The discussion of optimal p (sparsity parameter) on page 8 could be summarized in a more explicit way. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies a sketch-based approximation method of kernel regression. The following three items summarize the contributions of the authors. \n1. They derive a risk bound of the sketched kernel regression.\n2. They show that sketching by a sparse random matrix (sparse Gaussian or Rademacher) improves the time and space complexity while maintaining the theoretical guarantee.\n3. They derive an algorithm for training. ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Minor changes (not mandatory)\n1. It would be nice to discuss how tight the bounds of Theorem 4 are. For example, you can plot the theoretical bounds in Figure 1 (a) in addition to the empirical errors, allowing us to see the gap. \n2. You can summarize the conclusion of the optimal p~=0.7 in a theorem style so that readers easily digest it. ", "broader_impact_concerns": "No concern. ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "tWkLZFoQtC", "reviewer": "Reviewer_JjKC", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Contributions:\n- excess risk bounds for sketching of kernels with *generic Lipschitz loss* (Thm 2, Thm 4 for the multioutput case)\n- theoretical validity of sparsified Gaussian and Rademacher sketching techniques is proved\n- a technique to learn these kernels with an approximated feature map is proposed", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In the context of kernel methods, the paper studies a) learning with generic Lipschitz losses and b) faster computation of Gram matrix produced by so called *sparsified* sketching techniques.\nAs random projections methods, sketching techniques reduce the dimension of the problem while preserving useful information. However computing the resulting kernel matrices is either still costly, either results in poor quality sketches.\nBy introducing *p-sparsified sketches*, the authors allow for faster numerical solving of the subsequent optimization problems, as well as a reduced memory footprint, while preserving useful information from the full dataset.\n\nIn addition, the contribution goes beyond the squared loss by handling generic Lipschitz losses, which enables it to encompass the framework of robust regression (Huber loss), quantile regression with the pinball loss or the popular $\\epsilon$-insensitive losses.\nIn addition to standard scalar valued kernels, the results also apply to matrix-valued kernels, which play a role in multitask/multioutput learning, in the case of a *decomposable* kernel writing $K(\\cdot, \\cdot) = k(\\cdot, \\cdot) M$ for a fixed matrix $M$.\n\n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Questions:\n- P4 \"The existence and uniqueness of $\\delta_n^2$ is guaranteed for any unit ball: when what lives in the unit ball? And isn't existence clear since $\\mu_i$ are fixed, and the RHS grows like delta^2 while the LHS becomes constant for $\\delta \\geq \\mu_1$ and is strictly increasing before?\n- K-satisfiability of S with respect to ... instead of $\\delta^2_n$: where above is K-satisfiability proved with respect to $\\delta_n^2$, is it not assumed with respect to any constant $c$ in assumption 5? Clarifying this in assumption 5 (instead of saying below \"where $c$ is the constant of Ass. 5) may help.\n\nTypos:\n- P2, learning vector-valued... have been\n- P2 sketchin\n- P3, P4 is the notation arg inf instead of arg min standard?\n- why using the notation $\\mathcal{F}(\\mathcal{X}, \\mathbb{R})$ when above you have used $\\mathcal{Y}^\\mathcal{X}$?\n- In Definition 2 use $\\ll$ instead of $<<$\n- In Def 3 why using $\\tilde \\gamma$ instead of $\\hat \\gamma$, by analogy  with $\\hat \\alpha$?\n- $n^2$-identity matrix: is this the $n$ by $n$ identity matrix? If so why not using the introduced notation $I_n$?\n- P3 the excess of risk\n- P4 the estimators ... have bounded normS\n- An MVK : A MVK\n- P6: The 3rd term become\n- P17 We first prove first\n- P17 use of cite/citep/citet for Li et al, Bartlett etc.\n- \"let assume\" is not correct", "broader_impact_concerns": "n/a", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "WXNkRiCGML", "reviewer": "Reviewer_sd1L", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, I think this is a solid, well-written paper. To the best of my understanding, the theoretical results are not particularly ground-breaking, but having risk bounds for a wide variety of sketched ERM algorithms tied together by the satisfiability property is a nice result, even with the limitation to Lipschitz losses and bounded kernels; such restrictions can likely be (in principle) removed using some robustification techniques, although likely at great computational cost. As well, I cannot tell if the satisfiability guarantee for the authors $p$-sparse sketches is an easy extension of existing results for Gaussian sketches or not, but it is a natural extension and is nice to have.\n\nI am not really familiar with the related literature, but I felt that the authors were quite meticulous in their comparison of technical results here with those established in the existing literature.\n\nIn terms of weaknesses, my personal feeling is that more weight could have been placed on empirical analysis of the gains we get due to the \"decomposition trick,\" and how the statistical error changes across different datasets, all else equal. As $s$ gets larger, with $p$ fixed, it gets less and less likely that any column of $S$ will be completely filled with zeros, yet the computational boost (Fig 1(b)) doesn't change much. I found this a bit odd myself, but maybe I'm missing something. That said, the authors dedicate a lot of space to the scalar-->vector extension of their ERM analysis, a decision which I respect, but if the key point of interest in this paper is supposed to be the efficiency of $p$-sparse sketching, I think putting a bit more weight on related tests would be better aligned with the overall take-home message here.\n\nIf this were a review for an international conference, I think some reviewers would take issue that the ERM analysis is quite straightforward; it appears to be typical concentration plus porting the approximation error analysis of Yang et al 2017 under satisfiability. That said, this is TMLR, and as far as I can tell the results are valid, so the claims are solid. Plus, the multi-output extension has been formulated in a clean way, which I also think has value.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The broad context of this paper is that of designing and analyzing kernel-based machine learning methods, particularly for regression tasks, which are more computationally efficient than more naive implementations, while maintaining statistical guarantees (up to natural tradeoffs due to reduced precision).\n\nThe more refined context is that of \"sketching.\" The basic idea is to replace the Gram matrix that appears in standard formulations of kernel-based ERM with a matrix of smaller dimension (by replacing gram matrix $K$ with $SKS^{\\\\top}$), here sketch matrix $S$ is $s \\\\times s$ instead of $K$ which is $n \\times n$, where typically $s \\\\ll n$, and $n$ is the number of (iid) data points available for training. Arguably the simplest approach is to randomly sub-sample this data, but there are numerous other approaches. The first contribution of this work is a rather general analysis of ERM under sketching, with the critical technical starting point being the notion of \"satisfiability\", a property of the sketching matrix $S$. Roughly speaking, we want $S$ to preserve the \"important\" information in $K$ (preserving length of largest eigenvectors of $K$), while discarding the less important information. Their results take the form of excess risk bounds that hold for Lipschitz losses and bounded kernels, with ERM restricted to a ball around the true risk minimizer in RKHS, assuming that $S$ satisfies the \"satisfiability\" requirement. They handle both single output and multi-output regression tasks, but the overall flavor is the same.\n\nThe other main contribution is a method for generating sketch matrices (what they call $p$-sparsified sketches, of the Rademacher and Gaussian varieties), plus formal characterization of a \"good event\" (since the sketching is probabilistic) of high probability on which the $K$-satisfiability of the (random) $S$ generated by their method holds (Theorem 5). The notion of $p$-sparse sketches itself is not new here; the basic idea is simply to set a fixed probability (here, $1-p$) that each iid element of $S$ will be zero, otherwise it follows a typical (scaled) Gaussian/Rademacher distribution. According to the authors, the Rademacher variety has appeared previously (although with theoretical guarantees lacking), but the Gaussian variety is new here. By enforcing such sparsity, the authors show (Eqn 14) how the sketching can then be split into a fast sub-sampling sketch combined with a (small, post-sparsification) Gaussian sketch; this is only possible when full columns of $S$ are zero vectors, which is at least in principle possible here, whereas many existing sketching methods always have at least one non-zero element per column. They provide basic empirical evidence for the efficiency of their approach using both simulated and benchmark datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please see the comments above; none of the suggestions are critical to my recommendation, as I think the paper is already quite well-done as is, but I think would benefit the paper.\n\nBesides these points, here are a few minor points that I think the authors should consider, in no particular order of importance.\n\n- Why do the authors give the quadratic loss result in Thm 2? Bounding $\\\\mathcal{Y}$ and assuming a bounded kernel and finite diameter ball around the risk minimizer yields Lipschitzness for this function; is this special case for comparison with other works? Do there exist similar results (excess risk bounds) assuming the quadratic loss but without bounded $\\\\mathcal{Y}$? If so, I think this would be a natural point of comparison (such comparisons appear in more general regression contexts, outside of sketching/kernels).\n\n- A point on nomenclature; using the current naming, a $0$-sparsified sketch is completely sparse, which is somewhat unappealing aesthetically.\n\n- In the exposition following Eqn (14), the authors say that in their decomposition $S\\_{\\\\mathrm{SG}}$ is a \"sub-Gaussian sketch\". I don't see where this term comes from, and I'm not sure if it is meaningful here. Is this term used just because the elements of $S\\_{\\\\mathrm{SG}}$ are either Rademacher or Gaussian (assuming we are restricted to the confines of Defn 5), which are both sub-Gaussian? We could say the same thing about $S$ itself, so it seems weird to introduce the term here as if $S\\_{\\\\mathrm{SG}}$ is special in this sense. Maybe I'm missing something.\n\n- Assumption 5: I think it would be better to say \"$S$ is $K$-satisfiable with probability 1\" or something similar.\n\n- $\\\\tilde{f}\\_{s}$: this is a critical object in Thm 2, but it is buried in Defn 2; a link to this in Thm 2 would be nice.\n\n- Defn 2: replace $<<$ with $\\\\ll$.\n\n- Spelling/grammar: \"sketchin\" (p.2), \"kernels methods\" (p.2), ...", "broader_impact_concerns": "No relevant concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper introduces the concept of p-sparsified sketches, which enable quicker numerical solutions for subsequent optimization problems. This approach also reduces memory usage while retaining valuable information from the complete dataset. The paper offers a comprehensive theoretical analysis of p-sparsified sketches. Experimental results further illustrate that the proposed method surpasses existing techniques in performance.

The paper is well-written and has an important contribution to kernel methods. All reviewers have recommended its acceptance. I also vote for acceptance.
Decision: Accept (poster)