Marianne Clausel

Stephan Clémençon

Johan Suykens

Mastane, Eugène Alex Garcia

Pierre Ablin

Anas, Sholom Hamid Guillaume

Dimitri

Luc

À la mémoire de mon grand-père Jeannot. Résistant, Pompier, Pêcheur

Résumé

Les méthodes à noyaux sont au coeur de l'apprentissage statistique. Elles permettent de modéliser des fonctions à valeurs réelles dans des espaces de fonctions à fort potentiel représentatif, sur lesquels la minimisation de risques empiriques régularisés est possible et produit des estimateurs dont le comportement statistique est largement étudié. Lorsque les sorties ne sont plus réelles mais de plus grande dimension, les Espaces de Hilbert à Noyaux Reproduisants à valeurs vectorielles (vv-RKHSs) basés sur des Noyaux à Valeurs Opérateurs (OVKs) fournissent des espaces de fonctions similaires et permettent de traiter des problèmes tels que l'apprentissage multi-tâche, la prédiction structurée ou la regression à valeurs vectorielles.

Dans cette thèse, nous étudions différents problèmes liés à l'apprentissage de modèles à valeurs fonctionnelles, un cas que nous traitons en modélisant les sorties comme vivant dans un espace de Hilbert de dimension infinie. Ceci nous permet notamment d'exploiter la théorie associée aux vv-RKHSs précédemment évoquée, avec l'aide de noyaux à valeurs opérateurs agissant sur ces mêmes espaces. Ces problèmes requièrent des fonctions de perte dédiées, définies sur des espaces fonctionnels qui diffèrent des fonctions de perte usuelles utilisées en dimension finie. En particulier, nous nous intéressons aux fonctions de perte convexes pouvant s'exprimer sous la forme d'une intégrale sur l'espace de définition des fonctions de sortie. Ces pertes intégrales constituent une extension naturelle des fonctions de perte utilisées en dimension finie, et permettent d'aborder la recherche de fonctions cibles comme équivalent à la minimisation d'un risque dans les vv-RKHSs.

Le cadre d'apprentissage que nous proposons repose sur la minimisation de risques empiriques régularisés en présence de telles fonctions de perte. De par la nature fonctionnelle des sorties, les problèmes d'optimisation en résultant souffrent d'écueuils supplémentaires à ceux déjà présents en dimension finie. En premier lieu, et contrairement au cadre fini-dimensionnel, il n'est pas garanti que les estimateurs bénéficient d'une représentation de taille finie sur une base naturelle de l'espace de fonction considéré. De plus, même si une telle décomposition était trouvée, les pertes intégrales sont accompagnées de problèmes de calculabilité: en effet les intégrales (et leurs gradients) ne peuvent le plus souvent pas être exprimées sous forme analytique et doivent être estimées.

Une première contribution de cette thèse est de fournir un ensemble de techniques d'optimisation à même de répondre à ces deux problématiques, tout en s'assurant que les algorithmes proposés ne soient pas trop gourmands en ressources de calcul. En reposant sur un échantillonage des pertes intégrales, nous obtenons une représentation de dimension finie des estimateurs pour différents choix de régularisation dans les vv-RKHSs, rendant possible l'utilisation d'algorithmes de descente de gradient dont la convergence est bien étudiée. L'usage de la dualité lagrangienne vient compléter cette technique, en offrant un problème nouveau basé sur la transformée de Fenchel-Legendre des pertes intégrales, qui est explicitée pour quelques fonctions de perte utiles. Les difficultés liées à la représentation des variables duales sont traitées par l'usage de splines linéaires, qui conjuguées à l'approximation de l'opérateur intégral lié au noyaux rendent les problèmes d'optimisation solubles. Les problèmes de passages à l'échelle sont aussi traités par l'utilisation de noyaux approchés, dont les vv-RKHSs associés sont de dimension finie. Ceci permet de plus l'utilisation d'algorithmes de descente de gradient stochastique, les pertes intégrales étant interprétées comme des espérances. Enfin, nous proposons une architecture de modèle composée d'un réseau de neurone et d'une dernière couche à noyaux, qui rend possible l'apprentissage de représentations appropriées aux noyaux utiles dans les applications avec des données complexes comme les images.

Une seconde contribution de cette thèse est l'introduction d'une extension fonctionnelle originale du cadre multi-tâche appelée Apprentissage d'un Continuum de Tâches (ITL), qui permet de résoudre conjointement un continuum de tâches paramétrées, parmi lesquelles la régression quantile, la classification à coût assymétrique, ou l'estimation de niveaux de densité. Le cadre ITL est traité grâce aux techniques d'optimisation précédemment développées pour les pertes intégrales, avec quelques ajustements ayant pour but de gérer des contraintes souples sur la monotonicité des estimateurs dans le cas de la régression quantile notamment. Nous nous intéressons aussi aux capacités en généralisation des estimateurs, que nous contrôlons théoriquement en les étudiant sous l'angle de la stabilité uniforme. Enfin, cette approche fonctionnelle au cadre multi-tâche nous permet aussi de revisiter les problemes de transfert de style sous l'angle ITL, avec une application au transfert d'émotion.

Le dernier problème abordé dans cette thèse concerne la regression fonction-à-fonction robuste en présence de valeurs aberrantes. Nous substituons à la norme carrée des pertes convoluées, comme la perte ϵ-insensible ou la perte Huber, plus à même d'imposer des estimateurs parcimonieux ou robustes. Les problèmes d'optimisation sont traités dans le dual, où les synergies entre convolution infinitésimale et transformée de Fenchel-Legendre rendent possible le développement d'algorithmes de descente de gradient proximal. En particulier, certaines libertés dans la définition des pertes convoluées permettent l'émergence de différents types de parcimonie, ainsi que la résistance à différents types de valeurs aberrantes. Due to the increasing complexity of collected data, which can come from many sensors, embedded devices, IOT, or any other data acquisition pipeline, real world applications are in dire need of machine learning systems able to deal with sophisticated data. Especially, when the collected data correspond to the recording of the behaviour of a phenomenon through time or space for instance, there is an interest for considering the collection of datapoints as observations of a function. One can think for example to data describing the trajectory of a plane through time, the evolution of CO 2 levels in the atmosphere, or the transmission of biochemical signals inside the brain.

These real world scenarii have motivated the research in functional data analysis (FDA, Ramsay and Silverman 1997) whose goal is to provide a framework dedicated to data modeled as functions. The potential scope of FDA is enormous and diverse [START_REF] Ullah | Applications of functional data analysis: A systematic review[END_REF], and has attracted a great deal of attention [START_REF] Wang | Functional data analysis[END_REF].

In particular, early works focused on semi-functional regression where the explanatory variable is a function and the target variable is a scalar. The relationship between these two can then be modeled as linear [START_REF] Cardot | Spline estimators for the functional linear model[END_REF], nonlinear [START_REF] Cardot | Spline estimators for the functional linear model[END_REF] or hybrid [START_REF] Aneiros-Pérez | Semi-functional partial linear regression[END_REF], and many learning setting were investigated such as quantile regression [START_REF] Cardot | Quantile regression when the covariates are functions[END_REF], robust estimation [START_REF] Crambes | Robust nonparametric estimation for functional data[END_REF][START_REF] Azzedine | On robust nonparametric regression estimation for a functional regressor[END_REF] or variable selection [START_REF] Aneiros | Variable selection in semi-functional regression models[END_REF].

