Vivien Cabannes
email: vivien.cabannes@gmail.com

Loucas Pillaud-Vivien

Francis Bach

Alessandro Rudi

Overcoming the curse of dimensionality with Laplacian regularization in semi-supervised learning

As annotations of data can be scarce in large-scale practical problems, leveraging unlabelled examples is one of the most important aspects of machine learning. This is the aim of semi-supervised learning. To benefit from the access to unlabelled data, it is natural to diffuse smoothly knowledge of labelled data to unlabelled one. This induces to the use of Laplacian regularization. Yet, current implementations of Laplacian regularization suffer from several drawbacks, notably the well-known curse of dimensionality. In this paper, we provide a statistical analysis to overcome those issues, and unveil a large body of spectral filtering methods that exhibit desirable behaviors. They are implemented through (reproducing) kernel methods, for which we provide realistic computational guidelines in order to make our method usable with large amounts of data.

In the last decade, machine learning has been able to tackle amazingly complex tasks, which was mainly allowed by computational power to train large learning models on large annotated datasets. For instance, ImageNet is made of tens of millions of images, which have all been manually annotated by humans [START_REF] Deng | Imagenet: A large-scale hierarchical image database[END_REF]. The greediness in data annotation of such a current learning paradigm is a major limitation. In particular, when annotation of data demands in-depth expertise, relying on techniques that require zillions of labelled data is not viable. This motivates several research streams to overcome the need for annotations, such as self-supervised learning for images or natural language processing [START_REF] Devlin | BERT: pre-training of deep bidirectional transformers for language understanding[END_REF]. Aiming for generality, semi-supervised learning is the most classical one, assuming access to a vast amount of input data, but among which only a scarce percentage is labelled. To leverage the presence of unlabelled data, most semi-supervised techniques assume a form of low-density separation hypothesis, as detailed in the recent review of van Engelen and Hoos [START_REF] Van Engelen | A survey on semi-supervised learning[END_REF], and illustrated by state-of the-art models [START_REF] Berthelot | Mixmatch: A holistic approach to semi-supervised learning[END_REF][START_REF] Verma | Interpolation consistency training for semi-supervised learning[END_REF]. This hypothesis assumes that the function to learn from the data varies smoothly in highly populated regions of the input space, but might vary more strongly in scarcely populated areas, or that the decision frontiers between classes lie in regions with low-density. In such a setting, it is natural to enforce constraints on the variations of the function to learn. While semi-supervised learning is an important learning framework, it has not provided as much exciting realizations as one could have expected. This might be related to the fact that it is classically approached through graph-based Laplacian, a technique that does not scale well with the dimension of the input space [START_REF] Bengio | Label propagation and quadratic criterion[END_REF].

Paper organization. In Section 1, we motivate Laplacian regularization, and recall drawbacks of naive implementations. These limitations are overcome in Section 2 where we expose a theoretically principled path to derive well-behaved algorithms. More precisely, we unveil a vast class of estimates 35th Conference on Neural Information Processing Systems (NeurIPS 2021).

based on spectral filtering. We turn to implementation in Section 3 where we provide realistic guidelines to ensure scalability of the proposed algorithms. Statistical properties of our estimators are stated in Section 4.

Contributions. They are two folds. (i) Statistically, we explain that Laplacian regularization can be properly leveraged based on functional space considerations, and that those considerations can be turned into concrete implementations thanks to kernel methods. As a result, we provide consistent estimators that exhibit fast convergence rates under a low density separation hypothesis, and that, in particular, do not suffer from the curse of dimensionality. (ii) Computationally, we avoid dealing with large matrices of derivatives by providing a low-rank approximation that allows to deal with n γ log(n) × n γ log(n) matrices, with a parameter γ ∈ (0, 1] depending on the regularity of the problem, instead of n(d + 1) × n(d + 1) matrices, thus cutting down to O(log(n) 2 n 1+2γ d) the potential O(n 3 d 3 ) training cost.

Related work. Interplays between graph theory and machine learning were proven successful in the 2000s [START_REF] Smola | Kernels and regularization on graphs[END_REF]. The seminal paper of Zhu et al. [START_REF] Zhu | Semi-supervised learning using Gaussian fields and harmonic functions[END_REF] introduced graph-Laplacian as a transductive method in the context of semi-supervised learning. A smoothing variant was proposed by [START_REF] Zhou | Learning with local and global consistency[END_REF], which is coherent with the fact that enforcing constraints on labelled points leads to spikes [START_REF] Alaoui | Asymptotic behavior of p -based Laplacian regularization in semi-supervised learning[END_REF]. Interestingly, graph Laplacians do converge to diffusion operators linked with the weighted Laplace Beltrami operator [START_REF] Hein | Graph Laplacians and their convergence on random neighborhood graphs[END_REF][START_REF] García Trillos | Error estimates for spectral convergence of the graph Laplacian on random geometric graphs toward the Laplace-Beltrami operator[END_REF]. However, these local diffusion methods are known to suffer from the curse of dimensionality [START_REF] Bengio | Label propagation and quadratic criterion[END_REF]. That is, local averaging methods are intuitive learning methods that have been used for more than half a century [START_REF] Fix | Discriminatory analysis. nonparametric discrimination: Consistency properties[END_REF]. Yet, those methods do not scale well with the dimension of the input space [START_REF] Yang | Minimax nonparametric classification. I. Rates of convergence[END_REF]. This is related with the fact that to cover [0, 1] d , we need ε -d balls of radius ε. Interestingly, if the function to learn is m times differentiable with smooth partial derivatives, it is possible to leverage more information from function evaluations and overcome the curse of dimensionality when m d. This property is related to covering numbers (a.k.a. capacity) of Sobolev spaces [START_REF] Kolmogorov | ε-entropy and ε-capacity of sets in functional spaces[END_REF] and is leveraged by (reproducing) kernel methods [START_REF] Steinwart | Support Vector Machines[END_REF][START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF]. The crux of this paper is apply this fact to Laplacian regularization techniques. Note that derivative with reproducing kernel methods in machine learning have already been considered in different settings by [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF][START_REF] Rosasco | Nonparametric sparsity and regularization[END_REF][START_REF] Eriksson | Scaling Gaussian process regression with derivatives[END_REF].

Laplacian regularization

In this section, we introduce the notations and concepts related to the semi-supervised learning regression problem, noting that most of our results extend to any convex loss beyond least-squares. We motivate and describe Laplacian regularization that will allow us to leverage the low-density separation hypothesis. We explain statistical drawbacks usually linked with Laplacian regularization, and discuss on how to circumvent them.

In the following, we denote by X = R d the input space, Y = R the output space, and by ρ ∈ ∆ X ×Y the joint distribution on X × Y. For simplicity, we assume that ρ has compact support. In the following, we denote by ρ X the marginal of ρ over X , and by ρ| x the conditional distribution of Y given X = x. As usual, for p ∈ N * , L p (R d ) is the space of functions f such that f p is integrable. Moreover, we define usual Sobolev spaces: for s ∈ N, W s,p (R d ) stands for the space of functions whose weak derivatives of order s-th are in L p (R d ). When p = 2, they have a Hilbertian structure and we denote, H s (R d ) = W s,2 (R d ) these Hilbertian spaces. Ideally, we would like to retrieve the mapping g * : X → Y defined as

g * = arg min g∈L 2 (ρ X ) E (X,Y )∼ρ g(X) -Y 2 = arg min g∈L 2 (ρ X ) g -g ρ 2 L 2 (ρ X ) = g ρ , (1)
where g ρ : X → Y is defined as

g ρ (x) = E [Y | X = x].
In semi-supervised learning, we assume that we do not have access to ρ but we have access to n independent samples (X i ) i≤n ∼ ρ ⊗n X , among which we have n labels Y i ∼ ρ| Xi for i ≤ n , with n potentially much smaller than n. In other terms, we have n supervised pairs (X i , Y i ) i≤n , and n -n unsupervised samples (X i ) n <i≤n . While we restrict ourselves to real-valued regression for simplicity, our exposition indeed applies generically to partially supervised learning. In particular, it can be used off-the-shelve to complement the approaches of [START_REF] Cabannes | Structured prediction with partial labelling through the infimum loss[END_REF][START_REF] Cabannes | Disambiguation of weak supervision with exponential convergence rates[END_REF] as we detailed in Appendix A.

Semi-supervision setting

Baseline: = 0, = 1

Result: = 1, = 1/ Figure 1: Motivating example. (Left) We suppose given n = 2000 points in X = R 2 , represented as black dots, spanning 4 concentric circles. Among those points are n = 4 labelled points, with labels being either 1 represented in red, and -1 represented in blue. In this setting, it is natural to assume that g * should be constant on each circles, which can be encoded as ∇g * = 0 on supp ρX . (Middle) Kernel ridge regression estimate based on the labelled points with Gaussian kernel of bandwidth σ = .2r, r being the radius of the innermost circle. (Right) Laplacian regularization reconstruction. The reconstruction is based on approximate empirical risk minimization with p = n, which ensures a computational complexity of O(p 2 nd), instead of O(n 3 d 3 ) needed to recover the exact empirical risk minimizer [START_REF] Bengio | Label propagation and quadratic criterion[END_REF].

Diffusion operator L

In order to leverage unlabelled data, we will assume that g * varies smoothly on highly populated regions of X , and might vary highly on low density regions. For example, this is the case when data are clustered in well separated regions of space, and labels are constant on clusters. This is captured by the fact that the Dirichlet energy

X ∇g * (x) 2 ρ X (dx) = E X∼ρ X ∇g * (X) 2 =: L 1/2 g 2 L 2 (ρ X ) , (2)
is assumed to be small. Because the quadratic functional (2) will play a crucial role in our exposition, we define L as the self-adjoint operator on L 2 (ρ X ), extending the operator on H 1 (ρ X ) representing this functional. Under mild assumptions on ρ X , L -1 can be shown to be a compact operator, which we will assume in the following. In essence, we will assume that if we have a lot of unlabelled data and L 1/2 g can be well approximated for any function g, then we do not need a lot of labelled data to estimate correctly g * . To illustrate this, at one extreme, if we know that L 1/2 g * = 0, then g * is known to be constant on each connected component of ρ X so that, along with the knowledge of ρ X , only a few labelled points would be sufficient to recover perfectly g * . We illustrate those considerations on Figure 1.

Drawbacks of naive Laplacian regularization

Following the motivations presented precedently, it is natural to consider the regularized objective and solution defined, for λ > 0, as

g λ = arg min g∈H 1 (ρ X ) E (X,Y )∼ρ g(X) -Y 2 + λ E X∼ρ X ∇g(X) 2 R d = arg min g∈H 1 (ρ X ) g -g ρ 2 L 2 (ρ X ) + λ L 1/2 g 2 L 2 (ρ X ) = (I + λL) -1 g ρ . (3)
This regularization has nice properties. In particular, for small λ, it can be seen as a first order approximation of the heat equation solution e -λL g ρ , which represents the temperature profile at time t = λ, instantiated with the initial profile g ρ , and with ρ X modelling the thermal conductivity. It also has interpretations in term of random walk and Langevin diffusion [START_REF] Pillaud-Vivien | Statistical estimation of the poincaré constant and application to sampling multimodal distributions[END_REF][START_REF] Hamzi | Kernel-based approximation of the koopman generator and schrödinger operator[END_REF]. In a word, g λ is the diffusion of g ρ with respect to the density ρ X , which relates to the idea of diffusing labelled data with respect to the intrinsic geometry of the data, which is the idea captured by [START_REF] Zhu | Semi-supervised learning using Gaussian fields and harmonic functions[END_REF].

However, from a learning perspective, Eq. ( 3) is linked with the prior that g * belongs to H 1 (ρ X ), a prior that is not strong enough to overcome the curse of dimensionality as we saw in the related work section. Moreover, assuming we have enough unsupervised data to suppose known ρ X , and therefore L, Eq. ( 3) leads to the naive empirical estimate g (naive) ∈ arg min g:X →R

n i=1 g(X i ) -Y i 2 + n λ L 1/2 g 2 .
While the definition of g (naive) could seem like a great idea, in fact, such an estimate g (naive) is known to be mostly constant and spiking to interpolate the data (X i , Y i ) as soon as d > 2 [START_REF] Nadler | Statistical analysis of semi-supervised learning: The limit of infinite unlabelled data[END_REF]. This is to be related with the capacity of the space associated with the pseudo-norm L 1/2 g in L 2 . This capacity, related to H 1 , is too large for the Laplacian regularization term to constraint g (naive) in a meaningful way. In other terms, we need to regularize with stronger penalties.

Stronger regularization

In this subsection, we discuss techniques to overcome the issues encountered with g (naive) . Those techniques are based on functional space constraints or on spectral filtering techniques.

Functional spaces. A solution to overcome the capacity issue of H 1 in L 2 is to constrain the estimate of g * to belong to a smaller functional space. In the realm of graph Laplacian, [START_REF] Alaoui | Asymptotic behavior of p -based Laplacian regularization in semi-supervised learning[END_REF] proposed to solve this problem by considering the r-Laplacian regularization reading Ω r = X ∇g(X) r ρ(dx), with r > d. In essence, this restricts g to live in W 1,r (ρ X ) for r > d, and allows to avoid spikes associated with g (naive) . However considering high power of the gradient is likely to introduce instability (think that d is the potentially really big dimension of the input space), and from a learning perspective, the capacity of W 1,r , which compares to the one of H 2 , is still too big.

In this paper, we will rather keep the diffusion operator L, and add a second penalty to reduce the space in which we look for the solution. With G an Hilbert space of functions, we could look for, with µ > 0 a second regularization parameter

g λ,µ = arg min g:G∩H 1 (ρ X ) g -g ρ 2 L 2 (ρ X ) + λ L 1/2 g 2 L 2 (ρ X ) + λµ g 2 G . (4)
This formulation restricts g λ,µ to belong both to H 1 (ρ X ) (thanks to the term in λ) and G (thanks to the term in µ). In particular the resulting space H 1 (ρ X ) ∩ G to which g λ,µ belongs, has a smaller capacity in L 2 than the one of G in L 2 . In practice, we do not have access to ρ and ρ X but to (X i , Y i ) i≤n and (X i ) i≤n , and we might consider the empirical estimator defined through empirical risk minimization

g n ,n = arg min g∈G n -1 n i=1 g(X i ) -Y i 2 + λn -1 n i=1 ∇g(X i ) 2 + λµ g 2 G . (5)
For example, we could consider G to be the Sobolev space H m (dx). Note the difference between G linked with dx, the Lebesgue measure, that is known, and L linked with ρ X , the marginal of ρ over X , that is not known. In this setting, the regularization

L 1/2 g 2 + µ g 2 G reads X Dg(x) 2 ρ X (dx) + µ X m α=0 D α g(x)
2 dx. Because of the size of H m in L 2 , this allows for efficient approximation of g λ,µ based on empirical risk minimization. In particular, if n = +∞, we expect the minimizer (5) to converge toward g λ,µ at rates in L 2 scaling similarly to n -m/d in n . To complete the picture, depending on a prior on g ρ , g λ,µ might exhibit good convergence properties towards g ρ as λ and µ go to zero. This contrasts with the problem encountered with g (naive) . Those considerations are exactly what reproducing kernel Hilbert space will provide, additionally with a computationally friendly framework to perform the estimation. Note that quantity similar to g λ,µ were considered in [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF][START_REF] Rosasco | Nonparametric sparsity and regularization[END_REF]. Spectral filtering. Without looking for higher power-norm, [START_REF] Nadler | Statistical analysis of semi-supervised learning: The limit of infinite unlabelled data[END_REF] proposed to overcome the capacity issue by considering approximation of the operator L based on the graph-based technique provided by [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF][START_REF] Coifman | Diffusion maps[END_REF] and to reduce the search of g n on the space spanned by the first few eigenvectors of the Laplacian. In particular, on Figure 1, g * could be searched in the null space of L, that is, among functions that are constant on each connected component of supp ρ X . This technique exhibits two parts, the "unsupervised" estimation of L that will depend on the total number of data n, and the "supervised" search for g ρ on the first few eigenvectors of L that will depend on the number of labels n . While, at first sight, this technique seems to be completely different than Tikhonov regularization (4), it can be cast, along with gradient descent, into the same spectral filtering framework [START_REF] Lin | Optimal rates for spectral algorithms with least-squares regression over Hilbert spaces[END_REF]. This point of view enables the use of a wide range of techniques offered by spectral manipulations on the diffusion operator L. This paper is motivated by the fact that current well-grounded semi-supervised learning techniques are implemented based on graph-based Laplacian, which is a local averaging method that does not leverage smartly functional capacity. In particular, as recalled earlier, graph-based Laplacian is known to suffer from the curse of dimensionality, in the sense that the convergence of the empirical estimator L towards the L exhibits a rate of convergence of order O(n -1/d ) with d the dimension of the input space X [START_REF] Hein | Graph Laplacians and their convergence on random neighborhood graphs[END_REF]. In this work, we will bypass this curse of dimensionality by looking for g in a smooth universal reproducing kernel Hilbert space, which will lead to efficient empirical estimates.

Spectral Filtering with Kernel Laplacian

In this section, we approach Laplacian regularization from a functional analysis perspective. We first introduce kernel methods and derivatives in reproducing kernel Hilbert space (RKHS). We then translate the considerations provided in Section 1.3 in the realm of kernel methods.

Kernel methods and derivatives evaluation maps

In this subsection, we introduce kernel methods (see [START_REF] Aronszajn | Theory of reproducing kernels[END_REF][START_REF] Scholkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF][START_REF] Steinwart | Support Vector Machines[END_REF] for more details Consider (H, •, • H ) a reproducting kernel Hilbert space, that is a Hilbert space of functions from X to R such that the evaluation functionals L x : H → R; g → g(x) are continuous linear forms for any x ∈ X . Such forms can be represented by k x ∈ H such that, for any g ∈ H, L x (g) = k x , g H . A reproducing kernel Hilbert space can alternatively be defined from a symmetric positive semidefinite kernel k : X → X → R, that is a function such that for any n ∈ N and (x i ) i≤n ∈ X n the matrix (k(x i , x j )) i,j is symmetric positive semi-definite, by building (k x ) x∈X such that k(x, x ) = k x , k x H . From a learning perspective, it is useful to use the evaluation maps to rewrite

H = {g θ : x → k x , θ H | θ ∈ H}.
As such, kernel methods can be seen as "linear models" with features k x , allowing to parameterize large spaces of functions [START_REF] Micchelli | Universal kernels[END_REF]. In the following, we will differentiate θ seen as an element of H and g θ seen as its embedding in L 2 . To make this distinction formal, we define the embedding S : (H,

•, • H ) → (L 2 (ρ X ), •, • L2 ); θ → g θ , as well as its adjoint S : L 2 (ρ X ) → H.
Given a linear parametric model of functions g θ (x) = θ, k x H , it is possible to compute derivatives of g θ based on derivatives of the feature vector -think of H = R p and of k x = ϕ(x) as a feature vector with ϕ : R d → R p . For α ∈ N d , with |α| = i≤d α i , we have the following equality of partial derivatives, when k is 2 |α| times differentiable,

D α g θ (x) = θ, D α k x ,
where

D α = ∂ |α| (∂x 1 ) α1 (∂x 2 ) α2 • • • (∂x d ) α d .
Here and D α k x has to be understood as the partial derivative of the mapping of x ∈ X to k x ∈ H, which can be shown to belong to H [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF]. In the following, we assume that k is twice differentiable with continuous derivatives, and will make an extensive use of derivatives of the form ∂ i k x = ∂k x /∂x i for i ≤ d and x ∈ X . Note that, as well as we can describe completely the Hilbertian geometry of the space Span {k x | x ∈ X } through k(x, x ) = k x , k x , for x, x ∈ X , we can describe the Hilbertian geometry of Span {k

x | x ∈ X } + Span {∂ i k x | x ∈ X }, through ∂ 1,i k(x, x ) = ∂ i k x , k x H , and ∂ 1,i ∂ 2,j k(x, x ) = ∂ i k x , ∂ j k x H ,
where ∂ 1,i denotes the partial derivative with respect to the i-th coordinates of the first variable. This echoes to so-called "representer theorems".

Example 1 (Gaussian kernel). A classical kernel is the Gaussian kernel, also known as radial basis function, defined for σ > 0 as the following k, and satisfying, for i = j, the following equalities,

k(x, x ) = exp - x -x 2 2σ 2 , ∂ 1,i ∂ 2,j k(x, y) = - (x i -y i )(x j -y j ) σ 4 k(x, y), ∂ 1,i k(x, y) = - (x i -y i ) σ 2 k(x, y), ∂ 1,i ∂ 2,i k(x, y) = 1 σ 2 - (x i -y i ) 2 σ 4 k(x, y),
where x i designs the i-th coordinates of the vector x ∈ X = R d .

Eigen vector #1 ( 1 )

Eigen vector #4 ( 4 ) Eigen vector #8 ( 8 )

Figure 2: Few of the first generalized eigenvectors of ( Σ; L + µI) (with µ = 1/n). The first four eigenvectors correspond to constant functions on each circle, as shown with e1 and e4. The few eigenvectors after correspond to second harmonics localized on a single circle as shown with e8.

Tikhonov, spectral filtering and dimensionality reduction

Given the kernel k, its associated RKHS H and S the embedding of H in L 2 , we rewrite Eq. ( 4) under its "parameterized" version

g λ,µ = S arg min θ∈H Sθ -g ρ 2 L 2 (ρ X ) + λ L 1/2 Sθ 2 L 2 (ρ X ) + λµ θ 2 H . (4)
Do not hesitate to refer to Table 1 to keep track of notations. In the following, we will use that

L 1/2 Sθ 2 L 2 (ρ X ) + µ θ 2 H = (S LS + µI) 1/2 θ 2 H
. This equality explains why we consider µλ instead of µ in the last term. In the RKHS setting, the study of Eq. ( 4) unveils the three operators Σ, L, and I on H, (indeed g λ,µ = S arg min θ∈H {θ (Σ + λL + λµ)θ -2θ S g ρ }) where I is the identity, and, as we detail in Appendix C,

Σ = S S = E X∼ρ X [k X ⊗ k X ] ,
and

L = S LS = E X∼ρ X d i=1 ∂ j k X ⊗ ∂ j k X . (6)
Regularization and spectral filtering have been well-studied in the inverse-problem literature. In particular, the regularization Eq. ( 4) is known to be linked with the generalized singular value decomposition of [Σ; L + µI] (see, e.g., [START_REF] Edelman | The GSVD: where are the ellipses?, matrix trigonometry, and more[END_REF]), which is linked to the generalized eigenvalue decomposition of (Σ, L + µI) [START_REF] Golub | Matrix computations[END_REF]. We derive the following characterization of Eq. ( 4), whose proof is reported in Appendix D.

Proposition 1. Let (λ i,µ ) i∈N ∈ R N , (θ i,µ ) i∈N ∈ H N be the generalized eigenvalue decomposition of the pair (Σ, L + µI), that is (θ i,µ ) generating H and such that for any i, j ∈ N, Σθ i,µ = λ i,µ (L + µI)θ i,µ , and θ i,µ , (L + µI)θ j,µ = 1 i=j . Eq. (4) can be rewritten as

g λ,µ = i∈N ψ(λ i,µ )Sθ i,µ ⊗ Sθ i,µ g ρ = i∈N ψ(λ i,µ ) S g ρ , θ i,µ Sθ i,µ , (7)
with ψ : R + → R; x → (x + λ) -1 . Eq. (7) should be seen as a specific instance of spectral filtering based on a filter function ψ : R + → R.

Interestingly, the generalized eigenvalue decomposition of the pair (Σ, L+µI) was already considered by Pillaud-Vivien [START_REF] Pillaud-Vivien | Statistical estimation of the poincaré constant and application to sampling multimodal distributions[END_REF] to estimate the first eigenvalue of the Laplacian. Moreover, Pillaud-Vivien [START_REF] Pillaud-Vivien | Learning with Reproducing Kernel Hilbert Spaces: Stochastic Gradient Descent and Laplacian Estimation[END_REF] suggests to leverage this decomposition for dimensionality reduction based on the first eigenvectors of the Laplacian. As well as Eq. ( 4) contrasts with graph-based semi-supervised learning techniques, this dimensionality reduction technique contrasts with methods based on graph Laplacian provided by [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF][START_REF] Coifman | Diffusion maps[END_REF]. Remarkably, the semi-supervised learning algorithm that consists in using the unsupervised data to perform dimensionality reduction based on the Laplacian eigenvalue decomposition, before solving a small linear regression problem on the small resulting space, can be seen as a specific instance of spectral filtering, based on regularization by thresholding/cutting-off eigenvalue, which corresponds to ψ : x → x -1 1 x>λ for a given threshold λ > 0 in Eq. [START_REF] Cabannes | Structured prediction with partial labelling through the infimum loss[END_REF].

Implementation

In this section, we discuss on how to practically implement estimates for Eq. ( 7) based on empirical data (X i , Y i ) i≤n and (X i ) n <i≤n . We first review how we can approximate the integral operators of Eq. ( 6) based on data. We then discuss on how to implement our methods practically on a computer. We end this section by considering approximations that allow to cut down high computational costs associated with kernel methods involving derivatives.

Algorithm 1: Empirical estimates based on spectral filtering. Data:

(X i , Y i ) i≤n , (X i ) n <i≤n , a kernel k, a filter ψ, a regularizer µ Result: ĝp through c ∈ R p defining ĝp (x) = n i=1 c i k(x, X i ) = k x T a c Compute S n T a = (k(X i , X j )) i≤n,j≤p ∈ R n×p in O(pn) Compute Z n T a = (∂ 1,j k(X l , X i )) (j≤d,l≤n),i≤p ∈ R nd×p in O(pnd) Build T a ΣT a = n -1 (S n T a ) (S n T a ) in O(p 2 n) Build T a LT a = n -1 (Z n T a ) (Z n T a ) in O(p 2 nd) a Build T a T a = (k(X i , X j )) i,j≤p ∈ R p×p in O(1) as a partial copy of S n T a Get (λ i,µ , u i,µ ) i≤n the generalized eigenelements of (T a ΣT a , T a ( L + µI)T a ) in O(p 3 ) Get b = T a θ = (n -1 n i=1 Y i k(X i , X j )) j≤p ∈ R p in O(pn ) Return c = n i=1 ψ(λ i )u i u i b ∈ R p in O(p 3 ).
a Building this matrix can be avoided by using the generalized singular value decomposition rather than the generalized eigenvector decomposition. Implemented with Lapack, such a procedure will also requires O(p 2 nd) floating point operations, but with a smaller constant in the big O [START_REF] Golub | Matrix computations[END_REF].

Integral operators approximation

The classical empirical risk minimization in Eq. ( 5) can be understood as the plugging of the approximate distributions ρ = n -1 n i=1 δ Xi ⊗ δ Yi and ρX = n -1 n i=1 δ Xi instead of ρ and ρ X in Eq. ( 4). It can also be understood as the same replacement when dealing with integral operators, leading to the three following important quantities to rewrite Eq. [START_REF] Cabannes | Structured prediction with partial labelling through the infimum loss[END_REF],

Σ := n -1 n i=1 k Xi ⊗ k Xi , L := n -1 n i=1 d j=1 ∂ j k Xi ⊗ ∂ j k Xi , θ := S g ρ := n -1 n i=1 Y i k Xi . (8)
It should be noted that while considering n in the definition of Σ is natural from the spectral filtering perspective, to make it formally equivalent with the empirical risk minimization [START_REF] Bengio | Label propagation and quadratic criterion[END_REF], it should be replaced by n . Eq. ( 8) allows to rewrite Eq. ( 7) without relying on the knowledge of ρ, by considering ( λi,µ , θi,µ ) the generalized eigenvalue decomposition of ( Σ, L) and considering

ĝ = i∈N ψ( λi,µ ) S g ρ , θi,µ S θi,µ , (9)
We present the first eigenvectors (after plunging them in L 2 through S) of the generalized eigenvalue decomposition of ( Σ, L + µI) on Figure 2. The first eigenvectors allow to recover the null space of L. This explains clearly the behavior on the right of Figure 1. The canonical representation consists in taking p = n(d + 1) and considering for c ∈ R n(d+1) , the mapping [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF][START_REF] Rosasco | Nonparametric sparsity and regularization[END_REF]. This exact implementation implies dealing and finding the generalized eigen value decomposition of p × p matrices with p = n(d + 1), which leads to computational costs in O(n 3 d 3 ), which can be prohibitive. Two solutions are known to cut down prohibitive computational costs of kernel methods.

Matrix representation and approximation of operators

T c c = n i=1 c i0 k Xi + d j=1 c ij ∂ j k Xi
Both methods consist in looking for a space that can be parameterized by R p for a small p and that approximates well the space H n + H n,∂ ⊂ H. The first solution is provided by random features [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF]. It consists in approximating H with a space of small dimension p ∈ N, linked with an explicit representation ϕ : X → R p that approximate k(x, x ) k ϕ (x, x ) = ϕ(x), ϕ(x ) R p . In theory, it substitutes the kernel k by k ϕ . In practice, all computations can be done with the explicit feature ϕ.

Approximate solution. The second solution, which we are going to use in this work, consists in approximating H n + H n,∂ by H p = Span {k Xi } i≤p for p ≤ n. This method echoes the celebrated Nyström method [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF], as well as the Rayleigh-Ritz method for Sturm-Liouville problems. In essence, [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF] shows that, when considering subsampling based on leverage score, p = n γ log(n), with γ ∈ (0, 1] linked to the "size" of the RKHS and the regularity of the solution, is a good enough approximation, in the sense that it only downgrades the sample complexity by a constant factor. In theory, we know that the space H p will converge to H = Closure Span {k x } x∈supp ρ X as p goes to infinity. In practice, it means considering the approximation mapping T a : R p → H; c → p i=1 c i k Xi , and dealing with the p × p matrices T a ΣT a and T a LT a . It should be noted that the computation of T a LT a requires to multiply a p × nd matrix by its transpose. Overall, training this method can be done with O(p 2 nd) basic operations, and inference with this method can be done in O(p). The saving cost of this approximate method is huge: without compromising the precision of our estimator, we went from O(n 3 d 3 ) run time complexities to O(log(n) 2 n 1+2γ d) computations, with γ possibly very small. Similarly, the memory cost went from O(n 2 d 2 ) down to O(nd + n 2γ ). 1