More recently, the case of functional output regression (FOR) where both explanatory and target variables are functions has been under study. The FOR setting is more demanding than the semi-functional one, as it requires to model function-valued functions. The target variable can then belong to an infinite dimensional Hilbert space, a problem that needs to be addressed from a modeling point of view. This fundamental issue has lead to consider linear models [START_REF] Morris | Functional regression[END_REF], non-parametric modeling (Ferraty and Vieu, 2006), or kernel methods [START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel Hilbert spaces[END_REF]Kadri et [START_REF] Kadri | Operatorvalued kernels for learning from functional response data[END_REF][START_REF] Reimherr | Optimal prediction for additive function-on-function regression[END_REF]. The latter is of particular interest to us, and one goal of this thesis is to propose and study a general framework for function-valued regression beyond the square loss, and in the context of vector-valued RKHSs [START_REF] Pedrick | Theory of reproducing kernels for Hilbert spaces of vector-valued functions[END_REF].

Apart from functional data, we have noticed that it can be interesting to have predictive models with function-valued outputs even if the output observations do not directly correspond to a function. In particular, in a large number of ML problems, the task at hand depends on some hyperparameter varying continuously and solving the problem for a continuum of values presents many interests. A functional view on these problems can be seen as an extension of multi-task learning (MTL; [START_REF] Evgeniou | Regularized multi-task learning[END_REF], which is the second motivation for this thesis. MTL consists in learning jointly a finite number of machine learning tasks with a single model, the underlying assumption being that solving these tasks together will bring better results than solving them independently. One can think for example to tasks such as quantile regression (QR, [START_REF] Koenker | Regression quantiles[END_REF], anomaly detection using one class support vector machines (OC-SVM, Schölkopf et al. 2001b) or cost-sensitive classification (CSC, [START_REF] Elkan | The foundations of cost-sensitive learning[END_REF] where the task is characterized by some hyperparameter encoding information about the target function: the quantile level in QR, the fraction of outliers in the OCSVM, or the relative importance associated to false positive and false negatives in CSC. The links between MTL and vector-valued regression are fertile (Álvarez et al., 2012), and MTL can be tackled using kernels for vector-valued functions [START_REF] Micchelli | On learning vector-valued functions[END_REF].

In the context of learning function-valued functions, a natural problem to consider is the joint learning of infinitely many tasks. Such problem has been proposed in the seminal paper [START_REF] Takeuchi | Parametric task learning[END_REF] for tasks whose loss function depends piecewise linearly on the hyperparameter. Extending this framework to a larger class of tasks and models is one of the motivation of this thesis.

The organization of this chapter is as follows: in Section 1.1 we present the type of machine learning tasks tackled in this thesis, and raise related research questions in Section 1.2. Our contributions are summarized in Section 1.3, with associated publications and preprints in Section 1.4.

Machine Learning Tasks

We introduce below two examples of supervised learning tasks, namely the binary classification and least squares regression.

Binary classification:

The simplest example of a machine learning task is given by binary classification. In this setting, we are given two random variables (X, Y) ∈ X × {-1, 1}. X encodes the representation of the data, or features in the space X, and Y is the corresponding class. The goal of supervised classification is to build a classifier ĥ : X → {-1, 1} which attributes the right class to a realization of the random variable X, without knowing the probability distribution of (X, Y) and being only given a finite sample (x i , y i ) n i=1 independently drawn from the probability distribution of (X, Y), called the training sample [START_REF] Vapnik | An overview of statistical learning theory[END_REF].

Theoretically, the ideal classifier involves the optimization problem h † ∈ arg min h measurable E (X,Y) {-Y} (h(X)) .

(1.1)

MACHINE LEARNING TASKS

The loss is called the 0 -1 loss as it returns 0 whenever the prediction is correct and 1 when it fails. The best decision rule is given by the Bayes classifier

h † (x) :=    1 if P(Y = 1 | X = x) ≥ 1 2 0 otherwise.
However, given that the underlying distribution of (X, Y) is unknown, this rule cannot be applied in practice and one must trade Equation (1.1) against an empirical version based on the available training set (x i , y i ) n i=1 . Moreover, minimizing the 0 -1-loss is notoriously hard (can be NP-hard for many choices of H, Ben-David et al. 2003), so that a convex upper bound such as the hinge loss can be used instead. Combined with a margin based approach [START_REF] Vapnik | Statistical Learning Theory[END_REF], this results in ĥ ∈ arg min

h∈H 1 n n i=1 max (0, 1 -y i h(x i )) + Ω(h),
where H is a set of possible candidates denoted hypothesis space, and Ω(h) is a regularizer used to avoid overfitting. This leads to the celebrated support vector machines (SVM) problem formulation and its separating hyperplane interpretation when using kernel models [START_REF] Cortes | Support-vector networks[END_REF].

Least squares regression:

When the random variable Y is not binary but realvalued, inferring Y from X is called a regression task, the most well-known loss function for regression being the square loss. The goal is then to build a prediction function ĥ : X → R such that the residual Y -h † (X) is the smallest possible in expected square norm, having only access to a training dataset (x i , y i ) n i=1 . Minimizing the square loss has the benefit of estimating the conditional expectation of Y given X, as:

E[Y|X] = h † (X), h † ∈ arg min h measurable E (X,Y) (Y -h(X)) 2
The theoretical best prediction at point x is achieved by h † (x) = E Y|X=x Y . Again, this predictor is of little convenience in practice as one does not have access to the underlying distribution of (X, Y) and empirical versions based on the training dataset (x i , y i ) n i=1 must be considered, such as the regularized empirical risk minimization problem ĥ ∈ arg min

h∈H 1 n n i=1 (y i -h(x i )) 2 + Ω(h),
where H is hypothesis space and Ω(h) is a regularizer. These two simple examples fall under the regularized empirical risk minimization umbrella [START_REF] Tikhonov | Solution of Ill-posed Problems[END_REF][START_REF] Girosi | Regularization theory and neural networks architectures[END_REF], which consists in solving ĥ ∈ arg min

h∈H 1 n n i=1 ℓ(h(x i ), y i ) + Ω(h),
where ℓ is a loss function used to measure the discrepancy between a prediction h(x i ) and the ground truth y i , and Ω(h) is a regularizer. Ideally, an estimator ĥ should perform similarly on the training dataset (x i , y i ) n i=1 and on new data sampled from CHAPTER 1. MOTIVATION AND CONTRIBUTIONS P (X,Y) . The generalization capabilities of the estimator can be controlled by bounding E (X,Y) ℓ( ĥ(X), Y) -1 n n i=1 ℓ( ĥ(x i ), y i ) (1.2) in high probability [START_REF] Vapnik | Statistical Learning Theory[END_REF].

Parameterized Tasks

Parameterized tasks arise when the target function to estimate depends on a parameter that impacts the task. We give below three examples of parameterized tasks.

Cost-sensitive classification:

In the binary classification scenario, it may happen that the cost associated to making a mistake on one or the other class differs in real applications. For example in epidemiology, when designing a test to assess whether a person is contaminated by a virus or not, a false positive is of less impact than a false negative as the latter has dramatic consequence on the spread of the virus. It can then be of advantage to introduce an imbalanced coefficient that penalizes differently the two kinds of mistakes [START_REF] Bach | Considering cost asymmetry in learning classifiers[END_REF]. A parameter θ ∈ [-1, 1] is capable of encoding such asymmetry, and gives rise to the task ĥ ∈ arg min

h∈H 1 n n i=1 θ + 1 2 -{-1} (y i ) max (0, 1 -y i h(x i )) + Ω(h).
This way when θ is close to 1, mistakes made on ground truth samples with class -1 are almost not penalized, so that the classifier focuses on being right on the class 1, and vice versa.

Quantile regression:

In the regression setting, the conditional expectation may not suffice to give sufficiently good prediction of the output variable. One could think for example to data distributions where for all x ∈ X, P Y|X=x is bimodal and symmetric around 0, so that the conditional expectation E Y|X=x [Y] = 0 would be a poor choice of prediction as it would not be meaningful with respect to the underlying data distribution. In such cases, one can resort to estimating the conditional quantiles of Y given X. Given a quantile level θ ∈ (0, 1), it is defined as

q(x) := inf{u ∈ R : P(Y ≤ u | X = x) ≥ θ}.
Estimating the conditional quantile can be tackled through its variational formulation

q(x) = arg min u∈R E Y [ℓ(θ, u, Y) | X = x],
where ℓ is the so-called pinball loss [START_REF] Koenker | Regression quantiles[END_REF] ℓ(θ, u, y) = max (θ(y -u), (θ -1)(y -u)).

ERM based algorithms then allow to estimate q in various function spaces [START_REF] Steinwart | Estimating conditional quantiles with the help of the pinball loss[END_REF].

MACHINE LEARNING TASKS

Density level set estimation: In anomaly detection, the goal is to separate inliers, which are coherent observations, from outliers which are seen as anomalies in the data. Given an X-valued random variable X, one aims at building a prediction function ĥ : X → {-1, 1} used for deciding whether a new point is an inlier or an outlier. Given a desired fraction of outliers θ ∈ (0, 1), a natural way to achieve this goal is to estimate the minimum volume sets (MV-sets) associated to X [START_REF] Polonik | Minimum volume sets and generalized quantile processes[END_REF]:

MV(θ) = arg min G⊆X, G Borel Λ(G) | P X (G) ≥ θ ,
where Λ(G) is the Lebesgue measure of a Borel set G. MV-sets describe regions where X is most concentrated, and relate closely to the level sets of the density, in the sense that they describe the same sets but with a different parameterization (the level of the density versus the mass of the set with respect to P X ). Outlierness of a point x can then be characterized by its belonging to MV(θ). In particular, MV-sets can be estimated using one-class SVMs (Schölkopf et al., 2001b).

Learning Vector-Valued Functions

While aforementioned tasks focus on predicting a single value, complex machine learning systems often require to predict jointly several values. This scenario typically occurs in the regression setting when the output variable Y is R p -valued, or when when one's goal is to solve jointly multiple tasks [START_REF] Micchelli | On learning vector-valued functions[END_REF]. The latter is referred to as multi-task learning (MTL) and has attracted a great deal of attention in machine learning [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF] with far-reaching applications such as autonomous driving, climate science, or functional brain imaging.

In its most generic form, MTL proposes to solve jointly p tasks, where each task is encoded by a particular objective function, and observed dataset. Having access to p datasets (X j ) p j=1 where each X j := (x i,j , y i,j ) n j i=1 ∈ (X j × R) n j , the objective is then to build p estimators ( ĥj ) p j=1 such that each ĥj : X j → R performs well on task j. When the dataset is shared across the tasks, this simplifies to finding a vector-valued model ĥ : X → R p such that each coordinate of ĥ performs well on the corresponding task [START_REF] Micchelli | On learning vector-valued functions[END_REF]. This simpler setting is also referred to as multi-output learning (Álvarez et al., 2012) and is the setting considered in this thesis.

Compared to independent learning of each task, the advantage of MTL is to take advantage of the similarity between the tasks to ensure consistency between them. For example in quantile regression, the quantiles are by definition nondecreasing with respect to θ, and that should be reflected in the estimator that jointly estimates the conditional quantiles associated to some (θ j ) p j=1 . A simple way to jointly solve p tasks parameterized by (θ j ) p j=1 and loss functions (ℓ(θ j , •, •)) p j=1 is then to sum the corresponding loss functions:

ĥ ∈ arg min h∈H 1 np n i=1 p j=1 ℓ(θ j , h(x i ) j , y i ) + Ω(h),
where H ⊂ F(X, R p ) is a hypothesis space and Ω is some regularization term encoding similarities between tasks (Argyriou et al., 2008b,a;[START_REF] Baldassarre | Multi-output learning via spectral filtering[END_REF].

CHAPTER 1. MOTIVATION AND CONTRIBUTIONS

Relevant application-specific references include [START_REF] Sangnier | Data sparse nonparametric regression with ϵ-insensitive losses[END_REF] for quantile regression and [START_REF] Glazer | q-OCSVM: A q-quantile estimator for high-dimensional distributions[END_REF] for density level set estimation. We finally want to mention alternatives to summing the loss functions such as Pareto multi-task learning which make use of multiobjective optimization [START_REF] Lin | Pareto multi-task learning[END_REF].

Learning Function-Valued Functions

Learning function-valued functions is the goal of this thesis. Such scenario can be considered as an extension of the prior vector-valued case beyond Euclidean space R p to generalized Hilbert spaces.

Functional output regression:

The problem of learning function-valued function can appear in supervised learning when the output variable Y is a function itself. This happens for example in biomedical signal processing, epidemiology monitoring or climate science where the phenomena under study exhibit a functional nature, and the understanding of these phenomena depends on machine learning algorithms being able to reliably predict functional data [START_REF] Ramsay | Applied functional data analysis: methods and case studies[END_REF].

In this problem family the task is to regress to a functional output from a vectorial or functional input, a setting referred to as functional output regression (FOR). Assuming that Y takes values in a functional space L 2 [Θ, µ], where Θ is the domain of the realizations of Y and µ a probability measure on it, a natural way to tackle this problem is to estimate E[Y|X] = h † (X) where it is well-known that

h † = arg min h∈H E (X,Y) Y -h(X) 2 L 2 [Θ,µ]
,

with H being the set of all measurable functions from X to L 2 [Θ, µ]. As the underlying distribution is unknown, one can use regularized empirical risk minimization in suitable functional spaces to solve this problem, as proposed in the seminal works of [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF][START_REF] Kadri | Operatorvalued kernels for learning from functional response data[END_REF].

Jointly learning infinitely many tasks:

The problem of learning function-valued functions also emerges when considering the joint learning of infinitely many parameterized tasks. The goal is then to learn a mapping ĥ :

x → (θ → ĥ(x)(θ) ∈ R)
such that for all θ ∈ Θ, ĥ(•)(θ) performs well on the task θ. Here Θ is a space encoding the tasks, such as the quantile level in quantile regression (Θ = (0, 1)), the mass of the MV-sets in anomaly detection (Θ = (0, 1)) or the asymmetry coefficient in cost-sensitive classification (Θ = [-1, 1]).

The resulting optimization problem is

ĥ ∈ arg min h∈H 1 n n i=1 Θ ℓ(θ, h(x i )(θ), y i )dµ(θ) + Ω(h),
where µ is a probability measure describing the relative importance of each task, H ⊂ F(X, F(Θ, R)) is a suitable hypothesis space and Ω(h) is a regularizer. In their seminal work, [START_REF] Takeuchi | Parametric task learning[END_REF] introduce the parametric task learning framework able to handle a continuum of tasks when the loss function has a piecewise linear dependency to 1.2. RESEARCH QUESTIONS the θ parameter. By leveraging solution path techniques, they show that the associated optimization problem is tractable and yields a model whose dependency to θ is also piecewise linear.

We can notice that both problem families involve a class of loss functions that we denote integral losses, as the discrepancy between the function-valued prediction h(x) and the observed output y is computed by means of an integration of the local loss functions ℓ(θ, •, •) over the Θ space.

Research Questions

The goal of this thesis is to propose a general framework to learn function-valued functions, both to offer novel functional output regression tools and to extend multitask learning to a continuum of tasks. The outputs are supposed to be functions over a domain Θ, and we consider a family of losses that write as an integral over this domain. This leads to the formulation of the corresponding regularized empirical risk minimization problems, with the additional requirement of picking a suitable hypothesis space. This hypothesis space is chosen to be a vector-valued reproducible kernel Hilbert space (vv-RKHS), a generalization of RKHSs allowing function-valued models whose full introduction is postponed to Chapter 2. Solving these optimization problems is not straightforward, and raises various challenging questions:

• How can we solve problems involving integral losses while they are not computable analytically?

• How do we represent the models on a computer despite the vv-RKHS being an infinite-dimensional functional space?

• What compromises can be accepted on the hypothesis space to scale to larger datasets?

For the specific functional output regression problem, existing techniques in the literature mostly focus on the square loss as a measure of the discrepancy between predictions and ground truth, which is known to be sensitive to outliers in the data. Especially, in a learning scenario where the data is contaminated with erroneous values, square loss based estimators tend to perform poorly. We suggest the following lines of research:

• How can we go beyond the square loss in functional output regression? Can we take advantage of dual approaches to allow sparse or robust estimators?

In Chapter 5 we develop a learning framework able to jointly learn infinitely many parameterized tasks, referred to as infinite task learning (ITL). Interesting questions then include:

• Which type of regularization can be used? How can we induce consistency between the tasks?

• What kind of generalization guarantees can we get for the ITL estimator?

• Can we design kernels and appropriate optimization algorithms suited to handle complex structured data such as images?

• Can we extend the ITL framework to more involved learning tasks with higher dimensional Θ?

Contributions

We now list the contributions of this thesis to the aforementioned problems, followed by the description of the organization of the manuscript. In particular, Chapter 2 is not included in these contributions as it is a background chapter on convex optimization and (operator-valued) kernel methods used throughout the manuscript.

▶ Chapter 3 develops optimization algorithms to solve regularized empirical risk minimization problems in vector-valued RKHSs, in the presence of integral losses.

Closed-form solutions are obtained for the square loss based estimators in the fully and partially observed regime, extending known results to vectorial outputs. Primal optimization algorithms are proposed for general loss functions, relying either on some sampling of the integral loss that pertains the solution to a finitedimensional subspace of the vv-RKHS, or on stochastic optimization algorithms made possible by the use of random Fourier features. Finally, dual methods are explored, with compatibility conditions between the loss and the vv-RKHS allowing for different representation of the solution and dedicated (proximal) gradient descent algorithms.

▶ Chapter 4 investigates the functional output regression problem with a focus on robustness and sparsity. To ensure these properties, new loss functions are designed with building blocks being integral losses and infimal convolutions. In particular, these losses extend the classical Huber and ϵ-insensitive losses to the functional output case. The resulting optimization problems are solved using dual methods, shown to be well-suited to the nature of the new losses and enabling proximal gradient descent algorithms.

▶ Chapter 5 introduces the infinite task learning framework, able to jointly learn a continuum of parameterized tasks and exemplified in quantile regression, costsensitive classification and density level set estimation. Different choices of regularizers are proposed, corresponding to regularization in vv-RKHS norm or in mixed L 2 -RKHS norm, for which dedicated double representer theorems are stated. We analyse the generalization capabilities of the ITL estimator for the supervised setting using uniform stability. Finally, we show that hybrid models involving kernels and dedicated neural networks allow one to obtain more accurate predictions in the context of image processing.

▶ Chapter 6 casts a style transfer problem as an ITL problem with vectorial outputs. We consider the problem of emotion transfer for facial landmarks, for which a resolution in vv-RKHS is derived. We show that encoding emotions continuously in some embedding space makes sense for this problem and achieves good performance on two facial datasets, with the additional benefit over existing method of enjoying predictions over unseen emotions.

The algorithms developed in this thesis are gathered in the open source Python library torch_itl. It is designed as a high level python library exploiting pytorch's autodifferentiation capabilities and providing a scikit-like api. For now it supports quantile regres-1.4. PUBLICATIONS sion and vectorial infinite task learning (applied in the context of emotion transfer), but future release will include the rest of the proposed methods to enhance reproducibility.

Publications

These contributions have resulted in the following peer-reviewed publications and preprints (⋆ indicates equal contribution)

• R. Brault In this chapter, we present several mathematical tools that will be used throughout the manuscript. In Section 2.1 we start with a brief summary on optimization of convex functionals in general Hilbert spaces, including Fenchel-Legendre conjugates, parametric duality and proximal operators. In Section 2.2 we focus on operator-valued kernels (OVKs) for machine learning. These classes give rise to expressive spaces of functions called vector-valued reproducible kernel Hilbert spaces (vv-RKHSs), allowing to model function-valued functions and chosen as hypothesis spaces in all further applications.

Convex Optimization Tools

Ultimately, machine learning consists in the minimization of a function, in the sense that the output of a learning algorithm is the solution of a minimization problem. Among these problems, convex ones are notoriously easier to deal with than nonconvex ones. While nonconvex optimization has to deal with the existence of potentially different local minima, convex optimization problems enjoy the property that any local minimum is also a global minimum. In this thesis, we are concerned with convex optimization problems, over a potentially infinite-dimensional Hilbert space denoted by H. We refer to the monographs by e.g. [START_REF] Rockafellar | Convex analysis[END_REF]; [START_REF] Boyd | Convex optimization[END_REF]; [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF], and adopt the notations and terminology of the latter as the optimization framework they design is suited to minimization problems over general Hilbert spaces. We begin by defining a class of functions ubiquitous in convex optimization. Definition 2.1 (Proper, convex, lower semi-continuous functions). We denote by Γ 0 (H) the set of functions f : H → -∞, +∞ that are CHAPTER 2. BACKGROUND ON LEARNING WITH OPERATOR-VALUED KERNELS

1 2 • 2 H 1 2 • 2 H • p χ B p ⋆ 1 where 1 p + 1 p ⋆ = 1 ϵf ϵf ⋆ ( • ϵ ) for all ϵ > 0 f (• -y) f ⋆ + •, y H for all y ∈ H f □ g f ⋆ + g ⋆ • proper: dom(f ) := x ∈ H : f (x) < +∞ = ∅, • convex: ∀x, y ∈ H, ∀t ∈ [0, 1], f (tx + (1 -t)y) ≤ tf (x) + (1 -t)f (y), and
• lower semicontinuous: ∀x ∈ H, lim y→x f (y) ≥ f (x)
, where lim denotes limit inferior.

The Fenchel-Legendre conjugate of a function is an important notion for exploiting duality principles. In particular it plays an important role in convex machine learning as it is a key tool to solve variational problems.

Definition 2.2. The Fenchel-Legendre conjugate of a function f : H → -∞, +∞ is defined by ∀x ∈ H, f ⋆ (x) := sup y∈H x, y H -f (y).
(2.1)

The Fenchel-Legendre conjugate of a function f is always convex. It is also involutive on Γ 0 , meaning that (f ⋆ ) ⋆ = f for any f ∈ Γ 0 (H). We gather in Table 2.1 examples and properties of Fenchel-Legendre conjugates.

We now introduce the infimal convolution operator following [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF].

Definition 2.3 (Infimal Convolution). The infimal convolution of two functions f, g : H → -∞, +∞ is f □ g :   H → [-∞, +∞] x → inf x ′ ∈H f (x -x ′ ) + g(x ′ )   . (2.2)
The infimal convolution operator is commutative, which means that

f □ g = g □ f for all f, g : H →] -∞, +∞].
It allows to define a smooth approximation of a potentially non-smooth function f through the Moreau enveloppe defined as f □ 1 2γ • 2 H where γ > 0. One key property of the infimal convolution operator is that it behaves nicely under Fenchel-Legendre conjugation, as detailed in the following proposition.

Proposition 2.4 [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF], proposition 13.24). Let f, g : H → -∞, +∞ . Then

(f □ g) ⋆ = f ⋆ + g ⋆ .
We now define the proximal operator, used as a replacement for the classical gradient step in the presence of non-differentiable objective functions. Definition 2.5 (Proximal Operator, [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]). The proximal operator (or proximal map) is defined as

∀(f, x) ∈ Γ 0 × H, prox f (x) := arg min y∈H f (y) + 1 2 x -y 2 H . (2.3)
One advantage of working with functions in Γ 0 is that the proximal operator is always well-defined. Its computation is efficient for many practical losses thanks to the following proposition.

Proposition 2.6 (Moreau decomposition). Let f ∈ Γ 0 (H) and γ > 0. Then

Id = prox γf (•) + γ prox f ⋆ /γ (•/γ).
(2.4)

Example 2.7. A particular case of interest is the proximal operator associated to the indicator function of a closed convex set C ⊂ H. It holds that for all γ > 0,

∀h ∈ H, prox γχ C (•) (h) = Proj C (h)
where Proj C is the orthogonal projection on C in H.

The proximal operator allows to solve composite problems of the form

inf x∈H f (x) + g(x), (2.5)
where f, g ∈ Γ 0 (H) are such that f is Gâteaux differentiable with C-Lipschitz continuous gradient for some C > 0, while g is not.

Definition 2.8 (Gâteaux differentiability). Let f : H → -∞, +∞ be proper and x ∈ dom(f ). Let y ∈ H. The directional derivative of f in the direction y is

f ′ (x, y) = lim α↓0 f (x + αy) -f (x) α
provided that the limit exists. When f ′ (x, •) is linear in y and continuous, f is said to be Gâteaux differentiable at point x and there exist a unique vector ∇f

(x) ∈ H such that ∀y ∈ H, f ′ (x, y) = ∇f (x), y H .
While in the finite dimensional case, existence of C 1 partial derivatives ensures differentiablity of the function and existence of a gradient, this is no longer true in infinite dimension. Gâteaux differentiablity then helps to fill this void, and paves the way to first order optimization methods such as the proximal gradient descent algorithm, presented in Algorithm 2.1.

The subdifferential is a fundamental tool in convex analysis that generalizes the gradient to nonsmooth functions.

Definition 2.9 (Subdifferential).

Let f : H → -∞, +∞ . The subdifferential of f at x ∈ dom(f ) is: ∂f (x) := {u ∈ R d : ∀y ∈ H, f (y) ≥ f (x) + u, y -x H } , ( 2.6)
i.e. the set of "slopes" of all affine minorants of f which are exact at x.

CHAPTER 2. BACKGROUND ON LEARNING WITH OPERATOR-VALUED KERNELS

Algorithm 2.1 Proximal gradient descent for Problem 2.5 input : Lipschitz constant C, iteration number T init : x (0) for t = 1, . . . , T do

x (t) = prox 1 C g x (t-1) -1 C ∇f (x (t-1) ) return x (T )
If f is a convex Gâteaux differentiable function, the subdifferential at x ∈ dom f has only one element: ∇f (x). Subdifferentiability allows to generalize first-order optimality conditions to non-differentiable convex functions.

Proposition 2.10 (Fermat's rule). Let f ∈ Γ 0 (H). Then, for all x ∈ H:

x ∈ arg min x∈H f (x) ⇔ 0 ∈ ∂f (x).
(2.7)

We are now ready to state the strong duality theorem in the context of minimization under affine equality constraints, which will be used when exploiting duality in vv-RKHSs.

Definition 2.11 (Lagrangian). Let f ∈ Γ 0 (H), K be a Hilbert space, L ∈ L(H, K) and b ∈ K. We are interested in solving

inf x∈H f (x) + χ {0} (Lx -b) := P(x) (2.8)
that we refer to as the primal problem. The Lagrangian associated to Problem 2.8 is

L : H × K → -∞, +∞ (x, α) → f (x) + Lx -b, α K .
We then refer to sup

α∈K inf x∈H L(x, α) := D(α)
as the dual problem. Moreover, (x, α) is a saddle point if

P(x) = L(x, α) = D(α).
We can notice that P(x) = sup α∈K L(x, α) for ∀x ∈ H. Denoting by P * and D * the respective optimal values of the primal and dual problems, weak duality ensures that P * ≥ D * . It turns out that for a large class of problems, both values are equal, a setting refered to as strong duality.

Scalar-valued Kernels and RKHSs

Kernel methods stand as a cornerstone of machine learning. Not only they have played a preponderant role in extending linear models to nonlinear models but they come with a rigorous mathematical framework that eases the analysis of the learning algorithm (see [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF]; Berlinet and Thomas-Agnan (2004); [START_REF] Steinwart | Mercer's theorem on general domains: on the interaction between measures, kernels, and RKHSs[END_REF] for in depth presentation of kernel methods for machine learning).

Their development is strongly linked with convex optimization as practical algorithms make an extensive use of duality principles. They offer powerful modeling spaces, called reproducing kernel Hilbert spaces (RKHSs), which are Hilbert spaces of functions that enjoy a regularity property: convergence in the space implies pointwise convergence.

Equivalently, this can be summarized by the requirement of continuity for all functions evaluations.

Definition 2.13. Let X be any set and H ⊂ F(X, R) be a Hilbert space. H is said to be a reproducing kernel Hilbert space (RKHS) if and only if for all x ∈ X, the following evaluation mapping at x is continuous.

ev x : H → R h → h(x) .
Definition 2.13 is of limited practical interest. Indeed, we seldom have access first to a RKHS H, but rather construct it through the notion of reproducing kernel as presented below.

Definition 2.14. A (scalar-valued) kernel on X is a function k :

X × X → R such that 1. ∀(x, z) ∈ X 2 , k(x, z) = k(z, x), and 2. ∀n ∈ N * , (x i ) n i=1 ∈ X n , (α i ) n i=1 ∈ R n , n i=1 n j=1 α i α j k(x i , x j ) ≥ 0.
Remark 2.15. Throughout this thesis we consider kernels meant in the sense of positive definite functions while noting that the construction can be extended to indefinite kernels leading to Krein spaces [START_REF] Ong | Learning with non-positive kernels[END_REF][START_REF] Huang | Indefinite kernels in least squares support vector machines and principal component analysis[END_REF].

In the rest of this manuscript, we may with a slight abuse simply refer as kernel any scalar-valued kernel. Kernels can be constructed on sets X with very different structures including for instance permutations [START_REF] Jiao | The Kendall and Mallows kernels for permutations[END_REF], measures [START_REF] Cuturi | Semigroup kernels on measures[END_REF], graphs [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF], or time series [START_REF] Cuturi | A kernel for time series based on global alignments[END_REF]. In case of X = R d , we can mention the celebrated Gaussian kernel which reads

∀(x, z) ∈ R 2 , k(x, z) = e -γ x-z 2 2 , (γ > 0).
The following theorem originating from [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] states an equivalence between kernels and the existence of a feature map representation.

Theorem 2.16. A function k : X × X → R is a kernel if and only if there exists a Hilbert space V and a mapping ϕ :

V → H such that ∀(x, z) ∈ X 2 , k(x, z) = ϕ(x), ϕ(z) V .
(2.9)

The space V is then called a feature space and the mapping ϕ a feature map.

CHAPTER 2. BACKGROUND ON LEARNING WITH OPERATOR-VALUED KERNELS

The following theorem links kernels to RKHSs, saying that any kernel can be associated to a unique RKHS, for that reason kernels are sometimes called reproducing kernels.

Theorem 2.17. Let k be a kernel on X. Then there exists a unique Hilbert space of functions

H k ⊂ F(X, R) such that 1. ∀x ∈ X, k(•, x) ∈ H k , and 2. ∀(h, x) ∈ H k × X, h(x) = h, k(•, x) H k .
The mapping x → k( The above properties allow to define RKHSs on spaces X without any structural assumption. It turns out that endowing the input space with a metric structure and a probability measure leads to fertile connections with functional analysis and L 2 spaces.

Integral operators and RKHSs

We begin by introducing the notion of integral operator associated to a kernel on a compact space Θ.

Definition 2.18 (Integral Operator). Let Θ be a compact metric space endowed with a Borel probability measure µ, and k be a continuous kernel on Θ. The integral operator T k,µ is defined as

T k,µ :   L 2 [Θ, µ] → L 2 [Θ, µ] f → θ → Θ f (θ ′ )k(θ, θ ′ )dµ(θ ′ )   .
Remark 2.19. When there is no ambiguity on the considered measure, we omit the dependence in µ and abbreviate T k,µ as T k .

The continuity of k, combined with the compactness of Θ ensure that k is bounded on Θ, and thus T k is continuous and compact. Moreover, the symmetry of k implies that

T k is self-adjoint. Finally, it turns out that T k is positive (i.e. T k f, f L 2 [Θ,µ] ≥ 0 for ∀f ∈ L 2 [Θ, µ]).
The spectral theorem for self-adjoint compact operators guarantees the existence of an at most countable family of functions (ψ j ) j∈J forming an orthonormal system in

L 2 [Θ, µ] such that ∀f ∈ L 2 [Θ, µ], T k f = j∈J λ j f, ψ j ψ j , (2.10)
where (λ j ) j∈J are the family of positive eigenvalues ordered so that

λ 1 ≥ λ 2 ≥ • • • ≥ 0.
Additionally, the eigenvalues (ψ j ) j∈J can be chosen to be continuous since T k maps L 2 [Θ, µ] into the space C(Θ, R) of continuous functions. This provides a Mercer representation of the kernel k, as stated below.

Theorem 2.20 (Mercer). Let Θ be a compact metric space endowed with a Borel probability measure µ such that supp(µ) = Θ, and k be a continuous kernel on Θ. Let (λ j , ψ j ) j∈J be the eigendecomposition of T k with continuous eigenvectors. Then

∀(θ, θ ′ ) ∈ Θ 2 , k(θ, θ ′ ) = j∈J λ j ψ j (θ)ψ j (θ ′ ) (2.11)
where the convergence is uniform on Θ 2 and absolute for all (θ, θ ′ ) ∈ Θ 2 .

Remark 2.21. The compact assumption on Θ can be weakened to "H k is compactly embedded in L 2 [Θ, µ]" while still retaining most of Mercer's theorem, see Steinwart and Scovel (2012) for more details.

While the eigendecomposition of T k depends on µ, Equation (2.11) holds nevertheless.

When supp(µ) = Θ, a nice byproduct of Mercer's theorem is to provide a description of H k in terms of the decay rate of the Fourier coefficients in the basis (ψ j ) j∈J . Indeed, it holds that

H k =        f ∈ L 2 [Θ, µ] : j∈J f, ψ j L 2 [Θ,µ] 2 λ j < +∞        (2.12) so that h ∈ H k if and only if there exist f ∈ L 2 [Θ, µ] such that T k f = h.
Exhibiting the exact eigendecomposition of T k is a notoriously hard century-old problem for general k and µ [START_REF] Stone | Linear transformations in Hilbert space and their applications to analysis[END_REF][START_REF] Klus | Eigendecompositions of transfer operators in reproducing kernel Hilbert spaces[END_REF]. In some cases it can be carried out by solving a differential equation associated to the eigenvector problem, but this method can only be applied on an ad hoc basis.

λ j = 2 1 + c 2 j , ψ j : θ → c j cos(c j θ) + sin(c j θ) where (c j ) ∞ j=1 are solutions to the equation cot(c) = 1 2 c -1 c
where cot is the cotangent function.

If the exact eigendecomposition is impossible to get, one can resort to using approximate decomposition based on a sampling of the integral operator. ∼ µ. We consider the problem of finding a continuous eigenvector ψ of a sampled version of the integral operator T k with eigenvalue λ > 0:

∀θ ∈ Θ, 1 m m j=1 k(θ, θ j )ψ(θ j ) = λψ(θ).
(2.13)

In particular, by evalutating Equation (2.13) at points (θ j ) m j=1 , we get that (λm, ψ(θ j ) m j=1 ) is a pair of eigenvalue/eigenvector associated to the Gram matrix K ∈ M m (R) defined by K ij = k(θ i , θ j ) for all (i, j) ∈ [m] 2 . These can be computed using e.g. Singular Value Decomposition, and by backsubstitution in Equation (2.13) one gets an approximated eigenbasis of dimension at most m that can be used as a proxy in applications where the true eigenvectors of T k are required.

CHAPTER 2. BACKGROUND ON LEARNING WITH OPERATOR-VALUED KERNELS

The decay rate of the eigenvalues (λ j ) j∈J determines "how large" the RKHS is; the slower the decay the larger the space is. An approach to measuring the modeling capacity of the RKHS then consists in answering the question "Can functions in the RKHS approximate any continuous function for the uniform convergence ?". [START_REF] Gretton | A kernel statistical test of independence[END_REF][START_REF] Gretton | A kernel two-sample test[END_REF].

Kernel quadrature rules RKHSs can be used to approximate integrals, a topic of interest for this thesis. Given a measure µ on a compact space Θ ⊂ R d , the goal of quadrature rules is to find points (θ j ) m j=1 ∈ Θ and coefficients (η j ) m j=1 ∈ R m such that for some class of functions H the approximation

∀h ∈ H, Θ h(θ)dµ(θ) ≈ m j=1 η j h(θ j )
is reasonable. In particular, we want the approximation to have better convergence rates than the asymptotic O(m 1/2 ) obtained by Monte-Carlo methods where one uses

(θ j ) m j=1 i.i.d.
∼ µ and η j = 1 m . Taking H as the unit ball of some RKHS, this problem can be cast as the approximation of the mean embedding in the RKHS with well chosen elements [START_REF] Smola | A Hilbert space embedding for distributions[END_REF]. The resulting integration scheme (η j , θ j ) m j=1 is called a kernel quadrature rule and can enjoy asymptotic convergence rates faster than Monte-Carlo, for example of order O(m -r/d ) in a Sobolev space of order r [START_REF] Novak | Deterministic and stochastic error bounds in numerical analysis[END_REF]. Kernel quadrature rules enjoy rich connections with random feature expansions, among them the random Fourier features introduced later, which allows to express the optimal sampling measure of the (θ j ) m j=1 as a leverage score depending on the kernel and µ [START_REF] Bach | On the equivalence between kernel quadrature rules and random feature expansions[END_REF]. Finally we want to mention that it is possible to get convergence rates of the estimate given by kernel quadrature rules even when applied to functions outside the RKHS [START_REF] Kanagawa | Convergence analysis of deterministic kernel-based quadrature rules in misspecified settings[END_REF].

Empirical risk minimization in RKHSs

In the supervised learning setting, given a loss function ℓ : R 2 → R (taken proper, convex, lower semi-continuous with respect to its first argument) and a joint random variable (X, Y) ∈ X × R, an ideal prediction rule h ∈ F(X, R) minimizes R(h) := E (X,Y) ℓ(h(X), Y) over the set of admissible, i.e. measurable, prediction rules. Because the law of (X, Y) is unknown, one relies on i.i.d. samples (x i , y i ) n i=1 and trades the expectation for the empirical mean, to which a regularization term is added to prevent overfitting. The set of admissible prediction rule can be chosen to be a RKHS H k , resulting in the optimization problem

inf h∈H k 1 n n i=1 ℓ h(x i ), y i + λ 2 h 2 H k (λ > 0). ( 2
: X → R (called its signature) such that for all (x, z) ∈ X 2 , k(x, z) = k 0 (x-z).
The RFF model stems from Bochner's theorem [START_REF] Rudin | Fourier Analysis on Groups[END_REF], stating that any continuous, bounded and shift-invariant kernel k : R d × R d → R can be expressed as the Fourier transform of some finite Borel measure ρ k on R d , and the correspondence is CHAPTER 2. BACKGROUND ON LEARNING WITH OPERATOR-VALUED KERNELS one-to-one:

k(x, z) = R d cos ω, x -z dρ k (ω).
Given some integer m and (ω j ) m j=1 i.i.d. sampled from ρ k , the kernel k : X × X → R defined by

∀(x, z) ∈ X 2 , k(x, z) = 1 m m j=1 cos ( ω j , x -z ) (2.17)
is then a reasonable approximation for k (see the work of [START_REF] Sriperumbudur | Optimal rates for random Fourier features[END_REF] for optimal rates), and the associated RKHS H k can be used as hypothesis space in machine learning problems.

The feature map associated to the kernel is then φ : R d → R 2m defined for all x ∈ R d by

φ(x) = 1 √ m (cos (ω ⊤ 1 x), . . . , cos (ω ⊤ m x), sin (ω ⊤ 1 x), . . . , sin (ω ⊤ m x)) ⊤ .
(2.18)

Example 2.30 (RFF integral operator eigendecomposition). Let k be a RFF kernel on R d of type Equation (2.17) for some integer m and (ω j ) m j=1 i.i.d. sampled from ρ k . Let Θ be a compact subspace of R d , endowed with a probability measure µ. The eigendecomposition of T k is computable, as for any

f ∈ L 2 [Θ, µ], it holds that ∀θ ∈ Θ, T kf (θ) = Θ k(θ, θ ′ )f (θ ′ )dµ(θ ′ ) = 1 m m j=1 cos(ω ⊤ j θ) Θ cos(ω ⊤ j θ ′ )f (θ ′ )dµ(θ ′ ) + sin(ω ⊤ j θ) Θ sin(ω ⊤ j θ ′ )f (θ ′ )dµ(θ ′ ).
Plugging in some function h = n j=1 a j cos(ω ⊤ j θ) + b j sin(ω ⊤ j θ) and writing the eigenvector problem, we get that h is an eigenvector of T k if and only if (a 1 , . . . a m , b 1 , . . . b m ) ⊤ is an eigenvector of the matrix Ψ ∈ M 2m (R) encoding the scalar products in L 2 [Θ, µ] of the cosines and sines, defined ∀i, j ∈ [m] 2 by

Ψ i,j = Θ cos(ω ⊤ i θ) cos(ω ⊤ j θ)dµ(θ) Ψ i+m,j+m = Θ sin(ω ⊤ i θ) sin(ω ⊤ j θ)dµ(θ) Ψ i+m,j = Θ sin(ω ⊤ i θ) cos(ω ⊤ j θ)dµ(θ) Ψ i,j+m = Θ cos(ω ⊤ i θ) sin(ω ⊤ j θ)dµ(θ).
Moreover, they also share the same eigenvalues up to the 1 m factor, so that Singular Value Decomposition can be applied on Ψ to get the eigendecomposition of T k. Finally, notice that when µ is the Lebesgue measure, the coefficients of Ψ are computable in closed form, whereas they must be approximated for general µ.

RFF benefit from a large body of work, from pure computational aspects [START_REF] Yang | Quasi-Monte Carlo feature maps for shift-invariant kernels[END_REF][START_REF] Le | Fastfood -computing Hilbert space expansions in loglinear time[END_REF][START_REF] Zhang | Low-precision random Fourier features for memory-constrained kernel approximation[END_REF] to quantifying their effect in learning algorithms (Rudi and Rosasco, 2017;Li et al., 2019b).

Kernel Learning Though giving rise to expressive functional spaces, kernels are hard to tune in practice. The idea of kernel learning proposes to learn the kernel adapted to the task at hand. This can take many forms, including multiple kernel learning [START_REF] Sonnenburg | Large scale multiple kernel learning[END_REF] which proposes to learn a convex combination of kernels, kernel alignment [START_REF] Kandola | Optimizing kernel alignment over combinations of kernel[END_REF] which gives a measure of adequacy between the kernel and the data, spectral kernel learning where the goal is to learn a shift invariant kernel from its fourier spectral measure stemming from Bochner's theorem [START_REF] Oliva | Bayesian nonparametric kernel-learning[END_REF]Li et al., 2019a), and many others. In particular, we want to mention deep kernel learning, which consists in parameterizing a kernel by a neural architecture, and provides a way to train jointly the parameters of the estimator and the neural architecture. This idea has been introduced by [START_REF] Salakhutdinov | Using deep belief nets to learn covariance kernels for gaussian processes[END_REF], where authors use deep beliefs nets to learn covariance kernels for Gaussian processes. The pretraining of the network is carried out on a large, unsupervised dataset to learn relevant features, and the weights are fine-tuned for specific applications involving a smaller amount of data. The resulting kernel family is referred to as deep kernels (DK). Since then deep kernels have been exploited in various kernel-related tasks, such as maximum mean discrepancy for GAN discriminators [START_REF] Li | Towards deeper understanding of moment matching network[END_REF], density estimation via score matching in exponential families [START_REF] Wenliang | Learning deep kernels for exponential family densities[END_REF], semi supervised learning [START_REF] Jean | Semi-supervised deep kernel learning: Regression with unlabeled data by minimizing predictive variance[END_REF], non-parametric two-sample tests [START_REF] Liu | Learning deep kernels for non-parametric two-sample tests[END_REF], among others. The resulting optimization problems are costly, as learning deep kernels involves the joint minimization in the parameters of the neural architecture and in the coefficients of the representer theorem. Thus, most approaches focus on scalable approaches by choosing a kernel that brings computational advantages. We can mention [START_REF] Yang | Deep fried convnets[END_REF] for image classification with random Fourier features approximated using the FastFood method [START_REF] Le | Fastfood -computing Hilbert space expansions in loglinear time[END_REF], [START_REF] Wilson | Deep kernel learning[END_REF] for scaling Gaussian processes up to millions of points, [START_REF] Mehrkanoon | Scalable hybrid deep neural kernel networks[END_REF]; [START_REF] Mehrkanoon | Deep hybrid neural-kernel networks using random Fourier features[END_REF] for general hybrid RFF approaches. Deep kernel learning approaches allow to apply kernel methods to complex data such as images, as emphasized later in Section 5.4.4.

We now present the vv-RKHS extension to kernel methods used for modeling outputs which are not real-valued.

Operator-valued Kernels and vv-RKHSs

Vector-valued RKHSs (vv-RKHSs) are an extension of RKHSs allowing to model functions with outputs in a Hilbert space first developed in [START_REF] Pedrick | Theory of reproducing kernels for Hilbert spaces of vector-valued functions[END_REF]. While the definition and properties of vv-RKHSs are similar to the real output case (see Table 2.2 for a comparison between real and vector-valued cases), the kernels involved are not real-valued anymore but operator-valued. We refer to [START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces of integrable functions and Mercer theorem[END_REF][START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces and universality[END_REF] for an overview on the topic. In what follows, Y a separable Hilbert space. Definition 2.31. Let H ⊂ F(X, Y) be a Hilbert space. H is said to be a vector-valued reproducing kernel Hilbert space if and only if for all x ∈ X, the following evaluation mapping at x is continuous.

ev x : H → Y h → h(x) .
Similarly to the scalar case, vv-RKHSs can be constructed from OVKs.

k : X × X → R K : X × X → L(Y) positivity i,j∈[n] α i α j k(x i , x j ) ≥ 0 i,j∈[n] K(x i , x j )y i , y j Y ≥ 0 function space H k = Span{k(•, x) : x ∈ X} H K = Span{K(•, x)y : (x, y) ∈ X × Y} reproducing property h(x) = h, k(•, x) H k h(x) = K(•, x) # h feature map ϕ : X → V Φ : X → L(Y, V) parametrization h = ϕ(•), v V h = Φ(•) # v 1. ∀(x, z) ∈ X 2 , K(x, z) = K(z, x) ♯ , and 2. ∀n ∈ N * , (x i ) n i=1 ∈ X n , (y i ) n i=1 ∈ Y n , n i=1 n j=1 y i , K(x i , x j )y j ≥ 0.
OVKs can also be characterized by the existence of a feature map.

Theorem 2.33. A function K : X × X → L(Y) is an OVK if and only if there exist a Hilbert space V and a mapping Φ :

X → L(V, Y) such that ∀(x, z) ∈ X 2 , K(x, z) = Φ(x)Φ(z) ♯ .
We again call Φ a feature map and V a feature space. Given an OVK K and x ∈ X, K x : Y → F(X, Y) denotes the linear operator such that ∀y ∈ Y, ∀z ∈ X, K x y(z) = K(x, z)y.

The following theorem describes the vv-RKHS associated to an OVK.

Theorem 2.34. Let K be an OVK on X. Then there exists a unique Hilbert space of functions

H K ⊂ F(X, Y) such that 1. ∀(x, y) ∈ X × Y, K x y ∈ H K , and 2. ∀(h, x) ∈ H k × X, h(x) = K ♯ x h.
Remark 2.35. If H is a vv-RKHS in the sense of Definition 2.31 then its associated kernel is given by K(x, z) = ev x ev ♯ z .

We now define the family of separable kernels which are among the simplest OVKs to work with.

Definition 2.36. An operator-valued kernel K : X × X → L(Y) is called separable if there exists a scalar kernel k : X × X → R and an operator A ∈ L(Y) such that

∀(x, z) ∈ X 2 , K(x, z) = k(x, z)A.
Separable kernels are undoubtedly the most studied OVKs, for their simplicity and computational efficiency. Indeed the scalar product in the associated RKHS H K takes a straightforward form on elementary elements of the form K x y:

K x 1 y 1 , K x 2 y 2 H K = K(x 1 , x 2 )y 1 , y 2 Y = k(x 1 , x 2 ) Ay 1 , y 2 Y for all (x 1 , x 2 , y 1 , y 2 ) ∈ X 2 × Y 2 .
Other notables classes of OVKs include the curl-free and div-free matrix-valued kernels used for vector fields learning [START_REF] Macedo | Learning div-free and curl-free vector fields by matrix-valued kernels[END_REF], the transformable matrixvalued kernels for multi-view data [START_REF] Huusari | Multi-view metric learning in vector-valued kernel spaces[END_REF] and the recent entangled or partial trace OVKs deriving from quantum computing concepts [START_REF] Huusari | Entangled kernels -beyond separability[END_REF]. It is also worth mentioning that extensions of multiple kernel learning to OVKs have been investigated in [START_REF] Kadri | Multiple operator-valued kernel learning[END_REF].

Example 2.37. Let k X : X × X → R be a kernel, and A ∈ M s (R) for some integer s, such that A is a positive self-adjoint operator. Let

K : X × X → M s (R) (x, z) → k X (x, z)A .
Then K is an OVK on X and H K ⊂ F(X, R s ).

While in the scalar case the kernel associated to a RKHS is unique, in the operatorvalued case one may find OVKs acting on different output space that lead to the "same" space of function up to an isometry in the sense defined below.

Definition 2.38. Two Hilbert spaces H 1 , H 2 are said to be unitarily equivalent if there exist

W ∈ L(H 1 , H 2 ) such that W W ♯ = Id H 2 and W ♯ W = Id H 1 .
Example 2.39 [START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces and universality[END_REF], example 6). Let k X : X×X → R and k Θ : Θ×Θ → R be kernels. Let

K : X × X → L(H k Θ ) (x, z) → k X (x, z) Id H k Θ .
Then K is an OVK on X and the corresponding vv-RKHS

H K is unitarily equivalent to H k X ⊗ H k Θ .
Example 2.40 [START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces and universality[END_REF], example 7). Let k X : X×X → R and k Θ : Θ×Θ → R be kernels. Assume that Θ is a compact metric space endowed with a probability measure µ such that supp(µ) = Θ, and that k Θ is continuous. Let

K : X × X → L 2 [Θ, µ] (x, z) → k X (x, z)T k,µ .
Then K is an OVK on X and the corresponding vv-RKHS

H K is unitarily equivalent to H k X ⊗ H k Θ .
Remark 2.41. The modeling spaces from Example 2.39 and Example 2.40 are essentialy the same spaces of functions as the scalar RKHS associated to the product kernel k X ⊗ k Θ . What differs among these choices is whether we look at the elements of these spaces as functions (i) from X to

H k Θ , (ii) from X to L 2 [Θ, µ], or (iii) from X × Θ to R.
We leverage this flexibility in the design of optimization algorithms in Chapter 3.

Similarly to the scalar case, one can define the integral operator associated to an OVK K : Θ×Θ → L(Y) when H K ⊂ C(Θ, Y) and K is bounded in operator norm on supp(µ):

T K,µ :   L 2 [Θ, µ; Y] → L 2 [Θ, µ; Y] f → θ → Θ K(θ, θ ′ )f (θ ′ )dµ(θ ′ )   . ( 2
inf h∈H K 1 n n i=1 L(h(x i ), y i ) + λ 2 h 2 H K (λ > 0), (2.20)
where

K : X×X → L(Y)
ĥ = n i=1 K(•, x i )α i (2.21)
for some

(α i ) n i=1 ∈ Y n .
The benefit of Theorem 2.43 is to pertain the solution to Problem 2.20 to a specific subset of H K . The coefficients (α i ) n i=1 now belong to the Hilbert space Y, which can pose representation problems for infinite dimensional Y. This is discussed at length in Chapter 3.

Kernel methods have a rich history of exploiting parametric duality techniques to solve empirical risk minimization problem. In particular, these techniques translate well to OVK-based learning problems, and can be exploited to gain additional information about the coefficients (α i ) n i=1 . Below we present a result originating from [START_REF] Brouard | Input output kernel regression: supervised and semi-supervised structured output prediction with operator-valued kernels[END_REF] about the dualization of such problems with output in a general Hilbert space Y. We adopt the notation

I ℓ i : y ∈ Y → I ℓ (y, y i ) for any i ∈ [n].

Theorem 2.44 (Dualization, Brouard et al. (2016)). The solution of Problem 2.20 is given by

ĥ = 1 λn n i=1 K(•, x i )α i , (2.22) with (α i ) n i=1 ∈ Y n being the solution of the dual problem inf (α i ) n i=1 ∈Y n n i=1 I ⋆ ℓ i (-α i ) + 1 2λn n i,j=1 α i , K(x i , x j )α j Y .
(2.23)

Again, solving Problem 2.23 may not be straightforward and we refer to Chapter 3 for dedicated techniques in the context of integral losses.

Large scale learning

The vv-RKHS framework is even more challenging than the real-valued one from computational perspective. Approximated kernel schemes can however be applied to lighten it, as proposed by the operator random Fourier features (ORFF) methodology [START_REF] Brault | Infinite task learning in RKHSs[END_REF]. This approach holds for shift-invariant OVK, and relies on an operator-valued version of Bochner's theorem, stating that any shift-invariant Mercer OVK K : R d ×R d → R can be expressed as the Fourier transform of some finite operator-valued borelian measure Q K on R d , and the correspondence is one-to-one:

k(x, z) = R d cos ω, x -z dQ K (ω).
Under mild regularity assumptions verified for separable OVKs (see [START_REF] Brault | Random Fourier features for operatorvalued kernels[END_REF] for a precise statement), the operator-valued measure can be written dQ K (ω) = Q K (ω)dρ K (ω) for an operator-valued function Q K and bounded Borel measure ρ k . Given some integer m and (ω j ) m j=1 i.i.d. sampled from ρ K , the kernel K :

X × X → L(Y) defined by ∀(x, z) ∈ X 2 , K(x, z) = 1 m m j=1 cos ( ω j , x -z )Q K (ω j ) (2.24)
can be used as hypothesis space in machine learning problems. The feature map associated to the kernel is then Φ :

R d → L Y, Y 2m defined for all x ∈ R d by Φ(x) = 1 √ m m i=1 cos (ω ⊤ i x)B(ω j ) ♯ ⊕ sin (ω ⊤ m x)B(ω j ) ♯ (2.25) where ∀j ∈ [m], B(ω j )B(ω j ) ♯ = Q(ω j ).
In the case of a separable kernel K = k X A, dQ K = Adρ k so that it suffices to exhibit a decomposition BB ♯ = A for some B ∈ L(Y) and the feature map writes

Φ(x) = 1 √ m φX (x) ⊗ B ♯ . ( 2

Conclusion

In this chapter we provided a background in both convex optimization and vector-valued RKHSs theory, with focus on machine learning applications. In short, vv-RKHSs are extensions of RKHSs allowing to model functions with outputs in any Hilbert space, enabling the model of function-valued functions in particular. These spaces of functions are chosen as hypothesis spaces for the regularized empirical risk minimization problems raised in the rest of the thesis. In this chapter, we address the general question of solving regularized empirical risk minimization problems in vv-RKHSs with functional outputs, in the presence of integral losses. After introducing the family of considered problems in Section 3.1, we devote Section 3.2 to the study of the special case of the square loss function, for which closed-form solutions exist. In Section 3.3 we investigate ways to solve the problems based on the primal formulation, either by means of a double representer theorem or by the use of random features. We then delve into dual algorithms in Section 3.4, designing techniques involving adapted basis related to the operator-valued kernels used in practice, or exploiting linear splines bases to represent the dual variables, rephrasing the problem as the minimization of an approximated quadratic form under specific linear constraints. This chapter ends with a conclusion in Section 3.5.

Optimization Schemes for Learning with Integral Losses

Problem Formulation

Let (X, Y) be a pair of random variables taking values in a product space X×L 2 [Θ, µ; U]. Depending on the application considered, X will be R d or a functional space. The space Θ ⊂ R p is a compact domain endowed with a Borel probability measure µ, and U = R s for some integer s > 0. We begin by proposing a family of loss functions referred to as integral losses that involves an integration over the space Θ.

CHAPTER 3. OPTIMIZATION SCHEMES FOR LEARNING WITH INTEGRAL LOSSES

Integral Losses Losses on L 2 [Θ, µ; U] can be constructed from a family of ground losses on U parameterized by θ ∈ Θ. Let ℓ : Θ × U × U → R be such a family, and define a loss function I ℓ on L 2 [Θ, µ; U] 2 by integration with respect to µ:

I ℓ (f, g) = Θ ℓ(θ, f (θ), g(θ))dµ(θ), for f, g ∈ L 2 [Θ, µ; U]. (3.1)
The integrated loss I ℓ encodes the cost associated to the prediction of f given the ground truth g. As f and g are functions with domain Θ, the prediction error is measured pointwise (∀θ ∈ Θ) as ℓ(θ, f (θ), g(θ)); in I ℓ this error is aggregated on Θ with weighting µ (see Table 3.1 for a summarized view of all these parameters in different scenarii). We require that ℓ is a convex lower semi-continuous integrand; that is for all θ ∈ Θ, v ∈ U, ℓ(θ, •, v) is lower semi-continuous and convex. The integrand is also chosen to be proper:

∀(θ, u, v) ∈ Θ×U 2 , ℓ(θ, u, v) = -∞
and is not everywhere +∞. We assume in the sequel and refer to these requirements jointly as normal convex integrand. To avoid indefinite integral issues, by convention I ℓ (f, g) = +∞ if any of the positive or negative part of the integrand is infinite.

One of the challenges to tackle is that even computing I ℓ (f, g) in Equation (3.1) is not straightforward. Our goal is to minimize the risk (or in practice its empirical version) over a hypothesis class H (detailed below)

inf h∈H R(h) := E (X,Y) [I ℓ (h(X), Y)]. (3.2)
It is instructive to consider a few examples for Problem 3.2; the corresponding (X, Θ, U, ℓ) choices are summarized in Table 3.1.

• Functional output regression (FOR, Chapter 4): In this problem family the task is to regress to a functional output (Θ = [0, 1] is endowed with a probability measure

µ, U = R) from a vectorial (X = R d ) or functional input (X = L 2 [Θ 0 , µ 0 ]). A natural way to tackle this problem is to estimate E[Y|X] = h † (X)
where it is well-known that

h † = arg min h∈H E (X,Y) Y -h(X) 2 L 2 [Θ,µ]
, with H being the set of all measurable functions from

X to L 2 [Θ, µ]. As f 2 L 2 [Θ,µ] = Θ f 2 (θ) dµ(θ) for any f ∈ L 2 [Θ, µ], finding h † reduces to Problem 3.2 with ℓ(θ, u, v) = 1 2 (v -u) 2 .
• Joint quantile regression (JQR, Section 5.4): While classically quantiles are regressed for a single level [START_REF] Koenker | Regression quantiles[END_REF], one can consider the learning problem of simultaneously learning multiple quantiles [START_REF] Sangnier | Joint quantile regression in vector-valued RKHSs[END_REF] or that of the whole quantile function [START_REF] Brault | Infinite task learning in RKHSs[END_REF].

Particularly, let Θ = [0, 1], U = R and consider a pair of random variables (X, Y) ∈ R d × R or (X, Y) ∈ R d 1 ×d 2 × R. Here the r.v. Y is real valued, but is identified with a random variable of constant functions in L 2 [Θ, µ]. The θ-conditional quantile of Y|X for θ ∈ (0, 1) is defined as q(x)(θ) = inf{u ∈ R : P(Y ≤ u | X = x) ≥ θ} = arg min u∈R E Y [ℓ(θ, u, Y) | X = x], (3.3)
where in the variational description ℓ is the so-called pinball loss

ℓ(θ, u, v) = max (θ(v -u), (θ -1)(v -u)).
Task X Θ U ℓ(θ, u, v) QR (Section 5.4) R d [0,1] R max (θ(v -u), (θ -1)(v -u)) CSC (Section 5.5) R d [0, 1] R θ -{-1} (v) max(0, 1 -uv) FOR (Chapter 4) R d or L 2 [Θ 0 , µ 0 ] [0, 1] R 1 2 (u -v) 2 ET (Chapter 6) R d or R d 1 ×d 2 ⊂ R p X 1 2 u -v 2 U The variational form in Equation (3.3) implies that for any h ∈ H Θ E X [E Y [ℓ(θ, q(X)(θ), Y) | X]]dµ(θ) ≤ Θ E X [E Y [ℓ(θ, h(X)(θ), Y) | X]]dµ(θ).
This means that R(q) ≤ R(h) for any h ∈ H in Problem 3.2.

• Cost-sensitive classification (CSC, Section 5.5

): Let Θ = [0, 1], U = R and (X, Y) be random variables taking values in R d × {-1, 1}.
We denote c + (resp. c -) the cost associated to a false positive (resp. negative). The cost-sensitive hinge loss is defined as by ℓ hinge (c • Emotion transfer (ET, Chapter 6): Let Θ ⊂ R p be a compact set corresponding to an embedding space for emotions endowed with a probability µ measuring the importance of each emotion, and Y be a random variable taking its values in C(Θ, U)

+ , c -, y, y ′ ) = (c + {1} (y ′ ) + c -{-1} (y ′ )) max (0, 1 -yy ′ ) where (y, y ′ ) ∈ R × {-1, 1}. Given x ∈ X, a classification rule is obtained by minimizing E Y [ℓ(c + , c -, s, Y ) | X = x] for s ∈ R,
where

U = R d or U = R d 1 ×d 2 .
The random variable Y encodes the trajectory of a phenomenon (facial landmarks, face picture) with respect to the emotions. Let ϑ be a random variable on Θ with probability ν, independent from Y, and consider

X = Y(ϑ) the X-valued random variable (X = U). Estimating Y from X can be tackled by finding E[Y|X] = h † (X)
where

h † = arg min h∈H E (X,Y) Y -h(X) 2 L 2 [Θ,µ;U]
, with H being the set of all measurable function from

X to L 2 [Θ, µ; U]. Therefore, finding h † reduces to Problem 3.2 with ℓ(θ, u, v) = 1 2 v -u 2 U .
Hypothesis Space We make use of the modeling choices described in Example 2.37 and Example 2.39, given some kernels k X , k Θ , and output similarity encoding matrix A. Let H G be the vv-RKHS associated with the decomposable OVK G = k Θ A. Elements in H G model the Θ → U mapping. We then specify H K to be the vv-RKHS associated INTEGRAL LOSSES

to kernel K = k X Id H G . Elements h ∈ H K model the relation h : X → (Θ → U ∈H G ).
As seen in Remark 2.41, there is an equivalence between the views

K = k X Id H G and K = k X T G where T G is the integral operator associated to G defined in Equation (2.19).
Depending on the situation at hand, we may use one or the other representation, and denote by a unified notation Y the output space of functions in H K (either

H G or L 2 [Θ, µ; U]).
Choosing

H := H K in Problem 3.2 one arrives at inf h∈H K R(h) := E (X,Y) I ℓ (h(X), Y) . (3.4)
In practice, we do not have access to the true distribution of (X, Y), but we are given i.i.d. samples gathered in a dataset S := (x i , y i ) n i=1 . These samples define the empirical risk

R S (h) := 1 n n i=1 I ℓ (h(x i ), y i ),
and one can solve

inf h∈H K R S (h) + λ 2 h 2 H K (λ > 0). (3.5)
The term λ 2 h 2 H K acts as a regularizer. It can be interpreted two ways: from a statistical learning point of view, it is used to prevent over-fitting whereas from an optimization point of view it brings strong convexity to the problem and ensures the coercivity of the functional to minimize.

One of the advantages of working with vv-RKHSs is to have a fairly general representer theorem which ensures that the solution belongs to a specific subset of H K (see Theorem 2.43). Particularly, we know that Problem 3.5 admits a unique solution ĥ ∈ H K which can be written as ĥ

= n i=1 K(•, x i )α i (3.6)
for some

(α i ) n i=1 ∈ Y n . When Y is finite-dimensional,
this is enough to parameterize the solution and represent it on a computer. In our case, Y is a potentially infinite-dimensional functional space raises an additional challenge on the parameterization of the (α i ) n i=1 challenging. Therefore plugging Equation (3.6) into Problem 3.5 is not a viable choice to directly solve it, as we would be left with an optimization problem over an infinite-dimensional space Y n . Parametric duality can be exploited in the context of vv-RKHSs (see Theorem 2.44) and the dual problem takes the form inf where

(α i ) n i=1 ∈Y n n i=1 I ⋆ ℓ i (-α i ) + 1 2λn n i,j=1 α i , K(x i , x j )α j Y , (3.7)
I ℓ i : y ∈ Y → I ℓ (y, y i ) for any i ∈ [n].
Given that (α i ) n i=1 are solution to Problem 3.7, the resulting estimator is

ĥ = 1 λn n i=1 K(•, x i )α i , (3.8)
One can notice that Theorem 2.44 provides a representer expression, with the additional benefit of having access to a dual problem, which can be easier to solve than the direct plug-in of Equation (3.6) into Problem 3.5. The search for the solution in H K is again transferred to the search for optimal coefficients in Y n , which may be intractable in practice for infinite-dimensional Y.

Solving Problem 3.5 is thus a challenging task, and requires dedicated techniques. Two questions are raised here:

• How can we represent the coefficients (α i ) n i=1 so that the resulting estimator admits a finite-dimensional parametrization ?

• How can we learn the best estimator when we cannot even compute exactly each I ℓ (h(x i ), y i )?

These two questions are interlinked, and we propose in this chapter several algorithms to solve Problem 3.5. In Section 3.2 we focus on the particular case of the square loss, for which estimators enjoy closed-form solution. This is an extension of existing results on kernel ridge regression for real-valued output functions. In Section 3.3, we devise general optimization algorithms for any integral loss, using the primal version of Problem 3.5. In particular, we show that sampling the integral loss guarantees a finite parametrization of the estimator by means of a double representer theorem and allows at the same time the approximate computation of the loss function, enabling gradient descent algorithms. We then consider operator random Fourier features (ORFF) models that naturally enjoy finite parameterization and propose a stochastic gradient descent algorithm bypassing the need to compute each I ℓ (h(x i ), y i ). Finally, Section 3.4 is devoted to the exploration of dual algorithms. We investigate properties of the Fenchel-Legendre conjugate of integral losses on Y, and show that the flexibility of the vv-RKHS modeling is well-suited to their computation. We then propose to solve Problem 3.7 by leveraging appropriate representation bases for the dual variables, either based on the eigendecomposition of T G or using linear splines depending on the compatibility between the loss function and the output vv-RKHS H G .

The Special Square Loss Case

In the following, we derive closed-form expressions for the solution of the regularized empirical risk minimization with the square loss in two different settings: in case of modeling the observations (y i ) n i=1 as functions (Section 3.2.1) or to assume access to a sampled version of these (Section 3.2.2). Ultimately, even when modeling the outputs as functions they are represented as multi-dimensional vectors and relevant quantities have to be estimated from these vectors. This section corresponds to the extension of [START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel Hilbert spaces[END_REF][START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF][START_REF] Kadri | Operatorvalued kernels for learning from functional response data[END_REF] beyond real-valued output functions.

Functional Observation Case

In the case where ℓ(θ, u, v) = 1 2 u -v 2 U , one can notice that

I ℓ (f, g) = 1 2 f -g 2 L 2 [Θ,µ;U] for ∀f, g ∈ L 2 [Θ, µ; U].
Taking Y = L 2 [Θ, µ; U] and leveraging the view K = k X T G , leads to an easy-to-compute Fenchel-Legendre transform:

I ⋆ ℓ i (α) = 1 2 α 2 Y + α, y i Y
for all α ∈ Y. This is a direct consequence of the squared norm being the only fixed point of the Fenchel-Legendre transform [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF], as well as the fact that for any f :

Y → R and y ∈ Y, f (• -y) ⋆ = f ⋆ + •, y Y . Consequently, Problem 3.7 writes as inf (α i ) n i=1 ∈Y n n i=1 1 2 α i 2 Y -α i , y i Y + 1 2λn n i,j=1 k X (x i , x j ) α i , T G α j Y .
(3.9)

Before diving deeper into the solution, we introduce two useful notations: we denote by K X ∈ M n (R) the Gram matrix of points (x i ) n i=1 with kernel k X and by y = y i n i=1 ∈ Y n the aggregation of the observed outputs (y i ) n i=1 .

Proposition 3.1 (Closed-Form Functional Ridge). The solution ĥ of Problem 3.5 with the square loss is given by ĥ

= 1 λn n i=1 k X (•, x i )T G αi ,
where α := αi n i=1 ∈ Y n is the solution of the linear system

K X ⊗ T G + λn Id Y n α = nλy. (3.10)
Proof Using the K X and y notations, define

J(α) = 1 2 α Y n -α, y Y n + 1 2λn α, K X ⊗ T G α Y n , α ∈ Y n .
Using this notation Problem 3.9 corresponds to solving inf α∈Y n J(α). J is a (strongly) convex Gâteaux differentiable function (see Definition 2.8), and

∇J(α) = α -y + 1 λn K X ⊗ T G α.

THE SPECIAL SQUARE LOSS CASE

47

Setting the gradient to zero at optimum α ∈ Y n , one gets

α -y + 1 λn K X ⊗ T G α = 0 which translates to K X ⊗ T G + λn Id Y n α = λny.

Remark 3.2 (Relationship to Kadri et al. (2016)). The problem considered in Kadri et al. (2016) is a real-valued function-to-function regression problem using the squared norm in

L 2 [0, 1] (U = R, µ is the Lebesgue measure on Θ = [0, 1]).
The authors leverage a representer theorem associated to a directional derivative argument to obtain a closed-form solution which corresponds to the solution presented here with A = 1 ∈ R.

Solving Equation (3.10) can be carried out at the price of an operator inversion:

α = λn K X ⊗ T G + λn Id Y n -1 y. (3.11)
As a reminder, for the OVK G considered here we have

T G = A ⊗ T k Θ .
Computing the inverse operator in Equation (3.11) is intractable in practice for the large majority of output kernels k Θ because of the complexity of T k Θ . However, one can turn to spectral methods to get an approximate solution. This idea was first exploited in [START_REF] Kadri | Operatorvalued kernels for learning from functional response data[END_REF] in the setting U = R, Θ = [0, 1] and can be adapted to the more general case proposed here. It relies on the following observation: the eigen decomposition of K X ⊗ T G can be obtained from the eigen decompositions of K X and T G . Indeed, if (v, ψ) ∈ R n × Y are respectively eigenvectors of K X and T G with associated eigenvalues (σ, λ) then v ⊗ ψ is an eigenvector of K X ⊗ T G with eigenvalue σλ and all eigenvectors of K X ⊗ T G can be obtained this way.

Let us assume access to the eigen decomposition (λ j , ψ j ) m j=1 of rank m of the integral operator T G , similarly obtained using the eigen decomposition of A and T k Θ . A discussion about the eigen decomposition of

T k Θ is presented in Section 2.2.1. Let (σ i , v i ) n i=1 ∈ (R × R n
) n be the eigen decomposition of the matrix K X obtained by for instance using SVD. We propose to use

α = n i=1 m j=1 1 λ + σ i λ j y, v i ⊗ ψ j Y n v i ⊗ ψ j .
(3.12)

The cost associated to the eigen decomposition of K X (respectively A) scales as O(n 3 ) (respectively O(s 3 )). Finding the first m pairs of (eigenvalue, eigenvector) of T G can be hard and costly for general output kernels k Θ . One can then turn to using approximate eigen decomposition as presented in Example 2.23 or use the eigen decomposition of a Random Fourier Feature model from Example 2.30.

Remark 3.3. Equation (3.12) involves the scalar products y i , ψ j Y , which have to be estimated from the data.

We now investigate ways to solve Problem 3.5 in the square loss setting when one does not assume full access to the observed output functions (y i ) n i=1 referred to as the partial observation setting.

Partially Observed Case

In various cases, due to the functional nature of the (y i ) n i=1 , one cannot assume to have access to the whole functions (in the analogous signal sense, or being able to query at any θ ∈ Θ), and that makes the quantity R S impossible to compute. This will be the case in particular for the emotion transfer setting in Chapter 6. In this partially observed setting, we assume that each y i is observed at certain number (m ∈ N * ) of locations (θ ij ) m j=1 . Remark 3.4 (General Case). One could have considered that the number of observed locations varies with the numbering of the observation i, so that each y i is observed at locations (θ ij ) m i j=1 . To keep the notations simple, we consider that the number of locations observed per sample is fixed (m = m i for all i).

We moreover assume that for ∀i ∈

[n], (θ ij ) m j=1 i.i.d.
∼ µ, and define the sampled empirical risk as

R S (h) := 1 nm n i=1 m j=1 1 2 y i (θ ij ) -h(x i )(θ ij ) 2 U . (3.13)
. This choice ensures that R S (h) converges towards R(h) when the number of samples and locations grows towards infinity.

Remark 3.5 (Re-weighting Scheme). In the case where the locations are not sampled with probability µ, then there is no convergence of the sampled empirical risk towards the true risk. To remedy this bottleneck, one can consider a re-weighted sampled empirical risk of the form

R S (h) := 1 n n i=1 m j=1 η ij 2 y i (θ ij ) -h(x i )(θ ij ) 2 U ,
where (η ij ) n,m i,j=1 encodes a re-weighting scheme ensuring that the proposed sampled empirical risk converges towards the risk in Problem 3.4. The solution developed below can then be adapted to such settings.

We now introduce the problem based on the minimization of the sampled empirical risk:

inf h∈H K R S (h) + λ 2 h 2 H K . (3.14)
It is advantageous in this case to adopt the view K = k X Id H G as it helps in the derivation of the dual problem of Problem 3.14.

Proposition 3.6. The solution of Problem 3.14 is given by ĥ

(x)(θ) = 1 λnm n i=1 m j=1 k X (x, x i )k Θ (θ, θ j )Aα ij , ∀(x, θ) ∈ X × Θ, (3.15) with (α ij ) n,m i,j=1 ∈ U nm being the solution of the dual problem inf (α ij ) n,m i,j=1 ∈U nm n i=1 m j=1 1 2 α ij 2 U -y i (θ ij ), α ij U + 1 2λnm n i 1 ,i 2 =1 m j 1 ,j 2 =1 k X (x i 1 , x i 2 )k Θ (θ i 1 j 1 , θ i 2 j 2 ) α i 1 j 1 , Aα i 2 j 2 U .
(3.16)

Proof We introduce slack variables ξ := (ξ ij ) n,m i,j=1 ∈ U nm so that Problem 3.14 can be rewritten as

inf h∈H K ,ξ∈U nm n,m i,j=1 ℓ ij (ξ ij ) + λnm 2 h 2 H K s.t. ξ ij = h(x i )(θ ij ) ∀(i, j) ∈ [n] × [m]
where we adopt the notation

ℓ ij = • -y i (θ ij ) 2 U
. Denoting by α := (α ij ) n i,j=1 the dual variables, the Lagrangian associated to the problem is

L h, ξ, α = n,m i,j=1 ℓ ij (ξ ij ) + λnm 2 h 2 H K + α ij , ξ ij -h(x i )(θ ij ) U .
(3.17)

The dual function then writes as

g(α) = inf h∈H K ,ξ∈U nm L h, ξ, α = inf ξ∈U nm n,m i,j=1 ℓ ij (ξ ij ) + α ij , ξ ij U + inf h∈H K λnm 2 h H K - n,m i,j=1 α ij , h(x i )(θ ij ) U .
In the first term, we can recognize the Fenchel-Legendre conjugate of the squared norm:

inf ξ∈U nm n,m i,j=1 ℓ ij (ξ ij )+ α ij , ξ ij U = - n,m i,j=1 ℓ ⋆ ij (-α ij ) = -    n,m i,j=1 1 2 α ij 2 U -α ij , y i (θ ij ) U    .
For the second term, differentiating with respect to h yields the optimum

ĥ = 1 λnm n,m i,j=1 K(•, x i )G(•, θ ij )α ij
by using the reproducing property successively in H K and H G . It then holds that λnm 2 ĥ

H K - n,m i,j=1 α ij , ĥ(x i )(θ ij ) U = - 1 2λnm n i 1 ,i 2 =1 m j 1 ,j 2 =1 k X (x i 1 , x i 2 )k Θ (θ i 1 j 1 , θ i 2 j 2 ) α i 1 j 1 , Aα i 2 j 2 .
Finally the dual problem writes as

sup α∈U nm g(α) = -inf α∈U nm -g(α)
and by plugging in the expression of g(α) one arrives at Problem 3.16.

Proposition 3.6 allows to parameterize the solution to Problem 3.14 by a finite number of elements of the finite-dimensional space U. Thus, the solution can be represented on a computer, and solving Problem 3.16 boils down to the minimization of a quadratic form, solvable in closed-form.

Particularly, let us introduce three matrices encoding information about the dataset and dual variables. Let K ∈ M nm (R) be the Gram matrix associated to the problem, whose entries are defined as

K m(i 1 -1)+j 1 ,m(i 2 -1)+j 2 = k X (x i 1 , x i 2 )k Θ (θ i 1 j 1 , θ i 2 j 2 ), (i 1 , i 2 ) ∈ [n] 2 , (j 1 , j 2 ) ∈ [m] 2 . (3.18) Let Y ∈ M nm,s (R) be a matrix gathering observations ∀(i, j) ∈ [n] × [m], Y m(i-1)+j,: = y i (θ ij ) ⊤ . (3.19)
The dual variables

(α ij ) n,m i,j=1 ∈ U nm are collected in a matrix α ∈ M nm,s (R) such that ∀(i, j) ∈ [n] × [m], α m(i-1)+j,: = α ⊤ ij .
Using these notations, Problem 3.16 writes as

inf α∈Mnm,s(R) Tr 1 2 αα ⊤ -αY ⊤ + 1 2λnm KαAα ⊤ . (3.20)
It turns out that Problem 3.20 is akin to a Sylvester equation for which solutions can be computed in closed-form, as it is summarized in the following lemma. When A = Id s , the solution is analytic:

α = λnm K + λnm Id nm -1 Y. (3.22)
Proof The statement follows by setting the gradient of the objective function to zero. Indeed,

∇ α → 1 2 Tr αα ⊤ = α, ∇ α → Tr -αY ⊤ = -Y, ∇ α → 1 2λnm Tr KαAα ⊤ = 1 λnm KαA.
Remark 3.8. Lemma 3.7 writes as a generalization of the method developed in [START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel Hilbert spaces[END_REF]; [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF] to the vector-valued outputs case and non-regular grid.

If the matrix A is not the identity matrix, one can still solve Equation (3.21) by using dedicated solvers (El [START_REF] Guennouni | Block Krylov subspace methods for solving large Sylvester equations[END_REF]. We now move on to the general case of Problem 3.5 when the loss function no longer permits explicit closed-form solution.

SOLVING IN THE PRIMAL

51

Solving in the Primal

In this section, we consider the general case of any proper, convex lower semicontinuous integral loss and propose optimization algorithms based on the primal formulation of Problem 3.5, that reads as

inf h∈H K 1 n n i=1 I ℓ (h(x i ), y i ) + λ 2 h 2 H K .
The classical representer theorem from Theorem 2.43 states that there exist (α i ) n i=1 ∈ Y n such that the solution ĥ to Problem 3.5 takes the form

ĥ = n i=1 K(•, x i )α i .
However in our case, Y is a (potentially) infinite-dimensional Hilbert space, which raises the problem of being able to provide an expression of ĥ that can be explicitly calculated. Moreover, the computation of the loss function I ℓ is not straightforward, as it involves an integral over the Θ domain. To circumvent these difficulties, we propose below various optimization algorithms dedicated to solving Problem 3.5. Section 3.3.1 is devoted to solving an approximated problem based on a sampling of the integral loss I ℓ for which a finite-dimensional parametrization can be obtained by means of a double representer theorem, killing two birds with one stone. In Section 3.3.2, we explore the use of Random Fourier Features, with the double benefit of lowering the computational cost associated to the solution as well as benefiting by construction of a finite-dimensional representation of the solution.

Sampling Schemes and Representer Theorems

Sampling I ℓ appears as a natural approach to solve Problem 3.5. It provides the immediate benefit of making the loss function computable; the resulting sampled empirical risk is

R S (h) := 1 n n i=1 m j=1 η ij ℓ(θ ij , h(x i )(θ ij ), y i (θ ij )), (3.23)
where (θ ij , η ij ) n,m i,j=1 encodes information about a chosen sampling scheme. We are now interested in solving an approximated counterpart of Problem 3.5 based on using Equation (3.23) as a proxy for the empirical risk:

inf h∈H K R S (h) + λ 2 h 2 H K . (3.24)
Below we discuss various valid sampling schemes before stating the double representer theorem which allows to represent numerically the solution of Problem 3.24.

Monte-Carlo Sampling:

The Monte-Carlo (MC) method simply consists in sampling

i.i.d. random variables (θ j ) m j=1 i.i.d
∼ µ and estimating the integral I ℓ by its empirical counterpart

I ℓ (f, g) := 1 m m j=1 ℓ(θ j , f (θ j ), g(θ j )). INTEGRAL LOSSES
The MC methods have a probabilistic error scaling as O 1 √ m (with constants depending on the variance of the integrand), which is independent from the dimensionality of Θ. As R S consists in the sum of n terms involving I ℓ , a MC sampling scheme for R S is described by some

(θ ij ) n,m i,j=1 where ∀i ∈ [n], (θ ij ) m j=1 i.i.d. ∼ µ. By setting ∀(i, j) ∈ [n] × [m], η ij = 1
m we arrive at Equation (3.23). It is worth noting at this point that using m = 1 is an interesting choice: it amounts to a very loose approximation of each I ℓ (h(x i ), y i ) but as the number of samples grows it still holds that

lim n→∞ R S (h) = R(h).
Quasi Monte-Carlo Sampling: The idea behind quasi Monte-Carlo (QMC) sampling is to use a low discrepency sequence to estimate the integral. Such sequences provide successively finer uniform partitions of the compact space Θ and guarantee that the proportion of points falling in a subset B ⊆ Θ is close to µ(B). We give an example of such sequence called the Sobol sequence [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF], exemplified in case of Θ = [0, 1]. Start with θ 1 = 1 2 . Then build θ 2 = 1 4 , and θ 3 = 3 4 . Then θ 4 = 1 8 , θ 5 = 3 8 , and so on as to obtain a regular covering of Θ. Such sequence can be adapted to higher dimensional Θ.

Let η j = 1 m F -1 (θ j )
where µ is assumed to be absolutely continuous w.r.t. the Lebesgue measure and F is its associated c.d.f. , and (θ j ) m j=1 is a QMC sequence. We then define

I ℓ (f, g) = m j=1 η j ℓ(θ j , f (θ j ), g(θ j )).
The QMC methods has an approximation error scaling as O (log m) p m , which is faster that what is obtained using MC for many practical values of (m, p). However, the constants involve the Hardy-Kraus variation of the integrand, which is often hard to estimate. More details about this are given in Section 5.3. The QMC sampling scheme leads to Equation (3.23) with (η ij , θ ij ) n,m i,j=1 that are independent of i and simply consists in stacking the same locations and coefficients defined by µ and the Sobol sequence.

Kernel Quadrature Rules: Quadrature rules are suited for low-dimensional Θ. They are described in Section 2.2.1 and lead to an approximation

I ℓ (f, g) = m j=1 η j ℓ(θ j , f (θ j ), g(θ j )),
where (η j , θ j ) m j=1 are the set of weights and locations produced by the quadrature rule. By stacking n of them we recover Equation (3.23).

To solve Problem 3.24, we adopt the view K = k X Id H G , as it allows to apply a representer theorem on the coefficients (α i ) n i=1 ∈ Y n themselves, as expressed in the theorem below.

Theorem 3.9 (Double Representer). Problem 3.24 has a unique solution ĥ and it takes the form ĥ

(x)(θ) = n i=1 m j=1 k X (x, x i )k Θ (θ, θ ij )Aα ij , ∀(x, θ) ∈ X × Θ (3.25) for some coefficients αij ∈ U with i ∈ [n] and j ∈ [m]. Proof Let J(h) = R S (h) + λ 2 h 2 H K
denote the objective function to minimize. For all g ∈ H G , let K x g denote the function defined by (K x g)(t) = K(t, x)g ∀t ∈ X. Similarly, for all α ∈ U, G θ α stands for the function t → G(t, θ)α where t ∈ Θ. Let us take the finite-dimensional subspace

E = span K x i G θ ij α : i ∈ [t], j ∈ [m], α ∈ U .
The space H K can be decomposed as E and its orthogonal complement:

E ⊕E ⊥ = H K .
The existence of ĥ follows from the coercivity of J (i.e. J(h) → +∞ as h H K → +∞) which is the consequence of the quadratic regularizer and the lower boundedness of ℓ. Uniqueness comes from the strong convexity of the objective. Let us decompose ĥ = ĥE + ĥE ⊥ , and take any α ∈ U.

Then ∀(i, j) ∈ [n] × [m], ĥE ⊥ (x i )(θ ij ), α U (a) = ĥE ⊥ (x i ), G θ ij α H G (b) = ĥE ⊥ , K x i G θ ij α ∈E H K (c) = 0.
(a) follows from the reproducing property in H G , (b) is a consequence of the reproducing property in H K , and (c) comes from the decomposition

E ⊕ E ⊥ = H K . This means that ĥE ⊤ (x i )(θ ij ) = 0 ∀(i, j) ∈ [n] × [m], and hence R S ( ĥ) = R S ( ĥE ). Since λ ĥ 2 H K = λ ĥE 2 H K + ĥE ⊥ 2 H K ≥ λ ĥE 2 H K
we conclude that ĥE ⊤ = 0 and get that there exist coefficients αij ∈ U such that

ĥ = n i=1 m j=1 K x i G θ ij αij .
This evaluates for all (x, θ) ∈ X × Θ to ĥ

(x)(θ) = n i=1 m j=1 k X (x, x i )k Θ (θ, θ ij )Aα ij as claimed in Equation (3.25).
The benefit of Theorem 3.9 is to provide a finite-dimensional parameterization of the solution ĥ. Equation (3.25) then allows to write the predictions of the model in a compact form involving the Gram matrix of the problem and the similarity matrix A.

To that end, define K ∈ M nm (R) to be the Gram matrix whose entries are defined as

K m(i 1 -1)+j 1 ,m(i 2 -1)+j 2 = k X (x i 1 , x i 2 )k Θ (θ i 1 j 1 , θ i 2 j 2 ), (i 1 , i 2 ) ∈ [n] 2 , (j 1 , j 2 ) ∈ [m] 2 .
(3.26) Then, denoting by α ∈ M nm,s the matrix encoding information about the coefficients such that for all i, j ∈

[n] × [m], αi(m-1)+j,: = α⊤ ij , it holds that ĥ(x i )(θ ij ) = (K αA) i(m-1)+j,: , (i, j) ∈ [n] × [m].
(3.27)

Remark 3.10. Depending on the chosen sampling, the Gram matrices associated to the problem is structured differently. In particular, in the QMC and the quadrature rule schemes, the sampling locations are the same for all samples and the Gram matrix K has a tensorial structure K = K X ⊗ K Θ which can be used to speed up computations, and reduce the memory footprint. INTEGRAL LOSSES Associated to the fact that R S can be computed (the integration has been traded against a summation), we now have all ingredients needed to apply classical optimization tools to solve Problem 3.24. Indeed, defining

L(r) := 1 nm n i=1 m j=1 ℓ(θ ij , r m(i-1)+j,: , y i (θ ij )), r ∈ M nm,s
we can restate Problem 3.24 as

inf α∈Mnm,s L(KαA) + λ 2 Tr KαAα ⊤ . (3.28)
If one assumes access to the subgradients of ℓ, then by the separability of L it holds that,

∂L(r) m(i-1)+j,: = ∂ℓ(θ ij , •, y i (θ ij ))(r m(i-1)+j,: ), ∀r ∈ M nm,s , ∀(i, j) ∈ [n] × [m].
Problem 3.28 can be then be tackled by (sub)gradient descent methods making use of

∂ α → L(KαA) = K∂L(KαA)A, ∇ α → λ 2 Tr KαAα ⊤ = λKαA.
This allows to perform various gradient descent algorithms; the specific technique can be for instance vanilla (sub)gradient descent, accelerated gradient descent [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF], or quasi-Newton methods [START_REF] Zhu | Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization[END_REF], depending on the regularity properties of ℓ.

Random Features Based Learning

Random Fourier features (RFF) were initially introduced for shift-invariant scalar kernels [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF] before being extended to shift-invariant OVKs [START_REF] Brault | Random Fourier features for operatorvalued kernels[END_REF] under the operator random Fourier features (ORFF) framework. We refer to Section 2.2 for an overview on the topic. They allow to lower the computational burden associated to solving empirical risk minimization problems by providing finite-dimensional functional spaces whose functions approximate those in the original hypothesis space. In our case, applied to the OVK G, they also provide the additional benefit of bypassing the need to prove a representer theorem for the coefficients (α i ) n i=1 in Equation (3.6), since by construction the RKHSs associated to random feature maps are finite-dimensional.

We choose to work in output with an approximated OVK G obtained by applying the ORFF methodology. We call Φ : Θ → L(U, V) its feature map where V is the associated finite-dimensional feature space. As a reminder, V is of dimension 2ms where m is the chosen number of random features chosen to approximate the kernel k Θ . We refer to Equation (2.25) for the precise construction of this feature map, the key property exploited is that each function α ∈ H G can be written

α = Φ(•) ♯ v (3.29)
for some v ∈ V. This provides a finite-dimensional parametrization of the space H G , which we shall use in combination with Theorem 2.43. Indeed, we know that the solution of Problem 3.5 applied to the vv-RKHS associated to

K = k X Id H G writes as ĥ = n i=1 k X (•, x i )α i (3.
(v i ) n i=1 ∈V n n i=1 I ℓ    n j=1 k X (x i , x j ) Φ(•) ♯ v j , y i    + λ 2 n i,j=1 k X (x i , x j ) v i , v j V . (3.31)
This can be expressed in a more compact form by using tensor products. Indeed, by writing

v := (v i ) n i=1 ∈ V n one gets that n j=1 k X (x i , x j ) Φ(•) ♯ v j = K X ⊗ Φ(•) ♯ v i , ∀i ∈ [n]. (3.32)
We can then recognize in Problem 3.31 a general minimization problem of the form

min v∈V n E θ∼µ   1 n n i=1 ℓ θ, K X ⊗ Φ(θ) ♯ v i , y i (θ) + λ 2 Tr K X vv ⊤   . :=J(θ, v) (3.33)
These so-called stochastic optimization tasks are often encountered in machine learning, and can be solved using gradient methods based on stochastic approximations [START_REF] Bottou | Stochastic gradient learning in neural networks[END_REF].

We propose in Algorithm 3.1 a randomized optimization algorithm based on stochastic gradient descent to obtain v.

Algorithm 3.1 Stochastic Gradient Descent with ORFFs

input : ORFF feature map Φ, input Gram matrix K X , gradient steps

(γ t ) T -1 t=0 init : v (0) = 0 ∈ V n for epoch t from 0 to T -1 do // sampling step Sample θ ∼ µ // gradient step v (t+1) = v (t) -γ t ∂ v J(θ, v (t) ) return v (T )
Remark 3.11. This algorithm describes the vanilla stochastic gradient descent. Depending on the regularities of J, one can turn to more involved optimization algorithm employing e.g. Nesterov acceleration [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF] or any improvement over the vanilla version. The stepsizes (γ t ) T -1 t=0 are also to be chosen depending on the problem to ensure convergence of the algorithm, in our case the objective is strongly convex thanks to the regularizer so that a sufficiently small constant stepsize ensures convergence [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF].

Double ORFF:

If some application stems into large-scale learning, with n being so big that Algorithm 3.1 becomes impossible to employ, one can make use of the RFF trick on both kernels k X and G. This is possible when the input space X is the Euclidean space R d , which we assume in what follows. This amounts to working in the vv-RKHS associated to the OVK K = kX Id H G where G and kX are respectively ORFF and RFF INTEGRAL LOSSES approximations of kernels k X and G. This vv-RKHS is of finite dimension 4m X m Θ s where m X and m Θ are the number of random features used for approximating k X and k Θ . We denote by φX and Φ Θ the corresponding finite-dimensional feature maps, the dependence to each input space being made explicit for clarity. There is then a oneto-one correspondence between function in

H K and M 2m X ,2m Θ s (R): ∀h ∈ H K , ∃β ∈ M 2m X ,2m Θ s (R) such that h(x)(θ) = φX (x) ♯ β Φ Θ (θ), (x, θ) ∈ X × Θ. (3.34)
With this full parametrization of the model, Problem 3.5 writes as

min β∈M 2m X ,2m Θ s (R) E θ∼µ E i∼Unif([n]) ℓ θ, φX (x i ) ♯ β Φ Θ (θ), y i (θ) + λ 2 Tr ββ ⊤ . :=J i (θ, β) (3.
β (0) = 0 ∈ R 2m X ×2m Θ s for t from 0 to T -1 do // sampling step Sample i ∼ Unif([n]) Sample θ ∼ µ // gradient step β (t+1) = β (t) -γ t ∂ β J i (θ, β (t) )
return β (T )

Solving in the Dual

We now investigate dual algorithms for regularized empirical risk minimization in the presence of integral losses. As a reminder, Problem 3.7 writes as inf

(α i ) n i=1 ∈Y n n i=1 I ⋆ ℓ i (-α i ) + 1 2λn n i,j=1 α i , K(x i , x j )α j Y ,
We can see here that the objective function in Problem 3.7 is composed of two terms, the first term being related to the Fenchel-Legendre conjugate of the integral loss I ℓ and to the data, and the second term being a quadratic form involving the dual variables (α i ) n i=1 ∈ Y n and the OVK K.

SOLVING IN THE DUAL

57

The difficulties encountered while solving Problem 3.5 are still visible here: the search for the (α i ) n i=1 ∈ Y n cannot a priori be reduced to a finite-dimensional optimization problem. Also, it is unclear if the first term can be computed numerically, a challenge similar to what was observed in Problem 3.5 with integral losses. Finally, evaluating the quadratic form in the right term involves knowing how to compute each α i , K(x i , x j )α j Y which may not be straightforward.

To tackle these challenges, we first begin in Section 3.4.1 by providing general results on the Fenchel-Legendre conjugate of integral losses on vv-RKHSs. Then in Section 3.4.2 we propose a solution to Problem 3.7 based on expressing dual variables in a basis suitable to both the OVK K and the loss function. This is feasible only when there is some sort of compatibility between the loss function and the output vv-RKHS H G as detailed in Assumption 3. [START_REF] Vapnik | Statistical Learning Theory[END_REF]. Finally in Section 3.4.3 we develop optimization algorithms to be used when there is no such compatibility. The basis in which to express the dual variables is to be picked so as to make possible the use of proximal algorithms that depend on the loss function. The quadratic term in Problem 3.7 is then approximated to make the computations possible.

Fenchel-Legendre Conjugate of Integral Losses

To design dual optimization algorithm, it is crucial that we know how to compute the Fenchel-Legendre conjugate of an integral loss. It turns out that under certain assumptions, the Fenchel-Legendre conjugate of an integral loss is the integral of the Fenchel-Legendre conjugate, as was emphasized in the seminal paper [START_REF] Rockafellar | Conjugate duality and optimization[END_REF]. To that end, two key assumptions must be verified. The first one is that ℓ is a normal convex integrand, which we have assumed throughout this chapter, while the second one is the following: Assumption 3.12. The Hilbert space of functions Y is decomposable, which means that there exists an increasing sequence (Θ k ) k∈N of subsets of Θ with ∪ ∞ k=0 Θ k = Θ such that for ∀k ∈ N, for all bounded measurable functions f : Θ k → R, and for all g ∈ Y, the measurable function that coincides with f on Θ k and with g on Θ \ Θ k belongs to Y.

In a sense, Assumption 3.12 ensures that the space Y is big enough so that the Fenchel-Legendre conjugate of the integral loss is the integral of the Fenchel-Legendre conjugate, as detailed in the following proposition.

where ℓ ⋆ : Θ × U × U →] -∞, +∞] is the integrand defined for ∀(θ, u, v) ∈ Θ × U 2 as ℓ ⋆ (θ, u, v) = ℓ(θ, •, v) ⋆ (u).
It turns out that assuming Assumption 3.12 is strict in the sense that many commonlyused small spaces of functions do not satisfy it: linear models, polynomial functions, Gaussian RKHSs, and so on. However, L 2 [Θ, µ; U] for any probability measure µ verifies Assumption 3.12, which suggests to adopt the view Y = L 2 [Θ, µ; U] and the associated

OVK K = k X T G .
∈ H G sup g∈H G g, α H G -I ℓ (g, y) = I ℓ ⋆ (•, y) • T -1 G (α), (3.37)
where

T -1 G (α) is the only function in L 2 [Θ, µ; U] such that T G (T -1 G (α)) = α.
Proof We can remark that Im(T G ) = H G since supp µ = Θ. Moreover, A is injective, and universality of k Θ implies that T k Θ is injective (see Proposition 2.25), so that

T G = A ⊗ T k Θ is injective. Thus T -1 G (α) ∈ L 2 [Θ, µ; U] is well-defined for all α ∈ H G . Then for ∀g ∈ H G , g, α H G -I ℓ (g, y) = g, T G • T -1 G (α) H G -I ℓ (g, y) = g, T -1 G (α) L 2 [Θ,µ;U] -I ℓ (g, y)
where we used the fact that

T ♯ G is the canonical inclusion from H G to L 2 [Θ, µ; U] (Carmeli et al. (2010), Proposition 3). Since H G is dense in L 2 [Θ, µ; U] and g → g, T -1 G (α) L 2 [Θ,µ;U] -I ℓ (g,

Square Loss: When

ℓ(θ, u, v) = 1 2 u -v 2 U then ℓ ⋆ (θ, u, 0) = ℓ(θ, •, 0) ⋆ (u) = 1 2 u 2 U
and by integrating it along Θ with probability µ we recover what was used in Section 3.2.

The case with nonzero v can be deduced by translation properties of the Fenchel-Legendre conjugate (see Table 2.1).

Pinball loss:

This corresponds to U = R, Θ = [0, 1] and ℓ(θ, u, v) = max (θ(v -u), (θ -1)(v -u)).
In this case, for ∀θ ∈ Θ,

ℓ ⋆ (θ, u, 0) = ℓ(θ, •, 0) ⋆ (u) = χ {θ-1≤•≤θ} (-u). So that for ∀(θ, v) ∈ Θ × U, ℓ ⋆ (θ, u, v) = ℓ(θ, •, 0) ⋆ (u) -uv = χ {θ-1≤•≤θ} (-u) -uv
by using the translation properties of the Fenchel-Legendre conjugate. By integrating it, one gets for all (α, y) ∈ Y 2 :

I ℓ ⋆ (α, y) = χ C (-α) -α, y Y , where C = g ∈ Y : θ -1 ≤ g(θ) ≤ θ holds µ a.e .
Cost-sensitive hinge loss: In this case, U = R, Θ = [0, 1] and

ℓ(θ, u, v) = θ -{-1} (v) max(0, 1 -uv).
As we shall see in Chapter 5, in the application to cost-sensitive classification the data (y i ) n i=1 are constant functions being equal to ±1. We therefore only propose to compute

I ℓ (•, 1) ⋆ and I ℓ (•, -1) ⋆ . For all (θ, u) ∈ Θ × R, it holds that ℓ ⋆ (θ, u, 1) = ℓ(θ, •, 1) ⋆ (u) = χ [0,θ] (-u)
so by integration over Θ one arrives at

I ⋆ ℓ (α, 1) = χ C 1 (-u),
where C 1 = g ∈ Y : 0 ≤ g(θ) ≤ θ holds µ a.e . Similarly, we get

I ⋆ ℓ (α, -1) = χ C -1 (-u)
where

C -1 = g ∈ Y : 0 ≤ g(θ) ≤ 1 -θ holds µ a.e .
We now move on to the solution of Problem 3.7 when the problem benefits from compatibility between the integral operator T G and the loss functions I ℓ ⋆ i .

Integral Operator Eigenbasis Representation

As seen in Proposition 2.42, T G is a compact operator that admits an eigendecomposition

T G = j∈J λ j ψ j ψ ♯ j ,
where (λ j ) j∈J ∈ R |J| + is a non-increasing sequence of eigenvalues with limit 0 and (ψ j ) ∞ j=1 is an orthonormal family of L 2 [Θ, µ; U].

Remark 3.15. The eigenvalues (λ j ) j∈J are not to be confounded with the regularization parameter λ despite the similar notation.

While we have no guarantee that the sequence ( αi ) n i=1 has a finite expansion on this basis, we propose to perform the search of the dual variables (α i ) n i=1 in the subset of L 2 [Θ, µ; U] spanned by the first m eigenvectors of this basis, transforming the Problem 3.7 into a finite-dimensional optimization problem. The rationale behind this approach is that as j → ∞, the influence of a particular direction ψ j in the model will vanish, because

T G ψ j Y = λ j ψ j Y = λ j → j→∞ 0.
This approach can only be applied when there is compatibility between the loss function and the output RKHS H G . Indeed, for Problem 3.7 to be tractable, one needs to be able to compute the terms I ℓ i (α i ) based solely on the coefficients of the dual variables in the eigenbasis (ψ j ) m j=1 . We propose below an assumption that quantifies such behavior.

CHAPTER 3. OPTIMIZATION SCHEMES FOR LEARNING WITH INTEGRAL LOSSES

Assumption 3.16 (Basis Compatibility). The eigenbasis (ψ j ) m j=1 is compatible with the integral loss I ℓ in the sense that for all i ∈ [n], there exists

L i : R m → R such that ∀α ∈ R m , I ℓ ⋆ i   - m j=1 α j ψ j    = L i (α, y i ).
(3.38)

Remark 3.17. One major drawback of this approach is that the eigenbasis must be compatible with the computation of the Fenchel-Legendre conjugate, and Assumption 3.16 is seldom verified in practice. For example, it is verified in case of the squared loss, but not for the pinball or cost-sensitive hinge loss. However, the techniques developed here are useful beyond integral losses, in the context of convoluted losses. This will be exploited in Chapter 4 to gain sparsity using ϵ-insensitive losses or robustness through the Huber loss.

We remark here that there are multiple choices for obtaining the system (ψ j ) m j=1 , as presented in Section 2.2. In particular, the use of approximated eigenbasis Example 2.23 or Random Fourier Features Example 2.30 can be beneficial to the user. Working with a truncated basis is equivalent to solving Problem 3.7 in the vv-RKHS associated to kernel

K(x, z) := k X (x, z) m j=1 λ j ψ j ψ ♯ j .
Problem 3.7 can then be rewritten as

inf α∈Mn,m(R) n i=1 I ℓ ⋆ i   - m j=1 α ij ψ j    L i α i: , y i by (3.38) + 1 2λn
Tr K X αΛα ⊤ , (3.39)

where Λ = diag(λ j ) m j=1 encodes information about the eigenvalues of T G and K X is the input kernel Gram matrix over the data points (x i ) n i=1 . This optimization problem can then be tackled by different means depending on the regularity of the problem. We now explore in Section 3.4.3 optimization algorithms dedicated to solving Problem 3.7 when we cannot make Assumption 3.16.

Proximal Algorithms and Approximated Quadratic Forms

The solution of Problem 3.7 proposed in Section 3.4.2 relies on the strong Assumption 3. [START_REF] Vapnik | Statistical Learning Theory[END_REF]. In practice, this assumption is seldom verified and one cannot compute each I ℓ ⋆ i (-α i ) based on the coefficients of α i in an orthonormal basis. Thus, the approach consisting in choosing first a basis adapted to the quadratic form is flawed here, as the resulting problem is not amenable to optimization. We propose to work the other way around: choose a basis adapted to the I ⋆ ℓ i , and then find a way to compute the quadratic form, at the price of an approximation. Assumption 3.18. Let S be a finite-dimensional subspace of L 2 [Θ, µ; U], so that each function α ∈ S can be parameterized by a vector v ∈ R m . We say that S is compatible with

I ℓ if for all i ∈ [n], there exist R i : R m → R and Q i : R m → -∞, +∞ so that 1. for all α ∈ S parameterized by v ∈ R m , I ⋆ ℓ i (α) = R i (v) + Q i (v), 3.4. SOLVING IN THE DUAL 61
2. R i is differentiable, and 3. prox γQ i is computable in closed form for all γ > 0.

Assumption 3.18 is designed so that Problem 3.7 writes as a composite problem whose solution can be performed using proximal gradient descent algorithm (see Section 2.1 for reminders on these concepts). As an example, we consider the specific quantile regression scheme associated to the integral pinball loss. In this setting, Θ = [0, 1] and U = R are one-dimensional. The Fenchel-Legendre conjugate of the integral pinball loss computed in Section 3.4.1 reads

I ℓ ⋆ i (α) = -α, y i Y + χ C (-α), (α, i) ∈ Y × [n] with C = g ∈ Y : θ -1 ≤ g(θ) ≤ θ holds µ a.e . Moreover,
prox γχ C (•) (α) = Proj C (α)
for all γ > 0 and α ∈ Y. Thus, to satisfy Assumption 3.18 a suitable finite-dimensional representation space S must satisfy the two properties

1. for all (α, i) ∈ S × [n], if α is encoded by a vector v ∈ R m then R i (v) = -α, y i Y
is computable (or at least can be approximated)

2. for all α ∈ S, Proj C (α) can be computed in closed form.

Taking a closer look at this projection operator, we notice that

∀α ∈ Y, Proj C (α) : Θ → R θ → max(θ -1, min(α(θ), θ)).
(3.40)

We see here why the eigendecomposition of T G is not a suitable basis: Proj C (α) involve a pointwise projection of each α(θ) on the corresponding interval [θ -1, θ] so that Span{(ψ j ) m j=1 } may not stable with respect to this projection. This suggests to choose a representation for the dual variables that allows a pointwise control on the functions. We propose to use linear splines in this case, as their shape matches that of the boundaries of C.

Linear Splines A linear spline α ∈ Y is encoded by a set of ordered locations (θ j ) m j=1 ∈ Θ m and by a vector α := (α(θ j )) m j=1 ∈ R m . We denote by S m the set of linear splines obtained with those fixed locations. In between locations can be computed using linear interpolation:

∀θ ∈ [θ j , θ j+1 ], α(θ) = α(θ j ) + α(θ j+1 ) -α(θ j ) θ j+1 -θ j (θ -θ j ).
The first remark to be made about restricting dual variables to be linear splines is that linear splines can approximate any function in Y, so that it is a rich enough representation to work with. Moreover, although restricting the dual variables to be linear splines may appear a limitation, it turns out that in the estimator

ĥ = 1 λn n i=1 k X (•, x i )T G αi 62 CHAPTER 3. OPTIMIZATION SCHEMES FOR LEARNING WITH INTEGRAL LOSSES
the splines are smoothed by means of the integral operator T G , so that the model itself is not piecewise linear in θ, only the dual variables are. The drawback of this method is that one cannot compute exactly the quadratic terms α i , T G α j Y . We propose to approximate them as presented below.

Approximated Quadratic Forms: Let α, β be linear splines encoded by locations (θ j ) m j=1 ∈ Θ m distributed i.i.d. according to µ and the vectors of function evaluations α, β ∈ R m . Let K Θ be the Gram matrix associated to locations (θ j ) m j=1 and kernel k Θ . Then

T k Θ β (θ) = Θ k Θ (θ, θ ′ )β(θ)dµ(θ ′ ) ≈ 1 m m j=1 k Θ (θ, θ j )β(θ j ), θ ∈ Θ, and
α, T k Θ β Y = α, Θ k Θ (•, θ ′ )β(θ)dµ(θ ′ ) Y ≈ 1 m 2 m i,j=1 k Θ (θ i , θ j )α(θ i )β(θ j ) α ⊤ K Θ β
We thus propose to solve the following problem inf

(α i ) n i=1 ∈Sm n i=1 χ C (α i ) + n i=1 α i , y i Y + 1 2λnm 2 Tr K X αK Θ α ⊤ := J(α) , ( 3.41)
where given some splines (α i ) n i=1 ∈ S m , we denote by α ∈ M n,m (R) the matrix encoding row-wise the values associated to each spline. We present in Algorithm 3.3 a proximal gradient algorithm for solving Problem 3.41, that reads as a projected gradient descent: for each epoch we perform one step of gradient descent on J, followed by the proximal step that consists in the projection on the feasible set defined by C. Algorithm 3.3 Proximal Gradient Descent with Linear Splines for Quantile Regression

input : Splines locations (θ j ) m j=1 i.i.d. ∼ µ, Gram matrices K X , K Θ init : α = 0 ∈ M n,m (R) for epoch from 1 to T do // gradient step α = α -∇J(α)
// proximal step for row i from 1 to n do for column j from 1 to m do

α ij = Proj [θ j -1,θ j ] (α ij ) return α
Remark 3.19. While in the general case the α i , y i Y have to be estimated, as we will see in Chapter 5 the Infinite Task Learning framework uses constant functions as observed outputs (y i ) n i=1 . This makes the exact computation of these scalar products 3.5. CONCLUSION 63 possible when µ is the Lebesgue measure, as for a spline α ∈ S m α, y i Y = y i (0)

Θ α(θ)dθ = y i (0) m-1 j=1 θ j+1 θ j α(θ)dθ = y i (0) m-1 j=1 1 2 (α(θ j+1 ) + α(θ j ))(θ j+1 -θ j ).
This algorithm can be adapted to the integral cost-sensitive hinge loss with a slightly different projector suited to Section 3.4.1 and Section 3.4.1. The general idea of the algorithm can work in higher dimensions: provided that some set S allows a pointwise control on the dual variables and Assumption 3.18 is verified, the quadratic term can be approximated and a proximal gradient algorithm is possible.

Conclusion

In this chapter, we provide optimization algorithms to solve regularized empirical risk minimization problems in vv-RKHSs in the presence of integral losses. Beyond closedform solutions for the square loss, we address the question of representing the functional coefficients by several means. Primal approaches include the double representer theorem, which provides a workable parametrization of the solution to the price of an approximation on the integral loss, and random Fourier features amenable to optimization via stochastic algorithms. Dual techniques are also developed, and involve some compatibility conditions between the vv-RKHS and the loss for the representation of the dual variables to be viable. In particular, dual approaches pave the way to extend this framework beyond integral losses, by using the infimal convolution operator to design new loss functions. Finally, we want to emphasize the compatiblity of the developed algorithms with hybrid models involving the composition of a neural architecture and a kernel, allowing to tackle tasks with complex input data. This is exemplified in Section 5.4.4 for quantile regression with input being images. In this chapter, we propose a framework dedicated to regressing functional outputs with a focus on robustness. To do so, we extend the integral loss framework introduced in Chapter 3 to handle convoluted losses. This goal is achieved by exploiting duality principles which are well-suited to the nature of these loss functions. After a brief introduction in Section 4.1, we state the problem to solve in Section 4.2. In Section 4.3 we propose a methodology that handles different variants of the Huber loss, known to induce robustness to outliers. In Section 4.4 we study optimization problems based on ϵ-insensitive losses which are akin to bring sparsity to the estimator. Numerical experiments on real and synthetic datasets are gathered in Section 4.5, and conclusions are drawn in Section 4.6.

Introduction

Due to the increased availability of streaming data, learning to predict complex objects has attracted a great deal of attention in machine learning. Biomedical Signal Processing, Epidemiology Monitoring or Climate Science are examples of interdisciplinary research fields where the phenomena under study exhibit a functional nature, and the understanding of these phenomena depends on machine learning algorithms being able to reliably handle functional data. Functional data analysis (FDA, Ramsay and Silverman 1997; [START_REF] Wang | Functional data analysis[END_REF]) has been devoted to such setting, where one assumes that the measurements of the underlying phenomena come numerous enough so that modeling it with functions makes sense. In particular, Functional Output Regression (FOR) specializes to regression problems where the output variable is a function, the input variable being either finite-dimensional or of a functional nature themselves.

The simplest way to design an algorithm is then to model a linear dependency between the inputs and the outputs [START_REF] Morris | Functional regression[END_REF], at the price of being unable to cope with complex dependency within the data. To remedy this bottleneck, nonlinear approaches have been developed in recent years. In nonparametric statistics, [START_REF] Ferraty | Kernel regression with functional response[END_REF] propose a Banach-valued version of the Nadaraya-Watson estimator. Kernel methods have proven useful to tackle this problem as well, with works involving tri-variate regression problem in RKHSs [START_REF] Reimherr | Optimal prediction for additive function-on-function regression[END_REF], approximated kernel ridge regression (KRR) using orthonormal bases [START_REF] Oliva | Fast function to function regression[END_REF], and in the operator-valued kernel literature a function-valued KRR with double representer theorem [START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel Hilbert spaces[END_REF], solvers based on a sampling of the functional norm [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF] or purely functional methods relying on approximate inversion of integral operators [START_REF] Kadri | Operatorvalued kernels for learning from functional response data[END_REF]. We can also mention recent techniques using vv-RKHSs to learn finite-dimensional coefficients expressing the functional outputs with the help of a dictionary basis [START_REF] Bouche | Nonlinear functional output regression: A dictionary approach[END_REF].

Most of these applications focus on the square loss, which is known to induce estimates of the conditional expectation of the functional outputs given the input data. However, defective sensors or malicious attacks can lead to erroneous or contaminated measurements of the phenomena that result in functional outliers, and using the square loss has the major drawback of providing estimators which are sensitive to these outliers, producing unreliable prediction systems. In the scalar-valued case, variations of the square loss such as the Huber loss [START_REF] Huber | Robust estimation of a location parameter[END_REF] or ϵ-insensitive losses [START_REF] Lee | epsilon-SSVR: A smooth support vector machine for epsilon-insensitive regression[END_REF] have been introduced (among other techniques) as a way to mitigate such sensitivity. In the FDA setting, robustness aspects have been investigated using Bayesian methods [START_REF] Zhu | Robust, adaptive functional regression in functional mixed model framework[END_REF], trading the mean for the Banach-valued median [START_REF] Cadre | Convergent estimators for the l1-median of Banach valued random variable[END_REF], using bounded loss functions [START_REF] Maronna | Robust functional linear regression based on splines[END_REF], or leveraging principal component analysis [START_REF] Kalogridis | Robust functional regression based on principal components[END_REF]. We can also mention works involving the Huber loss function for the semi-functional setting (functional input variable, scalar output variable; Crambes In the operator-valued kernel literature, extension of ϵ-insensitive losses to the (finite) vector-valued regression setting has been proposed by [START_REF] Sangnier | Data sparse nonparametric regression with ϵ-insensitive losses[END_REF]. Using convex optimization tools such as the infimal convolution operator and parametric duality leads to efficient solvers and provide sparse estimators, an idea later exploited in [START_REF] Laforgue | Duality in RKHSs with infinite dimensional outputs: Application to robust losses[END_REF] where a generalization of this approach to infinite-dimensional outputs encompassing both the Huber and ϵ-insensitive losses is proposed. In this chapter, we extend the results from Laforgue et al. (2020) in the FOR setting and build a robust framework based on parametric duality in vv-RKHSs with functional outputs, adapting the dual optimization algorithms from Chapter 3 to the setting of convoluted losses.

Problem Setting

In the functional output regression (FOR) setting, the goal is to regress to a functional output Y taking its values in Y := L 2 [Θ, µ] where Θ := [0, 1] is endowed with a probability measure µ, from an input variable X that either takes its values in the Euclidean space X = R d or in the Hilbert space X = L 2 [Θ, µ] itself. A natural way to tackle this problem is to estimate E[Y|X] = h † (X) where it is well-known that

h † = arg min h∈H E (X,Y) Y -h(X) 2 Y (4.1)
with H being the set of all measurable functions from X to Y. The existence of h † is guaranteed by Theorem 10.2.2 from [START_REF] Dudley | Real Analysis and Probability[END_REF] as Y is a Polish space (i.e. a separable completely metrizable topological space).

In practice, the law of (X, Y) is unknown, and one has only information about it via i.i.d. samples (x i , y i ) n i=1 . Since working with the set of all measurable functions from X to Y is notoriously hard, we restrict the hypothesis space to be a vector-valued reproducing kernel Hilbert space (vv-RKHS) H K ⊂ F(X, Y) associated to some operator-valued kernel (OVK) K : X × X → L(Y). We assume that K is chosen to be a decomposable kernel of the form K = k X T k Θ with scalar-valued kernels k X : X × X → R and k Θ : Θ × Θ → R. We refer to Section 2.2.2 for more details about such functional spaces. We then arrive to the following regularized empirical risk minimization problem

inf h∈H K 1 n n i=1 y i -h(x i ) 2 Y + λ 2 h 2 H K , ( 4.2)
where λ > 0 is a hyperparameter controling the strength of the regularization in the vv-RKHS. Problem 4.2 has been introduced and solved in Kadri et al. ( 2016) with a closed-form estimator ĥ :

= n i=1 K(•, x i )α i , (4.3)
with optimal coefficients α := [ αi ] n i=1 ∈ Y n given by

α = K X ⊗ T k Θ + λn Id -1 y, (4.4)
where K X is the Gram matrix associated to the input data (x i ) n i=1 and kernel k X , and y = [y i ] n i=1 is the vector aggregating the functional responses (y i ) n i=1 . One main advantage of working with the square loss is to benefit from a closed-form solution. This is due to the fact that the gradient of the objective function linearly depends on the residuals y i -h(x i ), reducing Problem 4.2 to the solution of a linear system. This somehow innocent and positive property has dramatic consequences in the presence of contaminated data, as the estimator will be strongly influenced by a large shift in the data caused by some measurement error or malicious actor. One can then resort to using losses different than the square loss, that provide more robust estimators.

Remark 4.1. The notation here slightly differs from Chapter 3 where the loss functions were taking two arguments. The reason behind this was to encompass two kinds of residuals under the same umbrella: the regression case y -h(x) and the classification case yh(x). Since this chapter is devoted to regression, the loss function L in Problem 4.5 takes only one argument which is y -h(x).

Given a general loss function L : Y → R, we consider the problem

inf h∈H K 1 n n i=1 L(y i -h(x i )) + λ 2 h 2 H K (4.
(α i ) n i=1 ∈Y n n i=1 L ⋆ i (-α i ) + 1 2λn n i=1 n j=1 k X (x i , x j ) α i , T k Θ α j Y , (4.6)
where f ⋆ is the Fenchel-Legendre conjugate of any function f : Y → -∞, +∞ (see Definition 2.2) and by convention

L i = L(y i -•) ∀i ∈ [n].
The solution of Problem 4.5 then takes the form

ĥ = 1 λn n i=1 k X (•, x i )T k Θ αi
with (α i ) n i=1 being the solution of Problem 4.6. Among interesting losses, we can mention the Huber loss [START_REF] Huber | Robust estimation of a location parameter[END_REF] which is known to induce robust estimators, or ϵ-insensitive losses [START_REF] Lee | epsilon-SSVR: A smooth support vector machine for epsilon-insensitive regression[END_REF] leading to sparse solutions. A particularity of these losses is that they write as the infimal convolution of two loss functions, which can be used to our advantage in Problem 4.6. Indeed, constructing suitable loss functions L can be carried out using the infimal convolution operator (see Definition 2.3): given two loss functions f and g, taking L = f □ g yields a new loss function that can be used in Problem 4.5. In the duality framework, this type of loss functions is particularly appealing as the Fenchel-Legendre conjugate of L can be expressed as a simple addition L ⋆ = f ⋆ + g ⋆ in Problem 4.6. This way, to design a loss function of interest, we can for example take f = 1 2 || • || 2 Y and convolve it with a function g that brings an interesting property (sparsity, robustness) to the estimator resulting from solving Problem 4.6. This results in dual problems with objective function being the a smooth quadratic part, to which is added a non smooth term corresponding to g ⋆ . This idea was first developed in the operator-valued kernel literature in [START_REF] Sangnier | Data sparse nonparametric regression with ϵ-insensitive losses[END_REF], where it was applied to ϵ-insensitive losses in vv-RKHSs with finite-dimensional outputs and then extended to infinite-dimensional outputs in [START_REF] Laforgue | Duality in RKHSs with infinite dimensional outputs: Application to robust losses[END_REF].

We propose to deepen this approach in the FOR case, and introduce two families of problems that each bring special properties to the estimator. In Section 4.3, we develop optimization algorithms to enforce robust estimators by using different versions of the Huber loss. Then in Section 4.4, we apply this methodology to ϵ-insensitive loss functions, for which estimators benefit from sparsity properties. A summary of all these settings is given in Table 4.1. Finally, numerical experiments are presented in Section 4.5.2 to illustrate the effectiveness of the approach in different outlier scenarii.

ROBUST ESTIMATORS WITH HUBER LOSS

69

Robust Estimators with Huber Loss

The Huber loss has been introduced by Huber (1964) to provide a robust estimate of the mean of a random variable. In the real output case, it is defined as follows.

Definition 4.2. The real Huber loss with parameter κ > 0 is defined as

∀s ∈ R, ℓ H,κ (s) :=      1 2 s 2 if |s| ≤ κ κ |s| -κ 2 otherwise.
Due to its asymptotic behavior as κ| • |, the Huber loss is useful when the training data is heavy-tailed or contains outliers. It turns out that ℓ H,κ can be written in the more compact form

ℓ H,κ = 1 2 (•) 2 □ κ| • |. (4.7)
This allows to easily compute its Fenchel-Legendre transform as detailed in the following proposition.

Proposition 4.3. The Fenchel-Legendre transform of the real Huber loss ℓ H,κ is given by

∀s ∈ R, ℓ ⋆ H,κ (s) = 1 2 s 2 + χ [-κ,κ] (s). (4.8)
Proof This is a simple application of Proposition 2.4 to f = 1 2 (•) 2 and g = κ| • |, whose Fenchel-Legendre transforms can be found in Table 2

L H,κ := κ • Y □ 1 2 • 2 Y , (4.9)
or again:

∀y ∈ Y, L H,κ (y) =      1 2 y 2 Y if y Y ≤ κ κ y Y -κ 2
otherwise.

An illustration of the vectorial Huber loss is provided in Figure 4.1 for dimensions 1 and 2, in dimension 1 it coincides with the real Huber loss. Another way to extend the Huber loss to Y is to consider the integral loss associated to the real Huber loss (see Section 3.1 for more details about integral losses), as defined below.

Vectorial Huber Loss

Before proposing a solution to Problem 4.5 with the vectorial Huber loss, we quickly provide the expression of its Fenchel-Legendre transform.

Proposition 4.6. The Fenchel-Legendre transform of the vectorial Huber loss L H,κ is given by

∀y ∈ Y, L ⋆ H,κ (y) = 1 2 y 2 Y + χ Bκ (y) (4.11)
where B κ is the ball with radius κ according to norm • Y .

Proof This is again an application of Proposition 2.4 to f = 1 2 • 2 Y and g = κ • Y , whose Fenchel-Legendre transforms can be found in Table 2.1.

We are now ready to exemplify Problem 4.6 for the vectorial Huber loss scenario.

Proposition 4.7. The dual of Problem 4.5 in the vectorial Huber loss case writes as inf

(α i ) n i=1 ∈Y n n i=1 1 2 α i 2 Y -α i , y i Y + χ Bκ (α i ) + 1 2λn n i=1 n j=1 k X (x i , x j ) α i , T k Θ α j Y . (4.12)
Proof This is the direct instantiation of Problem 4.6 with Fenchel-Legendre transform given by Proposition 4.6, combined with the properties of the Fenchel-Legendre transform with respect to the translation operator (see Table 2.1).

We first notice that Problem 4.12 is similar to the dual problems tackled using integral losses in Section 3.4. Its solution can be performed using the techniques developed in Section 3.4.2. Indeed, Assumption 3.16 which assess of the compatibility between the loss function and the integral operator T k Θ is satisfied here: all the terms involved in Problem 4.12 can be computed provided that we possess a representation of the (α i ) n i=1 in an orthonormal basis. We thus represent the dual variables in the truncated eigenbasis (λ j , ψ j ) m j=1 associated to

T k Θ . We introduce the notation R ∈ M n,m (R) such that R ij = y i , ψ j Y and Λ = diag (λ j ) m j=1 ∈ M m (R).
Denoting by α ∈ M n,m (R) the matrix encoding the coefficients of the dual variables, Problem 4.12 can be rephrased as

inf α∈Mn,m(R) Tr 1 2 αα ⊤ -αR ⊤ + 1 2λn K X αΛα ⊤ := J(α) s.t. α 2,∞ ≤ κ. (4.13)
Algorithm 4.1 Projected Gradient Descent for Vectorial Huber input : Gram matrix K X , matrix of eigenvalues Λ, data scalar product matrix R, regularization parameter λ, Huber parameter κ, gradient step γ init :

α (0) = 0 ∈ R n×m or α (0) = λn K X ⊗ Λ + λn Id nm -1 R for epoch t from 0 to T -1 do // gradient step α (t+1) = α (t) -γ α (t) + 1 λn K X α (t) Λ -R // projection step for row i from 1 to n do α (t+1) i: = min    κ α (t+1) i: 2 , 1    α (t+1) i:
return α (T ) We recognize in Problem 4.13 a composite optimization problem (see Problem 2.5) that can be tackled using proximal gradient descent. Let γ > 0 and β ∈ R m , it holds that

prox γχ Bκ (•) (β) = Proj Bκ (β) = min κ/ β 2 , 1 β,
hence by the separability of α → n i=1 χ Bκ (α i: ) the proximal step is obtained by projecting each row of α onto B κ . The proximal gradient descent algorithm is then akin to a projected gradient descent because of the particular form of the proximal operator. The algorithm is summarized in Algorithm 4.1. The resulting estimator is then given by ĥ = 1 λn

n i=1 m j=1 λ j αij k X (•, x i )ψ j .
Remark 4.8. When κ is large, one recovers the unconstrained ridge regression problem from [START_REF] Kadri | Operatorvalued kernels for learning from functional response data[END_REF], whose practical solution is performed in the truncated eigenbasis (ψ j ) m j=1 (see Section 3.

for details). The solution is described by the closed-form expression

α = λn K X ⊗ Λ + λn Id nm -1 R. (4.14)

This can be used as an initialization point in Algorithm 4.1 if one has enough computational power. It does not correspond to a feasible dual point (except for large enough κ)

but after one epoch the projection step will ensure that the iterates satisfies α 2,∞ ≤ κ.

Remark 4.9. The gradient step γ to choose can be estimated from the parameters of the problem. Indeed, for guaranteed convergence, one must set γ < 2 C where C is the Lipschitz constant associated to the gradient of the objective function J. Here,

∇J(α) = α + 1 λn K X αΛ -R which is C-Lipschitz with C = 1 + 1 λn K X op λ 1 .

Integral Huber Loss

The solution of Problem 4.5 with the integral Huber loss follows similar steps, exploiting the corresponding Fenchel-Legendre transform whose expression is given below.

Proposition 4.10. The Fenchel-Legendre transform of the integral Huber loss I ℓ H,κ is given by

∀y ∈ Y, I ⋆ ℓ H,κ (y) = 1 2 y 2 Y + χ B ∞ κ (y), (4.15)
where

B ∞ κ = y ∈ Y : |y(θ)| ≤ κ , µ
(α i ) n i=1 ∈Y n n i=1 1 2 α i 2 Y -α i , y i Y + χ B ∞ κ (α i ) + 1 2λn n i=1 n j=1 k X (x i , x j ) α i , T k Θ α j Y . (4.16)
Proof This is the direct instantiation of Problem 4.6 with Fenchel-Legendre transform given by Proposition 4.10, combined with the properties of the Fenchel-Legendre transform with respect to the translation operator (see Table 2.1).

In contrast to what was obtained in Section 4.3.1, Assumption 3.16 is no longer satisfied, and we cannot perform the solution of Problem 4.16 in the eigenbasis associated to T k Θ . Indeed, given α i = m j=1 α ij ψ j it is not possible to evaluate χ B ∞ κ (α i ) based on the sole coefficients (α ij ) m j=1 , nor is it possible to easily project α i onto B ∞ κ . We can then turn to the method developed in Section 3.4.3 and represent the dual variables as linear splines, at the cost of approximating the quadratic term involving T k Θ . Indeed, a proximal gradient algorithm involves a projection on the feasible set B ∞ κ which has linear borders: it thus seems appropriate to use a linear spline basis for the representation of the dual variables (α i ) n i=1 . Given a set of locations (θ j ) m j=1 i.i.d.

∼ µ, we denote by S m the set of linear splines with locations (θ j ) m j=1 , and encode the dual variables (α i ) n i=1 by a matrix α ∈ M n,m (R) such that the rows of α correspond to the values of the dual variables at locations (θ j ) n j=1 . Finally, we approximate α i , y i Y ≈ 1 m m j=1 α i (θ j )y i (θ j ) and store the observed outputs y i (θ j ) in a matrix Y ∈ M n,m (R). Problem 4.16 then writes as

inf α∈Mn,m(R) Tr 1 2 αα ⊤ -αY ⊤ + 1 2λnm K X αK Θ α ⊤ s.t. α ∞ ≤ κ.
(4.17)

A proximal gradient descent algorithm to solve Problem 4.17 is presented in Algorithm 4.2, where projecting on B ∞ κ is equivalent to projecting each value α ij onto the interval [-κ, κ]. The resulting estimator is given by ∀x ∈ X, ĥ(

x) = 1 λnm n i=1 m j=1 αij k X (x, x i )k Θ (•, θ j ).
α (0) = 0 ∈ R n×m or α (0) = λnm K X ⊗ K Θ + λnm Id nm -1 Y for epoch t from 0 to T -1 do // gradient step α (t+1) = α (t) -γ α (t) + 1 λnm K X α (t) K Θ -Y // projection step for row i from 1 to n do for column j from 1 to m do 28 α (t+1) ij = sign α (t+1) ij min κ, α (t+1) ij return α (T )

Sparse Estimators with ϵ-Insensitive Losses

Let us recall the important notion of ϵ-insensitive losses. Learning with ϵ-insensitive losses forces the estimator to neglect small errors, preventing from overfitting and inducing some form of regularization, which will be highlighted later in the associated dual problems.

∀y ∈ Y, L ϵ (y) =      0 if y Y ≤ ϵ inf ∥d∥ Y ≤1
L(y -ϵd) otherwise.

In other terms, L ϵ (y) is the smallest value of L within the ball of radius ϵ centered at y.

Remark 4.13. The definition of L ϵ depends on the metric used to characterize B ϵ . We use the • Y norm as a natural choice, but one could envision variations of the ϵ-insensitive losses based on different metrics, leading to an alternative diverse family of dual problems.

In general, it is not possible to find an analytic expression for L ϵ . However, the ϵinsensitive version of the square loss enjoys a closed-form representation as detailed below.

Definition 4.14. We define the vectorial ϵ-ridge loss to be the ϵ-insensitive version of the square loss:

• 2 Y,ϵ := • 2 Y □ χ Bϵ (•) = max( • Y -ϵ, 0) 2 . (4.18)
The vectorial ϵ-ridge loss is illustrated in Figure 4.2 for dimension 1 and 2. Similarly to what was done with the Huber loss in Section 4.3, we can define an interesting loss on Y by integrating the local real ϵ-ridge losses. Having defined these two losses, we are ready to dive into the associated Problem 4.5. We start in Section 4.4.1 by proposing a solution method for the vectorial ϵ-insensitive ridge loss, before considering in Section 4.4.2 the integral ϵ-insensitive ridge loss function.

Vectorial ϵ-Ridge

The solution of Problem 4.6 with the vectorial ϵ-ridge is similar to its Huber counterpart presented in Section 4.3.1. We begin by giving the Fenchel-Legendre conjugate of • 2 Y,ϵ .

Proposition 4.16. The Fenchel-Legendre conjugate of the vectorial ϵ-insensitive ridge loss is given by

∀y ∈ Y, • 2 Y,ϵ ⋆ (y) = 1 2 y 2 Y + ϵ y Y . (4.20)
Proof This is again an application of Proposition 2.4 applied to f = 1 2 • 2 Y and g = χ Bϵ (•), whose Fenchel-Legendre conjugates can be found in Table 2.1.

We can now apply this to Problem 4.6 in the vectorial ϵ-insensitive ridge loss scenario as stated below.

Proposition 4.17. The dual to Problem 4.5 in the vectorial ϵ-insensitive ridge loss case writes as inf

(α i ) n i=1 ∈Y n n i=1 1 2 α i 2 Y -α i , y i Y +ϵ α i Y + 1 2λn n i=1 n j=1 k X (x i , x j ) α i , T k Θ α j Y . (4.21)
Proof This is the direct instantiation of Problem 4.6 with Fenchel-Legendre conjugate given by Proposition 4.16, combined with the properties of the Fenchel-Legendre conjugate with respect to the translation operator (see Table 2.1). Algorithm 4.3 Proximal Gradient Descent for Vectorial ϵ-Ridge input : Gram matrix K X , matrix of eigenvalues Λ, data scalar product matrix R, regularization parameter λ, ridge parameter ϵ, gradient step γ init :

α (0) = 0 ∈ R n×m or α (0) = λn K X ⊗ Λ + λn Id nm -1 R for epoch t from 0 to T -1 do // gradient step α (t+1) = α (t) -γ α (t) + 1 λn K X α (t) Λ -R // proximal step for row i from 1 to n do α (t+1) i: = BST α (t+1) i:
, γϵ return α (T ) We can notice that Assumption 3.16 is again satisfied here, and exploit the same methodology from Section 3.4.2 as in Section 4.3.1. Keeping the same notations, one can rephrase Problem 4.21 as

inf α∈Mn,m(R) Tr 1 2 αα ⊤ -αR ⊤ + 1 2λn K X αΛα ⊤ := J(α) +ϵ α 2,1 . (4.22)
Problem 4.22 consists in the minimization of the finite-dimensional quadratic form J with multi-task Lasso regularization ϵ • 2,1 . We can adapt the proximal gradient descent introduced in Algorithm 4.1 by choosing the suitable proximal operator which is the block soft thresholding (BST) operator, defined for all x ∈ R m by

BST(x, ϵ) := 1 - ϵ x 2 + x.
The algorithm is summarized in Algorithm 4.3. Remark 4.9 about the calibration of the gradient step γ is still valid here.

Integral ϵ-Ridge

The solution of Problem 4.5 with the integral ϵ-ridge loss follows closely the steps in Section 4.3.2. We begin by computing the associated Fenchel-Legendre conjugate.

Proposition 4.18. The Fenchel-Legendre conjugate of the integral ϵ-ridge loss is given by

∀y ∈ Y, I 1 2 (•) 2 ϵ ⋆ (y) = 1 2 y 2 Y + ϵ Θ |y(θ)|dµ(θ). (4.23)
Proof This is an application of Proposition 3.13 which states that

I ⋆ 1 2 (•) 2 ϵ = I 1 2 (•) 2 ϵ ⋆ .
We can now state Problem 4.6 for the integral ϵ-ridge loss scenario.

Proposition 4.19. The dual to Problem 4.5 in the integral ϵ-ridge loss case can be written as inf

(α i ) n i=1 ∈Y n n i=1 1 2 α i 2 Y -α i , y i Y +ϵ Θ |α i (θ)|dµ(θ)+ 1 2λn n i=1 n j=1 k X (x i , x j ) α i , T k Θ α j Y . (4.24)
Proof This is the direct instantiation of Problem 4.6 with Fenchel-Legendre conjugate given by Proposition 4.18, combined with the properties of the Fenchel-Legendre conjugate with respect to the translation operator (see Table 2.1).

In contrast to what was obtained in Section 4.4.1, Assumption 3.16 is no longer satisfied. This is due to the term ϵ Θ |α i (θ)|dµ(θ) in Problem 4.24 which cannot be computed from a representation of the dual variables in an orthonormal basis. We thus propose to rely on the splines methodology developed in Section 3.4.3, and begin by expressing the proximal operator further used in the proximal descent algorithm: as we recognize that ϵ

Θ |α i (θ)|dµ(θ) = ϵ α i 1
the proximal operator of the non-differentiable part can easily be computed using Moreau decomposition Proposition 2.6

∀γ > 0, prox γϵ∥•∥ 1 = Id -γ prox 1/γ ϵ • 1 ⋆ χ B ∞ ϵ • γ = Id -γ Proj B ∞ γϵ (•).
Thus linear splines are adapted to the problem in the sense of Assumption 3.18, as they provide pointwise control of the dual variables, suited to the projection on B ∞ γϵ . Using the same notations as in Section 4.3.2, an approximated Problem 4.24 writes as

inf α∈Mn,m(R) Tr 1 2 αα ⊤ -αY ⊤ + 1 2λnm K X αK Θ α ⊤ + ϵ α 1 . (4.25)
We recognize in Problem 4.25 the minimization of a quadratic form with Lasso penalty ϵ α 1 . We propose to solve it using a proximal gradient algorithm presented in Algorithm 4.4. The proximal step involves the soft-thresholding operator

ST(x, ϵ) := sign(x) |x| -ϵ + ,
which is known to induce sparsity in the iterates. The resulting estimator is given by ∀x ∈ X, ĥ

(x) = 1 λnm n i=1 m j=1 αij k X (x, x i )k Θ (•, θ j ).

Numerical Experiments

In this section, we present numerical experiments to illustrate the benefits of learning with convoluted losses. We consider the setting of function-to-function regression, where both the input variable X and the output variable Y are functions. The experiments relied on two datasets, one real and one synthetic.

α (0) = 0 ∈ R n×m or α (0) = λnm K X ⊗ K Θ + λnm Id nm -1 Y for epoch t from 0 to T -1 do // gradient step α (t+1) = α (t) -γ α (t) + 1 λnm K X α (t) K Θ -Y //

Real dataset

We choose to work with the Lip dataset from [START_REF] Ramsay | Functional data analysis[END_REF], which depicts the movement of the lip when a subject pronounces the word "bob" inside the phrase "say bob again". More precisely, the input observations (x i ) n i=1 are Electromyogram (EMG) recordings of a muscle that depresses the lower lip while the output observations (y i ) n i=1 describe the acceleration of the lower lip of the speaker, both seen as functions of the time. These functions are observed each millisecond, which forms a regular grid of size 641 considering the duration of the experiments; we rescaled this grid so that the domain of these functions is Θ := [0, 1]. The underlying goal is to map the EMG signals onto the movements of the lip, allowing for a better understanding of how the brain controls the diction. We have access to n = 32 such recordings, whose illustration can be found in Figure 4.3.

Synthetic dataset

To design a challenging task for our functional output regression framework, we adapted a synthetic dataset from [START_REF] Bouche | Nonlinear functional output regression: A dictionary approach[END_REF]. A pair of (input, output) functions is built as follows. We first draw r = 4 Gaussian processes (GPs) denoted by (GP in j ) r j=1 with mean 0 and Gaussian covariances of various bandwidths which will be used to model the input functions, and r other GPs to model the output functions, this time denoted by (GP out j ) r j=1 . Both sets of GPs were kept fixed, and to create each pair of samples (x i , y i ), we draw (a ij ) j∈[r] ∈ R r uniformly in [-2, 2] r , and

a ij GP in j (θ), y i (θ) = r j=1 a ij GP out j (θ).
Examples of such functions can be found in Figure 4.4.

Outliers generation

To illustrate the robustness properties brought by the use of the Huber and ϵ-insensitive losses, we propose to add two kinds of outliers to the aforementioned datasets. The first type of outliers, that we refer to as local outliers, model inconsistency in the data acquisition pipeline on a local level: we take a sample from the training set and add some noise to the measurements of the functions at m o number of locations ( θj ) mo j=1 chosen randomly. The second type of outliers consists in global outliers, that is functions whose shape do not match the ones of the training samples. In particular, we create them the following way: for the lip dataset, we consider the case of a data acquisition error where an input function x i is observed correctly, but associated to an output function ỹi = -s o y i where s o is the scale of the outlier, fixed during the experiment and y i is the normal observation. This is chosen to maximize the disalignment between inliers and outliers. For the synthetic dataset, we create new samples with different Gaussian processes in input and output, see Figure 4.5 for an illustration. We also use s o to refer to the scale of the outliers, in the sense that the coefficients

(a ij ) j∈[r] ∈ R r are drawn uniformly in [-s o , s o ] r .
Model Specification We use a decomposable kernel with respective input and output kernels being k X and k Θ . The former is chosen to be an integrated Gaussian kernel, i.e. k X (x, z) = Θ e -γ X [x(θ)-z(θ)] 2 . The bandwidth γ X is chosen via cross validation. For k Θ , we chose a Gaussian kernel with bandwidth γ Θ also chosen by cross validation. Since algorithms dealing with the vectorial Huber and ϵ-insensitive losses involve the eigendecomposition of T k Θ , in these cases we use a random Fourier features (RFF) approximation of the kernel for which the eigendecomposition of the associated integral operator is easy to get (see Example 2.30). The number of RFFs is set to 25 which provides functional spaces of dimension 50.

Huber Losses

In this section, we demonstrate the benefits of the two approaches developed in Section 4.3.1 and Section 4.3.2 in the global outlier scenario detailed above for both the real and synthetic dataset. We want to answer the following question:

• "Does the use of a Huber loss while training improves the mean square error at the testing phase when the training data is contaminated with global outliers?"

Real dataset Because of the limited amount of data provided for the Lip experiments, we are able to compute the leave-one-out (LOO) mean square error associated to the vectorial and integral Huber loss regressors. The setting is as follows: we augment the dataset with 4 global outliers (12.5%) generated by the procedure described above with various scale levels s o ∈ {0.5, 1, 2, 5}. With a slight abuse, the hyperparameters λ, γ X , γ Θ were computed once for each scale level using LOO cross validation for the ridge regressor, and kept fixed for the various settings obtained by trading the square loss against a Huber loss, thus results presented here are pessimistic w.r.t. the performances of the Huber loss estimators. For a given scale level, as κ grows, the constraint in the dual brought by the Huber loss becomes void, and we recover the ridge regression solution from [START_REF] Kadri | Operatorvalued kernels for learning from functional response data[END_REF] in the vectorial Huber case, and from [START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel Hilbert spaces[END_REF] in the integral Huber case. This phenomenon is illustrated in Figure 4.6, we observe that there exist a range of κ for which the LOO error (computed with the square loss) is better than the one of a model trained with a square loss itself, which suggests that in presence of outliers the dual constraint on the norm of the coefficients of the model helps to achieve robustness.

Synthetic dataset

For the synthetic dataset, the procedure is roughly the same, except that it is now cheap to have access to new test samples, so that we do not use LOO error but rather measure the efficiency of the estimator on a fresh test set drawn from the same distribution as the inliers. We consider 250 normal observations from the dataset illustrated in Figure 4.4, and augment the training dataset with 25 outliers generated as in Figure 4.5 with various scale levels s o ∈ {2, 4, 6, 10}. We then tune the hyperparameters by cross-validation on the gathered dataset, and report in Figure 4.7 mean square error obtained on a fresh test set generated from the inliers distribution. Numerically, there still exists a range of κ for which the Huber estimator performs better, but it is harder to see it with eye. This is certainly due to the fact that in this case it is harder to differentiate between the inliers and the outliers, whereas for the real dataset usecase the outliers were designed to maximize the confusion brought in to the estimator.

ϵ-insensitive Losses

In this section, we present numerical experiments to illustrate the effectiveness of the approaches developed in Section 4.4. In particular, we want to answer the following question:

• "At what price (in terms of test error) can we get sparse estimators?"

Here, the sparsity is to be understood in terms of the percentage of zero coefficients. As seen in the corresponding optimization problems, the vectorial and integral ϵ-insensitive losses induce sparsity as they respectively bring in Group-LASSO and LASSO penalization terms to the dual problem. When ϵ = 0, we recover the unconstrained ridge regression solution, and as ϵ grows, the optimal coefficients become sparser. We expect that as ϵ grows, there is a degradation of the performances of the model in terms of test mean square error. Intuitively, the ϵ-insensitive ridge loss does not penalize the residuals whose norm is smaller than ϵ, and we propose to experiment with these losses in the context of local outliers, when some amount of noise is artificially added to the data measurement points. We present an illustration of this behavior for the vectorial ϵ-ridge, the two estimators behaving similarly. Real Dataset Similarly to the Huber case, we used LOO error as a measure of performance and the same hyperparameters (λ, γ X , γ Θ ). For the vectorial ϵ-ridge case, we plot in Figure 4.8 the LOO error and the proportion of zero coefficients. Since the problem is akin to a group lasso, this proportion corresponds to the number of training samples used in the estimator, as the • 2,1 penalty forces some lines of the coefficient matrix to be zero. The local outliers were obtained by added randomly a Gaussian noise with standard deviation s o ∈ {0.2, 0.5, 0.8, 1.2} at 20% of the locations associated to 50% of the training data. As expected, the LOO error degrades when ϵ grows. We observe little influence of the scale of the outliers on this experiment.

Synthetic Dataset

We used test MSE on fresh test samples to assess the performance of the estimator. The local outliers were obtained by adding Gaussian noise with standard deviation s o ∈ {0.4, 1.2, 1.7} at 20% of the locations associated to 50% of the training data. The results are illustrated in Figure 4.9, we observe similarly to the real dataset case that as ϵ grows, the coefficients of the estimators grow sparser and the test MSE degrades.

Conclusion

In this chapter, we proposed a functional output regression framework based on duality in vv-RKHS that allows to go beyond the square loss case. The use of convoluted losses, particularly suited to duality thanks to their compatibility with the Fenchel-Legendre transform, produces a family of problems that read in the dual as the original regularized empirical risk minimization problem plus a term akin to a penalty, depending on the chosen loss. This term can enforce robustness, as seen with the family of the Huber losses, or sparsity when choosing the family of ϵ-insensitive losses. To obtain them, we convolve the square loss with a loss that is norm dependent, inviting to investigate in future work the choice of the norm used in the right term of the convolution. This would likely lead to different ways to enforce sparsity or robustness. In this chapter, we propose to exploit the general framework of learning with integral losses from Chapter 3 to adopt a new angle to multi-task learning. Multi-task learning consists in leveraging dependency across tasks to jointly solve multiple task with a single vector-valued model [START_REF] Micchelli | On learning vector-valued functions[END_REF]. We extend the idea of multiple tasks to a continuum of parameterized tasks, and present a principled way to jointly learn these tasks referred to as infinite task learning (ITL).

Infinite Task Learning

After a brief introduction in Section 5.1, we devote Section 5.2 to the definition of the ITL framework, before studying its generalization capabilities in Section 5.3. In Section 5.4 we apply the ITL methodology to quantile regression and demonstrate its efficiency on both synthetic and real-world datasets. Section 5.5 is devoted to the cost-sensitive classification problem, and finally Section 5.6 explores the unsupervised density level set estimation problem under the ITL angle. Conclusions are drawn in Section 5.7.

Introduction

Several fundamental problems in machine learning and statistics can be phrased as the minimization of a loss function described by a hyperparameter. The hyperparameter might capture numerous aspects of the problem: (i) the tolerance w.r.t. outliers as the ϵ-insensitivity in support vector regression (SVR; [START_REF] Drucker | Support vector regression machines[END_REF]), (ii) importance of smoothness or sparsity such as the weight of the l 2 -norm in Tikhonov regularization [START_REF] Tikhonov | Solution of Ill-posed Problems[END_REF], l 1 -norm in LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF], or more general structured-sparsity inducing norms [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF], (iii) density level set estimation (DLSE), see for example one-class support vector machines (OCSVM), (iv) confidence as exemplified by quantile regression (QR), or (v) importance of different decisions as implemented by cost-sensitive classification (CSC). In various cases including QR, CSC or DLSE, one is interested in solving the parameterized task for several hyperparameter values. Multi-task learning (MTL; Evgeniou and Pontil 2004) provides a principled way of benefiting from the relationship between similar tasks while preserving local properties of the algorithms: ν-property in DLSE [START_REF] Glazer | q-OCSVM: A q-quantile estimator for high-dimensional distributions[END_REF] or quantile property in QR [START_REF] Takeuchi | Nonparametric quantile estimation[END_REF].

A natural extension from the traditional multi-task setting is to provide a prediction tool being able to deal with any value of the hyperparameter. In their seminal work, [START_REF] Takeuchi | Parametric task learning[END_REF] extended multi-task learning by considering an infinite number of parameterized tasks in a framework called parametric task learning. Assuming that the loss is piecewise affine in the hyperparameter, the authors were able to get the whole solution path through parametric programming, relying on techniques developed by [START_REF] Hastie | The entire regularization path for the support vector machine[END_REF]. 1 In this chapter, we relax the affine model assumption on the tasks as well as the piece-wise linear assumption on the loss, and take a different angle. We propose infinite task learning (ITL) within the framework of function-valued function learning to handle a continuum number of parameterized tasks. For that purpose we leverage tools from operator-valued kernels and the associated vv-RKHS. The idea is that the output is a function over the hyperparameter space, modelled as an element of a scalar-valued RKHS. Properties of this kernel, e.g. continuity, give an explicit control on the relationship between tasks, and manipulating output functions in a RKHS also enables us to consider incorporate specific constraints on their nature. In the studied framework each task is described by a (scalar-valued) RKHS over the input space which is capable of dealing with nonlinearities. The resulting ITL formulation relying on vv-RKHS specifically encompasses existing multi-task approaches including joint quantile regression [START_REF] Sangnier | Joint quantile regression in vector-valued RKHSs[END_REF] or multi-task variants of density level set estimation [START_REF] Glazer | q-OCSVM: A q-quantile estimator for high-dimensional distributions[END_REF]) by encoding a continuum of tasks.

From Parameterized to Infinite Task Learning

In this section, we gradually define our goal by moving from single parameterized tasks (Section 5.2.1) to infinite task learning (ITL; Section 5.2.3) through multi-task learning (MTL; Section 5.2.2).

Learning Parameterized Tasks

A supervised parametrized task is defined as follows. Let (X, Y) ∈ X × R be a random variable with joint distribution P X,Y which is assumed to be fixed but unknown. We have access to n i.i.d. observations called training samples: S := (x i , y i ) n i=1 ∼ P ⊗n X,Y . Let Θ be the domain of hyperparameters, and ℓ : Θ × R × R → R be a loss function. Let H k X ⊂ F X; R denote our hypothesis class; it is assumed to be a RKHS with kernel k X : X × X → R. For a given θ, the goal is to estimate the minimizer of the expected risk

R(θ, h) := E X,Y [ℓ(θ, h(X), Y)]
(5.1) over H k X , using the training sample S. This task can be addressed by solving the regularized empirical risk minimization problem

min h∈H k X R S (θ, h) + Ω(h), (5.2)
where R S (θ, h) := 1 n n i=1 ℓ(θ, h(x i ), y i ) is the empirical risk and Ω : H k X → R is a regularizer. Below we give three examples.

Quantile regression (QR):

In this setting θ ∈ 0, 1 . For a given hyperparameter level θ, in QR the goal is to predict the θ-quantile of the real-valued output conditional distribution P Y|X . The task can be tackled using the pinball loss [START_REF] Koenker | Regression quantiles[END_REF]) defined in Equation (5.3) and illustrated in Figure 5.1.

ℓ(θ, h(x), y) = max(θ(y -h(x)), (θ -1)(y -h(x))),

(5.3)

Ω(h) = λ 2 h 2 H k X (λ > 0).

Cost-sensitive classification (CSC): Our next example considers binary classification (Y ∈ {-1, 1}

) where a (possibly) different cost is associated with each class, as it is often the case in medical diagnosis or default detection [START_REF] Japkowicz | The class imbalance problem: a systematic study[END_REF][START_REF] Elkan | The foundations of cost-sensitive learning[END_REF]. The sign of h ∈ H k X yields the estimated class and in cost-sensitive classification one takes

ℓ(θ, h(x), y) = θ + 1 2 -{-1} (y) max(0, 1 -yh(x)),
(5.4)

Ω(h) = λ 2 h 2 H k X λ > 0.
The θ ∈ [-1, 1] hyperparameter captures the trade-off between the importance of correctly classifying the samples having -1 and +1 labels. When θ is close to -1, the obtained h focuses on classifying well class -1, and vice-versa. Typically, it is desirable for a physician to choose a posteriori the value of the hyperparameter at which she wants to predict. Since this cost can rarely be considered to be fixed, this motivates to learn one model giving access to all hyperparameter values.

Density level sets estimation (DLSE):

Examples of parameterized tasks can also be found in the unsupervised setting. For instance in outlier detection, the goal is to separate outliers from inliers. A classical technique to tackle this task is the celebrated one-class support vector machines (OCSVM; [START_REF] Schölkopf | New support vector algorithms[END_REF]. OCSVM has a free parameter θ ∈ (0, 1], which can be proven to be an upper bound on the fraction of outliers. When using a Gaussian kernel with a bandwidth converging to zero, OCSVM consistently estimates density level sets [START_REF] Vert | Consistency and convergence rates of one-class SVMs and related algorithms[END_REF]. This unsupervised learning problem can be empirically described by the minimization of a regularized empirical risk R S (θ, t, h) + Ω(h), solved jointly over h ∈ H k X and t ∈ R with

ℓ(θ, h(x), t) = -t + 1 θ t -h(x) + , Ω(h) = 1 2 h 2 H k X .

Solving a Finite Number of Tasks as Multi-Task Learning

In all the aforementioned problems, one is rarely interested in the choice of a single hyperparameter value (θ) and associated risk R S (θ, •) , but rather in the joint solution of multiple tasks. The naive approach of solving the different tasks independently can easily lead to inconsistencies. A principled way of solving many parameterized tasks has been cast as a MTL problem [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF] which takes into account the similarities between tasks and helps providing consistent solutions. For example it is possible to encode the similarities of the different tasks in MTL through an explicit constraint function [START_REF] Ciliberto | Consistent multitask learning with nonlinear output relations[END_REF]. In the ITL approach, the similarity between tasks is designed in an implicit way through the use of a kernel on the hyperparameters. Moreover, in contrast to MTL, in our case the input space and the training samples are the same for each task; a task is specified by a value of the hyperparameter. This setting is sometimes referred to as multi-output learning (Álvarez et al., 2012).

Formally, assume that we have m tasks described by parameters (θ j ) m j=1 . The idea of multi-task learning is to minimize the sum of the local loss functions R S (θ j , •), i.e.

arg min

h∈H m j=1 R S (θ j , h j ) + Ω(h),
where the individual tasks are modelled by the real-valued h j functions, the overall R m -valued model is the vector-valued function h : x → (h 1 (x), . . . , h m (x)) belonging to some hypothesis space H, and Ω is a regularization term encoding similarities between tasks.

It is instructive to consider two concrete examples:

• In joint quantile regression one can use the regularizer to encourage that the predicted conditional quantile estimates for two similar quantile values are similar. This idea forms the basis of the approach proposed by Sangnier et al. (2016) who formulates the joint quantile regression problem in a vector-valued reproducing kernel Hilbert space with an appropriate decomposable kernel that encodes the links between the tasks. The obtained solution shows less quantile curve crossings compared to estimators not exploiting the dependencies of the tasks as well as an improved accuracy.

• A multi-task version of DLSE has been presented by [START_REF] Glazer | q-OCSVM: A q-quantile estimator for high-dimensional distributions[END_REF] with the goal of obtaining nested density level sets as θ grows. Similarly to joint quantile regression, it is crucial to take into account the similarities of the tasks in the joint model to efficiently solve this problem.

Towards Infinite Task Learning

In the following, we propose a novel framework called infinite task learning (ITL) in which we learn a function-valued function h ∈ H K where H K ⊂ F(X, F(Θ, R)) is a vv-RKHS. Our goal is to be able to handle new tasks after the learning phase and thus, not to be limited to given predefined values of the hyperparameter. Regarding this goal, our framework generalizes the parametric task learning (PTL) approach introduced by Takeuchi et al. ( 2013), by allowing a wider class of models and relaxing the hypothesis of piece-wise linearity of the loss function. Moreover a nice byproduct of this vv-RKHS based approach is that one can benefit from the functional point of view, design new regularizers and impose various constraints on the whole continuum of tasks, e.g. ,

• The continuity of the θ → h(x)(θ) function is a natural desirable property: for a given input x, the predictions on similar tasks should also be similar.

• Another example is to impose a shape constraint in QR: the conditional quantile should be increasing w.r.t. the hyperparameter θ. This requirement can be imposed through the functional view of the problem.

• In DLSE, to get nested level sets, one would want that for all x ∈ X, the decision

function θ → R + (h(x)(θ) -t(θ)
) changes its sign only once.

To keep the presentation simple, in the sequel we are going to focus on ITL in the supervised setting; we dedicate Section 5.6 to the unsupervised task of DSLE.

We introduce the integral loss function

I ℓ (h(x), y) := Θ ℓ(θ, h(x)(θ), y)dµ(θ), (5.5)
where ℓ : Θ × R × R → R is a loss function and µ is a probability measure on Θ which encodes the importance of the prediction at different hyperparameter values. Without prior information and for compact Θ, one may consider µ to be uniform. The true and empirical risks read then as

R(h) := E X,Y I ℓ (h(X), Y) , R S (h) := 1 n n i=1 I ℓ (h(x i ), y i ). (5.6)
Intuitively, minimizing the expectation of the integral over θ in a rich enough space corresponds to searching for a pointwise minimizer x → ĥ(x)(θ) of the parametrized tasks introduced in Equation (5.1) with, for instance, the implicit space constraint that θ → ĥ(x)(θ) is a continuous function for each input x. We show in Proposition 5.25 that this is precisely the case in QR.

Interestingly, the empirical counterpart of the true risk minimization can now be considered with a much richer family of penalty terms:

min h∈H K R S (h) + Ω(h).
(5.7)

Here, Ω(h) can be a weighted sum of various penalties

• imposed directly on h such as λ 2 h 2 H K
for λ > 0, in which case we recognize in Problem 5.7 the type of problems whose resolution has been discussed at length in Chapter 3 or 88

CHAPTER 5. INFINITE TASK LEARNING • integrated constraints on either θ → h(x)(θ) or x → h(x)(θ) such as X Ω 1 (h(x)(•))dP(x) or Θ Ω 2 (h(•)(θ))dµ(θ)
which allow the property enforced by Ω 1 or Ω 2 to hold pointwise on X or Θ respectively.

It is worthwhile to see a concrete example before turning to the statistical analysis in Section 5.3: in quantile regression, the monotonicity assumption of the θ → h(x)(θ) function can be encoded by choosing Ω 1 as

Ω 1 (f ) = λ nc Θ -f ′ (θ) + dµ(θ)
when the output RKHS H k Θ is populated with differentiable functions and λ nc > 0 controls the strength of the soft constraint.

Generalization Analysis through Uniform Stability

In this section, we study the generalization gap (or excess risk) associated to the ITL method. The goal here is to control with high probability the quantity E( ĥ) := R( ĥ) -R S ( ĥ),

(

where ĥ is the estimator resulting from the ITL method in the supervised setting:

ĥ = arg min h∈H K 1 n n i=1 I ℓ (h(x i ), y i ) + λ h 2 H K .
(5.9)

Remark 5.1. The reader can notice that the regularization parameter here is λ and not λ 2 as it helps readability of the generalization bound.

When E( ĥ) is small, the learned estimator performs comparably on unseen data drawn with probability P (X,Y) , whereas when it is large, the estimator overfits, meaning it performs poorly on new data. A myriad of techniques have been developed to study the generalization error [START_REF] Steinwart | Mercer's theorem on general domains: on the interaction between measures, kernels, and RKHSs[END_REF][START_REF] Zhou | The covering number in learning theory[END_REF]; Rudi and Rosasco, 2017), extending the seminal work of [START_REF] Vapnik | Principles of risk minimization for learning theory[END_REF]). In the particular kernel based learning, we can invoke the Rademacher complexity tool (Bartlett and Mendelson, 2002) which allows to bound sup h∈H K |E(h)| and benefits from a large body of work, with extensions to vector-valued learning [START_REF] Maurer | Bounds for vector-valued function estimation[END_REF]. However, these extensions rely on the assumption of trace class OVKs which is not verified for the kernel K = k X Id H k Θ . Moreover, bounding sup h∈H |E(h)| appears unnecessarily hard, since we are only interested in bounding E( ĥ).

These considerations lead us to work in the framework of uniform stability introduced in [START_REF] Bousquet | Stability and generalization[END_REF]. Uniform stability asks the question: "what happens when we learn an estimator with a slightly modified dataset ?" and answers by saying that when the two resulting estimators are close, then it is possible to bound the excess risk of the output of the algorithm.

Definition 5.2. Let S = x i , y i n i=1 be the training data. We call S i the training data S i = ((x 1 , y 1 ), . . . , (x i-1 , y i-1 ), (x i+1 , y i+1 ), . . . , (x n , y n )), i ∈ [n].

In the following, L : H K ×X×R is a loss function whose dependence w.r.t. to the model h is explicit, and the risks are

R(h) = E (X,Y) L(h, X, Y) , R S (h) = 1 n n i=1 L(h, x i , y i ).
Remark 5.3. To fit in this framework, we will choose either

L(h, x, y) = I ℓ (h(x), y) or L(h, x, y) = I ℓ (h(x), y) (5.10)
depending on whether we analyse the theoretical ITL estimator (see Section 5.3.1) or the one obtained from the practical solution with a sampled empirical risk and double representer theorem (see Section 5.3.2).

Definition 5.4. A learning algorithm mapping a dataset S to an estimator ĥS is said to be β-uniformly stable with respect to the loss function

L if for all n ∈ N * , i ∈ [n],
and

S training set, L( ĥS , •, •) -L( ĥS i , •, •) ∞ ≤ β.
(5.11) Assumption 5.5. There exists ξ ≥ 0 such that for all (x, y) ∈ X × Y and training set S, L( ĥS , x, y) ≤ ξ.

Having defined these concepts, we are now ready to state the generalization bound from [START_REF] Bousquet | Stability and generalization[END_REF]) that allows to quantify the excess risk of an estimator provided by a learning algorithm.

Proposition 5.6. [START_REF] Bousquet | Stability and generalization[END_REF] Let S → ĥS be a learning algorithm with uniform stability β with respect to a loss L satisfying Assumption 5.5. Then for all n ≥ 1 and δ ∈ (0, 1), with probability at least 1 -δ on the drawing of the samples, it holds that

E ĥS ≤ 2β + (4nβ + ξ) log (1/δ) n .
Proposition 5.6 ensures that whenever √ nβ → n→∞ 0, the true risk R( ĥS ) is controlled by its empirical version R S ( ĥS ), and the bounds are tight when β scales as 1 n . The bounds associated to uniform stability have recently benefited from sharper results proposed by Feldman and Vondrak (2018) and later [START_REF] Bousquet | Sharper bounds for uniformly stable algorithms[END_REF]. We present it below for completeness.

Proposition 5.7. [START_REF] Bousquet | Sharper bounds for uniformly stable algorithms[END_REF] Let S → ĥS be a learning algorithm with uniform stability β with respect to a loss L satisfying Assumption 5.5. Then for all n ≥ 1 and δ ∈ (0, 1), with probability at least 1 -δ on the drawing of the samples, it holds that

E ĥS ≤ c 1 β log 2 (n) log (e/δ) + c 2 ξ log (e/δ) n
where c 1 = 12 √ 2e and c 2 = 4e.

This bound shows both sub-Gaussian ( log e/δ term) and sub-exponential (log e/δ term) regimes and can be of great interest for β that decrease slower than 1 n (see [START_REF] Bousquet | Sharper bounds for uniformly stable algorithms[END_REF] for a discussion about this). In our case as we will see β scales as 1 n ; we present the results implied by Proposition 5.6 and note the analysis can be adapted to Proposition 5.7. This is the strength of the uniform stability framework, all we have to do is prove that the algorithms are β-stable, and we can profit from both bounds.

Uniform stability and vv-RKHS learning

We now exploit results from Kadri et al. (2016) (we can also mention the initial work of Audiffren and Kadri (2013)), who studied the uniform stability of learning algorithms based on empirical risk minimization with vv-RKHS norm regularization. There is a slight difference between their setting and ours, since they use losses defined for some y in the output space of the vv-RKHS, but this difference has no impact on the validity of the proofs. Their work relies on the following assumptions: Assumption 5.8. There exists κ > 0 such that for ∀x ∈ X, K(x, x) op ≤ κ 2 . Assumption 5.9. For all h 1 , h 2 ∈ H k Θ , the function

(x 1 , x 2 ) ∈ X × X → K(x 1 , x 2 )h 1 , h 2 H k Θ ∈ R is measurable.
We recall the choice of OVK we make: K(x, z) = k X (x, z)I H k Θ with (x, z) ∈ X × X, k X and k Θ are bounded scalar-valued kernels; in other words there exist (κ

X , κ Θ ) ∈ R 2 such that sup x∈X k X (x, x) ≤ κ 2 X and sup θ∈Θ k Θ (θ, θ) ≤ κ 2 Θ .
Remark 5.10. Assumptions 5.8, 5.9 are satisfied for our choice of kernel.

Assumption 5.11. The application (y, h, x) → L(y, h, x) is σ-admissible, i.e. convex with respect to h and Lipschitz continuous with respect to f (x), with σ as its Lipschitz constant.

We now state a proposition from Kadri et al. ( 2016) which gives β-stability guarantees for the output of the algorithm ĥS = arg min

h∈H K 1 n n i=1 L(h, x i , y i ) + λ h 2 H K .
(5.12) Proposition 5.12. [START_REF] Kadri | Operatorvalued kernels for learning from functional response data[END_REF] Under Assumption 5.8, Assumption 5.9, and Assumption 5.11, a learning algorithm that maps a training set S to the function ĥS defined in Equation (5.12) is β-stable with β = σ 2 κ 2 2λn .

Generalization Bound for Non-Approximated Scheme

In this section, we derive generalization bounds based on Proposition 5.12 for the non-approximated ITL scheme, in the sense that we analyse the theoretical estimator without considering the approximations induced by the practical optimization 5.3. GENERALIZATION ANALYSIS THROUGH UNIFORM STABILITY 91 algorithm:

ĥS = arg min h∈H K 1 n n i=1 L(h, x i , y i ) + λ h 2 H K (5.13) where L(h, x, y) = Θ ℓ(θ, h(x)(θ), y)dµ(θ).
(5.14)

Quantile regression

We recall that in this setting,

ℓ(θ, h(x)(θ), y) = max (θ(y -h(x)(θ)), (1 -θ)(y -h(x)(θ)))
and consider the associated L from Equation (5.14). Moreover, we will assume that |Y| is bounded almost surely by B ∈ R. We will therefore verify the hypothesis for y ∈ [-B, B] and not y ∈ R.

Lemma 5.13. In the case of the QR, the loss L is σ-admissible with σ = 2κ Θ .

Proof Let h 1 , h 2 ∈ H K and θ ∈ [0, 1]. ∀x, y ∈ X × R, it holds that ℓ(θ, h 1 (x)(θ), y) -ℓ(θ, h 2 (x)(θ), y) = (θ -t)(h 2 (x)(θ) -h 1 (x)(θ)) + (t -s)(y -h 1 (x)(θ)),
where s = y≤h 1 (x)(θ) and t = y≤h 2 (x)(θ) . We consider all possible cases for t and s :

• t = s = 0 : |(t -s)(y -h 1 (x)(θ))| ≤ |h 2 (x)(θ) -h 1 (x)(θ)|, • t = s = 1 : |(t -s)(y -h 1 (x)(θ))| ≤ |h 2 (x)(θ) -h 1 (x)(θ)|, • s = 1, t = 0 : |(t -s)(y -h 1 (x)(θ))| = |h 1 (x)(θ) -y| ≤ |h 1 (x)(θ) -h 2 (x)(θ)|, • s = 0, t = 1 : |(t -s)(y -h 1 (x)(θ))| = |y -h 1 (x)(θ)| ≤ |h 1 (x)(θ) -h 2 (x)(θ)|
because of the conditions on t, s.

Thus

ℓ(θ, h 1 (x)(θ), y) -ℓ(θ, h 2 (x)(θ), y) ≤ (θ + 1)|h 1 (x)(θ) -h 2 (x)(θ)| ≤ (θ + 1)κ Θ ||h 1 (x) -h 2 (x)|| H k Θ ,
where we used the Cauchy-Schwartz inequality applied to h 1 (x) -h 2 (x) and k Θ (•, θ) in H k Θ at the last line. By integrating this expression over the Θ space, we get that

L(h 1 , x, y) -L(h 2 , x, y) ≤ 1 0 (θ + 1)κ Θ ||h 1 (x) -h 2 (x)|| H k Θ dµ(θ) ≤ 2κ Θ ||h 1 (x) -h 2 (x)|| H k Θ
and L is σ-admissible with σ = 2κ Θ .

Lemma 5.14. Let S = ((x 1 , y 1 ), . . . , (x n , y n )) be a training set and λ > 0. Then for

all (x, θ) ∈ X × (0, 1), it holds that | ĥS (x)(θ)| ≤ κ X κ Θ B λ .
Proof Since ĥS is the output of our algorithm and 0 ∈ H K , it holds that

λ|| ĥS || 2 ≤ 1 n n i=1 1 0 ℓ(θ, 0, y i )dµ(θ) ≤ 1 n n i=1 1 0 max (θ, 1 -θ)|y i |dµ(θ) ≤ B. 92 CHAPTER 5. INFINITE TASK LEARNING Thus || ĥS || ≤ B λ . Moreover, ∀x, θ ∈ X × (0, 1), ĥS (x)(θ) = ĥS (x), k Θ (θ, •) H k Θ ≤ κ Θ ĥS (x) H k Θ by Cauchy-Schwartz inequality in H k Θ ≤ κ X κ Θ ĥS H K
by the bounded operator norm of K ♯ x which concludes the proof.

Lemma 5.15. Assumption 5.5 is satisfied for

ξ = 2 B + κ X κ Θ B λ .
Proof Let S = ((x 1 , y 1 ), . . . , (x n , y n )) be a training set and ĥS be the output of our algorithm. For all (x, y)

∈ X × [-B, B], it holds that L(y, ĥS , x) = 1 0 max (θ(y -ĥS (x)(θ)), (θ -1)(y -ĥS (x)(θ)))dµ(θ) ≤ 2 1 0 y -ĥS (x)(θ) dµ(θ) ≤ 2 1 0 |y| + ĥS (x)(θ) dµ(θ) ≤ 2   B + κ X κ Θ B λ   .
Proposition 5.16. The QR learning algorithm defined in Equation (5.13) is such that for all n ≥ 1 and δ ∈ (0, 1), with probability at least 1-δ on the drawing of the samples, it holds that

R ĥS ≤ R S ĥS + 4κ 2 X κ 2 Θ λn +    8κ 2 X κ 2 Θ λ + 2   B + κ X κ Θ B λ      log (1/δ) n . (5.15)
Proof This is a direct consequence of Proposition 5.12, Proposition 5.6, Lemma 5.13 and Lemma 5.15.

Cost-sensitive classification

In this setting, the loss ℓ is

ℓ(θ, h(x)(θ), y) = θ + 1 2 -{-1} (y) 1 -yh(x)(θ) +
and the associated L is given by Equation (5.14). One can verify (the way it was done for quantile regression) that the properties above still hold, but with constants

σ = κ Θ , β = κ 2 X κ 2 Θ 2λn , ξ = 1 + κ X κ Θ √ λ .
This allows to formulate the following bound.

Proposition 5.17. The CSC learning algorithm defined in Equation (5.13) is such that for all n ≥ 1 and δ ∈ (0, 1), with probability at least 1 -δ on the drawing of the samples, it holds that

R ĥS ≤ R S ĥS + κ 2 X κ 2 Θ λn + 2κ 2 X κ 2 Θ λ + 1 + κ X κ Θ √ λ log (1/δ) n .

Generalization Bound for Approximated Scheme

In this section, we derive generalization bounds based on Proposition 5.12 for the approximated ITL scheme ĥS = arg min

h∈H K 1 n n i=1 L(h, x i , y i ) + λ h 2 H K where L(h, x, y) = 1 m m j=1 ℓ(θ j , h(x)(θ j ), y).
(5. [START_REF] Vapnik | Statistical Learning Theory[END_REF] and (θ j ) m j=1 is a quasi Monte-Carlo sequence used to approximate the loss function as proposed in Section 3.3.1. We remind the reader of the statistical quantities:

R ĥS = E X,Y I ℓ ( ĥS (X), Y) , R S ĥS = 1 n n i=1 I ℓ ( ĥS (x i ), y i ) R ĥS = E X,Y    1 m m j=1 ℓ(θ j , ĥS (X)(θ j ), Y)    , R S ĥS = 1 nm n i=1 m j=1 ℓ(θ j , ĥS (x i )(θ j ), y i ).
The uniform stability framework then informs about the generalization gap between R and R S , and additional work is needed to control the quantity R -R S . In particular we make use of convergence properties of QMC integration which benefit from O log m m convergence rates for bounded variation functions.

Definition 5.18 (Hardy-Krause variation). Let Π be the set of subdivisions of the interval Θ = [0, 1]. A subdivision will be denoted by σ = (θ 1 = 0, θ 2 , . . . , θ p = 1) and f : Θ → R be a function. Using these notations, the Hardy-Krause variation of the function f is defined as

V (f ) = sup σ∈Π p-1 i=1 f (θ i+1 ) -f (θ i ) .
Remark 5.19. If f is continuous, V (f ) is also the limit as the mesh of σ goes to zero of the above quantity.

In the following, let f : θ → E X,Y ℓ(θ, ĥS (X)(θ), Y) . This function is of primary importance for our analysis, since in the quasi Monte-Carlo setting, the convergence properties of the integration scheme only hold if the function to integrate has finite Hardy-Krause variation.

CHAPTER 5. INFINITE TASK LEARNING

Quantile regression

The following lemma quantifies the Hardy-Krause variation of the function f defined above in the QR setting.

Lemma 5.20. Assume the boundeness of both scalar-valued kernels k X and k Θ . Assume moreover that k Θ is C 1 and that its partial derivatives are uniformly bounded by some constant C. Finally, assume that Y is almost surely bounded by some constant B. Then

V (f ) ≤ B + κ X κ Θ B λ + 2κ X 2BC λ .
(5.17)

Proof It holds that

V (f ) = sup σ∈Π p-1 i=1 f (θ i+1 ) -f (θ i ) = sup σ∈Π p-1 i=1 E X,Y ℓ(θ i+1 , ĥS (X)(θ i+1 ), Y) -E X,Y ℓ(θ i , ĥS (X)(θ i ), Y) = sup σ∈Π p-1 i=1 E X,Y ℓ(θ i+1 , ĥS (X)(θ i+1 ), Y) -ℓ(θ i , ĥS (X)(θ i ), Y) ≤ sup σ∈Π p-1 i=1 E X,Y ℓ(θ i+1 , ĥS (X)(θ i+1 ), Y) -ℓ(θ i , ĥS (X)(θ i ), Y) ≤ sup σ∈Π E X,Y   p-1 i=1 ℓ(θ i+1 , ĥS (X)(θ i+1 ), Y) -ℓ(θ i , ĥS (X)(θ i ), Y)   .
The supremum of the expectation is smaller than the expectation of each supremum, hence

V (f ) ≤ V (f x,y )dP X,Y , (5.18)
where f x,y : θ → ℓ(θ, y, ĥS (x)(θ)) is the local counterpart of the function f at point (x, y). To bound this quantity, let us first bound locally V (f x,y ). To that extent, we fix some (x, y) in the following. Since f x,y is continuous (because k Θ is C 1 ), then using Choquet (1969, Theorem 24.6), it holds that

V (f x,y ) = lim |σ|→0 p-1 i=1 f x,y (θ i+1 ) -f x,y (θ i ) .
Moreover since k ∈ C 1 and ∂k θ = (∂ 1 k)(•, θ) has a finite number of zeros for all θ ∈ Θ, one can assume that in the subdivision considered afterhand all the zeros (in θ) of the residuals y -ĥS (x)(θ) are present, so y -ĥS (x)(θ i+1 ) and y -ĥS (x)(θ i ) are always of the same sign. Indeed, if not, create a new, finer subdivision with this property and work with this one. Let us begin the proper calculation: let σ = (θ 1 , θ 2 , . . . , θ p ) be a subdivision of Θ, it holds that ∀i ∈ 1, . . . , p -1 :

f x,y (θ i+1 ) -f x,y (θ i ) = max (θ i+1 (y -ĥS (x)(θ i+1 )), (1 -θ i+1 )(y -ĥS (x)(θ i+1 ))) -max (θ i (y -ĥS (x)(θ i )), (1 -θ i+1 )(y -ĥS (x)(θ i ))) .
We now study the two possible outcomes for the residuals:

5.3. GENERALIZATION ANALYSIS THROUGH UNIFORM STABILITY 95 • If y -h(x)(θ i+1 ) ≥ 0 and y -h(x)(θ i ) ≥ 0 then f x,y (θ i+1 ) -f x,y (θ i ) = θ i+1 (y -ĥS (x)(θ i+1 )) -θ i (y -ĥS (x)(θ i )) = (θ i+1 -θ i )y + (θ i -θ i+1 ) ĥS (x)(θ i+1 ) + θ i ( ĥS (x)(θ i ) -ĥS (x)(θ i+1 )) ≤ (θ i+1 -θ i )y + (θ i -θ i+1 ) ĥS (x)(θ i+1 ) + θ i ( ĥS (x)(θ i ) -ĥS (x)(θ i+1 )) .
From Lemma 5.14, it holds that ĥS (x

)(θ i+1 ) ≤ κ X κ Θ B λ . Moreover, ĥS (x)(θ i ) -ĥS (x)(θ i+1 ) = ĥS (x), k Θ (θ i , •) -k Θ (θ i+1 , •) H k Θ ≤ ĥS (x) H k Θ k Θ (θ i , •) -k Θ (θ i+1 , •) H k Θ ≤ κ X B λ k Θ (θ i , θ i ) + k Θ (θ i+1 , θ i+1 ) -2k Θ (θ i+1 , θ i ) ≤ κ X B λ k Θ (θ i+1 , θ i+1 ) -k Θ (θ i+1 , θ i ) + k Θ (θ i , θ i ) -k Θ (θ i+1 , θ i ) .
Since k Θ is C 1 , with partial derivatives uniformly bounded by C, it holds that

k Θ (θ i+1 , θ i+1 ) -k Θ (θ i+1 , θ i ) ≤ C(θ i+1 -θ i ) and k Θ (θ i , θ i ) -k Θ (θ i+1 , θ i ) ≤ C(θ i+1 -θ i ) so ĥS (x)(θ i ) -ĥS (x)(θ i+1 ) ≤ κ X 2BC λ θ i+1 -θ i and overall f x,y (θ i+1 ) -f x,y (θ i ) ≤   B + κ X κ Θ B λ   (θ i+1 -θ i ) + κ X 2BC λ θ i+1 -θ i . • If y -h(x)(θ i+1 ) ≤ 0 and y -h(x)(θ i ) ≤ 0, then f x,y (θ i+1 ) -f x,y (θ i ) = (1 -θ i+1 )(y -ĥS (x)(θ i+1 )) -(1 -θ i )(y -ĥS (x)(θ i )) ≤ ĥS (x)(θ i ) -ĥS (x)(θ i+1 ) + (θ i+1 -θ i )y + (θ i -θ i+1 ) ĥS (x)(θ i+1 ) + θ i ( ĥS (x)(θ i ) -ĥS (x)(θ i+1 ))
CHAPTER 5. INFINITE TASK LEARNING so with similar arguments one gets

f x,y (θ i+1 ) -f x,y (θ i ) ≤   B + κ X κ Θ B λ   (θ i+1 -θ i ) + 2κ X 2BC λ θ i+1 -θ i .
(5. 19) Therefore, regardless of the sign of the residuals y -h(x)(θ i+1 ) and y -h(x)(θ i ), one gets Equation (5. 19). Since the square root function has Hardy-Kraus variation of 1 on the interval Θ = [0, 1], it holds that

sup σ∈Π p-1 i=1 |f x,y (θ i+1 ) -f x,y (θ i )| ≤ B + κ X κ Θ B λ + 2κ X 2BC λ .
Combining this with Equation (5.18) finally gives

V (f ) ≤ B + κ X κ Θ B λ + 2κ X 2BC λ .
We are now ready to compare the true risk R of the ITL estimator with its sampled version R based on the sampled integral loss I ℓ .

Lemma 5.21. Let R be the risk defined in Equation (5.6) for the quantile regression problem. Assume that (θ) m j=1 have been generated via the Sobol sequence and that k Θ is C 1 with its partial derivatives being uniformly bounded by some constant C. Then

R ĥS -R ĥS ≤   B + κ X κ Θ B λ + 2κ X 2BC λ   log(m) m .
(5.20)

Proof Let f : θ → E X,Y ℓ(θ, ĥS (X)(θ), Y) . It holds that |R( ĥS ) -R( ĥS )| ≤ V (f ) log(m) m
according to the Koksma-Hlawka inequality [START_REF] Morokoff | Quasi-Monte Carlo integration[END_REF], where V (f ) is the Hardy-Krause variation of f . Lemma 5.20 allows then to conclude.

When gathered together, Lemma 5.21 and Proposition 5.16 provide a bound on the true risk of the estimator w.r.t. its sampled empirical risk.

Proposition 5.22. Let ĥS be the quantile estimator resulting from Equation (5.13) with sampled loss from Equation (5. [START_REF] Vapnik | Statistical Learning Theory[END_REF]. Assume that Y is almost surely bounded, the kernels k X , k Θ are bounded and k Θ is C 1 with bounded derivatives. Then it holds that

R ĥS ≤ R S ĥS + O P X,Y 1 λ √ n + O log m √ λm . (5.21)
Proof This is a direct application of Lemma 5.21 and Proposition 5.16.

QUANTILE REGRESSION 97

Cost-sensitive classification The bound for the CSC case follows the same reasoning, and the Hardy-Kraus variation of the function f : θ → E X,Y ℓ(θ, ĥS (X)(θ), Y) is bounded by the following lemma.

Lemma 5.23. Assume the boundeness of both scalar-valued kernels k X and k Θ . Assume moreover that k Θ is C 1 and that its partial derivatives are uniformly bounded by some constant C. Then

V (f ) ≤ 1 + κ X κ Θ 1 λ + 2κ X 2C λ .
(5.22)

Proof The proof is similar to the one of Lemma 5.20, replacing the residual y -ĥS (x) by 1 -y ĥS (x) and performing the case distinction whether 1 -y ĥS (x) is positive of negative.

Quantile Regression

This section is dedicated to quantile regression when using ITL to learn the whole quantile function. As a reminder, given an X-valued random variable X and real-valued output random variable Y we define its conditional quantile function by

q(x)(θ) = inf u∈R t : P[Y ≤ u|X = x] = θ , (x, θ) ∈ X × (0, 1).
While the conditional quantile function is not defined for θ = 0 and θ = 1, it can be extended to these values by continuity whenever Y is a bounded random variable.

Remark 5.24. By definition, for all x ∈ X, the conditional quantile function q(x) is a nondecreasing function. However, estimates based on the ITL method do not necessarily satisfy this property, and it can be challenging to impose such shape constraint for the estimator. This will be discussed in Section 5.4.1.

It is well-known [START_REF] Koenker | Regression quantiles[END_REF] that for any quantile level θ ∈ (0, 1) and point x ∈ X, the quantile enjoys a variational formula

q(x)(θ) = arg min u∈R E Y|X=x ℓ(θ, u, Y) ,
where the pinball loss is defined by

ℓ(θ, u, v) : (u, v) ∈ R 2 → max (θ(v -u), (θ -1)(v -u)) ∈ R, (5.23)
and it is illustrated in Figure 5.1. The following proposition validates the introduction of the integral pinball loss I ℓ to tackle the estimation of the whole quantile function in one go.

Proposition 5.25. Let X, Y be two r.v. respectively taking values in X and R, with Y bounded almost surely. Let q : X → F([0, 1], R) be the associated conditional quantile function and µ be a positive probability measure on [0, 1] such that µ({0}) = µ({1}) = 0.

Then for ∀h ∈ F(X, F([0, 1], R))

R(h) -R(q) ≥ 0,
where R is the risk defined in Equation (5.6).

Proof

The proof is based on the one given [START_REF] Li | Quantile regression in reproducing kernel Hilbert spaces[END_REF] for a single quantile. Let f ∈ F(X, F([0, 1], R)), θ ∈ (0, 1) and (x, y) ∈ X × R. Let also

s =    1 if y ≤ f (x)(θ) 0 otherwise , t =    1 if y ≤ q(x)(θ) 0 otherwise . It holds that ℓ(θ, h(x)(θ), y) -ℓ(θ, q(x)(θ), y) = = θ(1 -s)(y -h(x)(θ)) + (θ -1)s(y -h(x)(θ)) -θ(1 -t)(y -q(x)(θ)) -(θ -1)t(y -q(x)(θ)) = θ(1 -t)(q(x)(θ) -h(x)(θ)) + θ((1 -t) -(1 -s))h(x)(θ) +(θ -1)t(q(x)(θ -h(x)(θ))) + (θ -1)(t -s)h(x)(θ) + (t -s)y = (θ -t)(q(x)(θ) -h(x)(θ)) + (t -s)(y -h(x)(θ)).
Then, notice that

E (θ -t)(q(X)(θ) -h(X)(θ)) = E E (θ -t)(q(X)(θ) -h(X)(θ)) |X = E E (θ -t)|X (q(X)(θ) -h(X)(θ))
by the tower rule for expectations, and since q is the true quantile function,

E t|X = E {Y≤q(X)(θ)} |X = P Y ≤ q(X)(θ)|X = θ, so E (θ -t)(q(X)(θ) -h(X)(θ)) = 0.
Moreover, (t -s) is negative when q(x)(θ) ≤ y ≤ h(x)(θ), positive when h(x)(θ) ≤ y ≤ q(x)(θ) and 0 otherwise, thus the quantity (t -s)(y -h(x)(θ)) is always positive. As a consequence,

R(h) -R(q) = 1 0 E ℓ(θ, h(X)(θ), Y) -ℓ(θ, q(X)(θ), Y) dµ(θ) ≥ 0
which concludes the proof.

Proposition 5.25 allows us to derive conditions under which the minimization of the risk above yields the true quantile function. Under the assumption that (i) q is continuous (as seen as a function of two variables), (ii) supp(µ) = [0, 1], then the minimization of the integrated pinball loss performed in the space of continuous functions yields the true quantile function on the support of P X,Y .

The ITL framework proposes to learn this quantile function in a vv-RKHS H K by solving a regularized empirical risk problem: ĥ = arg min

h∈H K 1 n n i=1 I ℓ (h(x i ), y i ) + λ 2 h 2 H K . (5.24)
This problem can then be solved using the optimization algorithm developed in Chapter 3, in particular by sampling the integral loss and benefiting from a double representer theorem (see Section 3.3.1) or by using dual algorithms (see Section 3.4.3). In the following, we further explore this problem: in Section 5.4.1 we study ways to enforce shape constraints for the estimate by means of a soft penalty added to the objective function in Problem 5.24. Later in Section 5.4.2 we propose to smooth the pinball loss, to the benefit of enjoying fast solvers exploiting the differentiability of the loss, then study in Section 5.4.3 the impact of the number of sampled locations m used to approximate the integral loss. Finally in Section 5.4.4 we apply the ITL framework to a quantile estimation problem where the input data are images, and propose to use a neural network to learn the input kernel.

Enforcing Shape Constraints

By definition, the quantile function is nondecreasing with respect to the parameter θ.

It is desirable that an estimator ĥ output by the ITL scheme possesses this property, so that it can be interpreted as a quantile function. Visually, this means that given two quantile levels (θ 1 , θ 2 ) ∈ Θ 2 the curves associated to ĥ(•)(θ 1 ) and ĥ(•)(θ 2 ) must not cross. This phenomenon, denoted as crossing quantiles, is a challenge to avoid in practice. While previous approaches in multi-output learning try to control the finite differences ĥ(x)(θ j+1 ) -ĥ(x)(θ j ) [START_REF] Sangnier | Joint quantile regression in vector-valued RKHSs[END_REF], having a functional model allows to consider the partial derivative of the estimator with respect to θ, which should be nonnegative to respect the nondecreasibility condition. The smoothness of the model can be ensured by the smoothness of the output kernel: following Ziemer (2012), whenever

k Θ ∈ C 2 (Θ 2 , R) then any f ∈ H k Θ ∈ C 1 (Θ, R).
Denoting by ∂ Θ h the partial derivative operator in H K with respect to θ, one would ideally want to solve ĥ = arg min

h∈H K 1 n n i=1 I ℓ (h(x i ), y i ) + λ 2 h 2 H K s.t. (∂ Θ h)(x)(θ) ≥ 0 ∀(x, θ) ∈ X × Θ.
(5.

(∂f )(θ) = f, (∂ 2 k Θ )(•, θ) H k Θ , (f, θ) ∈ H k Θ × Θ. ( 5
(∂ Θ h)(x)(θ) = h, K x (∂ 2 k Θ )(•, θ) H K , (h, x, θ) ∈ H K × X × Θ.
(5.27)

This type of constraint belongs to the family of hard constraints and in our case the functional constraint prevents a tractable optimization scheme. We can point to Aubin-Frankowski and Szabó 2020 for a recent work on enforcing hard shape constraints in kernel based learning. To mitigate this bottleneck, we propose to penalize the objective function in Problem 5.24 with a soft constraint based on the derivative of the model:

Ω nc (h) := λ nc X Θ -(∂ Θ h)(x)(θ) + dP X (x)dµ(θ).
Given that such a penalty is not computable analytically, we approximate it and define

Ω nc (h) := λ nc 1 n nc m nc nnc i=1 mnc j=1 -(∂ Θ h)(x i )( θj ) + .
where (x i ) nnc i=1 and ( θj ) mnc j=1 form a grid at which the penalty is computed. They can be chosen different from the real data (x i ) n i=1 and sampled locations (θ j ) m j=1 used in the double representer theorem. We then arrive at the following optimization problem ĥ = arg min

h∈H K 1 n n i=1 I ℓ (h(x i ), y i ) + λ 2 h 2 H K + Ω nc (h).
(5.28)

Is is to be noticed that the resulting estimator does not necessarily prevents crossing quantile, hence the term soft constraints. The double representer expression from Theorem 3.9 can then be adapted to encompass this type of penalty, as emphasized in the theorem below.

Theorem 5.27. Problem 5.28 admits a unique solution ĥ ∈ H K , and there exist (α ij ) n,m i,j=1 ∈ R nm and ( βij ) nnc,mnc i,j=1

∈ R nncmnc such that for all (x, θ) ∈ X × Θ, ĥ(x

)(θ) = n,m i,j=1 αij k X (x, x i )k Θ (θ, θ j ) + nnc,mnc i,j=1 βij k X (x, xi )∂ 2 k Θ (θ, θj ).
(5.29)

Proof

The proof is similar to the one of Theorem 3.9. Define

E 1 = K x i k Θ (•, θ j ) n,m i,j=1 : (i, j) ∈ [n] × [m] , E 2 =    K xi ∂ 2 k Θ (•, θj ) nnc,mnc i,j=1 : (i, j) ∈ [n nc ] × [m nc ]    and E = Span {E 1 ∪ E 2 } ⊂ H K . Let J(h) = 1 n n i=1 I ℓ (h(x i ), y i ) + λ 2 h 2 H K
+ Ω(h). J is coercive and bounded below, so that there exists a minimizer ĥ ∈ H which is unique by strong convexity of J. Since E is a finite-dimensional subspace of H K , it is closed in H K , and

E ⊕ E ⊥ = H K .
By decomposing ĥ = ĥE + ĥE ⊥ and applying the reproducing property, we get that

J( ĥE + ĥE ⊥ ) = J( ĥE ) + λ 2 ĥE ⊥ 2 H K
, which in turns implies that ĥE ⊥ = 0 and validates the decomposition in Equation (5.29).

The benefit of Theorem 5.27 is to provide a workable finite-dimensional parameterization of ĥ. The size of the representation is nm + n nc m nc , in particular if we use the data points (x i ) n i=1 and sampled locations (θ j ) m j=1 in the computation of Ω, then the size of the representation is 2nm. The resulting finite dimensional problem can then be tackled by first order optimization methods such as subgradient descent.

The efficiency of the non-crossing penalty is illustrated in Figure 5.2 on a synthetic dataset from [START_REF] Sangnier | Joint quantile regression in vector-valued RKHSs[END_REF]. This dataset consists in a sine curve at 1Hz modulated by a sine envelope at 1/3Hz and mean 1, distorted with a Gaussian noise of mean 0 and a linearly decreasing standard deviation from 1.2 at x = 0 to 0.2 at x = 1.5.

Here n = 40 data samples and m = 20 sampling locations have been generated, and we used the same anchors for the non-crossing penalty. Many crossings are visible on the right plot, while they are almost not noticeable on the left plot, using the non-crossing penalty.

Concerning our real-world examples, to study the efficiency of the proposed scheme in quantile regression the following experimental protocol was applied. We worked with 20 UCI datasets, each dataset was splitted randomly into a training set (70%) and a test set (30%). We optimized the hyperparameters by minimizing a 5-folds cross validation with a Bayesian optimizer2 Once the hyperparameters were obtained, a new regressor was learned on the whole training set using the optimized hyperparameters. We report the value of the pinball loss and the crossing loss on the test set for three methods: our technique is called infinite quantile regression (IQR), we refer to Sangnier et al. (2016)'s approach as joint quantile regression (JQR), and independent learning (abbreviated as IND-QR) represents a further baseline. We repeated 20 simulations (different random training-test splits); the results are also compared using a Mann-Whitney-Wilcoxon test. A summary is provided in Table 5.1.

Notice that while JQR is tailored to predict finite many quantiles, our IQR method estimates the whole quantile function hence solves a more challenging task. Despite the more difficult problem solved, as Table 5.1 suggest the performance in terms of pinball loss of IQR is comparable to that of the state-of-the-art JQR on all the twenty studied benchmarks, except for the 'crabs' and 'cpus' datasets (pval < 0.25%). In addition, when considering the non-crossing penaltycone can observe that IQR outperforms the IND-QR baseline on eleven datasets (pval < 0.25%) and JQR on two datasets. This illustrates the efficiency of the constraint based on the continuum scheme.

Smoothing the Loss Function

There are several ways to solve the non-smooth optimization problems associated to the QR task. One could proceed for example by duality-as presented in Section 3.4.3-, or apply sub-gradient descent techniques (which often converge quite slowly). In our experiments we used the L-BFGS-B [START_REF] Zhu | Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization[END_REF]) optimization scheme which is widely popular in machine learning, with non-smooth extensions [START_REF] Skajaa | Limited memory BFGS for nonsmooth optimization[END_REF][START_REF] Keskar | A limited-memory quasi-Newton algorithm for boundconstrained non-smooth optimization[END_REF]. The technique requires only evaluation of objective function along Notice that the non-crossing penalty does not provide crossings to occur in the regions where there is no points to enforce the penalty (e.g. x ∈ [0.13, 0.35]). This phenomenon is alleviated by the regularity of the model. with its gradient, which can be computed automatically using reverse mode automatic differentiation (as in [START_REF] Abadi | Tensorflow: Large-scale machine learning on heterogeneous distributed systems[END_REF]). To benefit from the available fast smooth implementations [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF][START_REF] Fei | Parallel L-BFGS-B algorithm on GPU[END_REF], we smoothed the non-differentiable part of the loss function. To do that, remark that the pinball loss writes as

ℓ(θ, u, v) = θ -R -(v -u) |v -u|
and replace the absolute value by a smoothed approximation given by its Moreau enveloppe [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF]. The resulting κ-smoothed (κ > 0) absolute value is as follows:

ψ κ 1 (r) := |•| □ 1 2κ |•| 2 (r) =    1 2κ r 2 if |r| ≤ κ |r| -κ 2 otherwise,
Minimizing the κ-smoothed pinball loss

ℓ κ (θ, h(x), y) = θ -R -(y -h(x)) ψ κ 1 (y -h(x))
yields the quantiles when κ → 0, the expectiles as κ → +∞. The intermediate values are known as M-quantiles [START_REF] Breckling | M-quantiles[END_REF]. The influence of the κ parameter is illustrated in Figure 5.3. For this experiment, 10000 samples have been generated from the sine wave dataset described in Section 5.4.1, and the model have been trained on 100 quantiles generated from a Gauss-Legendre quadrature. When κ is small the quantiles are recovered (the dashed lines on the right plot match the theoretical quantiles in plain lines). It took 225s (258 iteration, and 289 function evaluations) to train for κ = 1 • 10 1 , 1313s for κ = 1 • 10 -1 (1438 iterations and 1571 function evaluations), 931s for κ = 1e -3 (1169 iterations and 1271 function evaluations) and 879s for κ = 0 (1125 iterations and 1207 function evaluations). We used a GPU Tensorflow implementation and run the experiments in float64 on a computer equipped with a GTX 1070, and intel i7 7700 and 16GB of DRAM.

Influence of the Number of Sampled Locations

In this section, we quickly present an additional experiment on the number of sampled locations m used to approximate the integral loss I ℓ .

In the experiment presented in Figure 5.4, on the sine synthetic benchmark, we draw n = 1000 training points and study the impact of increasing m on the quality of the quantiles at θ ∈ {0.05, 0.25, 0.5, 0.75, 0.95}. We notice that when m ≥ 34 ≈ √ 1000 there is little benefit to draw more m samples: the quantile curves do not change on the n test = 2000 test points.

Deep Kernel Learning for Quantile Regression on Images

In this section, we apply the IQR method in the problem of estimating the age of fishes (Y) from images (X). This task is critical in sustainable management of fish resources with the goal of avoiding overfishing. The age of a fish can be inferred from its otholiths, which are calcium carbonate structures growing behind their brains. As the fish gets older, the otholith becomes larger and more convoluted, as illustrated in Figure 5.5. The estimation problem is rather challenging, and typically involves experts who carefully examine the growth zones of the otholiths. This process is costly, especially when we realize that the number of otholiths analysed yearly around the world scales in millions [START_REF] Campana | Otoliths, increments, and elements: keys to a comprehensive understanding of fish populations?[END_REF]. Automating the task is of fundamental interest, and has been addressed recently using convolutional neural networks (CNNs) in [START_REF] Moen | Automatic interpretation of otoliths using deep learning[END_REF], yielding satisfactory predictions and interpretable results [START_REF] Ordonez | Explaining decisions of deep neural networks used for fish age prediction[END_REF].

Beyond the single real-valued prediction of the age of the fish, we propose to learn the conditional quantile function of Y (age of the fish) given X (photography of its otholith). In the ITL framework, this requires to define a kernel on the input space, which can be difficult. As CNN are known to learn relevant features stored in their last layer [START_REF] Goodfellow | Deep learning[END_REF], a first idea would be to use a pretrained CNN to extract features for X, and then apply the ITL framework with an input kernel on the resulting feature space, denoted V in what follows. This amounts to using a kernel k X (x, y) = k V (ϕ ω (x), ϕ ω (y)),

(5.30) where k V is a kernel on the feature space, and ϕ ω : X → V is the neural architecture used for extracting features parameterized by its weights ω. In a deep kernel learning spirit, we propose here to learn jointly the parameters of the neural architecture and that of the estimator. To do so, we consider an input kernel of the form Equation (5.30), where ϕ ω is an Inception v3 architecture [START_REF] Szegedy | Rethinking the inception architecture for computer vision[END_REF], and k V is a Gaussian random Fourier feature (RFF) kernel. We also choose k Θ to be a Gaussian RFF kernel, which allows to exploit Algorithm 3.2 for the optimization. To train jointly on both parameters of the kernel and the estimator, at each epoch we make a gradient step on the kernel parameters, starting from pre-trained weights used in Ordonez et al. (2020) for predicting fish ages. Many tricks are needed to make the training smooth, such as rescaling images, feature augmentation, and so on; we refer to Ordonez et al. (2020) for details. While visualization is a bit hard for this problem, we illustrate in Figure 5.6 a few quantile functions associated to random samples from the dataset. We can see that the quantile functions are nondecreasing with respect to θ, which is at in agreement with theoretical properties of quantiles.

Cost-Sensitive Classification

In this section we present some numerical illustration to the CSC application of the ITL scheme. The optimization is performed using a L-BFGS solver on the finite dimensional problem provided by the double representer expression from Theorem 3.9, associated to a smoothing of the imbalanced Hinge loss obtained by its Moreau enveloppe, similarly to the quantile regression case developed in Section 5.4.2.

We used the Iris UCI dataset with 4 attributes and 150 samples, the two synthetic scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] datasets Two-Moons (noise=0.4) and Circles (noise=0.1) with both 2 attributes and 1000 samples. As detailed in section 5.2.3, CSC on a continuum Θ = [-1, 1] that we call infinite cost-sensitive classification (ICSC) can be tackled by our proposed technique. In this case, the hyperparameter θ controls the tradeoff between the importance of the correct classification with labels -1 and +1. When θ = -1, class -1 is emphasized; the probability of correctly classified instances with this label (called specificity) is desired to be 1. Similarly, for θ = +1, the probability of correct classification of samples with label +1 (called sensitivity) is ideally 1. We chose k X to be a Gaussian kernel with bandwidth σ X = (2γ X ) (-1/2) the median of the Euclidean pairwise distances of the input points (Jaakkola et al., 1999). The output kernel k Θ is also a Gaussian kernel with bandwidth γ Θ = 5. We used m = 20 for all datasets in the double representer theorem associated to the sobol sequence.

As a baseline we trained independently 3 CSC classifiers with θ ∈ {-0.9, 0, 0.9}. We repeated 50 times a random 50 -50% train-test split of the dataset and report the average test error and standard deviation (in terms of sensitivity and specificity).

Our results are illustrated in Table 5.2. For θ = -0.9, both independent and joint learners give the desired 100% specificity; the joint CSC scheme however has significantly higher sensitivity value (15% vs 0%) on the dataset Circles. Similar conclusion holds for the θ = +0.9 extreme: the ideal sensitivity is reached by both techniques, but the joint learning scheme performs better in terms of specificity (0% vs 12%) on the dataset Circles.

The results from Table 5.2 are highlighted in Figure 5.7 and Figure 5.8 where we present a vizualization of the independent vs continuum learning problems.

Density Level Set Estimation

In this section, we adapt the ITL framework to an unsupervised learning problem, the DLSE carried out by means of OCSVM. As a reminder, given some data (x i ) n i=1 following an unknown law P X and a hyperparameter value θ ∈ (0, 1], the classical OCSVM in a scalar RKHS H k X can be formulated as ĥ, t = arg min

h∈H k X ,t∈R 1 nθ n i=1 max (0, t -h(x i )) -t + 1 2 h 2 H k X .
(5.31)

Given the resulting estimator ĥ, t ∈ H k X × R, the rule d(x) := sign ĥ(x) -t is then a reasonable choice for deciding whether a new point x is an inlier (d(x) = 1) or an outlier (d(x) = -1). In particular, it is known (Schölkopf et al., 2001b) that if the distribution P X does not contain discrete components, and the kernel k X is analytic and non-constant, then this approach is theoretically grounded as the decision function d will separate the input space X in two sets (inliers and outliers) with respective mass 1 -θ and θ, a property refered to as the θ-property.

Our objective is to generalize this algorithm to learning an OCSVM jointly for the continuum θ ∈ Θ := (0, 1].

Ω(h) = 1 2 Θ h(•)(θ) 2 H k X dµ(θ), 5
K ,t∈H k b 1 n n i=1 I ℓ (h(x i ), t) + Ω(h) + λ 2 t 2 H k b , (5.32)
where λ > 0 is a regularization parameter, and I ℓ is the integral loss associated to

ℓ(θ, u, v) = 1 θ max (0, v -u) -v
and to a probability measure µ. The decision function is then able to predict the outlierness of a new data point x for all θ ∈ (0, 1] by d(x)(θ) = sign ĥ(x)(θ) -t(θ) .

Remark 5.28. The regularizer Ω(h) is uniformly weaker than a vv-RKHS norm based regularizer

1 2 h 2 H K : ∀θ ∈ Θ, h(•)(θ) 2 H k X = h(•)(θ), h(•)(θ) H k X ≤ k Θ (θ, θ) h, h H K so that after integration Ω(h) ≤ 1 2 Θ k(θ, θ)dµ(θ) h 2 H K .
This prevents the direct use of a representer theorem, as the regularizer does not write as a nondecreasing function of h H K .

Remark 5.29. Due to the non-integrability of 1 θ at a neighborhood of 0, the integral in I ℓ is not assured to converge to a finite value. To alleviate this bottleneck, one can choose the measure µ appropriately and consider either a measure with support [ϵ, 1] for some small ϵ > 0 or a measure with full support that ensures proper convergence of I ℓ (e.g. the measure of density θ → 2θ).

We now devote the following part to solving Problem 5.32, that we refer to as infinite one-class SVM (IOCSVM). In Section 5.6.1, we propose a representer theorem for a sampled version of Problem 5.32 that is amenable to optimization and emphasize the benefit of this approach with numerical experiments in Section 5.6.2.

Representer Theorem for Mixed Regularization

One difficulty of working with Problem 5.32 is that the classical representer theorem from Theorem 2.43 does not apply. To get a finite parametrization of the estimator, we will solve a sampled version of the problem, and prove a double representer theorem for ĥ. Given a sampling scheme (η j , θ j ) m j=1 ∈ (R×Θ) m using for instance quasi Monte-Carlo or quadratures (see Section 3.3.1), we introduce the notation .

I ℓ (h(x), t) := m j=1 η j ℓ(θ j , h(x i )(θ j ), t(θ j )), Ω(h) := 1 2 m j=1 η j h(•)(θ j ) 2 H k X (5.
(θ, θ ′ ) → k Θ (θ, θ ′ ) Id H k X . Then for all m ∈ N * and (θ j ) m j=1 ∈ Θ m , + m j=1 Im(K θ j ) ⊕ ∩ m j=1 Ker(K ♯ θ j ) = H K , ( 5
Proof The statement boils down to proving that

V := + m j=1 Im(K θ j ) is closed in H K , since it is straightforward by double inclusion that V ⊥ = ∩ m j=1 Ker K ♯ θ j . Let (e j ) l j=1 be an orthonormal basis of Span {(k Θ (•, θ j )) m j=1 } ⊂ H k Θ , with l = m potentially. Such basis can be obtained by applying the Gram-Schmidt orthonormalization method to (k Θ (•, θ j )) m j=1 . Then, V = Span {e j ⊗ f, 1 ≤ j ≤ k, f ∈ H k X }. Notice also that for ∀(j 1 , j 2 ) ∈ [l] × [l], ∀f, g ∈ H k X , e j 1 ⊗ f, e j 2 ⊗ g H K = e j 1 , e j 2 H k Θ f, g H k X .
(5.36)

Let (h n ) n∈N * be a sequence in V converging to some h ∈ H K . By definition, one can find sequences (f 1,n ) n∈N * , . . . , (f l,n ) n∈N * ∈ H k X such that for ∀n ∈ N * , h n = l j=1 e j ⊗ f j,n . Let p, q ∈ N * . It holds that, using the orthonormal property of (e j ) l j=1 and Equation (5.36),

h p -h q 2 H K = l j=1 e j (f j,p -f j,q ) 2 H K = l j=1 f j,p -f j,q 2 H k X . Since (h n ) n∈N * converges in H K , it
is a Cauchy sequence, thus so are the sequences (f j,n ) n∈N * . But H k X is a complete space, so these sequences converge in H k X , and by denoting f j = lim n→∞ f j,n , one gets h = l j=1 e j ⊗ f j . Therefore h ∈ V, V is closed and the orthogonal decomposition in Equation (5.35) holds.

Lemma 5.31. Let k X , k Θ be two scalar-valued kernels, with k Θ a strictly positive definite kernel (i.e. any Gram matrix associated to k Θ is strictly positive), and

K : (θ, θ ′ ) → k Θ (θ, θ ′ )I H k X . Let also m ∈ N * , (θ j ) m j=1 ∈ Θ m , and V = + m j=1 Im(K θ j ) . Then I : V → R defined as I(h) := m j=1 h(θ j ) 2 H k X
is γ-strongly convex where γ is the smallest eigenvalue of the Gram matrix K Θ associated to (θ j ) m j=1 and k Θ . In particular, I is coercive.

Proof Notice that I is the quadratic form on V associated to the linear mapping L ∈ L(H K ) defined by L := m j=1 K θ j K ♯ θ j , in the sense that for all h ∈ V, I

(h) = h, Lh H K . Indeed, ∀h ∈ V, I(h) = m j=1 K ♯ θ j h, K ♯ θ j h H k X = m j=1 h, K θ j K ♯ θ j h H K = h, Lh H K . Let h ∈ V, by definition there exist (f j ) m j=1 ∈ H m k X such that h = m j=1 K θ j f j . We begin by computing h 2 H K : h, h H K = m i=1 K θ i f i , m j=1 K θ j f j H K = m i=1 m j=1 K θ i f i , K θ j f j H K = m i=1 m j=1 K ♯ θ j K θ i f i , f j H k X = m i=1 m j=1 k Θ (θ i , θ j ) f i , f j H k X = K Θ , K f F , = Tr K Θ K f ,
where we introduced the notation K Θ ∈ M m (R) for the Gram matrix over the points (θ j ) m j=1 with the kernel k Θ and K f ∈ M m (R) for the Gram matrix over the (f j ) m j=1 with the linear kernel. We compute I(h) similarly, starting by

Lh = m j=1 K θ j K ♯ θ j   m i=1 K θ i f i   = m i=1 m j=1 k Θ (θ i , θ j )K θ j f i ,
and following by

h, Lh H K = m l=1 K θ l f l , m i=1 m j=1 k Θ (θ i , θ j )K θ j f i H K = m l=1 m i=1 m j=1 k Θ (θ i , θ j ) K θ l f l , K θ j f i H K = m l=1 m i=1 m j=1 k Θ (θ i , θ j ) K ♯ θ j K θ l f l , f i H k X = m l=1 m i=1 m j=1 k Θ (θ i , θ j )k Θ (θ j , θ l ) f l , f i H k X = K Θ , K f K Θ F = Tr K Θ K Θ K f .
By the strictly positive assumption on k Θ , K Θ has a minimum eigenvalue γ > 0. Then K Θ -γ Id m and K Θ K f are positive symetric matrices and

Tr K Θ -γ Id m K Θ K f ≥ 0 (5.37)
which can be rewritten as

Tr K Θ K Θ K f ≥ γ Tr K Θ K f or I(h) ≥ γ h 2 H K .
Finally, we recognize in the left term of Equation (5.37) a convex function of h, so

I -γ • 2 H K is convex
(x, θ) ∈ X × Θ ĥ(x)(θ) = n i=1 m j=1 αij k X (x, x i )k Θ (θ, θ j ), t(θ) = m j=1 βj k b (θ, θ j ) for some (α ij ) n,m i,j=1 ∈ R n×m and ( βj ) m j=1 ∈ R m . Moreover the solution is unique on V × H k b with V defined in Lemma 5.31. Proof Introduce K ′ : (θ, θ ′ ) ∈ Θ × Θ → k Θ (θ, θ ′ )I H k X ∈ L(H k X ).
The spaces H K and H K ′ are unitarily equivalent by means of the feature operator W :

H K ′ → H K such that ∀(f, x, θ) ∈ H K ′ × X × Θ, W f (x)(θ) = f (θ)(x)
. This means that they are essentially the same spaces of functions, whose elements are seen either as mappings

X → (Θ → R) or Θ → (X → R). Define J :   H K × H k b → R (h, t) → 1 n n i=1 I ℓ (h(x i ), t) + Ω(h) + λ 2 t 2 H k b   . Let V = W + m j=1 Im(K ′ θ j ) . Since W is an isometry, thanks to Equation (5.35), it holds that V ⊕ V ⊥ = H K . Let (h, t) ∈ H K × H k b , there exists unique h V ⊥ ∈ V ⊥ , h V ∈ V such that h = h V + h V ⊥ . Notice that J(h, t) = J(h V + h V ⊥ , t) = J(h V , t) since ∀j ∈ [m], ∀x ∈ X, h V ⊥ (x)(θ j ) = W -1 h V ⊥ (θ j )(x) = 0.
Moreover, J is bounded by below so its infimum is well-defined, and inf

(h,t)∈H K ×H k b J(h, t) = inf (h,t)∈V×H k b J(h, t).
Finally, notice that J is coercive on V × H k b endowed with the sum of the norms in V and H k b (which makes it a Hilbert space): if

(h n , t n ) n∈N * ∈ V × H k b is such that h n H K + t n H k b → n→∞ +∞, then either ( h n H K ) n∈N or ( t n H k b
) n∈N has to diverge :

• If t n H k b → n→∞ +∞, since t n (θ j ) = t n , k b (•, θ j ) H k b ≤ k b (θ j , θ j ) t n H k b ≤ κ b t n H k b (∀j ∈ [m]), then J(h n , t n ) ≥ λ 2 t n 2 H k b - m j=1 w j t(θ j ) → n→∞ +∞. 5.7. CONCLUSION 113 • If h n H K → n→∞ +∞, according to Lemma 5.31, J(h n , t n ) → n→∞
+∞ as long as all η j are strictly positive, which is verified for the quasi Monte-Carlo scheme.

Thus J is coercive, so Bauschke et al. (2011, Proposition 11.15) allows to conclude that J has a minimizer ( ĥ, t) on V × H k b . Then, in the same fashion as Theorem 3.9, define

U 1 = Span {(K(•, x i )k Θ (•, θ j )) n,m i,j=1 } ⊂ V and U 2 = Span {(k b (•, θ j )) m j=1 } ⊂ H k b ,
and use the reproducing property to show that ( ĥ, t) ∈ U 1 × U 2 , so there exist (α ij ) n,m i,j=1 and (β j ) m j=1 such that for ∀(x, θ) ∈ X × Θ, ĥ

(x)(θ) = n i=1 m j=1 αij k X (x, x i )k Θ (θ, θ j ), t(θ) = m j=1 βj k b (θ, θ j ).
Finally, using the γ notation from Lemma 5.31, J is min (λ, γ)-strongly convex on V × H k b which ensures that the resulting ( ĥ, t) is unique.

Theorem 5.33 guarantees a finite-dimensional expansion for the estimator ( ĥ, t), which can be used as a plug-in in Problem 5.34. The next section highlights the performance of this approach on some synthetic and real-world datasets.

Numerical Experiments

To assess the quality of the estimated model by IOCSVM, we illustrate the θ-property [START_REF] Schölkopf | New support vector algorithms[END_REF]: the proportion of inliers has to be approximately 1 -θ (∀θ ∈ (0, 1)). For the studied datasets (Wilt, Spambase) we used the raw inputs without applying any preprocessing. After experimenting with a few input kernels, we settled for the exponentiated χ 2 kernel k X (x, z) := exp -γ X d k=1 (x k -z k ) 2 /(x k + z k ) with bandwidth γ X = 0.25. A Gauss-Legendre quadrature rule provided the integral approximation in Equation (5.33), with m = 100 samples. We chose the Gaussian kernel for k Θ ; its bandwidth parameter γ Θ was the 0.2-quantile of the pairwise Euclidean distances between the θ j 's obtained via the quadrature rule. The margin (bias) kernel was k b = k Θ . As it can be seen in Fig. 5.9, the θ-property holds for the estimate which illustrates the efficiency of the proposed continuum approach for density level-set estimation. In the scalar case, the OCSVM algorithm with the Gaussian kernel in input estimates the level sets of the density of X. Combined with the θ-property, this means that the IOCSVM estimator should be able to estimate these level sets as a function of θ. We show in Figure 5.10 that this is the case on a simple Gaussian mixture underlying probability P X . The strength of the approach is then to be able to provide a score of outlierness for a new point x ∈ X:

s(x) := sup θ∈Θ θ : d(x)(θ) ≥ 0 .
This is well-suited to anomaly detection problems as it indicates whether a point x ∈ X is a rare event (tail of the distribution, s(x) close to 0) or a common event (s(x) close to 1).

Conclusion

In this chapter we proposed infinite task learning (ITL), a novel nonparametric framework aiming at jointly solving parametrized tasks for a continuum of hyperparameters.

Introduction

Recent years have witnessed an increasing attention around style transfer problems (Gatys et [START_REF] Fu | Style transfer in text: Exploration and evaluation[END_REF] as well as audio signal processing [START_REF] Grinstein | Audio style transfer[END_REF] where objects at hand are contents in which style is inherently part of their perception. Style transfer is one of the key components of data augmentation [START_REF] Mikołajczyk | Data augmentation for improving deep learning in image classification problem[END_REF] as a means to artificially generate meaningful additional data for the training of deep neural networks. Besides, it has also been shown to be useful for counterbalancing bias in data by producing stylized contents with a well-chosen style (see for instance [START_REF] Geirhos | ImageNet-trained CNNs are biased towards texture; increasing shape bias improves accuracy and robustness[END_REF]) in image recognition. More broadly, style transfer fits into the wide paradigm of parametric modeling, where a system, a process or a signal can be controlled by its parameter value. Adopting this perspective, style transfer-like applications can also be found in digital twinning [START_REF] Tao | Digital twin in industry: State-of-the-art[END_REF][START_REF] Barricelli | A survey on digital twin: Definitions, characteristics, applications, and design implications[END_REF][START_REF] Lim | A state-of-the-art survey of digital twin: techniques, engineering product lifecycle management and business innovation perspectives[END_REF], a field of growing interest in health and industry.

In this chapter, we propose a novel principled approach for style transfer, exemplified in the context of emotion transfer of face images. Given a set of emotions, classical emotion transfer refers to the task of transforming face images according to these target emotions. The pioneering works in emotion transfer include that of [START_REF] Blanz | A morphable model for the synthesis of 3D faces[END_REF] who proposed a morphable 3D face model whose parameters could be modified for facial attribute editing. [START_REF] Susskind | Generating facial expressions with deep belief nets[END_REF] designed a deep belief net for facial expression generation using action unit (AU) annotations.

More recently, extensions of generative adversarial networks (GANs, [START_REF] Goodfellow | Generative adversarial nets[END_REF]) have proven to be particularly powerful for tackling image-to-image translation problems [START_REF] Zhu | Unpaired image-to-image translation using cycle-consistent adversarial networks[END_REF]. Several works have addressed emotion transfer for facial images by conditioning GANs on a variety of guiding information ranging from discrete emotion labels to photos and videos. In particular, StarGAN [START_REF] Choi | StarGAN: Unified generative adversarial networks for multi-domain image-to-image translation[END_REF]) is conditioned on discrete expression labels for face synthesis. ExprGAN [START_REF] Ding | ExprGAN: Facial expression editing with controllable expression intensity[END_REF] proposes synthesis with the ability to control expression intensity through a controller module conditioned on discrete labels. Other GAN-based approaches make use of additional information such as AU labels [START_REF] Pumarola | GANimation: Anatomically-aware facial animation from a single image[END_REF], target landmarks [START_REF] Qiao | Geometry-contrastive GAN for facial expression transfer[END_REF], fiducial points [START_REF] Song | Geometry guided adversarial facial expression synthesis[END_REF] and photos/videos [START_REF] Geng | Warp-guided GANs for singlephoto facial animation[END_REF]. While GANs have achieved high quality image synthesis, they come with some pitfalls: they are particularly difficult to train and require large amounts of training data.

Unlike previous approaches, we adopt a functional point of view: given some person, we assume that the full range of the emotional faces can be modelled as a continuous function from emotions to images. This view exploits the geometry of the representation of emotions [START_REF] Russell | A circumplex model of affect[END_REF], assuming that one can pass a facial image "continuously" from one emotion to an other. We then propose to address the problem of emotion transfer by learning an image-to-function model able to predict for a given facial input image represented by its landmarks [START_REF] Tautkute | I know how you feel: Emotion recognition with facial landmarks[END_REF], the continuous function that maps an emotion to the image transformed by this emotion. This function-valued regression approach relies on a the infinite task learning (ITL) technique developed in Chapter 5. ITL enlarges the scope of multi-task learning (Evgeniou and Pontil, 2004; [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF]) by learning to solve simultaneously a set of tasks parametrized by a continuous parameter. While strongly linked to other parametric learning methods such the one proposed by [START_REF] Takeuchi | Nonparametric quantile estimation[END_REF], the approach differs from previous works by leveraging the use of operator-valued kernels and vector-valued reproducing kernel Hilbert spaces (vv-RKHS; Pedrick 1957; [START_REF] Micchelli | On learning vector-valued functions[END_REF][START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces of integrable functions and Mercer theorem[END_REF]) to model function-valued functions with vectorial outputs. vv-RKHSs have proven to be relevant in solving supervised learning tasks such as multiple quantile regression [START_REF] Sangnier | Joint quantile regression in vector-valued RKHSs[END_REF] or unsupervised problems like anomaly detection (Schölkopf et al., 2001b). A common property of these works is that the output to be predicted is a real-valued function of a real parameter.

To solve the emotion transfer problem, we present an extension of ITL, vector ITL (or shortly vITL) which involves functional outputs with vectorial representation of the faces and the emotions, showing that the approach is still easily controllable by the choice of appropriate kernels guaranteeing continuity and smoothness. In particular, the functional point of view by the inherent regularization induced by the kernel makes the approach suitable even for limited and partially observed emotional images. We demonstrate the efficiency of the vITL approach in a series of numerical experiments showing that it can achieve state-of-the-art performance on two benchmark datasets.

The chapter is structured as follows. We formulate the problem and introduce the vITL framework in Section 6.2. Section 6.3 is dedicated to the underlying optimization problem. Numerical experiments conducted on two benchmarks of the domain are presented in Section 6.4. Discussion and future work conclude the chapter in Section 6.5.

Problem Setting

In this section we define our problem. Our aim is to design a system capable of transferring emotions: having access to the face image of a given person our goal is to convert his/her face to a specified target emotion. In other words, the system should implement a mapping of the form (face, emotion) → face.

(6.1)

In order to tackle this task, one requires a representation of the emotions, and similarly that of the faces. The classical categorical description of emotions deals with the classes 'happy', 'sad', 'angry', 'surprised', 'disgusted', 'fearful'. The valence-arousal model [START_REF] Russell | A circumplex model of affect[END_REF]) embeds these categories into the 2-dimensional space. The resulting representation of the emotions are points θ ∈ R 2 , each coordinate of these vectors encoding the valence (pleasure to displeasure) and arousal (high to low) associated to the emotions. This is the emotion representation we use while noting that there are alternative encodings in higher dimension (Θ ⊂ R p , p ≥ 2; The resulting mapping from Equation (6.1) is illustrated in Figure 6.1: starting from a neutral face and the target happy one can traverse to the happy face; from the happy face, given the target emotion surprise one can get to the surprised face.

In an ideal world, for each person, one would have access to a trajectory z mapping each emotion θ ∈ Θ to the corresponding landmark locations x ∈ X; this function z : Θ → X can be taken for instance to be the element of L 2 [Θ, µ; X], the space of R d -valued square-integrable function w.r.t. to a measure µ. The probability measure µ allows capturing the frequency of the individual emotions. In practice, one has realizations (z i ) i∈ [n] , each z i corresponds to a single person possible appearing multiple times. The . 1 In order to capture the relation from Equation (6.1) one can rely on a hypothesis space H with elements h : X → (Θ → X). (6.

2)

The value h(x)(θ) represents the landmark prediction from face x and target emotion θ.

We consider two tasks for emotion transfer:

• Single emotional input: In the first problem, the assumption is that all the faces appearing as the input in Equation (6.1) come from a fixed emotion θ0 ∈ Θ.

The data which can be used to learn the mapping h consists of t = n triplets2

x i = z i θ0 ∈ X, Y i = z i θi,j

=:y i,j j∈[m] ∈ X m , (θ i,j ) j∈[m] = θi,j j∈[m] ∈ Θ m , i ∈ [t].
To measure the quality of the reconstruction using a function h, one can consider a convex loss ℓ : X × X → R on the landmark space. The resulting objective function to minimize is

R S (h) := 1 tm i∈[t] j∈[m]
ℓ(h(x i )(θ i,j ), y i,j ). (

The risk R S (h) captures how well the function h reconstructs on average the landmarks y i,j when applied to the input landmark locations x i .

• Joint emotional input: In this problem, the faces appearing as input in Equation (6.1) can arise from any emotion. The observations consist of triplets

x m(i-1)+l = z i θi,l ∈ X, Y m(i-1)+l = z i θi,j which we detail below. A general reminder on these spaces of functions is presented in Chapter 2. The vector (R d )-valued capability is beneficial to handle the Θ → X = R d mapping; the associated R d -valued RKHS H G is uniquely determined by a matrixvalued kernel G : Θ × Θ → R d×d = L(X). Similarly, in Equation (6.4) the X → H G mapping is modelled by a vv-RKHS H K corresponding to an operator-valued kernel K : X × X → L(H G ). We chose K and G to be separable kernels of the form

G(θ, θ ′ ) = k Θ (θ, θ ′ )A, K(x, x ′ ) = k X (x, x ′ ) Id H G (6.5)
with a scalar-valued kernel k X : X × X → R and k Θ : Θ × Θ → R, and symmetric, positive definite matrix A ∈ R d×d . This choice corresponds to the intuition that for similar input landmarks and target emotions, the predicted output landmarks should also be similar, as measured by k X , k Θ and A, respectively. More precisely, smoothness (analytic property) of the emotion-to-landmark output function can be induced for instance by choosing a Gaussian kernel k Θ (θ, θ ′ ) = exp(-γ θ -θ ′ 2 2 ) with γ > 0. The matrix A when chosen as A = Id d corresponds to independent landmarks coordinates while other choices encode prior knowledge about the dependency among the landmarks coordinates (Álvarez et al., 2012). Similarly, the smoothness of function h can be driven by the choice of a Gaussian kernel over X while the identity operator on H G is the simplest choice to cope with functional outputs. By denoting the norm in H K as

• H K , the final objective function is

min h∈H K R S (h) + λ 2 h 2 H K (6.6)
with a regularization parameter λ > 0 which balances between the data-fitting term (R S (h)) and smoothness ( h

2 H K
). We refer to Problem 6.6 as vector-valued infinite task learning (vITL).

Remark: This problem is a natural adaptation of the ITL framework [START_REF] Brault | Infinite task learning in RKHSs[END_REF] learning with operator-valued kernels mappings of the form X → (Θ → Y) where Y is a subset of R; here Y = X. An other difference is µ: in ITL this probability measure is designed to approximate integrals via quadrature rule, in vITL it captures the observation mechanism.

Optimization

This section is dedicated to the solution of Problem 6.6 which is an optimization problem over functions (h ∈ H K ). General tools to handle problems such as Problem 6.6 have been proposed in Chapter 3, in particular throughout this chapter we consider the squared loss ℓ(x, x ′ ) = 1 2 x -x ′ 2 2 ; in this case we recognize the partially observed scenario with square loss from Section 3.2.2. Thus, Lemma 3.7 resumes solving Problem 6.6 to solving a Sylvester equation, as it holds that the solution ĥ ∈ H K writes ĥ(x)(θ) = t i=1 m j=1 k X (x, x i )k Θ (θ, θ i,j )Aα i,j , ∀(x, θ) ∈ X × Θ (6.7)

for some (α i,j ) • Regularization in vv-RKHS: Using the notations above, for any h ∈ H K parameterized by a matrix α, it holds that h 2 H K = Tr KαAα ⊤ . Given two matrices A 1 , A 2 and associated vv-RKHSs H K 1 and H K 2 , if A 1 and A 2 are invertible then any function in H K 1 parameterized by α also belongs to H K 2 (and vice versa), within which it is parameterized by αA -1 2 A 1 . This means that the two spaces contain the same functions, but their norms are different.

Experiments

In this section we demonstrate the efficiency of the proposed vITL technique in emotion transfer. We first introduce the two benchmark datasets we used in our experiments and give details about data representation and choice of the hypothesis space in Section 6.4.1. Then, in Section 6.4.2, we provide a quantitative performance assessment of the vITL approach (in mean squared error and classification accuracy sense) with a comparison to the state-of-the-art StarGAN method. Section 6.4.3 is dedicated to investigation of the role of A (see Equation (6.5)) and the robustness of the approach w.r.t. missing obervation. These two sets of experiments (Section 6.4.2 and Section 6.4.3) are augmented with a qualitative analysis (Section 6.4.4). The code written for all these experiments is available on GitHub.

Experimental Setup

We used the following two popular face datasets for evaluation. In our experiments, we used frontal images and seven emotions from each of these datasets. An edge map illustration of landmarks for different emotions is shown in Figure 6.2.

At this point, it is worth recalling that we are learning a function-valued function, h : X → (Θ → X) using a vv-RKHS as our hypothesis class (see Section 6.2). In the following we detail the choices made concerning the representation of the landmarks in X, that of the emotions in Θ, and in the kernel design k X , k Θ and A.

Landmark representation, pre-processing: We applied the following pre-processing steps to get the landmark representations which form the input of the algorithms. To extract 68 landmark points for all the facial images, we used the standard dlib library. The estimator is based on dlib's implementation of [START_REF] Kazemi | One millisecond face alignment with an ensemble of regression trees[END_REF], trained on the iBUG 300-W face landmark dataset. Each landmark is represented by its 2D location. The alignment of the faces was carried out by the Python library imutils.

The method ensures that faces across all identities and emotions are vertical, centered and of similar sizes. In essence, this is implemented through an affine transformation computed after drawing a line segment between the estimated eye centers. Each image was resized to the size 128 × 128. The landmark points computed in the step above were transformed through the same affine transformation. These two preprocessing steps gave rise to the aligned, scaled and vectorized landmarks x ∈ R 136=2×68 .

for 10 epochs with a batch size of 16.

We report the mean and standard deviation of the aforementioned metrics over ten 90%-10% train-test splits of the data. The test set for each split is constructed by removing 10% of the identities from the data. For each split, the best γ X , γ Θ and λ values were determined by 6-fold and 10-fold cross-validation on KDEF and RaFD, respectively.

Baseline: We used the popular StarGAN [START_REF] Choi | StarGAN: Unified generative adversarial networks for multi-domain image-to-image translation[END_REF] system as our baseline.

Other GAN-based studies use additional information and are not directly comparable to our setting. For fair comparison, the generator G and discriminator D were modified to be fully-connected networks that take vectorized landmarks as input. In particular, G was an encoder-decoder architecture where the target emotion, represented as a 2D emotion encoding as for our case, was appended at the bottleneck layer. It contained approximately one million parameters, which was chosen to be comparable with the number of coefficients in vITL (839, 664 = 126×7×7×136 for KDEF). ReLU activation function was used in all layers except before bottleneck in G and before penultimate layers of both G and D. We used their default parameter values in the code.3 Experiments over each split of KDEF and RaFD were run for 50K and 25K iterations, respectively.

MSE results:

The test MSE for the compared systems is summarized in Table 6.1. As the table shows, the vITL technique outperforms StarGAN on both datasets. One can observe low reconstruction cost for vITL in both the single and the joint emotional input case. Interestingly, a performance gain is obtained with vITL-joint on the RaFD data in MSE sense. We hypothesize that this is due to the joint model benefiting from input landmarks for other emotions in the small data regime (only 67 samples per emotion for RaFD). Despite our best efforts, we found it quite difficult to train StarGAN reliably and the diversity of its outputs was low.

Classification results:

The emotion classification accuracies are available in Table 6.2. The classification results clearly demonstrate the improved performance and the higher quality of the generated emotion of vITL over StarGAN; the latter also produces predictions with visible face distortions as it is illustrated in Section 6.4.4. To provide further insight into the classification performance we also show the confusion matrices for the joint vITL model on a particular split of KDEF and RaFD datasets in Figure 6.4. For both the datasets, the classes 'happy' and 'surprised' are easiest to detect. Some confusions arise between the classes 'neutral' vs 'sad' and 'fearful' vs 'surprised'. Such mistakes are expected when only using landmark locations for recognizing emotions.

Analysis of Additional Properties of vITL

This section is dedicated to the effect of the choice of A (in kernel G) and to the robustness of vITL w.r.t. missing obervation.

Influence of A in the matrix-valued kernel G: Here, we illustrate the effect of matrix A (see Equation (6.5)) on the vITL estimator and show that a good choice of A can lead to lower dimensional models, while preserving the quality of the prediction. The choice of A is built on the knowledge that the empirical covariance matrices of the output training data contains structural information that can be exploited with vv-RKHS [START_REF] Kadri | A generalized kernel approach to structured output learning[END_REF]. In order to investigate this possibility, we performed the These experiments and qualitative results demonstrate the efficiency of the vITL approach in emotion transfer.

Conclusion

In this chapter we introduced a novel approach to style transfer based on functionvalued regression, and exemplified it on the problem of emotion transfer. The proposed vector-valued infinite task learning (vITL) framework relies on operator-valued kernels. vITL (i) is capable of encoding and controlling continuous style spaces, (ii) benefit from a representer theorem for efficient computation, and (iii) facilitates regularity control via the choice of the underlying kernels. The framework can be extended in several directions. Other losses can be leveraged to produce outlier-robust or sparse models in the spirit of Chapter 4. Instead of being chosen prior to learning, the input kernel could be learned using deep architectures [START_REF] Mehrkanoon | Deep hybrid neural-kernel networks using random Fourier features[END_REF][START_REF] Liu | Learning deep kernels for non-parametric two-sample tests[END_REF] opening the door to a wide range of applications.

Conclusion and Perspectives

In this thesis, we have introduced a general learning framework based on vector-valued reproducible kernel Hilbert spaces and integral losses, akin to estimate predictive models with functional outputs. Due to the complexity of the task, certain approximations are required for the resulting optimization problems to be tractable. These can be applied either on the loss side, gaining a workable representation of the solution by means of a double representer theorem or on the hypothesis space by restricting the coefficients to live in a finite-dimensional subspace of the vv-RKHS, chosen to be well-suited to the problem. These approximations enable primal or dual algorithm using variations of gradient descent. Among them, the use of random Fourier features allows to tackle large-scale learning settings. This learning framework is proved to be useful to extend functional output regression to more involved losses such as robust losses. Moreover, scenarii where the observed training outputs are not functional but for which the task to be solved is parameterized can also benefit from the functional point of view, resulting in the infinite task learning framework.

We exemplified the proposed framework in quantile regression, cost-sensitive classification, density level set estimation and emotion transfer with non scalar outputs and parameters. In these applications, the space of parameters Θ is relatively simple; future work could include learning scenarii where Θ is a complex structured set, over which the design of the output kernel is crucial. One could extend for instance quantile regression to Hilbert-valued output random variables that result in Θ being the unit ball of a certain Hilbert space.

Functional output regression in vv-RKHSs was tackled by leverage the duality view, allowing to design losses through infimal convolution whose associated estimator shows robustness or sparsity. This idea could be pushed further by considering different norms to convolute with the loss function, giving rise to a large family of dual problems that might influence the behavior of the estimator in interesting ways.

From a statistical learning point of view, many questions are still to be investigated. In particular, while we have studied the generalization capabilities of the ITL estimator, the convergence of the risk associated to the estimator towards the Bayes risk remains an open question. The analysis of such risk could be performed in the context of random Fourier features, possibly allowing for the derivation of optimality conditions on the sampling measure to achieve good performance with the fewer features possible.

From a modeling point of view, we have evocated the combination of neural networks with kernels to deal with objects over which choosing a good kernel is notoriously hard like images or documents. These hybrid approaches seem promising and could be further developed, as learning the kernel brings nonconvexity to the optimization problems and requires dedicated solvers. Moreover, the class of vv-RKHSs used in this thesis rely on separable kernels which can be less appropriate in certain applications as they induce a tensorial structure in the Gram matrix. Going beyond such kernels requests to start anew from the problem formulation as we would lose the computability advantages brought by the separability.

Conclusion and perspectives

Concerning further applications, we are particularly interested in applying the developed framework to the study of dynamical systems. Imagining applications where the output variable is a trajectory having to satisfy a certain differential equation, choosing appropriate output vv-RKHS would ensure by design a control over the properties of the model, and using the reproducing property for derivatives in a RKHS paves the way to tractable optimization problems.

Finally, we want to mention the potential of ITL for meta-modeling and meta-learning [START_REF] Ton | Noise contrastive meta-learning for conditional density estimation using kernel mean embeddings[END_REF], which could be deepened in further works.

  Learning Tasks . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.1.1 Parameterized Tasks . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.1.2 Learning Vector-Valued Functions . . . . . . . . . . . . . . . . . 19 1.1.3 Learning Function-Valued Functions . . . . . . . . . . . . . . . . 20 1.2 Research Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 1.4 Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Example 2. 22 (

 22 Laplacian kernel eigendecomposition, Kadri et al. (2016)). Let Θ = [0, 1] and µ be the Lebesgue measure. Let k(θ, θ ′ ) = e -|θ-θ ′ | . The eigendecomposition of T k is given by

Example 2 .

 2 23 (Approximate eigendecomposition, Hoegaerts et al. 2005). Let k be a continuous kernel on a compact metric space Θ endowed with a measure µ. Let m > 0 and (θ j ) m j=1 i.i.d.

  and predicting sign(s) (Bach et al., 2006). Defining the asymmetry coefficient θ = c + c + +c -and restraining the values of (c + , c -) to the line c + +c -= 1 the training can be performed jointly for all values of θ ∈ [0, 1] by considering (3.2) with ℓ(θ, u, v) = θ -{-1} (v) max(0, 1 -uv), µ encoding the user importance to different values of θ, and identifying the random variable Y ∈ {-1, 1} with a random variable of constant functions in L 2 [Θ, µ] (Brault et al., 2019).

Lemma 3. 7 (

 7 Optimization task for α). The solution α of Problem 3.20 is the solution of the following linear equation K αA + λnm α = λnmY.(3.21)

Proposition 3 .

 3 13 (Fenchel-Legendre conjugate of Integral Losses, Rockafellar (1976)). Let ℓ : Θ × U × U → -∞, +∞ be a normal convex integrand. If Y satisfies Assumption 3.12 then for all y ∈ Y, I ℓ (•, y) ⋆ = I ℓ ⋆ (•, y). (3.36)

  y) is continuous, this results in Equation (3.37). The result in Proposition 3.14 endorses the modeling choice Y = L 2 [Θ, µ], since it avoids the difficulty brought in by the term involving T -1G and additional assumptions. One can then directly work with dual variables (α i ) n i=1 ∈ L 2 [Θ, µ; U] and solve Problem 3.7. Before delving into this solution, we list below the Fenchel-Legendre conjugates of a various integral losses used in practice.

  et al. 2008; Shin and Lee 2016; Qingguo 2017; Boente et al. 2020).

Definition 4. 5 .

 5 The integral Huber loss of parameter κ > 0 is given by ∀y ∈ Y, I ℓ H,κ (y) := Θ ℓ H,κ (y(θ))dµ(θ). (4.10) In Section 4.3.1, we propose a solution of Problem 4.6 with the vectorial Huber loss, before considering the use of the integral Huber loss in Section 4.3.2. Both losses add a different term in Problem 4.6 that constrains the dual variables, thus resulting in robustness of the estimator to outliers.

Figure 4

 4 Figure 4.1 -Standard square loss and Huber loss in 1 and 2 dimensions (κ = 0.8).

4. 4 .

 4 SPARSE ESTIMATORS WITH ϵ-INSENSITIVE LOSSES 73 Algorithm 4.2 Projected Gradient Descent for Integral Huber input : Gram matrices K X , K Θ , data matrix Y, regularization parameter λ, Huber parameter κ, gradient step γ init :

  As proposed in Sangnier et al. (2017); Laforgue et al. (2020) we extend them from R p to any Hilbert space Y. Definition 4.12. Let L : Y → R + be a convex loss such that L(0) = 0, and ϵ > 0. The ϵ-insensitive version of L, denoted L ϵ , is defined by L ϵ (y) = (L □ χ Bϵ )(y), or again:

Figure 4

 4 Figure 4.2 -Standard and ϵ-insensitive versions of the square loss in 1 and 2 dimensions (ϵ = 1.5).

Algorithm 4. 4

 4 Proximal Gradient Descent for Integral ϵ-ridge input : Gram matrices K X , K Θ , data matrix Y, regularization parameter λ, Huber parameter κ, gradient step γ init :

proximal step for row i from 1 to n do for column j from 1 to m do 39 αFigure 4

 394 Figure 4.3 -Lip dataset.

Figure 4

 4 Figure 4.4 -Synthetic dataset. 50 samples are shown.

Figure 4

 4 Figure 4.5 -Synthetic dataset global outliers. 15 samples are shown.

Figure 4 . 6 -

 46 Figure 4.6 -Leave-one-out (LOO) generalization error when learning with Huber loss and contaminated real data for several scales of outliers s o ∈ {0.5, 1, 2, 5}. Left: vectorial Huber loss. Right: integral Huber loss. For large enough Huber loss parameter κ the estimates (solid) coincide with that of the ridge solution (dashed).

Figure 4

 4 Figure 4.7 -Test MSE when learning with Huber loss and contaminated synthetic data for several scales of outliers s o ∈ {2, 4, 6, 10}. Left: vectorial Huber loss. Right: integral Huber loss.

  (a) LOO error. (b) Sparsity of the coefficients.

Figure 4

 4 Figure 4.8 -LOO error and sparsity when learning with vectorial ϵ-ridge loss and contaminated real data for several scales of local outliers s o ∈ {0.2, 0.5, 0.8, 1.2}. Left: Test MSE. Right: Proportion of zero coefficients.

98 CHAPTERFigure 5 . 1 -

 9851 Figure 5.1 -Pinball loss for θ = 0.8.

θ

  Figure 5.2 -Impact of crossing penalty on toy data. Left plot: strong non-crossing penalty (λ nc = 10). Right plot: no non-crossing penalty (λ nc = 0). The plots show 100 quantiles of the continuum learned, linearly spaced between 0 (blue) and 1 (red).Notice that the non-crossing penalty does not provide crossings to occur in the regions where there is no points to enforce the penalty (e.g. x ∈ [0.13, 0.35]). This phenomenon is alleviated by the regularity of the model.

Figure 5

 5 Figure 5.3 -Impact of the smoothing of the pinball loss for different values of κ.

Figure 5

 5 Figure 5.4 -Impact of the number of hyperparameters sampled.

Figure 5

 5 Figure 5.5 -Fish bone images of otholiths at different ages. Dataset from the Norwegian Marine Data Center.

Figure 5 . 6 -

 56 Figure 5.6 -Estimated quantiles of a few otholiths samples chosen randomly in the dataset. The colorbar on the right indicates the ground truth age of the fish, and on the left is plotted its estimated quantile as a function of θ. We can see that the quantiles are nondecreasing thanks to the non crossing penalty.

(a)

  Figure 5.7 -Illustration of ICSC on toy datasets used inTable 5.2.

Figure 6 .

 6 Figure 6.1 -Illustration of emotion transfer.

  Vemulapalli and Agarwala 2019) to which the presented framework can be naturally adapted. Throughout this work faces are represented by landmark points. Landmarks have been proved to be a useful representation in facial recognition (Saragih et al., 2009; Scherhag et al., 2018; Zhang et al., 2015), 3D facial reconstruction and sentiment analysis. Tautkute et al. (2018) have shown that emotions can be accurately recognized by detecting changes in the localization of the landmarks. Given M number of landmarks on the face, this means a description x ∈ X := R 2M , and hence a d := 2M -dimensional representation.

120CHAPTER 6 .

 6 EMOTION TRANSFERtrajectories are observable at finite many emotions θi,j j∈[m]

=

  :y m(i-1)+l,j j∈[m] ∈ X m (θ m(i-1)+l,j ) j∈[m] = θi,j j∈[m] ∈ Θ m , where (i, l) ∈ [n] × [m]and the number of pairs is t = nm. Having defined this dataset one can optimize the same objective in Equation (6.3) as before. Particularly, this means that the pair (i, l) plays the role of index i of the previous case. The(θ i,j ) i,j∈[t]×[m] is an extended version of the θi,j i,j∈[t]×[m]to match the indices going from 1 to t in Equation (6.3).We leverage the flexible class of vector-valued reproducing kernel Hilbert spaces (vv-RKHS;[START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces and universality[END_REF]) for the hypothesis class schematically illustrated in Equation (6.2). Learning within vv-RKHS has been shown to be relevant for tackling 6.3. OPTIMIZATION 121 function-valued regression[START_REF] Kadri | Operatorvalued kernels for learning from functional response data[END_REF]. The construction follows the structure h : X → (Θ → X)

Figure 6 . 2 -

 62 Figure 6.2 -Illustration of the landmark edge maps for different emotions and both datasets.

  al., 2010; Ferraty et al., CHAPTER 1. MOTIVATION AND CONTRIBUTIONS 2011; Oliva et al., 2015;

  ⋆ , A. Lambert ⋆ , Z. Szabó, M. Sangnier and F. d'Alché-Buc. Infinite Task Learning in RKHSs. In International Conference on Artificial Intelligence and Statistics (AISTATS), pages 1294-1302, 2019.

	Background on Learning with
	Operator-Valued Kernels

• A. Lambert ⋆ , R. Brault ⋆ , Z.

Szabó

, M. Sangnier and F. d'Alché-Buc. A Functional Extension of Multi-Output Learning In International Conference on Machine Learning (ICML): Adaptive & Multitask Learning workshop (AMTL), 2019. • P. Laforgue, A. Lambert, L. Brogat-Motte, and F. d'Alché-Buc. Duality in RKHSs with Infinite Dimensional Outputs: Application to Robust Losses In International Conference on Machine Learning (ICML), pages 5598-5607, 2020. • A. Lambert ⋆ , S. Parekh ⋆ , Z. Szabó, and F. d'Alché-Buc. Emotion Transfer Using Vector-Valued Infinite Task Learning Technical report, 2021. (https: //arxiv.org/abs/2102.05075). (under review) 2 "Kernels are dead! Long live kernels!" Contents 2.1 Convex Optimization Tools . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2 Kernel Methods for Machine Learning . . . . . . . . . . . . . . . . . . . 28 2.2.1 Scalar-valued Kernels and RKHSs . . . . . . . . . . . . . . . . . 29 2.2.2 Operator-valued Kernels and vv-RKHSs . . . . . . . . . . . . . . 35 2.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

Table 2

 2

.1 -Useful Fenchel-Legendre conjugates, for any f, g : H → -∞, +∞ . Function Fenchel-Legendre conjugate

  •, x) is called the canonical feature map, with V in Equation (2.9) being H k itself. If one gets a RKHS H from Definition 2.13, then the existence for all x ∈ X of a function k(•, x) ∈ H stems from Riesz representation theorem, and taking x → k(•, x) as a feature map in Equation (2.9) ensures that H is the RKHS associated to the kernel k.

Regularized Empirical risk minimization in vv-RKHSs Let

			CHAPTER 2. BACKGROUND ON LEARNING WITH
			OPERATOR-VALUED KERNELS
	Such object can be defined for the larger class of Mercer kernels, which are OVKs such
	that x → K(x, x)	op	is locally bounded, and H K ⊂ C(X, Y) (see Carmeli et al. (2010)
	for more details about this).
	Proposition 2.42 (Carmeli et al. (2010)). Let Θ be a compact metric space endowed
	with a Borel probability measure µ, and K : Θ × Θ → L(Y) be an OVK. Assume that
	H K ⊂ C(Θ, Y) and that K is bounded in operator norm on supp(µ). If K(θ, θ) is a
	compact operator for all θ ∈ Θ, then T K,µ is a compact operator.
	Proposition 2.42 allows to perform the eigendecomposition of T K,µ , and we get a similar
	characterization of H K than the one in Equation (2.12). We refer the interested reader
	to Carmeli et al. (2010).
			(X, Y) ∈ X × Y be
	random variables, we assume access to i.i.d. samples (x i , y i ) n i=1 . Given a proper, convex
	lower-semicontinuous loss function L : Y × Y → R (with respect to its first argument),
	we consider the learning problem
			.19)

  Contents 3.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2 The Special Square Loss Case . . . . . . . . . . . . . . . . . . . . . . . . 46

3.2.1 Functional Observation Case . . . . . . . . . . . . . . . . . . . . 46 3.2.2 Partially Observed Case . . . . . . . . . . . . . . . . . . . . . . . 48 3.3 Solving in the Primal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3.1 Sampling Schemes and Representer Theorems . . . . . . . . . . . 51 3.3.2 Random Features Based Learning . . . . . . . . . . . . . . . . . 54 3.4 Solving in the Dual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4.1 Fenchel-Legendre Conjugate of Integral Losses . . . . . . . . . . 57 3.4.2 Integral Operator Eigenbasis Representation . . . . . . . . . . . 59 3.4.3 Proximal Algorithms and Approximated Quadratic Forms . . . . 60 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

Table 3

 3

.1 -Considered problem family: examples. QR: quantile regression. CSC: cost-sensitive classification. FOR: functional output regression. ET: emotion transfer.

Table 3 .

 3 2 -Summary of proposed algorithms. GD: (sub) gradient descent. SGD: stochastic (sub) gradient descent. PGD: proximal gradient descent.

	Section type	Y	parameterization	loss	algorithm
	3.2.1	closed form L 2 [Θ, µ; U] eigenbasis of T G	square loss	analytic
	3.2.2	closed form H G	double representer square loss	analytic
	3.3.1	primal	H G	double representer sampled	GD
	3.3.2	primal	H G	ORFF	any	SGD
	3.4.2	dual	L 2 [Θ, µ; U] eigenbasis of T G	Assumption 3.16 GD
	3.4.3	dual	L 2 [Θ, µ; U] linear splines	Assumption 3.18 PGD

  INTEGRAL LOSSESWe can however wonder what would happen if we tried to compute I ℓ (•, y) ⋆ with Y = H G . The following lemma provides a response with additional assumptions of universality of kernel G and continuity of I ℓ . Let k Θ be a universal kernel, and A ∈ M s (R) an invertible matrix. Let H G be the vv-RKHS associated to the kernel G = k Θ A. Assume that µ has full support, in other words supp µ = Θ. Assume finally that I ℓ is continuous. Then for all y

	Proposition 3.14.

  4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.2 Problem Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.3 Robust Estimators with Huber Loss . . . . . . . . . . . . . . . . . . . . 69 4.3.1 Vectorial Huber Loss . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.3.2 Integral Huber Loss . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.4 Sparse Estimators with ϵ-Insensitive Losses . . . . . . . . . . . . . . . . 73 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

4.4.1 Vectorial ϵ-Ridge . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.4.2 Integral ϵ-Ridge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.5 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.5.1 Huber Losses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.5.2 ϵ-insensitive Losses . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.6

  .1.

	Extending the Huber loss to higher dimension can be done multiple ways. A first one
	would consist in replacing 1 2 (•) by 1 2 • 2

Y and κ| • | by κ • Y in Equation (4.7), which we refer to as the vectorial Huber loss. Definition 4.4. The vectorial Huber loss of parameter κ > 0 is given by

  a.e. is the ball of radius κ for the ∞-norm. The dual to Problem 4.5 in the integral Huber loss case writes inf

	Proof This is an application of Proposition 3.13 which states that I ⋆ ℓ H,κ = I ℓ ⋆ H,κ combined with Proposition 4.3 that gives an expression for ℓ ⋆ H,κ .
	We can now present Problem 4.6 for the integral Huber loss scenario.
	Proposition 4.11.

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.2 From Parameterized to Infinite Task Learning . . . . . . . . . . . . . . . 84 5.2.1 Learning Parameterized Tasks . . . . . . . . . . . . . . . . . . . 85 5.2.2 Solving a Finite Number of Tasks as Multi-Task Learning . . . . 86 5.2.3 Towards Infinite Task Learning . . . . . . . . . . . . . . . . . . . 87 5.3 Generalization Analysis through Uniform Stability . . . . . . . . . . . . 88 5.3.1 Generalization Bound for Non-Approximated Scheme . . . . . . 90 5.3.2 Generalization Bound for Approximated Scheme . . . . . . . . . 93 5.4 Quantile Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.4.1 Enforcing Shape Constraints . . . . . . . . . . . . . . . . . . . . 99 5.4.2 Smoothing the Loss Function . . . . . . . . . . . . . . . . . . . . 101 5.4.3 Influence of the Number of Sampled Locations . . . . . . . . . . 104 5.4.4 Deep Kernel Learning for Quantile Regression on Images . . . . 104 5.5 Cost-Sensitive Classification . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.6 Density Level Set Estimation . . . . . . . . . . . . . . . . . . . . . . . . 108 5.6.1 Representer Theorem for Mixed Regularization . . . . . . . . . . 109 5.6.2 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . 113 5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

	Contents
	5.1

  For a kernel k Θ (θ 1 , θ 2 ) of two variables, assumed to be C 2 , we denote ∂ 1 k Θ the partial derivative of k with respect to θ 1 and ∂ 2 k Θ the partial derivative of k with respect to θ 2 . From[START_REF] Zhou | Derivative reproducing properties for kernel methods in learning theory[END_REF], for all θ ∈ Θ, (∂ 2 k Θ )(•, θ) ∈ H k Θ and the reproducing property in H k Θ reads as

	25)
	Remark 5.26. There is a slight notation overloading with the subdifferential operator
	∂ introduced in Section 2.1. In this chapter we are dealing with differentiable functions
	and identify ∂f (θ) = {f ′ (θ)} with f ′ (θ)'?.

Table 5 .

 5 1 -Quantile regression on 20 UCI datasets. Reported: 100×value of the pinball loss, 100×crossing loss (smaller is better). pval: outcome of the Mann-Whitney-Wilcoxon test of JQR vs IQR and Independent vs IQR. Boldface: significant values w.r.t. IQR.

	dataset	JQR				IND-QR		IQR
	(pinball	pval)	(cross	pval)	(pinball	pval)	(cross	pval)	pinball	cross

Table 5 .

 5 2 -ICSC vs Independent (IND)-CSC. Higher is better.

	5.5. COST-SENSITIVE CLASSIFICATION				107
	Dataset	Method	θ = -0.9	θ = 0	θ = +0.9
			sensitivity specificity sensitivity specificity sensitivity specificity
	Two-Moons	IND ICSC	0.3 ± 0.05 0.32 ± 0.05	0.99 ± 0.01 0.99 ± 0.01	0.83 ± 0.03 0.84 ± 0.03	0.86 ± 0.03 0.87 ± 0.03	0.99 ± 0 1 ± 0	0.32 ± 0.06 0.36 ± 0.04
	Circles	IND ICSC	0 ± 0 0.15 ± 0.05	1 ± 0 1 ± 0	0.82 ± 0.02 0.82 ± 0.02	0.84 ± 0.03 0.84 ± 0.03	1 ± 0 1 ± 0	0 ± 0 0.12 ± 0.05
	Iris	IND ICSC	0.88 ± 0.08 0.89 ± 0.08	0.94 ± 0.06 0.94 ± 0.05	0.94 ± 0.05 0.94 ± 0.06	0.92 ± 0.06 0.92 ± 0.05	0.97 ± 0.05 0.97 ± 0.04	0.87 ± 0.06 0.90 ± 0.05

  To that end, we use a function-valued model h ∈ H K where K = k X Id H k Θ , and model t as a function belonging to an RKHS H k b where k

b : Θ × Θ → R is a scalar-valued kernel on Θ. We want each ĥ(•)(θ), t(θ) to be (close to) a solution of Problem 5.31. For this reason we introduce the L 2 -RKHS mixed regularizer

  It turns out that the estimator resulting from Problem 5.34 enjoy a representer expression. To prove it, we need the following two lemmas. They make use of the equivalence of views between H K and H k X ⊗ H k Θ (see Remark 2.41). Let k X and k Θ be two scalar-valued kernels on X and Θ, and K :

	110	CHAPTER 5. INFINITE TASK LEARNING
	and propose to solve
	ĥ, t = arg min h∈H K ,t∈H k b	1 n	n i=1	I ℓ (h(x i ), t) + Ω(h) +	λ 2	t	2 H k b	.	(5.34)
	Lemma 5.30.
									33)

  .35) where + m j=1 Im(K θ j ) denotes the Minkowski sum of the sets Im(K θ j )

j∈

[m]

  Let k X , k Θ , k b be three scalar-valued kernels, with k Θ a strictly positive definite kernel, k b bounded, and K = k X Id H k Θ . Then Problem 5.34 admits a solution ĥ, t ∈ H K × H k b which writes for all

	Remark 5.32. Lemma 5.31 highlights the difference between the regularizers Ω(h),
	Ω(h) and 1 2 h	2 H K	. While the latter is uniformly stronger than the two others, on
	the specific subspace V, Ω(h) is γ-strongly convex. As m grows, the strong convexity
	constant γ goes to zero, and when m = +∞, we recover Ω(h) which does not bring
	strong convexity to the problem and is thus strictly weaker than 1 2 h	2 H K	.
	We are now ready to prove a double representer theorem for the solution of Prob-
	lem 5.34.
	Theorem 5.33.

which proves the γ-strong convexity of I on V.

  al., 2016; Wynen et al., 2018; Jing et al., 2020) in machine learning. In a nutshell, style transfer refers to the transformation of an object according to a target style. It has found numerous applications in computer vision (Ulyanov et al., 2016; Choi et al., 2018; Puy and Pérez, 2019; Yao et al., 2020), natural language processing

  i,j∈[t]×[m] ∈ U nt which can be gathered in a matrix α ∈ M nt,d (R) so that α is solution to K αA + tmλ α = Y (6.8) where the Gram matrix K = [k i,j ] i,j∈[tm] ∈ M tm (R), and the matrix consisting of all the observations Y = [Y i ] i∈[tm] ∈ M tm,d (R) are defined ask m(i 1 -1)+j 1 ,m(i 2 -1)+j 2 := k X (x i 1 , x i 2 )k Θ (θ i 1 ,j 1 , θ i 2 ,j 2 ), (i 1 , j 1 ), (i 2 , j 2 ) ∈ [t] × [m], Y m(i-1)+j := y ⊤ i,j , (i, j) ∈ [t] × [m].Computational complexity: In case of A = Id d , the complexity of the closed form solution is O (tm) 3 . If all the samples are observed at the same locations(θ i,j ) i,j∈[t]×[n] , i.e. θ i,j = θ l,j for ∀(i, l, j) ∈ [t] × [t] × [m], then the Gram matrix K has a tensorial structure K = K X ⊗ K Θ with K X = [k X (x i , x j )] i,j∈[t] ∈ R t×t and K Θ = [k Θ (θ 1,i , θ 1,j )] i,j∈[m] ∈ R m×m .In this case, the computational complexity reduces to O t 3 + m 3 . If additional scaling is required one can leverage recent dedicated kernel ridge regression solvers(Rudi et al., 2017;[START_REF] Meanti | Kernel methods through the roof: handling billions of points efficiently[END_REF]). If A is not identity but invertible, then multiplying Equation (6.8) with A -1 gives K α + tmλ αA -1 = YA -1 which is a Sylvester equation for which efficient custom solvers exist (El[START_REF] Guennouni | Block Krylov subspace methods for solving large Sylvester equations[END_REF]. If A is not invertible, one can use Singular Value Decomposition to reduce the dimensionality of α and fall back on the invertible case.

	Remarks:
	•

Alternative optimization techniques to deal with countable or continuous hyperparameter spaces could include semi-infinite(Stein,

2012) or bi-level programming[START_REF] Wen | Linear bi-level programming problemsa review[END_REF].

We used a Gaussian process model and minimized the expected improvement. The optimizer was initialized using 27 samples from a Sobol sequence and ran for 50 iterations.

(a) Theoretical level sets of the Gaussian mixture, scaled by the value of the density. (b) Estimated level sets of the Gaussian mixture, scaled by 1mass of the sets Figure 5.10 -Proof of concept of ∞-OCSVM for DLSE on a simple Gaussian mixture.

To keep the notation simple, we assume that m is the same for all the zi-s.

In this case θi,j is a literal copy of θi,j which helps to get a unified formulation with the joint emotional input setting.

The code is available at https://github.com/yunjey/stargan.

Remerciements

We provided excess risk guarantees for the studied ITL scheme, and demonstrated its practical efficiency and flexibility in various tasks including cost-sensitive classification, quantile regression and density level set estimation.

While in all aforementioned tasks Θ and Y are one-dimensional, the ITL framework can be adapted to more complex tasks such as emotion transfer tasks as emphasized in Chapter 6.

Emotion representation:

We represented emotion labels as points in the 2D valencearousal space (VA, [START_REF] Russell | A circumplex model of affect[END_REF]. Particularly, we used a manually annotated part of the large-scale AffectNet database [START_REF] Mollahosseini | AffectNet: A database for facial expression, valence, and arousal computing in the wild[END_REF]. For all samples of a particular emotion in the AffectNet data, we computed the centroid (data mean) of the valence and arousal values. The resulting ℓ 2 -normalized 2D vectors constituted our emotion representation as depicted in Figure 6.3. The normalization is akin to assuming that the modeled emotions are of the same intensity. In our experiments, the emotion 'neutral' was represented by the origin. Such an emotion embedding allowed us to take into account prior knowledge about the angular proximity of emotions in the VA space, while keeping the representation simple and interpretable for post-hoc manipulations.

Kernel design:

We took the kernels k X , k Θ to be Gaussian on the landmark representation space and the emotion representation space, with respective bandwidth γ X and γ Θ . A was assumed to be Id d unless specified otherwise.

Quantitative Performance Assessment

In this section we provide a quantitative assessment of the proposed vITL approach.

Performance measures:

We applied two metrics to quantify the performance of the compared systems, namely the test mean squared error (MSE) and emotion classification accuracy. The classification accuracy can be thought of as an indirect evaluation. To compute this measure, for each dataset we trained a ResNet-18 classifier to recognize emotions from ground-truth landmark edge maps (as depicted in Figure 6.2). The trained network was then used to compute classification accuracy over the predictions at test time. To rigorously evaluate outputs for each split of the data, we used a classifier trained on RaFD to evaluate KDEF predictions and vice-versa; this also allowed us to make the problem more challenging. The ResNet-18 network was appropriately modified to take grayscale images as input. During training, we used random horizontal flipping and cropping between 90-100% of the original image size to augment the data. All the images were finally resized to 224 × 224 and fed to the network. The network was trained from scratch using the stochastic gradient descent optimizer with learning rate and momentum set to 0.001 and 0.9, respectively. The training was carried out CHAPTER 6. EMOTION TRANSFER

and train a vITL system with the resulting A. While in this case A is no more invertible, each coefficient αi,j from Equation (6.7) belongs to the r-dimensional subspace of R d generated by the eigenvectors associated to the r largest eigenvalues of Y ⊤ Y. This makes a reparameterization possible and leads to a decrease in the size of the model, going from t × m × d parameters to t × m × r. We report in Figure 6.5 the resulting test MSE performance (mean ± standard deviation) obtained from 10 different splits, and empirically observe that r = 20 suffices to preserve the optimal performances of the model.

Learning under a missing obervation regime:

To assess the robustness of vITL w.r.t. missing data, we considered a random mask (η i,j ) i∈[n],j∈[m] ∈ {0, 1} n×m ; a sample z i (θ i,j ) was used for learning only when η i,j = 1. Thus, the percentage of missing data was p :=

The experiment was repeated for 10 splits of the dataset, and on each split we averaged the results using 4 different random masks (η i,j ) i∈[n],j∈ [m] . The resulting test MSE of the predictor as a function of p is summarized in Figure 6.6. As it can be seen, the vITL approach is quite stable in the presence of missing data on both datasets.

Qualitative Analysis

In this section we show example outputs produced by vITL in the context of discrete and continuous emotion generation. While the former is the classical task of synthesis given input landmarks and target emotion label, the latter serves to demonstrate a key benefit of our approach, which is the ability to synthesize meaningful outputs while continuously traversing the emotion embedding space.

Discrete emotion generation:

In Figure 6.7 and Figure 6.8 we show qualitative results for generating landmarks using discrete emotion labels present in the datasets. For vITL, not only are the emotions recognizable, but landmarks on the face boundary are reasonably well-synthesized and other parts of the face visibly less distorted when compared to StarGAN. The identity in terms of the face shape is also better preserved.

Continuous emotion generation:

Starting from neutral emotion, continuous generation in the radial direction is illustrated in Figure 6.9. The landmarks vary smoothly and conform to the expected intensity variation in each emotion on increasing the radius of the vector in VA space. We also show in Figure 6.10 the capability to generate intermediate emotions by changing the angular position, in this case from 'happy' to 'surprised'. For a more fine-grained video illustration traversing from 'happy' to 'sad' along the circle, see the GitHub repository. In this thesis, we introduce an original functional extension of multi-output learning called Infinite Task Learning (ITL), that allows to jointly solve an infinite number of parameterized tasks, including for instance quantile regression, cost-sensitive classification and density level set estimation. We propose a learning framework based on convex integral losses that encompasses the ITL problem and functionvalued regression. Optimization schemes dedicated to solving the associated regularized empirical risk minimization pro-blems are designed. By sampling the integral losses, we derive finite-dimensional representation of the solution under several choices of regularizers or shape constraints penalties, while keeping theoretical guarantees over their generalization capabilities. We also employ dualization techniques with the benefit of bringing desirable properties such as robustness or sparsity to the estimators thanks to the use of convoluted losses. Scalability issues are addressed by deriving optimization algorithms in the the context of approximated OVKs whose corresponding vv-RKHSs are of finite dimension. The use of trainable deep architectures composed by a neural network followed by a shallow kernel layer is also investigated as a way to learn the kernel used in practice on complex data such as images.

We apply these techniques to various ITL problems and to robust function-to-function regression, that are tackled in the presence of outliers. We also cast style transfer problems as a vectorial output ITL problem and demonstrate its efficiency in emotion transfer.

Institut Polytechnique de Paris 91120 Palaiseau, France