Statistical analysis

In this section, we are interested in quantifying the risk of the learnt mapping ĝ. We study it through the generalization bound, which consists in obtaining a bound on the averaged excess risk

E data ĝ -g ρ 2 L 2 .
In particular, we want to answer the following points. 1. How, and under which assumptions, Laplacian regularization boost learning? 2. How the excess of risk relates to the number of labelled and unlabelled data?

In terms of priors, we want to leverage a low-density separation hypothesis. In particular, we can suppose that when diffusing g ρ with e -tL we stay close to g ρ , or that g ρ is supported on a finite dimensional space of functions on which L 1/2 g (which measures the variation of g) is small. Both those assumptions can be made formal by assuming the g ρ is supported by the first eigenvectors of the diffusion operator L.

Assumption 1 (Source condition). g ρ is supported on a finite dimensional space that is left stable by the diffusion operator L. In other terms, if (e i ) ∈ (L 2 ) N are the eigenvectors of L, there exists r ∈ N, such that g ρ ∈ Span {e i } i≤r .

We will also assume that the diffusion operator L can be well approximated by the RKHS associated with k. In practice, under mild assumptions, c.f. Appendix C, the eigenvectors of the Laplacian are known to be regular, in particular to belong to H m for m ∈ N bigger than d. As such, many classical kernels would allow to verify the following assumption.

Assumption 2 (Approximation condition). The eigenvectors (e i ) of L belongs to the RKHS H.

We add one technical assumptions regarding the eigenvalue decay of the operator Σ compared to the operator L, with denoting the Löwner order (i.e., for A and B symmetric,

A B if B -A is positive semi-definite).

Assumption 3 (Eigenvalue decay).

There exists a ∈ [0, 1] and c > 0 such that L cΣ α .

Note that, in our setting, L is compact and bounded and Assumption 3 is always satisfied with a = 0. For translation-invariant kernel, such as Gaussian or Laplace kernels, based on considerations linking eigenvalue decay of operators with functional space capacities [START_REF] Steinwart | Support Vector Machines[END_REF], under mild assumptions, we can take a > 1 -2/d. We discuss all assumptions in more details in Appendix C.

To study the consistency of our algorithms, we can reuse the extensive literature on kernel ridge regression [START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF][START_REF] Lin | Optimal rates for spectral algorithms with least-squares regression over Hilbert spaces[END_REF]. This literature body provides an extensive picture on convergence rates relying d+4 log(n) as suggested by graph-based theoretical results [START_REF] Hein | Graph Laplacians and their convergence on random neighborhood graphs[END_REF]. We report classification error as a function of the number of samples n. The error is averaged over 50 trials, with errorbars representing standard deviations. We fixed the ratio n /n to one tenth, and generated the data according to two Gaussians in dimension d = 10 with unit variance and whose centers are at distance δ = 3 of each other (similar to the setting of [START_REF] Castelli | On the exponential value of labeled samples[END_REF][START_REF] Lelarge | Asymptotic Bayes risk for Gaussian mixture in a semisupervised setting[END_REF]). Our method discovers the structure of the data much faster than graph-based Laplacian (to get a 20% error we need 40 points, while graph-based need 700). (Right) Time to perform training with graph-based Laplacian in orange, with Algorithm 1 in blue (with the specification of the left figure), and with the naive representation in R n(d+1) of the empirical minimizer Eq. ( 5) in green. When dealing with 1000 points, our algorithm, as well as graph-based Laplacian, can be computed in about one tenth of a second on a 2 GHz processor, while the naive kernel implementation requires 10 seconds. We show in Appendix B that this cut in costs is not associated with a loss in performance.

on various filters and assumptions of capacity, a.k.a effective dimension, and source conditions. Our setting is slightly different and showcases two specificities: (i) the eigenelements (λ i,µ , θ i,µ ) are dependent of µ; (ii) the low-rank approximation in Algorithm 1 is specific to settings with derivatives. We end our exposition with the following convergence result, proven in Appendix E. Note that the dependency of p in n can be improved based on subsampling techniques that leverage expressiveness of the different (k Xi ) [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF]. Moreover, universal consistency results could also be provided when the RKHS is dense in H 1 , as well as convergence rates for other filters and laxer assumptions which we discuss in Appendix E (in particular, the source condition can be relaxed by considering the biggest q ∈ (0, 1] such that g ∈ im L q ). Theorem 1 (Convergence rates). Under Assumptions 1, 2 and 3, for n , n ∈ N, when considering the spectral filtering Algorithm 1 with ψ λ : x → (x + λ) -1 , there exists a constant C independent of n, n , λ, µ and p such that the estimate ĝp defined in Algorithm 1 verifies

E Dn ĝp -g ρ 2 L 2 ≤ C λ 2 + λµ + σ 2 n -1 + n -2 + n -1 λµ + log(p) p + λ log(p) a p a , (10)
with σ 2 is a variance parameter that relates to the variance of the variable Y (I + λL) -1 δ X , inheriting its randomness from (X, Y ) ∼ ρ. In particular, when the ratio r = n /n is fixed, with the regularization scheme

λ n = λ 0 n -1/4 , µ n = µ 0 n -1/4
, for any λ 0 > 0 and µ 0 > 0, and the subsampling scheme p n = p 0 n s log(n) for any p 0 > 0 and with s = max( 1 /2, 1 /4a), there exists a constant C independent of n and n such that the excess of risk verifies

E Dn ĝp -g ρ 2 L 2 ≤ C (n -1/2 + σ 2 n -1/2 ). (11)
Theorem 1 answers the two questions asked at the beginning of this section. In particular, it characterizes the dependency of the need for labelled data to a variance parameter linked with the diffusion of observations (X i , Y i ) based on the density ρ X through the operator L. Intuitively, if the index set I ⊂ {1, 2, • • • , n} of data (X i ) i∈I we labelled does not change the profile of the diffusion solution ĝ, then we do not need that much labelled data -as this is the case on Figure 1.

Finally, Theorem 1 is remarkable in that it exhibits no dependency to the dimension of X in the power of n and n . This contrasts with graph-based Laplacian methods that do not scale well with the input space dimensionality [START_REF] Bengio | Label propagation and quadratic criterion[END_REF][START_REF] Hein | Graph Laplacians and their convergence on random neighborhood graphs[END_REF]. Indeed, Figure 3 shows the superiority of our method over graph-based Laplacian in dimension d = 10, with a mixture of Gaussians. We provide details as well as additional experiments in Appendix B.

Conclusion

Diffusing information or enforcing regularity through penalties on derivatives are natural ideas to tackle many machine learning problems. Those ideas can be captured informally with graph-based techniques and finite element differences, or captured more formally with the diffusion operator we introduced in this work. This formalization allowed us to shed lights on Laplacian regularization techniques based on statistical learning considerations. In order to make our method usable in practice, we provided strong computational guidelines to cut down prohibitive cost associated with a naive implementation of our methods. In particular, we were able to develop computationally efficient semi-supervised techniques that do not suffer from the curse of dimensionality.

This work paves the way to many extensions beyond semi-supervised learning. For example, in Appendix A, we describe its usefulness to the partial supervised learning problem, where minimizing the Dirichlet energy provide a learning principle, in order to bypass the restrictive non-ambiguity assumption usually made in this setup [START_REF] Cour | Learning from partial labels[END_REF][START_REF] Liu | Learnability of the superset label learning problem[END_REF][START_REF] Cabannes | Structured prediction with partial labelling through the infimum loss[END_REF][START_REF] Cabannes | Disambiguation of weak supervision with exponential convergence rates[END_REF]. Moreover, in the context of active learning, retaking the strategy of Karzand and Nowak [START_REF] Karzand | Maximin active learning in overparameterized model classes[END_REF], this energy provides a computationally-effective, theoretically-grounded, data-dependent score to select the next point to query. As such, follow-ups would be of interest to see how this introductory theoretical paper makes its way into the world of concrete applications.

Ethical considerations

This work aims at advancing our understanding of weakly supervised learning. Weakly supervised learning enrolls in the quest of an automated artificial intelligence, free from the need of human supervision. Automation, which is at the basis of computer science [START_REF] Turing | Computing machinery and intelligence[END_REF], is known to increase productivity at a reduced human labor cost, and is associated with several political/societal issues. In term of concrete applications, reducing the need for annotations is especially useful when humans reproduce biases when annotating data, or when the lack of output annotation restricts the outreach of a method (e.g., learning to translate languages by collecting input/output pairs based on books already translated by humans can hardly be applied to languages with few written resources).

A Extensions: Least-square surrogate and partially supervised learning

In this section, we first show how our work can be extended to generic semi-supervised learning problem, beyond real-valued regression. This first extension is based on the least-square surrogate introduced by Ciliberto et al. [START_REF] Ciliberto | A general framework for consistent structured prediction with implicit loss embeddings[END_REF] for structured prediction problems. We later show how our work can be extended to generic partially-supervised learning. This second extension is based on the work of Cabannes et al. [START_REF] Cabannes | Structured prediction with partial labelling through the infimum loss[END_REF].

A.

f * = arg min f :X →Y R(f ), with R(f ) = E (X,Y )∼ρ [ (f (X), Y )] . (12)
Ciliberto et al. [START_REF] Ciliberto | A general framework for consistent structured prediction with implicit loss embeddings[END_REF] showed that as soon as can be decomposed through two mappings ϕ : Y → H Y and ψ : Y → H Y with H Y a Hilbert space as (y, z) = ϕ(y), ψ(z) H Y , it is possible to leverage the least-square regression by considering the surrogate problem

g * ∈ arg min g:X →H Y E (X,Y )∼ρ g(X) -ϕ(Y ) 2 H Y . (13)
This surrogate problem relates to the original one through the decoding d that relates a surrogate estimate g : X → H Y to an estimate of the original problem f :

X → Y as f = d(g) defined through, for x ∈ supp ρ X , f (x) = arg min z∈Y ψ(z), g(x) H Y . (14)
In the real-valued regression case, presented precedently, our estimates for g n can all be written as g n (x) = n i=1 β i (x)Y i , where β i (x) is a function of the (X i ) i≤n , involving the kernel k and its derivatives. Those estimates can be cast to vector-valued regression by considering coordinateswise regression2 , which leads to g n (x) = n i=1 β i (x)ϕ(Y i ), and to the original estimates, for any x ∈ supp ρ X ,

f n (x) ∈ arg min z∈Z n i=1 β i (x) (z, Y i ). ( 15
)
The behavior of f n being independent of the decomposition (ϕ, ψ) of was referred to as the loss trick. In particular, Ciliberto et al. [START_REF] Ciliberto | A general framework for consistent structured prediction with implicit loss embeddings[END_REF] showed that convergence rates derived between g n -g * L 2

does not change if we consider g : X → R or g : X → H Y and that those rates can be cast directly as convergence rates between R(f n ) and R(f * ) with f n = d(g n ) defined by Eq.( 14). Moreover, when Y is a discrete output space, it is possible to get much better generalization bound on R(f n ) -R * by introducing geometrical considerations regarding g * and decision frontier between classes [START_REF] Cabannes | Fast rates in structured prediction[END_REF].

Example 2 (Binary classification). This framework aims at generalizing well known surrogate considerations in the case of the binary classification. Binary classification corresponds to Y = {-1, 1}, the 0 -1 loss. In this setting, H Y = R, ϕ : Y → R; y → y, and ψ = -ϕ. This definition verifies (y, z) = .5 -.5ϕ(y) ϕ(z) ϕ(y) ψ(z). This corresponds to the usual least-square surrogate, which is R S

(g) = E[ g(X) -Y 2 ], g(x) = E [Y | X = x] and f = sign g.
Beyond least-squares. Considering a least-square surrogate assumes that retrieving g * (13) is the way to solve the original problem [START_REF] Chen | Deep neural tangent kernel and Laplace kernel have the same RKHS[END_REF] and that the low-density separation hypothesis can be expressed as Assumption 1 being verified by g * . We would like to point out that the low-density separation could be expressed under a much weaker form, which is that there exists g such that f * = d(g) ( 14) and g verifies Assumption 1. In particular, the cluster assumption [START_REF] Rigollet | Generalization error bounds in semi-supervised classification under the cluster assumption[END_REF] could be understood as assuming that g = ϕ(f * ), the trivial embedding of f * in H Y , is constant on clusters, with means that g belongs to the kernel of the Laplacian operator L. Yet,

g * : x → E[ϕ(Y )|X = x],
which depends on the labelling noise, could be really non-smooth, even under the cluster assumption.

A.2 Partially supervised learning

Partial supervision is a popular instance of weak supervision, which generalizes semi-supervised learning. It has been known under the name of partial labelling [START_REF] Cour | Learning from partial labels[END_REF], superset learning [START_REF] Liu | Learnability of the superset label learning problem[END_REF], as well as learning with partial label [START_REF] Grandvalet | Logistic regression for partial labels[END_REF], with partial annotation [START_REF] Lou | Structured learning from partial annotations[END_REF], with candidate labeling set [START_REF] Luo | Learning from candidate labeling sets[END_REF] or with multiple label [START_REF] Jin | Learning with multiple labels[END_REF]. It encompasses many problems such as "classification with partial labels" [START_REF] Nguyen | Classification with partial labels[END_REF][START_REF] Cour | Learning from partial labels[END_REF], "multilabelling with missing labels" [START_REF] Yu | Large-scale multi-label learning with missing labels[END_REF], "ranking with partial ordering" [START_REF] Hüllermeier | Label ranking by learning pairwise preferences[END_REF], "regression with censored data" [START_REF] Tobin | Estimation of relationships for limited dependent variables[END_REF], "segmentation with pixel annotation" [START_REF] Verbeek | Scene Segmentation with CRFs Learned from Partially Labeled Images[END_REF][START_REF] Papandreou | Weakly-and semisupervised learning of a deep convolutional network for semantic image segmentation[END_REF], as well as instances of "action retrieval", especially on instructional videos [START_REF] Alayrac | Unsupervised learning from narrated instruction videos[END_REF][START_REF] Miech | Howto100m: Learning a text-video embedding by watching hundred million narrated video clips[END_REF].

It consists, for a given input x, in not observing its label y ∈ Y, but observing a set of potential labels s ∈ 2 Y that contains the labels (y ∈ s). Typically, if Y is the space S m of orderings between m items (e.g. movies on a streaming website), for a given input x (e.g. some feature vectors characterizing a user) s might be specified by a partial ordering that the true label y should satisfy (e.g. the user prefers romantic movies over action movies).

In this setting, it is natural to create consensus between the different sets giving information on (y|x), which has been formalized mathematically by the infimum loss (z, s) ∈ Y ×2 Y → inf y∈s (z, y) ∈ R for : Y × Y → R a specified loss on the underlying fully supervised learning problem. This leads, for τ ∈ ∆ X ×2 Y encoding the distribution generating samples (X, S), to the formulation

f * ∈ F = arg min f :X →Y E (X,S)∼τ [inf Y ∈S (f (X), Y )] .
To study this problem, a non-ambiguity assumption is usually made [START_REF] Cour | Learning from partial labels[END_REF][START_REF] Luo | Learning from candidate labeling sets[END_REF][START_REF] Liu | Learnability of the superset label learning problem[END_REF][START_REF] Cabannes | Structured prediction with partial labelling through the infimum loss[END_REF][START_REF] Cabannes | Disambiguation of weak supervision with exponential convergence rates[END_REF]. This is a very strong assumption to ensure that F is, in essence, a singleton. Highly adequate to this setting, the Laplacian regularization allows to relax this assumption, assuming that F can be big, but that we can discriminate between function in F by looking for the smoothest one in the sense defined by the Laplacian penalty. Moreover, the loss trick (15) allows to endow, in a off-the-shelf fashion, the recent work of Cabannes et al. [START_REF] Cabannes | Structured prediction with partial labelling through the infimum loss[END_REF][START_REF] Cabannes | Disambiguation of weak supervision with exponential convergence rates[END_REF] on the partial supervised learning problem with our considerations on Laplacian regularization.

B Experiments B.1 Low-rank approximation

Cutting computation cost thanks to low-rank approximation, as we did by going from the naive exact empirical risk minimizer ĝ Eq. ( 5) to the smart implementation ĝp Algorithm 1, is associated with a trade-off between computational versus statistical performance. This trade-off can be studied theoretically thanks to Theorem 1, which shows that under mild assumptions, considering p = n 1/2 log(n) does not lead to any loss in performance, in the sense that the convergence rates in n, the number of samples, are only changed by a constant factor. We show on Figure 4 that in the setting of Figure 3, our low-rank approximation is not associated with a loss in performance. Actually low-rank The fact that the error of the graph-based method stalls around one, is due to the amplitude of the estimate being very small, which is coherent with behaviors described in [START_REF] Nadler | Statistical analysis of semi-supervised learning: The limit of infinite unlabelled data[END_REF].

Trainset Kernelized estimate

Graph-based estimate approximation can even be beneficial as it tends to lower the risk for overfitting as discussed by Rudi et al. [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF].

B.2 Comparison with graph-based Laplacian

One the main goal of this paper is to make people drop graph-based Laplacian methods and adopt our "kernelized" technique. As such, we would like to discuss in more detail our comparison with graph-based Laplacian. In particular, we will discuss how and why we choose the hyperparameters and the setting of Figure 3.

The setting of Figure 3 is the one of Figure 5, we considered two Gaussians with unit variance and whose centers are at distance δ = 3 of each other. We chose Gaussians distributions as it is a well-understood setting. We chose δ = 3 so that there is an mild overlap between the two distributions.

For the bandwidth parameter, we considered σ n = Cn -1 d+4 log(n) as this is known to be the optimal bandwidth for graph Laplacian [START_REF] Hein | Graph Laplacians and their convergence on random neighborhood graphs[END_REF]. We chose C = 1 as this leads to σ n of the order of δ. We chose λ = 1 to enforce Laplacian regularization and µ n = 1/n, as this is a classical regularization parameter in RKHS. We did not cross-validate parameters in order to be fair with graph-Laplacian that do not have as much parameters as our kernel method. We compute the error in a transductive setting, retaking the exact problem and algorithm of Zhu et al. [START_REF] Zhu | Semi-supervised learning using Gaussian fields and harmonic functions[END_REF]. We choose d = 10, as we know that this is a good dimension parameter in order to illustrate the curse of dimensionality phenomenon without needing too much data. 3

( Σ + ) -1 ★ ( ˆ + ) -1 ★

B.3 Usefulness of Laplacian regularization

It is natural to ask about the relevance of Laplacian regularization. To give convergence results, we have used Assumptions 1 and 2, which imply that g * belongs to the RKHS H, and we got convergence rates in n 1/2 l , which is not better than the rates we could get with pure kernel ridge regression. In particular, our algorithm can be split between an unsupervised part that learn the penalty L 1/2 g 2 L 2 (ρ X ) and a supervised part, that solve the problem of estimating g λ from few labels (X i , Y i ) given the penalty associated to L. But the same method can be used for pure kernel ridge regression: unsupervised data could be leveraged to learn the covariance matrix Σ (6), and supervised data could be used to get S g ρ to converge towards S g ρ . The same analysis would yield the same type of convergence rates. Yet the parameter σ appearing in Theorem 1 would not be linked with the variance of Y (I + λL + λµK -1 ) -1 δ X but with the variance of Y (I + µK -1 ) -1 δ X . This is a key fact, the geometry of the covariance operator Σ is not supposed to be that relevant to the problem, while the one of L is. We illustrate this fact on Figure 6.

C Central Operators

The paper makes an intensive use of operators. This section aims at providing details and intuitions on those operators, in order to help the reader. In particular, we discuss on Assumptions 1 and 2 and we prove the equality in Eq. ( 6).

C.1 The diffusion operator

In this subsection, we discuss on the diffusion operator, and recall its basic properties.

The diffusion operator is a well-known operator in the realm of partial differential equation. Let us assume that ρ X admit a smooth density ρ X (dx) = p(x)dx, say p ∈ C 2 (R d ), that cancels outside a domain Ω ⊂ R d . Then the diffusion operator L can be explicitly written, for g twice differentiable, as

Lg(x) = -∆g(x) + 1 p(x)
∇p(x), ∇g(x) . This follows from the fact that for f once and g twice differentiable, using Stokes theorem,

f, Lg L 2 (ρ X ) = ∇f, ∇g L 2 (ρ X ) = ∇f, p∇g L 2 (dx) = X div(f p∇g) dx -f, div(p∇g) L 2 (dx) = -f, div(p∇g) L 2 (dx) = -f, p -1 div(p∇g) L 2 (ρ X ) = -f, div ∇g + p -1 (∇p).∇g L 2 (ρ X ) .
Note that when the distribution is uniform on Ω, the diffusion operator is exactly the opposite usual Laplacian operator ∆. As for the Laplacian case, it can be shown that under mild assumption on p, whose smoothness properties directly translates to the smoothness properties of the boundary of Ω, that the diffusion operator L has a compact resolvent (that is, for λ / ∈ spec(L), (L + λI) -1 is compact). This is a standard result implied by a standard version of the famous Rellich-Kondrachov compactness embedding theorem:

H 2 (Ω) is compactly injected in L 2 (Ω) whenever Ω is a bounded open with C 2 -boundaries.
In such a setting, we can consider the eigenvalue decomposition of L -1 , that is, (λ i , e i ) ∈ (R + ×L 2 ) N , with (e i ) i∈N an orthonormal basis of L 2 and (e i ) i≤dim ker L generating the null space of L, with the convention λ i = M for i ≤ dim ker L, with M an abstraction representing +∞, and (λ i ) decreasing towards zero afterwards. This decomposition reads

L -1 = i∈N λ i e i ⊗ e i . ( 16
)
Note that the fact that all the (λ i ) are positive, is due to the fact that L -1 is the inverse of a positive self-adjoint operator. As a consequence, the diffusion operator has discrete spectrum, and can be written as

L = i∈N λ -1 i e i ⊗ e i . (17)
In such a setting, the kernel-free Tikhonov regularization Eq. (3) reads

g λ = i∈N ψ(λ i ) g ρ , λ 1/2 i e i L 2 λ 1/2 i e i , (18)
with ψ : x → (x + λ) -1 , and the convention M ψ(M ) = 1.

C.2 Regularity of the eigen vectors of the diffusion operator

In this subsection, we discuss on the regularity assumed in Assumption 2.

Introducing the kernel k and its associated RKHS H is useful when the eigen vectors of L can be well approximated by functions in H. In applications, people tends to go for kernels that are translationinvariant, which implied that the RKHS H is made of smooth functions, could it be analytical functions (for the Gaussian kernel) or functions in H m (for Sobolev kernels). As a consequence, we should investigate on the regularity of those eigen vectors. Indeed, if ρ derives from a Gibbs potential, that is ρ(dx) = e -V (x) dx, the eigen vectors of L can be shown to inherit from the smoothness of the potential V [START_REF] Pillaud-Vivien | Learning with Reproducing Kernel Hilbert Spaces: Stochastic Gradient Descent and Laplacian Estimation[END_REF]. For example if V belongs to H m , and H m ⊂ H, we expect (e i ) to belongs to H, thus verifying Assumption 2.

Counter-example and beyond Assumption 2. Note that if ρ has several connected components of non-empty interiors, the null space of L is made of functions that are constants on each connected components of supp ρ X . Those functions are not analytic. In such a setting, the Gaussian kernel is not sufficient for Assumption 2 to hold, and one should favor kernel associated with richer functional space such as the Laplace kernel or the neural tangent kernel [START_REF] Chen | Deep neural tangent kernel and Laplace kernel have the same RKHS[END_REF]. However, as illustrated by Figure 2, e i not belonging to H does not mean that e i can not be well approximated by H. Indeed it is well known that the approximation power of H for e i can be measure in the biggest power p such that e i ∈ im K p [START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF], where K = SS * . Assumption 2 corresponds to p = 1/2, but it should be seen as a specific instance of more generic approximation conditions.

Handling constants in RKHS. Finally, note also that many RKHS do not contain constant functions, and therefore might not contain the constant function e 0 (although we are only looking for equality in the support of ρ X ), however this specific point with e 0 can easily be circumvent either by assuming that g ρ has zero mean, either by centering the covariance matrices Σ and Σ [START_REF] Pillaud-Vivien | Learning with Reproducing Kernel Hilbert Spaces: Stochastic Gradient Descent and Laplacian Estimation[END_REF]. This relates with the usual technique for SVM consisting in adding a unpenalized bias [START_REF] Steinwart | Support Vector Machines[END_REF].

C.3 Low-density separation

In this subsection, we discuss on how Assumption 1 relates to the idea of low-density separation.

Low-variation intuition. The low-density separation supposes that the variations of g * take place in region with low-density, so that L 1/2 g * / g * is small. As such, using Courant-Fischer principle, Assumption 1 can be reformulated as g * belonging to the space

Span {e i } i≤r = arg min F ⊂L 2 ; dim F =r max g∈F L -1/2 f 2 L 2 f 2 L 2 .
In other terms, Assumption 1 can be restated as g * belonging to a finite dimensionnal space that minimizes a measure of variation given by the Dirichlet energy.

To tell the story differently, suppose that we are in a classification setting, i.e. Y ∈ {-1, 1}, and that the supp ρ X is connected. Then we know that the null space of L is made of constant functions. Then the first eigen vector e 2 of L is a function that is orthogonal to constants. Hence e 2 is a function that changes its sign and that is "balanced" in the sense that

E[e 2 ] = 0 -i.e. if e 2 (x) = E µ [Y | X = x]
for some measure µ, we have E µ [Y ] = 0, meaning that classes are "balanced". Moreover, in order to minimize L 1/2 e 2 , the variations of e 2 should take place in low-density regions of X .

Diffusion intuition. Finally, as L is a diffusion operator, we also have an interpretation of Assumption 1 is term of diffusion. Consider (λ i , e i ) the eigen elements of Eq. ( 17). The diffusion of g ρ according the density ρ X can be written as, for t ∈ R,

g t = e -tL g ρ = i∈N e -tλ -1 i g ρ , e i e i .
This diffusion will cut off the high frequencies of g ρ that corresponds to g ρ , e i for big i, and big λ -1 i . Indeed, the difference between the diffusion and the original g ρ can be measured as

g t -g ρ 2 L 2 = i∈N (e -tλ -1 i -1) 2 g ρ , e i 2 = i∈N t 2 λ -2 i g ρ , e i 2 + o(t 2 λ -2 i ).
So that assuming that g ρ is supported on few of the first eigen vectors of L, can be rephrased as saying that the diffusion of g ρ does not modify it too much.

The variance σ . Theorem 1 shows that the need for labels depends on the variance parameter σ 2 .

It is natural to wonder on how this parameters relates to the low-density hypothesis. As we discussed, this parameter is linked to the variance of Z = Y (I + λL) -1 δ X . We can separate the variability of this variable due to X and the variability due to

Y Z = Z X + Z Y , with Z X = (I + λL) -1 g ρ (X)δ X , Z Y = (I + λL) -1 (Y -E[Y | X])δ X .
As such we see that this variance depends on the structure of the density ρ X with the variance of (I + λL) -1 δ X , and the labelling noise with the variance of (Y | X). The low-density separation does not tell us anything about the level of noise in Y or the diffusion structure linked with ρ X , but additional hypothesis could be made to characterize those.

C.4 Kernel operators

In this subsection, we define formally the operators S and Σ.

We now turn towards operators linked with the Hilbert space H. Recall that for k : X → X → R a kernel, H is defined the closure of the span of the elements k x under the scalar product k

x , k x = k(x, x ). In particular, k x 2 H = k(x, x). H parameterize a vast class of function in R X through the mapping S : H → R X θ → ( k x , θ ) x∈X .
Under mild assupmtions, S maps H to a space of function belongs to L 2 . Proposition 2. When x → k(x, x) belongs to L 1 (ρ X ), S is a continuous mapping from H to L 2 (ρ X ). This is particularly the case when ρ X has compact support and k is continuous.

Proof. Consider θ ∈ H, we have Sθ 2 L 2 = X k x , θ 2 ρ X (dx) ≤ X k x , θ 2 H ρ X (dx) ≤ X k x 2 H θ 2 H ρ X (dx) = θ 2 H X k(x, x)ρ X (dx) = θ 2 H x → k(x, x) L 1 .
Moreover, when ρ X has compact support and k is continuous, k is bounded on the support of ρ X therefore x → k(x, x) belongs to L 1 .

As a continuous operator from the Hilbert space H to the Hilbert space L 2 , S is naturally associated with many linear structure. In particular its adjoint S , but also the self-adjoint operators K = SS and Σ = S S. Proposition 3. The adjoint of S is defined as

S : L 2 → H g → X g(x)k x ρ X (dx) = E X∼ρ X [g(X)k X ].
To S is associated the kernel self-adjoint operator on

L 2 K := SS : L 2 → L 2 g → (x → X k(x, x )g(x )ρ X (dx )),
as well as the (not-centered) covariance on H, Σ :

= S S = E X∼ρ X [k X ⊗ k X ].
Proof. We shall prove the equality defining those operators. Consider θ ∈ H and g ∈ L 2 , we have

S g, θ H = g, Sθ L 2 = E X∼ρ X [g(X) k X , θ H ] = E X∼ρ X [g(X)k X ], θ H . We also have, for x ∈ X , (SS g)(x) = k x , E X∼ρ X [g(X)k X ] H = E X∼ρ X [g(X) k x , k X H ] = E X∼ρ X [g(X)k(X, x)]. Finally, we have S Sθ = E X∼ρ X [Sθ(X)k X ] = E X∼ρ X [ θ, k X H k X ] = E X∼ρ X [k X ⊗ k X ]θ.
This provides the last of all the equalities stated above.

The functional space H. In the main text, we have written everything in term of θ, highlighting the parametric nature of kernel methods. This made it easier to dissociate the norm on functions derived from H and the one derived from L 2 or H 1 . In literature, people tends to keep everything in term of functions g θ = Sθ without even mentioning the dependency in θ. Such a setting consists in considering directly the RKHS H whose scalar product is defined for g, g ∈ (ker K) ⊥ by g, g H = g, K -1 g L 2 .

C.5 Derivative operators

In this subsection, we discuss on derivative in RKHS and we define formally the operator L.

As well as evaluation maps can be represented in H, under mild assumptions, derivative evaluation maps can benefited of such a property. Indeed, for g θ = Sθ, x ∈ X and u ∈ B X (0, 1) a unit vector, we have

∂ u g θ (x) = lim t→0 g θ (x + tu) -g θ (x) t = lim t→0 θ, k x+tu H -g, k x H t = lim t→0 θ, k x+tu -k x t H
As a linear combination of elements in H, the difference quotient evaluation map t -1 (k x+tu -k x ) belongs to H and has a norm

k x+tu -k x t 2 H = k(x + tu, x + tu) -2k(x + tu, x) + k(x, x) t 2 .
In order for the limit when t goes to zero to belong to H, we see the importance of k to be twice differentiable. This limit ∂ u k x , whose existence is proven formally by Zhou [START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF], provides a derivative evaluation map in the sense that

∂ u g θ (x) = θ, ∂ u k x H .
From this equality, we derive that

∂ 1i k(x, x ) = k x , ∂ i k x , and recursively that ∂ i k x , ∂ j k x = ∂ 1i ∂ 2j k(x, x ).
Similarly to the operator S, we can introduce the operators Z i for i ∈ [1, d], defined as

Z i : H → R X θ → ( ∂ i k x , θ H ) i≤d .
Once again, under mild assumptions, im Z i inherit from an Hilbertian structure.

Proposition 4. When x → ∂ 1i ∂ 2i k(x, x) belongs to L 1 (ρ X ), Z i is a continuous mapping from H to L 2 (ρ X )
. This is particularly the case when ρ X has compact support and k is twice differentiable with continuous derivatives.

Proof. Consider θ ∈ H, similarly to before, we have

Zθ 2 L 2 = X ∂ i k x , θ 2 ρ X (dx) ≤ θ 2 H X ∂ i k x 2 H ρ X (dx) = θ 2 H x → ∂ 1,i ∂ 2,i k(x, x) L 1 .
Moreover, when ρ X has compact support and

∂ 1,i ∂ 2,i k is continuous, ∂ 1,i ∂ 2,i k is bounded on the support of ρ X therefore x → ∂ 1,i ∂ 2,i k belongs to L 1 .
Among the linear operator that can be build from Z i , in the theoretical part of this paper, we are mainly interested in Z i Z i . In the empirical part however, we might be interested in Z i Z j as well as Z i S as it might appear in Algorithm 1 (where Z n has to be understood as the empirical version of

Z = [Z 1 ; • • • ; Z d ]).
Proposition 5. The energy Dirichlet on H can be represented through the operator

S LS = d i=1 Z i Z i = d i=1 E X∼ρ X [∂ i k X ⊗ ∂ i k X ].
Proof. Let θ ∈ H and g θ = Sθ, we have

g θ , Lg θ L 2 = θ, S LSθ H = E X∼ρ X ∇g θ (X) 2 = d i=1 E X∼ρ X (∂ i g θ (X)) 2 = d i=1 E X∼ρ X ∂ i k X , θ 2 H = d i=1 Z i θ 2 L 2 = θ, d i=1 Z i Z i θ H = d i=1 E X∼ρ X [ θ, (∂ i k X ⊗ ∂ i k X )θ H ] = θ, d i=1 E X∼ρ X [∂ i k X ⊗ ∂ i k X ] θ H .
Since the three operators are self-adjoint and they all represent the same quadratic form, they are equals.

C.6 Relation between Σ and L

In this subsection, we discuss on the relation between Σ and L and show that we can expect the existence of a ∈ (1 -2/d, 1] and c > 0 such that L cΣ a .

Informal capacity considerations. We want to compare Σ and L, as L cΣ a with the biggest a possible. This depend on how fast the eigen values are decreasing, which is linked to the entropy numbers of those two compact operators. Those entropy numbers are linked with the capacity of the functional spaces g

∈ L 2 K -1/2 g L 2 < ∞ and g ∈ L 2 K -1/2 L -1/2 g L 2 < ∞ .
The first space is the reproducing kernel Hilbert space linked with k, the second space is, roughly speaking, a space of function whose integral belongs to the first space. As such, if the first space is H m , the second is H m-1 , and we can consider a = (m -1)/m. Because we are considering kernel, we have m > d/2 (this to make sure that the evaluation functionals L X : f → f (x) are continuous), so that a > 1 -2/d. Without trying to make those "algebraic" considerations more formal, we will give an example on the torus.

Translation-invariant kernel and Fourier transform. Consider L 2 ([0, 1] d , dx) the space of periodic functions in dimension d, square integrable against the Lebesgue measure on [0, 1] d . For simplicity, we will suppose that ρ X is the Lebesgue measure on [0, 1] d . Consider a translation invariant kernel k(x, y) = q(x -y) for q : R d → R that is one periodic. In this setting, the operator K, operating on L 2 , is the convolution by q, that is

K : L 2 → L 2 g → q * g , hence Kg = qĝ.
Where we have used the fact that convolutions can be represented by a product in Fourier. Note that, from Böchner theorem, we know that k being positive definite implies that the Fourier transform of q exists and is not negative. Let us define the Fourier coefficient and the inverse Fourier transform as

∀ ω ∈ Z d , ĝ(ω) = [0,1] d g(x)e -2iπω x dx, and ∀ x ∈ [0, 1] d , g(x) = ω∈Z d e 2iπω x ĝ(ω).
K being a convolution operator, it is diagonalizable with eigen elements (q(ω), x → e 2iπω x ) ω∈Z d .

From this, we can explicit many of our abstract operators. First of all, using Perceval's theorem,

g 2 H = g, K -1 g L 2 = ω∈Z d |ĝ(ω)| 2 q(ω) .
Hence we can parametrize H with (θ ω ) ω∈Z d ∈ C Z d and the 2 -metric, where θ ω = ĝ(ω)/ q(ω) and

(Sθ)(x) = g θ (x) = ω∈Z d e 2iπω x q(ω)θ ω .
Note that this is not the usual parameterization of H by elements θ ∈ H as (C Z d , 2 ) is not a space of functions. However, such a parametrization of H does not change any of the precedent algebraic considerations on the operators S, Σ, K, and L.

Diffusion operator and Fourier transform. As well as convolution operators are well represented in Fourier, derivation operators are. Indeed, when g is regular, we have

L 1/2 g 2 L 2 = ∇g 2 L 2 = d j=1 ∂ j g 2 L 2 = d j=1 ω∈Z d ω 2 j |ĝ(ω)| 2 .
As a consequence, using the expression of Sθ, we have

Σθ = ω∈Z d q(ω)θ ω , while Lθ = ω∈Z d ω 2 2 q(ω)θ ω , where ω 2 2 = d j=1 ω 2 i .
With this parameterization, the eigen elements of Σ are (q(ω), δ ω ) ω∈Z d while the one of L are

( ω 2 2 q(ω), δ ω ) ω∈Z d .
Eigen value decay comparison. Hence, having L cΣ a is equivalent to having ω 2 2 q(ω) ≤ cq(ω) a . Now suppose that the decay of q is governed by

c 1 (1 + σ -1 ω 2 2 ) -m ≤ q(ω) ≤ c 2 (1 + σ -1 ω 2 2 ) -m , for two constants c 1 , c 2 > 0.
In particular, this is verified for Matérn kernels, corresponding to the fractional Sobolev space H m , and for the Laplace kernel with m = (d + 1)/2, which reads k(x, y) = exp(-σ -1 x -y ). The Gaussian kernel could be seen as m = +∞ as it has exponential decay. With such a decay we have, assuming without restrictions that we are in one dimension

ω 2 q(ω) ≤ c 2 ω 2 (1 + σ -1 ω 2 ) -m ≤ c 2 σ(1 + σ -1 ω 2 ) -(m-1) ≤ c m m-1 1 c 2 σ q(ω) m-1 m .
In other terms, we can consider c = c m m-1 1 c 2 σ and a = (m -1)/m. Assuming that q is squareintegrable, so is q, which implies that 2m > d. As a consequence, we do have a > 1 -2/d. Note that this reasoning could be extended to the case where ρ X has a density against the Lebesgue measure, that is bounded above and below away from zero.

D Spectral decomposition

In this section, we recall facts on spectral regularization, before proving Proposition 1 and extending it to the case µ = 0.

D.1 Generalized singular value with matrices

Generalized singular value decomposition. Let A ∈ R m1×n and B ∈ R m2×n be two matrices. There exists U ∈ R m1×m1 , V ∈ R m2×m2 two orthogonal matrices, c ∈ R m1×r and s ∈ R m2×r be two 1-diagonal matrices such that c c + s s = I r , and H ∈ R n×r a non-singular matrix such that

A = U cH -1 , B = V sH -1 .
To be more precise c is such that only entries c ii = cos(θ i ) for i < min(r, m 2 ) are non-zeros and s such that only entries s m1-i,r-i = sin(θ r-i ) for i < min(r, m 1 ) are non-zeros, with θ i ∈ [-π/2, π/2] an angle. Here, c stands for cosine, s for sinus and r for rank.

Link with generalized eigenvalue problem. As well as the singular value of A is linked with the eigen value of A A, the generalized singular value decomposition of [A; B] is linked with the generalized eigenvalue problem linked with (A A, B B). Indeed, we have

A A = H -c cH -1 , B B = H -s sH -1 .
In particular, with (e i ) the canonical basis of R r , and h i the i-th column of H, we get

H A Ah i = cos(θ i ) 2 e i = tan(θ i ) -2 sin(θ i ) 2 e i = tan(θ i ) -2 H B Bh i . From which we deduce that, since im A ∪ im B ⊂ im H , A Ah i = tan(θ i ) -2 B Bh i , h j B Bh i = sin(θ i ) 2 1 i=j . So if we denote by f i = |sin(θ i )| -1 h i and λ i = tan(θ i ) -2
, assuming λ i = 0 for all i ≤ r (which corresponds to ker B ⊂ ker A), (λ i ) i≤r , (f i ) i≤r provide the generalized eigenvalue decomposition of (A A, B B) in the sense that

A Af i = λ i B Bf i , f j B Bf i = 1 i=j , f j A Af i = λ i 1 i=j .

D.2 Tikhonov regularization

Define the Tikhonov regularization

x λ = arg min x∈R n Ax -b 2 + λ Bx 2 .
When this problem is well-defined, the solution is defined as

x λ = (A A + λB B) † A b.
With the generalized singular value decomposition of A and B, we have

A A + λB B = H -γ λ H -1 , with γ λ = c c + λs s.
Using the fact that A b = H -c U b, we get

x λ = Hγ -1 λ c U b = r i=1 cos(θ i ) cos(θ i ) 2 + λ sin(θ i ) 2 h i ⊗ u i b.
Now, we would like to replace c ii , s ii , h i and u i with quantities that depend on λ i , f i and A.

To do so recall that AH = U c, therefore cos(θ i )u i = Ah i , and recall that h i = sin(θ i )f i and λ i = cos(θ i ) 2 / sin(θ i ) 2 . Inputting those equality in the last expression of x λ we get

x λ = r i=1 sin(θ i ) 2 cos(θ i ) 2 + λ sin(θ i ) 2 f i ⊗ Af i b = r i=1 1 λ i + λ f i ⊗ Af i b.
Finally,

b λ = Ax λ = r i=1 ψ(λ i ) Af i , b Af i , where ψ(x) = 1 x + λ .

D.3 Extension to operators

To end the proof of Proposition 1, we should prove that we can apply the generalized eigenvalue decomposition to operators. We will only prove that it is possible for (Σ, L + µ) based on simple considerations.

Proposition 6. When k is continuous and supp ρ X is bounded, Σ is a compact operator.

Proof. We have Σ = E[k X ⊗ k X ] and k x = k(x, x). Since k is continuous and supp ρ X is compact, for x ∈ supp ρ X , k(x, x) is bounded. Hence Σ is a nuclear operator, hence trace class and compact.

Proposition 7. When k is twice differentiable with continuous derivative, and supp ρ X is compact, L is a compact operator. As a consequence, L has a compact spectrum, and has a pseudo inverse that we will denote, with a slight abuse of notation, by L -1 .

Proof. The proof is similar to the one showing that Σ is compact, based on the fact that L =

d i=1 E[∂ i k X ⊗ ∂ i k X ], and
∂ i k X 2 = ∂ 1i ∂ 2i k(x, x). Proposition 8. When Σ is compact, for all µ > 0, (L + µ) -1/2 Σ(L + µ) -1/2 is a compact operator.
Proof. The proof is straightforward

Tr((L+µ) -1/2 Σ(L+µ) -1/2 ) = Tr(Σ(L+µ) -1 ) ≤ (L + µ) -1 op Tr(Σ) ≤ µ -1
Tr(Σ) < +∞. Therefore the operator is trace class, hence compact. Proposition 9. For any µ > 0, the generalized eigen value decomposition of (Σ, L + µ) as defined in Proposition 1 exists.

Proof. Using the spectral theorem, since (L + µ) -1/2 Σ(L + µ) -1/2 is positive self-adjoint compact operator, there exists (ξ i ) a basis of H and (λ i ) ∈ R + a decreasing sequence (note that ker(L + mu) = ker Σ = {0}), such that

(L + µ) -1/2 Σ(L + µ) -1/2 = i∈N λ i ξ i ⊗ ξ i . Taking θ i = (L + µ) -1/2 ξ i , we get Σθ i = λ i Lθ i . Because (ξ i ) generates H, and (L + µ) -1/2 is bijective (since L is compact, (L + µ) -1 is coercive), ((L + µ) -1/2 ξ i ) generates H.
Proposition 1 follows from prior discussion on Tikhonov regularization extended to infinite summations.

D.4 The case µ = 0

When µ = 0, Eq. ( 7) should be seen as the rewritting of Eq. ( 18) based on the RKHS G = im S. This can only be done when the eigen vectors of L appearing in Eq. ( 17) belongs to G = im S = im K 1/2 , which is exactly what Assumption 2 provides. In such a setting, we can find (θ i ) ∈ H N to write λ 1/2 i e i = Sθ i as soon as λ i = 0 (write M e i = Sθ i for M an abstraction representing +∞ when λ i = 0, handling the potential fact that ker B ⊂ ker A), we get θ i S * Sθ j = λ i 1 i=j , and

Lθ i = λ -1
i Σθ i , and we can extend Proposition 1 to the case µ = 0, with

g λ = i∈N ψ(λ i ) S g ρ , θ i Sθ i , (19)
where we handle the null space of L with the equality M ψ(M ) = 1, verified by M our abstraction representing +∞, so that ψ(M ) S g ρ , θ i Sθ i = g ρ , e i e i as soon as λ i = 0.

Beyond Assumption 2. Assumption 2 could be made generic by considering the biggest (p i ) ∈ R N + such that K -pi e i belongs to L 2 , and rewriting Eq. ( 19) under the form

g λ = i∈N ψ(λ i ) (S 0 K pi ) g ρ , θ i S 0 K pi θ i , with θ i = λ -1/2 i
S -1 0 K -pi θ i and S 0 = K -1/2 S the isomorphism between H and L 2 (assuming that S is dense in L 2 ). Such an assumption would describe all situations from no assumption (p i = 0 for all i), Assumption 2 (p i = 1/2 for all i) to even more optimistic assumptions (p i ≥ 1 for all i).

E Consistency analysis

This section is devoted to the proof of Theorem 1. The proof is based on Eqs. ( 7) and [START_REF] Eriksson | Scaling Gaussian process regression with derivatives[END_REF], and splits the error of g ρ -ĝp 2 L 2 into several components linked with how well we approximate S g ρ , and how well we approximate the eigenvalue decomposition (λ i , θ i ) of (Σ, L).

E.1 Sketch and understanding of the proof

In this subsection, we explain how work the proofs for consistency theorems such as Theorem 1.

Let us define the mapping G : H × C → L 2 with C the set of pairs of self-adjoint operators on H that admit a generalized eigen value decomposition, as

G(θ, (A, B)) = i∈N ψ(λ i ) θ, θ i Sθ i with (λ i , θ i ) ∈ GEVD(A, B). (20)
G(θ, (A, B)) ∈ L 2 corresponds to writing θ ∈ H in the basis associated with the generalized eigen value decomposition (GEVD) of (A, B).

Proposition 10. Under Assumptions 1 and 2, and with ψ defined in Theorem 1 g λ = G(S g ρ , (Σ, L)), and ĝp = G( Ŝ g ρ , (P ΣP, P LP + µP )),

with P the projection matrix from H to Span {k Xi } i≤p .

Proof. This is a direct application of Assumptions 1, 2, Eq. ( 19) and Algorithm 1.

The main point of the proof is to relate g ρ to ĝp . To do so, we will use several functions in L 2 generated by G. We detail our steps in Table 2. Table 2 gives a first answer to the two questions asked in the opening of Section 4. The number of unlabelled data controls the convergence of the operators (P ΣP, P LP + µP ) towards (Σ, L + µ). The number of labelled data controls the convergence of the vector S g ρ towards S g ρ . Priors on the structure of the problem, such as source and approximation conditions, control the convergence of the bias estimate g λ,µ towards g ρ . Furthermore, a more precise study reveals that the concentration of operators are related to efficient dimension [START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF] and are accelerated by capacity assumptions on the functional space whose norm is g = (L + K -1 ) 1/2 g L 2 , and that the concentration of the vector S g ρ is accelerated by assumptions on moments of the variable Y (I + λL) -1 δ X (inheriting randomness from (X, Y ) ∼ ρ).

Low-rank approximation [45] ↓ ĝ S g ρ ( Σ, L + µ)
Concentration for self-adjoint operators [START_REF] Minsker | On some extensions of Bernstein's inequality for self-adjoint operators[END_REF] ↓

g n S g ρ (Σ, L + µ)
Concentration for vector in Hilbert space [START_REF] Vladimirovich | On an infinite-dimensional version of S. N. Bernstein's inequalities[END_REF] ↓

g λ,µ S g ρ (Σ, L + µ)
Bias controlled with source condition [START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF][START_REF] Lin | Optimal rates for spectral algorithms with least-squares regression over Hilbert spaces[END_REF] ↓

g λ S g ρ (Σ, L)
↓ Bias controlled with source condition [START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF][START_REF] Lin | Optimal rates for spectral algorithms with least-squares regression over Hilbert spaces[END_REF] g ρ Control of biases. We begin our study in a downward fashion regarding Table 2. Indeed, for Tikhonov regularization Eq. ( 3), we can show that for q ∈ [0, 1],

g λ -g ρ L 2 ≤ λ q L q g ρ L 2 .
Meaning that if we have the source condition g ρ ∈ im L q (which is a condition on how fast ( g ρ , e i ) i∈N decreases compared to (λ i ) i∈N for (λ i , e i ) the eigen value decomposition of L -1 ), the rates of convergence of this term when n goes to infinity is controlled by the regularization scheme λ q n . In a similar fashion to the kernel-free bias above, for (q i ) ∈ (0, 1) N , we can have

g λ,µ -g λ 2 L 2 ≤ 2 i∈N λ 2qi µ 2qi λ i λ + λ i 2 | e i , g ρ | 2 K -qi e i 2 L 2 .
This shows explicitly the usefulness of controlling at the same time how g ρ is supported on the eigen spaces of L and how the eigen vectors are well approximated by the RKHS H, which can be read in the value of (q i ) such that all e i ∈ im K qi .

Vector concentration. Let us now switch to concentration of S g ρ = n -1 n i=1 Y i k Xi towards S g ρ = E (X,Y )∼ρ [Y k X ], it will allow to control g n -g λ,µ 2 L 2 with the notations appearing in Table 2. Note how this convergence should be measured in term of the reconstruction error i∈N ψ(λ i,µ ) S g ρ -S g ρ , θ i,µ Sθ i,µ L 2 . This error might behave in a must better fashion than the L 2 error between SS g ρ and S S g ρ . In particular, on Figure 1, we can consider ψ(λ i,µ ) = 0 for i > 4, and we might have Y k X , θ i,µ = Y 1 x∈Ci , for i ≤ 4 and (X, Y ) ∈ X × Y, where C i is the i-th innermost circle. In this setting, when all four (Y | X ∈ C i ) are deterministic, we only need one labelled point per circle to clear the reconstruction error. Based on concentration results in Hilbert space, when |Y | is bounded by a constant c Y , and x → k(x, x) by a constant κ 2 , we have, with D n ∼ ρ ⊗n the dataset generating the labelled data

E Dn g n -g λ,µ 2
L 2 ≤ 2σ 2 (µλn ) -1 + 4 9 c 2 Y κ 2 (µλn 2 ) -1 .
where σ 2 ≤ c 2 Y Tr(Σ) is a variance parameter to relate to the variance of Y (I + λL) -1 δ X (where the randomness is inherited from (X, Y ) ∼ ρ). The fact that the need for labelled data depends on the variance of (X, Y ) after being diffused through L is coherent with the results obtained by Lelarge and Miolane [START_REF] Lelarge | Asymptotic Bayes risk for Gaussian mixture in a semisupervised setting[END_REF] in the specific case of a mixture of two Gaussians.

Basis concentration. We are left with the comparison of g n , which is the filtering of S g ρ with the operators (Σ, L + µ), and ĝp , which is the filtering of the same vector with the operators (P ΣP, P LP + µP ). As the number of samples grows towards infinity, we know that (P ΣP, P LP + µP ) will converge in operator norm towards (Σ, L + µ). Yet, how to leverage this property to quantify the convergence of ĝp towards g n ? Let us write (λ i , θ i ) = GEVD(Σ, L + µ), and

(λ i , θ i ) = GEVD(P ΣP, P LP + µP ), we have ĝp -g n L 2 = i∈N ψ(λ i ) θ i , θρ Sθ i -ψ(λ i ) θ i , θρ Sθ i with θρ = S g ρ .
The generic study of this quantity requires to control eigenspaces one by one. Note that we expect convergence of eigenspaces to depend on gaps between eigenvalues. However, when considering Tikhonov regularization, this quantity can be written under a simpler form. In particular, the concentration of operators is controlled, up to few leftovers, through the quantity

(Σ + λL + µλ) -1/2 ((Σ -Σ) + λ(L -L))(Σ + λL + µλ) -1/2
op where • op designs the operator norm. In this setting, the low-rank approximation is controlled through (Σ + λL) 1/2 (I -P ) op , and when L cΣ α , this term can be controlled by Σ 1/2 (I -P ) op + λ 1/2 Σ 1/2 (I -P ) α op which can be controlled based on the work of Rudi et al. [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF].

E.2 Risk decomposition

In this subsection, we decompose the risk appearing in Theorem 1.

E.2.1 Control of biases

We begin by splitting the error g ρ -ĝp L 2 between a bias term due to the regularization parameters and a variance term due to the data. With the notation of Table 2,

g ρ -ĝp L 2 ≤ g ρ -g λ L 2 + g λ -g λ,µ L 2 + g λ,µ -ĝp L 2 . ( 21
)
We will control the first two terms here, and the last term in the following subsections.

Proposition 11 (Bias in λ). Under Assumption 1

g λ -g ρ L 2 ≤ λ Lg ρ L 2 . ( 22
)
Proof. Based on the definition of g λ = (I + λL) -1 g ρ , we have g λ -g ρ = ((I + λL) -1 -I)g ρ = -λLg ρ .

Because g ρ is supported on the first eigen vectors of the Laplacian (Assumption 1), we have g ρ = r i=1 g ρ , e i e i , with e i the eigen vector of L appearing in Eq. ( 17), and

Lg ρ 2 L 2 = r i=1 λ -1 i g ρ , e i e i 2 L 2 = r i=1 λ -2 i g ρ , e i 2 ≤ λ -2 r g ρ 2 L 2 < +∞.
This ends the proof of this proposition.

Proposition 12 (Bias in µ). Under Assumptions 1 and 2, we have

g λ,µ -g λ 2 L 2 ≤ λµc 2 a g ρ 2 L 2 , with c 2 a = r i=1 K -1/2 e i 2 L 2 = r i=1 e i 2 H . (23)
Proof. Before diving into the proof, recall that the RKHS norm penalization can be written as

g H = K -1/2 g L 2 . Using the fact that A -1 -B -1 = A -1 (B -A)B -1 , we have g λ,µ -g λ = ((I + λL + µλK -1 ) -1 -(I + λL) -1 )g ρ = -(I + λL + µλK -1 ) -1 λµK -1 (I + λL) -1 g ρ = -(λµ) 1/2 (I + λL + µλK -1 ) -1/2 (I + λL + µλK -1 ) -1/2 • • • × (λµK -1 ) 1/2 K -1/2 (I + λL) -1 g ρ .
As a consequence,

g λ,µ -g λ 2 L 2 ≤ λµ K -1/2 (I + λL) -1 g ρ 2 L 2
, where we used the fact that I + λL + µλK -1 I, so that (I + λL + µλK -1 ) -1/2 op ≤ 1 (with • op the operator norm), and that

(I + λL + µλK -1 ) -1/2 (λµK -1 ) 1/2 2 op = λµ K -1/2 (I + λL + µλK -1 ) -1 K -1/2 op = λµ (K + λK 1/2 LK 1/2 + µλ) -1 op ≤ 1.
We continue the proof with

K -1/2 (I + λL) -1 g ρ = r i=1 λ i λ + λ i g ρ , e i K -1/2 e i ≤ r i=1 λ i λ + λ i | g ρ , e i | K -1/2 e i ≤ r i=1 | g ρ , e i | K -1/2 e i L 2 ≤ g ρ L 2   i≤r K -1/2 e i 2 L 2   1/2 .
Putting all the pieces together ends the proof.

E.2.2 Vector concentration

We are left with the study of the variance ĝp -g λ,µ . To ease derivations, we denote C = Σ + λL, Ĉ = Σ + λ L, θ ρ = S g ρ , θρ = S g ρ and P the projection from H to Span {k Xi } i≤p . We have, for Tikhonov regularization ĝp -g λ,µ L 2 = S P (P ĈP + λµ)

-1 P θρ -(C + λµ) -1 θ ρ L 2 = Σ 1/2 P (P ĈP + λµ) -1 P θρ -(C + λµ) -1 θ ρ H .
We begin by isolating the dependency to labelled data ĝp -g λ,µ L 2 ≤ Σ 1/2 P (P ĈP + λµ) -1 ( θρ -θ ρ )

H • • • + Σ 1/2 P (P ĈP + λµ) -1 θ ρ -(C + λµ) -1 θ ρ H . (24)
We will control the first term here, and the second term in the following subsection.

Lemma 13 (Vector term). When (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1/2 op ≤ 1/2, we have

Σ 1/2 P (P ĈP + λµ) -1 P ( θρ -θ ρ ) H ≤ 2 (C + λµ) -1/2 ( θρ -θ ρ ) H . (25)
Proof. We begin with the splitting

Σ 1/2 P (P ĈP + λµ) -1 P ( θρ -θ ρ ) H ≤ Σ 1/2 P (P ĈP + λµ) -1 P (C + λµ) 1/2 op • • • × (C + λµ) -1/2 ( θρ -θ ρ ) H .
The first term will concentrate towards a matrix smaller than identity, while the second term concentrates towards zero. We can make those considerations more formal. Following basic properties with the Löwner order on operators, we have

(C + λµ) -1/2 (C -Ĉ)(C + λµ) -1/2 t ⇒ Ĉ (1 -t)C -tλµ ⇒ P ĈP (1 -t)P CP -tλµP (1 -t)P CP -tλµ ⇒ P ĈP + λµ (1 -t)(P CP + λµ) ⇒ (P ĈP + λµ) -1 (1 -t) -1 (P CP + λµ) -1
⇒ (C + λµ) 1/2 P (P ĈP + λµ) -1 P (C + λµ) 1/2

(1 -t) -1 (C + λµ) 1/2 P (P CP + λµ)

-1 P (C + λµ) 1/2 (1 -t) -1 ,
where we have used the fact that the last operator is a projection. As a consequence, for any t ∈ (0, 1), we have

(C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1/2 op ≤ t ⇒ (C + λµ) 1/2 P (P ĈP + λµ) -1 (C + λµ) 1/2 op ≤ (1 -t) -1 . ⇒ Σ 1/2 P (P ĈP + λµ) -1 P (C + λµ) 1/2 op ≤ (1 -t) -1 .
Where the last implication follows from the fact that C + λµ = Σ + λL + λµ Σ.

E.2.3 Basis concentration

We are left with the study of the basis concentration with the number of unlabelled data.

Lemma 14 (Basis term).

When (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1/2 op ≤ 1/2, we have Σ 1/2 (P (P ĈP + λµ) -1 -(C + λµ) -1 )θ ρ H ≤ 3 C 1/2 (I -P ) op g λ H + 2 (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1 θ ρ H . (26)
Notice that Assumptions 1 and 2 imply g λ H ≤ c a g ρ L 2 < +∞.

Proof. First of all, using that A

-1 -B -1 = A -1 (B -A)B -1 , notice that Σ 1/2 (P (P ĈP + λµ) -1 -(C + λµ) -1 )θ ρ H = Σ 1/2 P (P ĈP + λµ) -1 P (C -ĈP )(C + λµ) -1 θ ρ -Σ 1/2 (I -P )(C + λµ) -1 θ ρ H ≤ Σ 1/2 P (P ĈP + λµ) -1 P ( ĈP -C)(C + λµ) -1 θ ρ H + Σ 1/2 (I -P )(C + λµ) -1 θ ρ H ≤ Σ 1/2 P (P ĈP + λµ) -1 P (C + λµ) 1/2 op (C + λµ) -1/2 P ( ĈP -C)(C + λµ) -1 θ ρ H • • • + Σ 1/2 (I -P ) op (C + λµ) -1 θ ρ H .
Because Σ Σ + λL = C, we have Σ 1/2 (I -P ) op ≤ C 1/2 (I -P ) op , and we also have

(C + λµ) -1 θ ρ H ≤ C -1 θ ρ H = K -1/2 SC -1 θ ρ H = K -1/2 g λ L 2 = g λ H .
Recall, that, for any t ∈ (0, 1), we have already shown that

(C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1/2 op ≤ t ⇒ Σ 1/2 P (P ĈP + λµ) -1 P (C + λµ) 1/2 op ≤ (1 -t) -1 .
We are left with one last term to work on

(C + λµ) -1/2 P ( ĈP -C)(C + λµ) -1 θ ρ H ≤ (C + λµ) -1/2 P ( Ĉ -C)P (C + λµ) -1 θ ρ H • • • + (C + λµ) -1/2 C(I -P )(C + λµ) -1 θ ρ H .
We control the first term with the fact for A, B, C three self-adjoint operators and x a vector we have

AP BP Cx = AP BP Cx ⊗ xCP BP A 1/2
op , and that

P Cx ⊗ xCP Cx ⊗ xC ⇒ P BP Cx ⊗ xCP BP BP Cx ⊗ xCP B BCx ⊗ xCB ⇒ AP BP Cx ⊗ xCP BP A ABCx ⊗ xCBA, so that (C + λµ) -1/2 P ( Ĉ -C)P (C + λµ) -1 θ ρ H ≤ (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1 θ ρ H .
We control the second term with

(C + λµ) -1/2 C(I -P )(C + λµ) -1 θ ρ ≤ (C + λµ) -1/2 C 1/2 C 1/2 (I -P ) (C + λµ) -1 θ ρ .
Using that (C + λµ) -1/2 C 1/2 I, we can add up everything to get the lemma.

For the part concerning g λ H , notice that

g λ H = K -1/2 g λ L 2 = r i=1 λ i λ i + λ g ρ , e i K -1/2 e i L 2 ≤ r i=1 |g ρ | e i K -1/2 e i ≤ g ρ L 2 d i=1 K -1/2 e i 2 L 2 1/2 = c a g ρ L 2 ,
with c a defined as before.

E.2.4 Conclusion

Based on the last subsections, we have proved the following proposition.

Proposition 15 (Risk decomposition). Under the Assumptions 1 and 2, when

(C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1/2 ≤ 1/2, ĝp -g ρ 2 L 2 ≤ 4λ 2 Lg ρ 2 L 2 + 4λµc 2 a g ρ 2 L 2 + 8 (C + λµ) -1/2 ( θρ -θ ρ ) 2 H • • • + 12c 2 a C 1/2 (I -P ) 2 op g ρ 2 L 2 + 8 (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1 θ ρ 2 H . (27)
We are left with the quantification of the different convergences when the number of labelled and unlabelled data grows towards infinity. We will quantify those convergences based on concentration inequalities.

E.3 Probabilistic inequalities

In this subsection, we bound each term appearing in Eq. ( 27) based on concentration inequalities.

E.3.1 Vector concentration

The concentration of θρ = S g ρ towards θ ρ = S g ρ is controlled through Bernstein inequality.

Theorem 2 (Concentration in Hilbert space [START_REF] Vladimirovich | On an infinite-dimensional version of S. N. Bernstein's inequalities[END_REF]). Let denote by A a Hilbert space and by (ξ i ) a sequence of independent random vectors in A such that E[ξ i ] = 0, that are bounded by a constant M , with finite variance

σ 2 = E[ n i=1 ξ i 2 ]. For any t > 0, P( n i=1 ξ i ≥ t) ≤ 2 exp - t 2 2σ 2 + 2tM/3 .
Proposition 16 (Vector concentration). When |Y | is bounded by a constant c Y , and x → k(x, x) by a constant κ 2 , we have, with D n ∼ ρ ⊗n the dataset generating the labelled data

P Dn (C + λµ) -1/2 ( θρ -θ ρ ) H ≥ t ≤ 2 exp - n t 2 2σ 2 (µλ) -1 + 2tc Y (µλ) -1/2 κ/3 , (28)
where σ 2 ≤ c 2 Y Tr(Σ) is a variance parameter to relate with the variance of Y (I + λL) -1 δ X (where the randomness is inherited from (X, Y ) ∼ ρ).

Proof. Recall that

(C + λµ) -1/2 ( θρ -θ ρ ) = (Σ + λL + λµ) -1/2 (n -1 n i=1 Y i k Xi -E ρ [Y k X ])
We want to apply Bernstein inequality to the vector ξ i = (Σ + λL + µλ) -1/2 Y i k Xi , after centering it. Let us denote by c Y a bound on |Y |, c Y ∈ R exists since we have supposed ρ of compact support. We have

σ 2 = E[ n i=1 ξ i -E[ξ i ] 2 ] = n E[ ξ -E[ξ i ] 2 ] ≤ n E[ ξ 2 ] = n E (X,Y )∼ρ Y 2 k X , (Σ + λL + µλ) -1 k X ≤ n c 2 Y E X∼ρ X k X , (Σ + λL + µλ) -1 k X = n c 2 Y Tr (Σ + λL + µλ) -1 Σ ≤ n c 2 Y Tr(Σ) (Σ + λL + µλ) -1 op ≤ n c 2 Y Tr(Σ)(µλ) -1
. Note that we have proceed with a generic upper bound, but we expect this variance, which is related to the variance of Y (I + λL + λµK -1 ) -1 δ X to be potentially much smaller -if we remove the term in P the vector concentration is the concentration of the vector S(S * S + λS LS

+ λµ) -1 Y k X K 1/2 (K + λK 1/2 LK 1/2 + λµ) -1 K 1/2 S -Y k X = (I + λL + λµK -1 ) -1 Y S -k X (I + λL + λµK -1 ) -1 Y δ X .
To capture this fact, we will write σ 2 ≤ n σ 2 (µλ) -1 , with σ = c Y Tr(Σ) 1/2 in our analysis, but potentially much smaller under refined hypothesis and in practice. Similarly to the bound on the variance, we have

ξ -E[ξ] ≤ ξ = (Σ + λL + µλ) -1/2 Y i k Xi ≤ (µλ) -1/2 c Y κ,
with κ an upper bound on k(x, x) 1/2 for x ∈ supp ρ X . As a consequence, applying Bernstein concentration inequality, we get, for any t > 0,

P Dn ( n -1 n i=1 ξ i -E[ξ i ] ≥ t) ≤ 2 exp - n t 2 2σ 2 (µλ) -1 + 2tc Y (µλ) -1/2 κ/3 .
This ends the proof.

E.3.2 Operator concentration

The convergence of Ĉ towards C is controlled with Bernstein inequality for self-adjoint operators.

Theorem 3 (Bernstein inequality for self-adjoint [START_REF] Minsker | On some extensions of Bernstein's inequality for self-adjoint operators[END_REF]). Let A be a separable Hilbert space, and (ξ i ) a sequence of independent random self-adjoint operators operators on A Assume that (ξ i ) are bounded by M ∈ R, in the sense that, almost everywhere, ξ op < M , and have a finite variance

σ 2 = n i=1 E[ξ 2 i ] op . For any t > 0, P   n i=1 (ξ i -E[ξ i ]) op > t   ≤ 2 1 + 6 σ 2 + M t/3 t 2 Tr n i=1 E[ξ 2 i ] n i=1 E[ξ 2 i ] op exp - t 2 2σ 2 + 2tM/3 . Proposition 17 (Operator concentration). When x → k(x, x) is bounded by κ 2 a nd x → ∂ 1,j ∂ 2,j k(x, x
) is bounded by κ 2 j , we have

P Dn (C + µλ) -1 2 (C -Ĉ)(C + µλ) -1 2 op > 1/2 ≤ 2 + 56 κ 2 + λ d i=1 κ 2 j λµn • • • × (1 + λµ C -1 op ) κ 2 + λ d j=1 κ 2 j λµ exp   - λµn 10 κ 2 + λ d j=1 κ 2 j   . (29)
Proof. We want to apply the precedent concentration inequality to

ξ i = (Σ + λL + λµ) -1/2 (k Xi ⊗ k Xi + λ d j=1 ∂ j k Xi ⊗ ∂ j k Xi )(Σ + λL + λµ) -1/2 ,
since we have, based on the fact that C = Σ+λL and that Σ

= E[k X ⊗k X ] and L = E[ n j=1 ∂ j k X ⊗ ∂ j k X ], (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1/2 op = n -1 n i=1 ξ i -E[ξ i ] op .
We bound ξ with

ξ op = (C + µλ) -1 2   k X ⊗ k X + λ d j=1 ∂ j k X ⊗ ∂ j k X   (C + µλ) -1 2 op ≤ Tr   (C + µλ) -1 2   k X ⊗ k X + λ d j=1 ∂ j k X ⊗ ∂ j k X   (C + µλ) -1 2   = Tr (C + µλ) -1 2 k X ⊗ k X (C + µλ) -1 2 • • • + λ d j=1 Tr (C + µλ) -1 2 ∂ j k X ⊗ ∂ j k X (C + µλ) -1 2 = (C + µλ) -1 2 k X 2 H + λ d j=1 (C + µλ) -1 2 ∂ j k X 2 H ≤ (λµ) -1   κ 2 + λ d j=1 κ 2 j   .
With κ 2 an upper bound on the kernel k and κ 2 j an upper bound on ∂ 1,j ∂ 2,j k. For the variance we have, using Löwner order,

E[ξ 2 ] sup X∈X ξ(X) op E[ξ] (λµ) -1   κ 2 + λ d j=1 κ 2 j   E[ξ] = (λµ) -1   κ 2 + λ d j=1 κ 2 j   (C + λ) -1 C (λµ) -1   κ 2 + λ d j=1 κ 2 j   .
Therefore, we get for any t > 0,

P Dn (C + µλ) -1 2 (C -Ĉ)(C + µλ) -1 2 op > t ≤ 2 1 + 6 (κ 2 + λ d i=1 κ 2 j )(1 + t/3) λµnt 2 C op + λµ C op Tr (C + λ) -1 C • • • × exp   - nt 2 2(λµ) -1 κ 2 + λ d j=1 κ 2 j (1 + t/3)   .
Remark that

Tr (C + µλ) -1 C ≤ (C + µλ) -1 op Tr(C) ≤ (λµ) -1 (κ 2 + λ d j=1 κ 2 j ).
Taking t = 1/2 ends the lemma.

E.3.3 Basis concentration

Similarly we could control (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1 θ ρ H by using concentration of self-adjoint, yet this will lead to laxer bounds, than using concentration on vectors.

Proposition 18 (Basis concentration). When x → k(x, x) is bounded by κ 2 , x → ∂ 1,j ∂ 2,j k(x, x) is bounded by κ 2 j , with Assumptions 1 and 2, we have

P Dn (C + µλ) -1 2 (C -Ĉ)(C + µλ) -1 θ ρ H > t ≤ 2 exp - µλnt 2 2c 1 (c 1 + λ 1/2 µ 1/2 t/3) . ( 30
)
with c 1 = (κ 2 + λ d i=1 κ 2 j )c a g ρ L 2 .
Proof. We want to apply Bernstein concentration inequality to the vectors

ξ i = (C + µλ) -1/2   k Xi ⊗ k Xi + λ d j=1 ∂ j k Xi ⊗ ∂ j k Xi   (C + λµ) -1 θ ρ , since (C + µλ) -1 2 (C -Ĉ)(C + µλ) -1 θ ρ H = n -1 n i=1 ξ i -E[ξ i ] H .
We bound ξ, reusing prior derivations, with

ξ i H = (C + µλ) -1/2   k Xi ⊗ k Xi + λ d j=1 ∂ j k Xi ⊗ ∂ j k Xi   (C + λµ) -1 θ ρ H ≤ (C + µλ) -1/2 op   k Xi ⊗ k Xi + λ d j=1 ∂ j k Xi ⊗ ∂ j k Xi   op (C + µλ) -1 θ ρ H . ≤ (µλ) -1/2 (κ 2 + λ d i=1 κ 2 j )c a g ρ L 2 .
For the variance, we have, similarly to prior derivations,

E[ ξ 2 ] ≤ sup X∈X k X ⊗ k X + λ d j=1 ∂ j k X ⊗ ∂ j k X op (C + λµ) -1 θ ρ 2 • • • × E   (C + µλ) -1   k X ⊗ k X + λ d j=1 ∂ j k Xi ⊗ ∂ j k X   op   ≤ κ 2 + λ d i=1 κ 2 i c 2 a g ρ 2 L 2 • • • E   (C + µλ) -1 k X ⊗ k X op + λ d j=1 (C + µλ) -1 ∂ j k Xi ⊗ ∂ j k X op   = κ 2 + λ d i=1 κ 2 i c 2 a g ρ 2 L 2 Tr (C + µλ) -1 C ≤ (λµ) -1 κ 2 + λ d i=1 κ 2 i 2 c 2 a g ρ 2 L 2 .
As a consequence, using Bernstein inequality,

P n -1 n i=1 ξ i -E[ξ i ] > t ≤ 2 exp - µλnt 2 2c 1 (c 1 + λ 1/2 µ 1/2 t/3) ,
with c 1 = (κ 2 + λ d i=1 κ 2 j )c a g ρ L 2 . Note that we have bound naively the variable ξ and its variance, but we have shown how appears sup X∈X (C + λµ) -1 k X +λ d i=1 (C + λµ) -1 ∂ i k X and Tr((C + λµ) -1 C), which under interpolation and capacity assumptions but be controlled in a better fashion.

E.3.4 Low-rank approximation

We now switch to Nyström approximation.

Proposition 19 (Low-rank approximation). When x → k(x, x) is bounded by κ2 , for any p ∈ N and t > 0, we have P Dp (I -P )Σ where we used the fact that (I -P ) a = (I -P ) and that A s B s ≤ AB s for s ∈ [0, 1] and A, B positive self-adjoint.

E.4 Averaged excess of risk -Ending the proof

Based on the precedent excess of risk decomposition, and precedent concentration inequalities, we have all the elements to derive convergence rates of our algorithm. We will enunciate this convergence in term of the averaged excess of risk of E Dn ĝp -g ρ To control the conditional expectation, we use that, when X is positive

E [X | A] P (A) = E[X] -E [X | c A] P( c A) ≤ E[X]
. This ends the proof.

Based on deviation inequalities, we can control expectations based on the equality, for X positive, E[X] = +∞ 0 P(X > t) dt. Lemma 22. In the setting of the paper,

E Dn (C + λµ) -1/2 ( θρ -θ ρ ) 2 H ≤ 8σ 2 (n µλ) -1 + 8c 2 Y κ 2 (n 2 µλ) -1 . ( 33
)
Proof. First, recall that P (C + λµ) This is the result stated in the lemma.

Lemma 23. In the setting of the paper, We conclude the proof similarly to the precedent lemma.

E Dn (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1 θ ρ 2 H ≤ 8(κ 2 + λ∂κ 2 ) 2 c 2 a g ρ 2 L 2 • • • × (µλn) -1 + (µλn 2 ) -1 , (34)
Proof. Once again, this result comes from integration of the tail bound obtained on C 1/2 (I -P )

2 op

through the one we have on Σ 1/2 (I -P ) Once again this is optimized for cb = log(p).

Remark 25 (Leverage scores). Out of simplicity, we only present low rank approximation with random subsampling. Yet, we can improve the result by considering subsampling based on leverage scores. If we consider the Gaussian kernel, Sk x ∈ L 2 can be thought as a function that is a little bump around x ∈ X . In essence, subsampling based on leverage scores, consists in representing the solution on a subsampled sequence (k Xi ) i∈I where the X i are far from one another so that the bump functions (Sk Xi ) can approximate a maximum of functions. [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF] shows that with leverage scores, we can take p = (µλ) γ log(n), with γ linked with the capacity of the RKHS linked with the kernel k.

If we add all derivations, we have derived the following theorem.

Figure 3 :

 3 Figure 3: (Left) Comparison between our kernelized Laplacian method (Tikhonov regularization version with λ = 1, µn = n -1 , p = 50) and graph-based Laplacian based on the same Gaussian kernel with bandwidth σn = n -1d+4 log(n) as suggested by graph-based theoretical results[START_REF] Hein | Graph Laplacians and their convergence on random neighborhood graphs[END_REF]. We report classification error as a function of the number of samples n. The error is averaged over 50 trials, with errorbars representing standard deviations. We fixed the ratio n /n to one tenth, and generated the data according to two Gaussians in dimension d = 10 with unit variance and whose centers are at distance δ = 3 of each other (similar to the setting of[START_REF] Castelli | On the exponential value of labeled samples[END_REF][START_REF] Lelarge | Asymptotic Bayes risk for Gaussian mixture in a semisupervised setting[END_REF]). Our method discovers the structure of the data much faster than graph-based Laplacian (to get a 20% error we need 40 points, while graph-based need 700). (Right) Time to perform training with graph-based Laplacian in orange, with Algorithm 1 in blue (with the specification of the left figure), and with the naive representation in R n(d+1) of the empirical minimizer Eq. (5) in green. When dealing with 1000 points, our algorithm, as well as graph-based Laplacian, can be computed in about one tenth of a second on a 2 GHz processor, while the naive kernel implementation requires 10 seconds. We show in Appendix B that this cut in costs is not associated with a loss in performance.

Figure 4 :

 4 Figure4: Cut in computation cost are not associated with a loss in performance. The estimate ĝp Algorithm 1 (in blue), based on low-rank approximation that cut computation cuts performs as well as the exact computation of ĝ Eq. (5). (Left) Classification error in the setting of Figure3. (Right) Regression error in the same setting. The fact that the error of the graph-based method stalls around one, is due to the amplitude of the estimate being very small, which is coherent with behaviors described in[START_REF] Nadler | Statistical analysis of semi-supervised learning: The limit of infinite unlabelled data[END_REF].

Figure 5 :

 5 Figure 5: Setting of Figure 3 with n = 1000. (Left) Training set. We represent a cut of X ⊂ R d according to the two first coordinates {(x1, x2) | (x1, x2, • • • , x d ) ∈ X }. We have two Gaussians distribution with unit variance, one centered at x = (0, 0, • • • 0) and the other one centered at x = (3, 0, • • • , 0). One of the Gaussian distribution is associated to the blue class, the other one with the orange class. We consider n = 1000 unlabelled points, represented by small points, colored according to their classes, and n = 100 labelled points, represented in colour with black edges. (Middle) Reconstruction with our kernelized Laplacian methods. Our method uncovers correctly the structure of the problem, allowing to make a quite optimal reconstruction. The optimal decision frontier being illustrated by the grey line ∂X. (Right) Reconstruction with graph-Laplacian. The graph-Laplacian diffuses information too far away from what it should, leading to many incorrect guesses.

Figure 6 :

 6 Figure6: Usefulness of Laplacian regularization. We illustrate the reconstruction based on our spectral filtering techniques based on the sole use of the covariance matrix Σ on the left, and on the sole use of the Laplacian matrix L on the right. We see that the covariance matrix does not capture the geometry of the problem, which contrasts with the use of Laplacian regularization.

Lemma 20 .. ( 31 ) 2 op+ 2 op+ 2 op+

 2031222 When L ≤ c d Σ a , we have (I -P )C 1/2 2 op ≤ (I -P )Σ 1/2 2 op + c d λ (I -P )Σ 1/2 2a op Proof. This follows from the fact that C 1/2 (I -P ) 2 op = (I -P )C(I -P ) op = (I -P )(Σ + λL)(I -P ) op ≤ (I -P )Σ(I -P ) op + λ (I -P )L(I -P ) op ≤ (I -P )Σ(I -P ) op + λc d (I -P )Σ a (I -P ) op = (I -P )Σ 1/2 λc d (I -P )Σ a/2 2 op = (I -P )Σ 1/2 λc d (I -P ) a Σ a/2 2 op ≤ (I -P )Σ 1/2 λc d (I -P )Σ 1/2 2a op ,

Lemma 21 . 2 L 2 2 • • • + 4λ 2 ρ 2 L 2 + 4λµc 2 a g ρ 2 L 2 • • • + 8 E 2 H + 12c 2 a g ρ 2 L 2 E 2 op• • • + 8 Eθ ρ 2 H. ( 32 ) 2 ,E Dn ĝp -g ρ 2 L 2 ≤ sup Dn ĝp -g ρ 2 P 2 L 2 ≤ 4c 2 Y 2 L 2 2 L 2 P(A) + 4λµc 2 a g ρ 2 L 2 P 2 H+ 12c 2 a g ρ 2 L 2 E+ 8 Eθ ρ 2 HA

 2122222222822228232222222222222222282 Under Assumptions 1 and 2,E Dn ĝp -g ρ ≤ 4c 2 Y P (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1/2 ≤ 1/Lg Dn (C + λµ) -1/2 ( θρ -θ ρ ) Dn C 1/2 (I -P ) Dn (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1Proof. We proceed using the fact thatE[X] = E[X | c A] P( c A) + E[X | A] P(A) ≤ sup X P( c A) + E[X | A] P(A), with A = D n (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1/2 ≤ 1/( c A) + E Dn ĝp -g ρ 2 A P(A).When Y is bounded by c Y , because g ρ is a convex combination of Y , we know that g ρ L 2 ≤ c Y , as a consequence, we can clip ĝp to [-c Y , c Y ], which will only improve the estimation of g ρ , as a consequence, we can consider the clipping estimate for which we have sup Dn ĝp -g ρ . Regarding the second part, we have already decomposed the risk under the eventA = D n (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1/2 ≤ 1/2 .As a consequence, we haveE Dn ĝp -g ρ ≤ 4c 2 Y P( c A) + 4λ 2 Lg ρ (A) • • • + 8 E Dn (C + λµ) -1/2( θρ -θ ρ ) Dn C 1/2 (I -P ) Dn (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1 P(A).

with ∂κ 2

 2 Let us denote byA the quantity (C + λµ) -1/2 ( Ĉ -C)(C + λµ) -1 θ ρ H , and ∂κ 2 = d i=1 κ 2 j . Recall that P (A > t) ≤ 2 exp -µλnt 2 2c 1 (c 1 + λ 1/2 µ 1/2 t/3)

Lemma 24 .•

 24 Under Assumption 3, E Dn C 1/2 (I -P ) 2a c d λ log(p) a p a

2 opΣ 1 / 2 (I -P ) 2 op

 2122 and the fact that C 1/2 (I -P ) + c d λ Σ 1/2 (I -P ) 2a op . For any a, b > 0, we have E Dn Σ 1/2 (I -P )

Theorem 4 . 2 + 2 L 2 + 4λµc 2 a g ρ 2 L 2 +

 422222 Under Assumptions 1, 2 and 3,E Dn ĝ -g ρ λ∂κ 2 ) • • • + 4λ 2 Lg ρ 64σ 2 (n µλ) -1 + 57c

  H n + H n,∂ ⊂ H, with H n = Span {k Xi | i ≤ n} and H n,∂ = Span {∂ j k Xi | i ≤ n, j ≤ d}.This subspace of H is of dimension at most n(d + 1) and if T : R p → H n + H n,∂ (with p ≤ n(d + 1)) parameterizes H n + H n,∂ , our problem can be cast in R p by considering the p × p matrices T ΣT and T ( L + µI)T instead of the operators Σ and L + µI.

	Currently, we are dealing with operators ( Σ, L) and vectors (e.g., θ) in the Hilbert space H. It is natural to wonder on how to represent this on a computer. The answer is the object of representer theorems (see Theorem 1 of [61]), and consists in noticing that all the objects introduced are actually defined in, or operate on,

  Those considerations are related to an open problem in machine learning, which is that we do not know what is the best statistical way (and the best surrogate problem) to solve the fully supervised binary classification problem [see e.g. 59]. However, many points introduced in the work could be retaken with other surrogate, could it be SVM (which leads to g

* = ϕ(f * ), with g * minimizing the Hinge loss), softmax regression (used in deep learning) or others.

Table 1 :

 1 Notations

	Symbol	Description
	(X i ) i≤n (Y i ) i≤n l ρ g ρ λ, µ g λ , g λ,µ ĝ ĝp H k S S Σ = S S K = SS L L = S LS g θ λ i e i	n samples of input data n l labels Distribution of (X, Y ) Function to learn (1) Regularization parameters Biased estimates (3, 4) Empirical estimate (5) Empirical estimate with low-rank approximation (Algo. 1) Reproducing kernel Hilbert space Reproducing kernel Embedding of H in L 2 Adjoint of S, operating from L 2 to H Covariance operator on H Equivalent of Σ on L 2 Diffusion operator (a.k.a Laplacian) Restriction of the diffusion operator to H Generic element in L 2 Generic element in H Generic eigen value Generic eigen vector in L 2

Table 2 :

 2 Steps in the consistency analysis

	Estimate Vector	Property of convergence	Basis
	ĝp	S g ρ

(P ΣP, P LP + µP )

  Proof. Reusing Proposition 3 of Rudi et al.[START_REF] Rudi | Less is more: Nyström computational regularization[END_REF], for any γ > 0, we have, with P the projection on Span {k Xi } i≤p andΣ = p -1 p i=1 k Xi ⊗ k Xi , Dp γ Σ 1/2 ( Σ + γ) -1 Σ 1/2 -γt -1). Taking γ = t/2, this term is simplified as which is the object of this proposition.

			1/2 2	> t ≤ 2 +	116κ 2 tp	(2 + t Σ	-1 op )	κ 2 t	exp -	pt 10κ 2 ,
	(I -P )Σ 1/2 2	≤ γ ( Σ + γ) -1/2 Σ 1/2 2	≤ γ Σ 1/2 ( Σ + γ) -1 Σ 1/2	.
										op	op
	As a consequence, skipping derivations that can be retaken from our precedent proofs,
	P Dp	(I -P )Σ 1/2 2	> t ≤ inf			op	> t
		≤ inf
		≤ inf γ>0	2 + 56	κ 2 γp	(1 + γ Σ	-1 op )	κ 2 γ	exp -	pγu 2 2κ 2 (1 + u/3)	,
	with u = (1 P Dp	(I -P )Σ 1/2 2	> t ≤ 2 + 116	κ 2 tp	(2 + t Σ	-1 op )	κ 2 t	exp -	pt 10κ 2 ,

γ>0 P γ>0 P Dp (Σ + γ) -1/2 ( Σ -Σ)(Σ + γ) -1/2 op > (1 -γt -1 )

  -1/2 ( θρ -θ ρ ) (µλ) -1 + 2tc Y (λµ) -1/2 κ/3 ≤ 2 exp -n t 2 2 max 2σ 2 (µλ) -1 , 2tc Y (λµ) -1/2 κ/3As a consequenceE (C + λµ) -1/2 ( θρ -θ ρ ) P (C + λµ) -1/2 ( θρ -θ ρ )

				2	+∞	2
					=	> t dt
				H	0	H
	≤ 2 exp -	n µλt 4σ 2	dt + 2 exp -	3n µ 1/2 λ 1/2 t 1/2 4c Y κ	dt.
	= 8σ 2 (n µλ) -1 +	64c 2 Y κ 2 9	(n 2 µλ) -1 .
			H	> t ≤ 2 exp -2σ 2 ≤ 2 exp -n µλt 2 4σ 2	n t 2 + 2 exp -	3n µ 1/2 λ 1/2 t 4c Y κ	.

  2 Σ This last quantity is optimized for ab = log(p), which leads to the first part of lemma. Similarly E Dn Σ 1/2 (I -P )

							2a op	=		0	∞	P Dn	Σ 1/2 (I -P )	op 2a	> t dt
	=	0	∞	P Dn		Σ 1/2 (I -P )	op 2	> t 1/a dt
	≤	0	∞	min 1, 2κ 2 Σ		-1 op	1 +	58κ 2 t 1/a p	1 +	2κ 2 t 1/a exp -	pt 1/a 10κ 2	dt
	=	10 a κ 2a ac a p a	0	∞	min u a-1 , 2κ 2 Σ	-1 op	1 +	58 10cu	1 +	p 5cu	1 u 1-a exp (-cu) du
	≤	10 a κ 2a ac a p a	b a a	+		b	∞	2κ 2 Σ	-1 op	1 +	6 cu	1 +	p 5cu	1 u 1-a exp (-cu) du
	≤	10 a κ 2a p a	(cb) a + 2κ 2 Σ	-1 op	1 +	6 cb	1 +	p 5cb	1 (cb) 1-a exp (-cb) .
														-1 op	1 +	58κ 2 tp	1 +	2κ 2 t	exp -	pt 10κ 2	dt
			=	10κ 2 a p	0	∞	min 1, 2κ 2 Σ	-1 op	1 +	58 10au	1 +	p 5au	exp (-au) du
			≤	10κ 2 a p	b +	b	∞	2κ 2 Σ	-1 op	1 +	6 au	1 +	p 5au	exp (-au) du
			≤	10κ 2 p		ab + 2κ 2 Σ	-1 op	1 +	6 ab	1 +	p 5ab	exp (-ab) .

  2 Y κ 2 (n 2 µλ) -1 • • • + 64(κ 2 + λ∂κ 2 ) 2 c 2 a g ρ 2 L 2 (µλn) -1 + 57(κ 2 + λ∂κ 2 ) 2 c 2

a g ρ 2 L 2 (µλn 2 ) -1

  where c Y is an upper bound on Y , κ 2 is an upper bound on x → k(x, x), ∂κ 2 = → ∂ 1i ∂ 2i ∂k xi , c d and a the constants appearing in Assumption 3, c a a constant such that g H ≤ c a g L 2 and σ 2 ≤ c 2 Y κ 2 a variance parameter linked with the variance of Y (I + λL) -1 δ X . Theorem 1 is a corollary of this theorem.

									(36)
	• • • + 12c 2 a g ρ	2 L 2	10κ 2 log(p) p	+	10 a κ 2a c d λ log(p) a p a
	• • • × 1 +	2κ 2 Σ op log(p)	1 +	6 log(p)	1 p	+	1 5 log(p)	.
	bound on x								d i=1 κ 2 i with κ 2 i a

Our code is available online at https://github.com/VivienCabannes/partial_labelling.

To parameterize functions g from X to HY , we can parameterized independently each coordinates g, ei H Y , for (ei) a basis of HY , by the space G -note that it is possible to generalize real-valued kernel to parameterize coordinates in a joint fashion[START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF]. The coordinate-wise parameterization corresponds to the tensorization H = HY ⊗ H and to the parametric space G = {x → Θkx | Θ ∈ H } of functions from X to HY . G naturally inherits of the Hilbertian structure of H , itself inherited from the structure of H and HY .

Note that our consistency result Theorem 1 describes a convergence regime that applies to a vast class of problems. Such a regime usually takes place after a certain number of data (depending on the value of the constant C). Before entering this regime, describing the error of our algorithm would require more precise analysis specific to each problem instance, eventually involving tools from random matrix theory.

L 2 .

Acknowledgments and Disclosure of Funding

This work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). We also acknowledge support of the European Research Council (grants SEQUOIA 724063 and REAL 94790).