Maxime Sangnier
email: maxime.sangnier@upmc.fr

Olivier Fercoq
email: olivier.fercoq@telecom-paristech.fr

Yung-Kyun Noh

Min-Ling Zhang

Data sparse nonparametric regression with -insensitive losses

Keywords: Quantile regression, Expectile regression, Operator-valued kernel

Leveraging the celebrated support vector regression (SVR) method, we propose a unifying framework in order to deliver regression machines in reproducing kernel Hilbert spaces (RKHSs) with data sparsity. The central point is a new definition of -insensitivity, valid for many regression losses (including quantile and expectile regression) and their multivariate extensions. We show that the dual optimization problem to empirical risk minimization with -insensitivity involves a data sparse regularization. We also provide an analysis of the excess of risk as well as a randomized coordinate descent algorithm for solving the dual. Numerical experiments validate our approach.

Introduction

In supervised learning, sparsity is a common feature [START_REF] Hastie | Statistical Learning with Sparsity: The Lasso and Generalizations[END_REF]. When it is put on the table, the first reaction is certainly to invoke feature selection [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF]. In many domains such as finance and biology, feature selection aims at discovering the covariates that really explain the outcomes (hopefully in a causal manner). However, when dealing with nonparametric modeling, sparsity may also occur in the orthogonal direction: data, which support ultimately the desired estimator. Data sparsity is mainly valuable for two reasons: first, controlling the number of training data used in the model can prevent from overfitting. Second, data sparsity implies less data to store and less time needed for prediction as well as potentially less time needed for training. The main evidence of this phenomenon lies in the celebrated support vector machines (SVMs) [START_REF] Boser | A Training Algorithm for Optimal Margin Classifiers[END_REF][START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]Mohri et al., 2012). The principle of margin maximization results in selecting only a few data points that support the whole prediction function, leaving the others sink into oblivion.

Data sparsity is particularly known for classification thanks to SVMs but is less predominant for regression. Indeed, while classifying relies on the few points close to the frontier (they are likely to be misclassified), regression is supported by points for which the outcome is far from the conditional expectation. Those points, that are likely to be poorly characterized by the prediction, gather most of the training points. Thus, by nature, regressors tend to be less data sparse than classifiers. In order to obtain the valued data sparsity for a wide class of regression problems, we consider regression relying on empirical risk minimization [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF] with non-linear methods based on RKHSs for uni and multidimensional outcomes [START_REF] Micchelli | On Learning Vector-Valued Functions[END_REF].

The companion of SVM for regression, SVR [START_REF] Drucker | Support Vector Regression Machines[END_REF], is probably the most relevant representative of data sparse regression methods in the frequentist domain (in the Bayesian literature, relevance vector machines can deliver sparse regressors [START_REF] Tipping | Sparse Bayesian Learning and the Relevance Vector Machine[END_REF]). SVR makes use of the celebrated -insensitive loss ξ ∈ R → max (0, |ξ| -), where > 0, in order to produce data sparsity. In practice, residues close to 0 are not penalized, and consequently corresponding points are not considered in the prediction function. Later, [START_REF] Park | Quantile regression with an epsilon-insensitive loss in a reproducing kernel Hilbert space[END_REF] extended this -insensitive loss to quantile regression (QR) [START_REF] Koenker | Quantile Regression[END_REF] in an ad hoc manner.

These approaches are in deep contrast compared to kernel ridge regression (using the loss ξ ∈ R → 1 2 ξ 2 ), that does not benefit from data sparsity. To remedy that lack of sparsity, [START_REF] Lim | Operator-valued kernel-based vector autoregressive models for network inference[END_REF] leveraged the representer theorem [START_REF] Schölkopf | A Generalized Representer Theorem[END_REF] for ridge regression and penalized the coefficients associated to each training point. The same procedure was used by [START_REF] Zhang | On Quantile Regression in Reproducing Kernel Hilbert Spaces with the Data Sparsity Constraint[END_REF] for QR but with a sparse constraint instead of a sparse regularization. Contrarily to -insensitive losses, the latter approaches constrain the hypotheses while keeping the vanilla regression losses.

The previous works are either restricted to unidimensional losses [START_REF] Drucker | Support Vector Regression Machines[END_REF][START_REF] Park | Quantile regression with an epsilon-insensitive loss in a reproducing kernel Hilbert space[END_REF][START_REF] Zhang | On Quantile Regression in Reproducing Kernel Hilbert Spaces with the Data Sparsity Constraint[END_REF] or compelled to assume a finite representation of the optimal regressor [START_REF] Lim | Operator-valued kernel-based vector autoregressive models for network inference[END_REF][START_REF] Zhang | On Quantile Regression in Reproducing Kernel Hilbert Spaces with the Data Sparsity Constraint[END_REF]. Our contribution in this context is to provide a novel and unifying view of -insensitive losses for data sparse kernel regression, that is theoretically well founded (Section 3). In the spirit of SVR, we show that, for a large number of regression problems (uni or multivariate, quantile or traditional regression), data sparsity can be obtained by properly defining an -insensitive variant of the loss in the data-fitting term which turns to be a data sparse regularization in the dual. This framework contains SVR as a special case but also a new loss for QR (that is different from the one introduced by [START_REF] Park | Quantile regression with an epsilon-insensitive loss in a reproducing kernel Hilbert space[END_REF]) and expectile regression [START_REF] Newey | Asymmetric Least Squares Estimation and Testing[END_REF], as well as their multivariate extensions (Section 4) in the context of vector-valued RKHSs. We also provide an analysis of the excess of risk (Section 5). Special attention is paid to the dual form of the optimization problem and a randomized primal-dual block coordinate descent algorithm for estimating our data sparse regressor (Section 6) is derived. Numerical experiments conclude the discussion (Section 7).

Framework

Notations In the whole paper, we denote vectors with bold-face letters. Moreover, τ ∈ (0, 1) is a quantile level, τ ∈ (0, 1) p is a vector of quantile levels, n and p are two positive integers (samples size and dimension), [n] is the range of integers between 1 and n, α l ∈ R n denotes the l th row vector of any α = (α i ) 1≤i≤n ∈ (R p ) n , I • is the indicator function (set to 1 when the condition is fulfilled and 0 otherwise), 1 denotes the all-ones vector, χ • is the characteristic function (set to 0 when the condition is met and to ∞ otherwise), (respectively ≺) denotes the pointwise (respectively strict) inequality, diag is the operator mapping a vector to a diagonal matrix, ∂ψ is the subdifferential of any function ψ (when it exists), prox ψ denotes its proximal operator [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF], proj 1 is the projector onto the span of 1.

General framework Let X be a nonempty input space and R p be the output Hilbert space (for a given integer p). Let also H be a Hilbert space of hypotheses h : X → R p and {(x i , y i )} 1≤i≤n be a training sample of n couples (x i , y i ) ∈ X × R p . We consider the setting of regularized empirical risk minimization for regression, based on a real-valued convex loss function : R p → R. This is to provide a minimizer to the optimization problem:

minimize h∈H, b∈R p λ 2 h 2 H + 1 n n i=1 (y i -(h(x i ) + b)) , (P1)
where λ > 0 is a trade-off parameter and (Mohri et al., 2012).

Let, for any x ∈ X , E x : h ∈ H → h(x) be the evaluation map and E * x be its adjoint operator. A common theorem in nonparametric regression, known as the representer theorem [START_REF] Schölkopf | A Generalized Representer Theorem[END_REF], states that any estimator based on (P1) is a finite expansion of local functions supported by the training examples. More formally, the functional part of any solution ( ĥ, b) of Problem (P1) admits a representation of the form ĥ = n i=1 E * x α i , where

α i ∈ R p for all i ∈ [n].
A proof of this statement based on Karush-Kuhn-Tucker (KKT) conditions can be found in Supplementary B, but many proofs already exist for kernel methods in the scalar (p = 1) and the vectorial (p > 1) cases [START_REF] Schölkopf | A Generalized Representer Theorem[END_REF][START_REF] Brouard | Input Output Kernel Regression: Supervised and Semi-Supervised Structured Output Prediction with Operator-Valued Kernels[END_REF]. The representer theorem states that ĥ can be parametrized by a finite number of α i ∈ R p (which boil down to be Lagrange multipliers), such that each α i is associated to a training point x i . Therefore, as soon as α i = 0, x i does not appear in ĥ (that is, x i is not a support vector). When many α i = 0, we say that ĥ is data sparse. This feature is very attractive since it makes possible to get rid of useless training points when computing ĥ(x), and speeds up the prediction task. The next section introduces a way to get such a sparsity.

-insensitive losses and the associated minimization problem

In the whole discussion, we will assume that the loss function is convex, lower semicontinuous and has a unique minimum at 0 with null value. Then, we develop the idea that data sparsity can be obtained by substituting the original loss function by a slightly different version , where is a kind of soft-thresholding of . More formally, let be a positive parameter and define the -insensitive the loss by:

∀ξ ∈ R p : (ξ) =    0 if ξ 2 ≤ inf d∈R p : d 2 =1 (ξ -d) otherwise. Name Loss (ξ) -insensitive loss (ξ) Conjugate loss (α) Unidimensional Absolute |ξ| max(0, |ξ| -) χ -1≤α≤1 (ξ ∈ R, α ∈ R) Squared 1 2 ξ 2 1 2 (max (0, |ξ| -)) 2 1 2 α 2 Pinball τ -I ξ<0 ξ    0 if |ξ| ≤ τ (ξ -) if ξ ≥ (τ -1)(ξ + ) if ξ ≤ - χ τ -1≤α≤τ
Squared pinball Put another way, the new loss value associated to a residue ξ is set to 0 if the magnitude of the residue is sufficiently small, or to the smallest loss value in a neighborhood centered at the original residue. Let us remark that, by the previous assumptions, is convex (see Supplementary B). To illustrate this definition, Table 1 as well as Figure 1 provide several examples of -insensitive losses, particularly for regression, QR and expectile regression. Before developing further these examples in the next section, we explain why using this -insensitive loss helps in getting data sparsity for the regressor.

1 2 τ -I ξ<0 ξ 2      0 if |ξ| ≤ τ 2 (ξ -) 2 if ξ ≥ 1-τ 2 (ξ + ) 2 if ξ ≤ - 1 2 |τ -I α<0 | -1 α 2 Multidimensional 1 -norm ξ 1 No closed-form found χ -1 α 1 (ξ ∈ R p , α ∈ R p ) 2 -norm 1 2 ξ 2 2 1 2 ξ -min( ξ 2 , ) ξ ξ 2 2 2 1 2 α 2 2 Multiple pinball p j=1 τ j -I ξ j <0 ξ j No closed-form found χ τ -1 α τ Multiple squared pinball 1 2 p j=1 τ j -I ξ j <0 ξ 2 j No closed-form found 1 2 p j=1 τ j -I α j <0 -1 α 2 j Table 1: Examples of -insensitive losses.
Proposition 1 (Dual optimization problem) A dual to the learning problem

minimize h∈H, b∈R p λ 2 h 2 H + 1 n n i=1 (y i -(h(x i ) + b)) (P2) is minimize ∀i∈[n], α i ∈R p 1 n n i=1 (α i ) + 1 2λn 2 n i,j=1 α i , E x i E * x j α j 2 - 1 n n i=1 α i , y i 2 + n n i=1 α i 2 s. t. n i=1 α i = 0, (P3) where : α ∈ R p → sup ξ∈R p α, ξ 2 -(ξ)
is the Fenchel-Legendre transform of . In addition, let ( αi ) i∈[n] be solutions to Problem (P3). Then, the function ĥ = 1 λn n i=1 E *

x αi ( ĥ ∈ H) is solution to Problem (P2) for a given intercept b ∈ R p .

We see that turning into an -insensitive loss leads to blending a sparsity regularization in the learning process: the 1 / 2 -norm n i=1 α i 2 , which is known to induce sparsity in the vectors α i (Bach et al., 2012). Even though a proof of this proposition is proposed in Supplementary B, we can grasp this phenomenon by remarking that turns out to be an infimal convolution [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF], noted , of two functions. Thus, its Fenchel-Legendre transform is the sum of Fenchel-Legendre transforms of each contribution [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]

: = χ • 2 ≤ = + χ • 2 ≤ = + • 2 .
The next section discusses several applications of -insensitive losses.

Examples of -insensitive losses

To illustrate the link between -insensitive losses and data sparsity, we focus on RKHSs H, based either on a scalar-valued kernel k : X ×X → R [START_REF] Steinwart | Support Vector Machines[END_REF] or on a matrix-valued kernel (MVK) K : X × X → R p×p [START_REF] Micchelli | On Learning Vector-Valued Functions[END_REF]. Therefore, the base functions

E * x α i , appearing in Proposition 1 become x ∈ X → α i k(x , x) ∈ R (α i ∈ R) and x ∈ X → K(x , x)α i ∈ R p (α i ∈ R p )
, respectively for the uni and multidimensional case.

Least absolute deviation

The notion of -insensitive loss is very well illustrated in SVR: let us study the problem of least absolute deviation, that is, solving (P1), where there is a single output (h : X → R,

y i ∈ R, b ∈ R, α i ∈ R)
and the loss to be minimized is (ξ) = |ξ|. In this situation, it is easy to check that boils down to: (ξ) = max(0, |ξ| -), which is the well known -insensitive loss considered in SVR [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]. Now, since (α) = 0 when -1 ≤ α ≤ 1 and +∞ otherwise, Problem (P3) becomes (up to normalizing the objective function by n):

minimize ∀i∈[n], α i ∈R 1 2λn n i,j=1 α i α j k(x i , x j ) - n i=1 α i y i + n i=1 |α i | s. t. ∀i ∈ [n], -1 ≤ α i ≤ 1, n i=1 α i = 0. Using that |α| = inf α + ,α -≥0 α=α + -α - {α + + α -}, we can introduce auxiliary variables α + i ≥ 0 and α - i ≥ 0, decouple α i in α + i -α - i and replace |α i | by α + i + α - i .
This turns the previous optimization problem into:

minimize ∀i∈[n], α + i , α - i ∈R 1 2λn n i,j=1 (α + i -α - i )(α + j -α - j )k(x i , x j ) - n i=1 (α + i -α - i )y i + n i=1 (α + i + α - i ) s. t. ∀i ∈ [n], 0 ≤ α + i ≤ 1, 0 ≤ α - i ≤ 1, n i=1 α + i -α - i = 0,
which is the well known dual of SVR. In this perspective, SVR is a data sparse version of least absolute deviation, for which the sparsity inducing norm appears in the dual.

Least mean squares

Besides SVR, Kernel Ridge Regression is also well known among kernel methods. In particular, for the general case of multivariate regression (y i ∈ R p ), we aim at minimizing the squared loss (

•) = 1 2 • 2 2 . Referring to Table 1, Problem (P3) then becomes: minimize ∀i∈[n], α i ∈R p 1 2 n i,j=1 α i , I p + 1 λn K(x i , x j ) α j 2 - n i=1 α i , y i 2 + n i=1 α i 2 s. t. n i=1 α i = 0, ( P4
)
where I p is the p-dimensional identity matrix. This is essentially similar to the dual of kernel ridge regression but with an extra 1 / 2 -norm on Lagrange multipliers.

As far as we know, data sparsity for multivariate kernel ridge regression has been first introduced by [START_REF] Lim | Operator-valued kernel-based vector autoregressive models for network inference[END_REF]. In this work, the authors assume first that the optimal predictor ĥ has a finite representation ĥ = 1 λn n i=1 K(•, x i )c i without relying on a representer theorem and regularize the vectors of weights c i associated to each point thanks to an 1 / 2 -norm. They also dismiss the intercept b for simplicity. The learning problem, as presented in [START_REF] Lim | Operator-valued kernel-based vector autoregressive models for network inference[END_REF] (with normalization chosen at purpose), is:

minimize h∈H, ∀i∈[n], c i ∈R p λ 2 h 2 H + 1 2n n i=1 y i -h(x i ) 2 2 + λn 2 n i=1 c i 2 s. t. h = 1 λn n i=1 K(•, x i )c i .
For the sake of readability, we write c = (c 1 , . . . , c n ) the vector of all weights and c 1 / 2 = n i=1 c i 2 . Let also K = (K(x i , x j )) 1≤i,j≤n be the symmetric positive semidefinite (PSD) kernel matrix. Fermat's rule indicates that a solution ĉ of the previous learning problem must verify K y -

1 λn Kĉ -ĉ ∈ ∂( • 1 / 2 )(ĉ).
To compare Lim et al.'s approach to the framework proposed in this paper, let us dismiss the intercept b, which means omitting the dual constraint n i=1 α i = 0 in (P4), and let α = ( α 1 , . . . , α n ) be a solution to (P4). Then Fermat's rule gives:

y -1 λn K α -α ∈ ∂( • 1 / 2 )( α).
It appears that both approaches induce data sparsity thanks to a structured norm and are equivalent up to a normalization by K.

Quantile regression

A slight generalization of least absolute deviation consists in changing the slope of the absolute loss (•) = | • |, asymmetrically around 0. For a quantile level τ ∈ (0, 1), the so called pinball loss is defined by: Such a loss is used to estimate a conditional quantile of the output random variable, instead of a conditional median such as in standard regression [START_REF] Koenker | Quantile Regression[END_REF].

∀ξ ∈ R, ρ τ (ξ) = (τ -I ξ<0 ) ξ = max (τ ξ, (τ -1)ξ) .
In order to improve prediction, some works propose to estimate and to predict simultaneously several quantiles, which is referred to as joint QR [START_REF] Takeuchi | Parametric Task Learning[END_REF]Sangnier et al., 2016). In this context, close to multivariate regression, the output space is R p and the loss is (ξ) = p j=1 ρ τ j (ξ j ), where τ is a vector of p quantile levels τ j ∈ (0, 1). Even though we do not have a closed-form expression for the corresponding -loss, Figure 2 provides graphical representations of .

Besides the loss described above, the multiple QR framework comes with y i ∈ R p , α i ∈ R p , b ∈ R p and h : X → R p , such that the j th component of the prediction value, f j (x) = h j (x) + b j , estimates the τ j -conditional quantile of the output random variable.

In a manner very similar to least absolute deviation, the Fenchel-Legendre transformation of the loss of interest is (α) = 0 when τ -1 α τ and +∞ otherwise. Therefore, Problem (P3) becomes:

minimize ∀i∈[n], α i ∈R p 1 2λn n i,j=1 α i , K(x i , x j )α j 2 - n i=1 α i , y i 2 + n i=1 α i 2 s. t. ∀i ∈ [n], τ -1 α i τ , n i=1 α i = 0, (P5)
where the sparsity inducing term is the 1 / 2 -norm n i=1 α i 2 (Bach et al., 2012). Very recently, [START_REF] Park | Quantile regression with an epsilon-insensitive loss in a reproducing kernel Hilbert space[END_REF] also proposed an -insensitive loss for single QR (that is, in the unidimensional case, when estimating only a single τ -conditional quantile), based on a generalization of the one used for SVR. The proposition made by [START_REF] Park | Quantile regression with an epsilon-insensitive loss in a reproducing kernel Hilbert space[END_REF] is:

PK (ξ) = max (0, ρ τ (ξ) -) = 0 if ρ τ (ξ) ≤ ρ τ (ξ) - otherwise.
That is different from the -loss proposed in this paper, as illustrated in Figure 1. As a comparison, the single QR dual in [START_REF] Park | Quantile regression with an epsilon-insensitive loss in a reproducing kernel Hilbert space[END_REF]) is similar to (P5) but with max α i τ , α i τ -1 instead of |α i |. Let us remark that, because of the asymmetry, it is difficult to get a sparsity intuition for this modified 1 -penalization.

The advantage of our approach over [START_REF] Park | Quantile regression with an epsilon-insensitive loss in a reproducing kernel Hilbert space[END_REF] is to provide a unifying view of -losses in the multidimensional setting, which matches SVR as a special case. Moreover, the possibility to extend Park and Kim's loss to joint QR is unclear and could not guarantee data sparsity: as explained in Remark 1 (Supplementary A), this goal is achieved thanks to the 1 / 2 -norm in the dual, which is equivalent to our proposition of -loss.

Expectile regression

As for QR, it is worth considering weighting the squared loss asymmetrically around 0. The resulting framework is called expectile regression [START_REF] Newey | Asymmetric Least Squares Estimation and Testing[END_REF], by analogy with quantiles and conditional expectation estimation with the squared loss. Despite the fact that the expectile regression estimator has no statistical interpretation, it can be seen as a smooth (and less robust) version of a conditional quantile. This explains that this topic, at least in the unidimensional context, attracts a growing interest in the learning community [START_REF] Farooq | An SVM-like Approach for Expectile Regression[END_REF][START_REF] Yang | Flexible Expectile Regression in Reproducing Kernel Hilbert Space[END_REF][START_REF] Farooq | Learning Rates for Kernel-Based Expectile Regression[END_REF].

Following the guiding principle of QR (Sangnier et al., 2016), we can imagine learning simultaneously several expectiles in a joint expectile regression framework. As far as we know, this multivariate setting did not appear in the literature yet. Therefore, we consider the multiple squared pinball loss:

∀ξ ∈ R, (ξ) = p j=1 ψ τ j (ξ j ), with ψ τ (ξ) = 1 2 |τ -I ξ<0 | ξ 2 .
Similarly to the squared loss, the Fenchel-Legendre transformation of is :

α ∈ R p → 1 2 α ∆(α)α
, where ∆(α) is the p×p diagonal matrix whose entries are τ j -

I α j <0 -1 1≤j≤p .
Therefore, a dual optimization problem of sparse joint expectile regression is:

minimize

∀i∈[n], α i ∈R p 1 2 n i,j=1 α i , ∆(α i ) + 1 λn K(x i , x j ) α j 2 - n i=1 α i , y i 2 + n i=1 α i 2 s. t. n i=1 α i = 0.
We see that the only difference compared to multivariate regression is the anisotropic ridge ∆(α i ). Moreover, compared to other kernelized approaches for expectiles estimation such as [START_REF] Farooq | An SVM-like Approach for Expectile Regression[END_REF][START_REF] Yang | Flexible Expectile Regression in Reproducing Kernel Hilbert Space[END_REF], this setting goes beyond both in its multivariate feature as well as its data sparsity.

Statistical guarantees

In this section, we focus on the excess of risk for learning with the proposed -insensitive loss.

As it is usual in statistical learning, we consider the constrained version of Problem (P2), which is legitimated by convexity of [START_REF] Tikhonov | Solutions of ill-posed problems[END_REF]:

minimize h∈H : h H ≤c b∈R p 1 n n i=1 (y i -(h(x i ) + b)) , ( P6
)
where c is a positive parameter linked to λ. This change of optimization problem is motivated by the fact that any solution to (P2) is also solution to (P6) for a well chosen parameter c.

Let (X, Y ) be the couple of random variables of interest, following an unknown joint distribution. We assume being provided with n independent and identically distributed (iid ) copies of (X, Y ), denoted (X i , Y i ) 1≤i≤n . Let us now fix an intercept b ∈ R p . F = {h(•) + b : h ∈ H, h H ≤ c} is a class of hypotheses. We denote R n (F) the Rademacher average of the class F ( Bartlett and Mendelson, 2002)

. Let f † ∈ arg min f ∈F E [ (Y -f (X))]
be the target function, where the expectation is computed jointly on X and Y . For instance, in the context of single QR, if we assume that the conditional quantile lives in F, then the target f † can be this conditional quantile (up to uniqueness). Moreover, let f ∈ arg min f ∈F

1 n n i=1 (Y i -f (X i )) be the empirical estimator of f † .
Theorem 2 (Generalization) Let us assume that the loss is L-Lipschitz (L > 0) and that residues are bounded: ∀f ∈ F, Y -f (X)

2 ≤ M almost surely (where M > 0). Then, for any δ ∈ (0, 1], with probability at least 1 -δ over the random sample

(X i , Y i ) 1≤i≤n : E Y -f (X) -E Y -f † (X) ≤ 2 √ 2LR n (F) + 2LM log(2/δ) 2n + L .
A proof of this theorem is given in Supplementary B. The bound is similar to what we usually observe [START_REF] Boucheron | Theory of Classification: a Survey of Some Recent Advances[END_REF]Mohri et al., 2012), but suffers from an extra term, which is linear in . The latter embodies the bias induced by soft-thresholding the original loss to induce data sparsity.

A primal-dual training algorithm

Description

For traditional kernel methods (real or vector-valued ones), pairwise coordinate descent as well as deterministic and randomized coordinate descents are popular and efficient training algorithms [START_REF] Platt | Fast training of support vector machines using sequential minimal optimization[END_REF][START_REF] Shalev-Shwartz | Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization[END_REF][START_REF] Minh | A Unifying Framework in Vector-valued Reproducing Kernel Hilbert Spaces for Manifold Regularization and Co-Regularized Multi-view Learning[END_REF]. However, few of current algorithms are able to handle multiple non-differentiable contributions in the objective value (such as the ones introduced by ) and multiple linear constraints (coming from considering an intercept in our regressor, which is mandatory for QR (Takeuchi et al., 2006;Sangnier et al., 2016)). For these reasons, we propose to use a randomized primaldual coordinate descent (PDCD) technique, introduced by [START_REF] Fercoq | A Coordinate Descent Primal-Dual Algorithm with Large Step Size and Possibly Non Separable Functions[END_REF] and utterly workable for the problem at hand. Moreover, PDCD has been proved favorably competitive with several state-of-the-art approaches [START_REF] Fercoq | A Coordinate Descent Primal-Dual Algorithm with Large Step Size and Possibly Non Separable Functions[END_REF].

The learning problem we are interested in is:

minimize

∀i∈[n], α i ∈R p s(α 1 , . . . , α n ) + n i=1 (α i ) differentiable + n i=1 α i 2 not differentiable + χ n i=1 α i =0 not differentiable , ( P7
)
where s(α 1 , . . . ,

α n ) = 1 2λn n i,j=1 α i , K(x i , x j )α j 2 -1 n n i=1 α i , y i 2 is a quadratic (differentiable) function.
Depending on the kind of regression (see in Table 1), the mapping may be differentiable or not. In any case, the objective function in (P7) can be written as the summation of three components: one differentiable and two non-differentiable. Given this three-term decomposition, PDCD (see Algorithm 1) dualizes the second nondifferentiable component and deploys a randomized block coordinate descent, in which each iteration involves the proximal operators of the first non-differential function and of the Fenchel-Legendre transformation of the second one. We can see in Algorithm 1 that PDCD uses dual variables θ ∈ (R p ) n (which are updated during the descent) and has two sets of parameters ν ∈ R n and µ ∈ R n , that verify ∀i ∈ [n]: µ i < 1 λ max,i +ν i , where λ max,i is the largest eigenvalue of K(x i , x i ). In practice, we keep the same parameters as in [START_REF] Fercoq | A Coordinate Descent Primal-Dual Algorithm with Large Step Size and Possibly Non Separable Functions[END_REF]: ν i = 10λ max,i and µ i equal to 0.95 times the bound.

Algorithm 1 Primal-Dual Coordinate Descent. Initialize α i , θ i ∈ R p (∀i ∈ [n]). repeat Choose i ∈ [n] uniformly at random. Set θ l ← proj 1 θ l + diag(ν)α l for all l ∈ [p]. Set d i ← ∇ α i s(α 1 , . . . , α n ) + 2θ i -θ i . Set α i ← prox µ i ( • 2 + ) (α i -µ i d i ). Update coordinate i: α i ← α i , θ i ← θ i ,
and keep other coordinates unchanged. until convergence PDCD is a workable algorithm for non-differentiable as soon as the proximal operator prox • 2 + can be computed (when is differentiable, it can be moved to s and there is only prox • 2 , which is easy to compute (Bach et al., 2012)). For the situations of interest in this paper, non-differentiable appear for multiple QR and take the form of the characteristic function of box constraints (see Table 1). Nevertheless, as far as we know, computing the proximal operator of the sum of the 2 -norm and box constraints has not been done yet. Therefore, the following section provides a way to compute it for box constraints intersecting both the negative and the positive orthants (proofs are in Supplementary B.4).

Proximal operator for multiple quantile regression

Let a and b be two vectors from R p with positive entries, and λ > 0. From now on, let us denote [•] b -a : R p → R p the clip operator, defined by:

∀y ∈ R p , ∀j ∈ [p] : [y] b -a j =    y j if -a j < y j < b j b j if y j ≥ b j -a j if y j ≤ -a j .
Lemma 3 Let y ∈ R p such that y 2 ≥ λ be a fixed vector and µ ∈

[0, 1]. The equation   1 + λ [µy] b -a 2    µ = 1 (1)
is well defined and has a solution µ ∈ [0, 1].

Proposition 4 ∀y ∈ R p , prox λ • 2 +χ -a • b (y) = [µy] b -a ,
where µ = 0 if y 2 ≤ λ and solution to Equation (1) otherwise.

Proposition 4 states that prox λ • 2 +χ -a • b can be computed as the composition of a clip and a scaling operator, for which the scaling factor can be easily obtained by a bisection of a Newton-Raphson method. As a straightforward consequence of Proposition 4, we can express the proximal operator of the sum of an 1 / 2 -norm and box constraints as the concatenation of prox λ • 2 +χ -a • b (see Supplementary B for a formal statement), which involves as many scaling factors as groups (n in this case).

Active set for multiple quantile regression

Learning sparse models gains a lot in identifying early the active points and optimizing only over them. This is a way to speed up learning, that is particularly topical [START_REF] Ndiaye | GAP Safe screening rules for sparse multi-task and multi-class models[END_REF][START_REF] Shibagaki | Simultaneous Safe Screening of Features and Samples in Doubly Sparse Modeling[END_REF][START_REF] Ndiaye | GAP Safe Screening Rules for Sparse-Group Lasso[END_REF]. Active points may have different meanings depending on the context (sparse regression, SVR, SVM) but the unifying concept is to detect optimizing variables for which the optimal value cannot be figured out beforehand.

In the context of QR (but this also holds true for SVR and SVM [START_REF] Shibagaki | Simultaneous Safe Screening of Features and Samples in Doubly Sparse Modeling[END_REF]), active points correspond to every data point x i for which the dual vector α i is neither null nor on the border of the box constraint. These active points can be identified thanks to optimality conditions.

Let f = ĥ + b be an optimal solution to (P2) and ( αi ) 1≤i≤n be optimal dual vectors to (P3). In the general setting, primal feasibility and stationarity from KKT conditions yield:

∀i ∈ [n] : y i -f (x i ) ∈ ∂ ( ) ( αi ), (2)
where, in the case of QR, =

• 2 + χ τ -1 • τ . Therefore, ∀α ∈ R p : ∂ ( ) (α) =            {d : d ∈ R p , d 2 ≤ } if α = 0 α 2 α if α = 0, τ -1 ≺ α ≺ τ d ∈ ∂ (χ τ -1 • τ ) (α) : d, α α 2 2 ≥ otherwise.
Consequently, optimality condition in Equation ( 2) indicates that:

1. y i -f (x i ) 2 < =⇒ αi = 0; 2. y i -f (x i ) 2 > =⇒ ∀j ∈ [p], ( αi ) j = τ j or ( αi ) j = τ j -1.
Therefore, for each situation (1. or 2.), if both conditions are fulfilled for the current estimates f and α i , the corresponding dual vector α i is put aside (at least temporally) and is not updated until optimality conditions are violated. In section 7.1, we will show that this strategy dramatically speed up the learning process when is large enough.

Numerical experiments

This section presents two numerical experiments in the context of multiple QR with data sparsity. The first experiment deals with training time. We compare an implementation of Algorithm 1 with an off-the-shelf solver and study the impact of the active set strategy.

The second experiment analyses the effect of on quantile prediction and on the number of support vectors (data points x i for which α i = 0). Following Sangnier et al. (2016), we use a matrix-valued kernel of the form K(x, x ) = k(x, x )B, where B = (exp(-γ(τ i -τ j ) 2 )) 1≤i,j≤p and k(x, x ) = exp(-xx 2 2 /2σ 2 ). In the first experiment (simulated data), γ = 1 and σ 2 = 1, while in the second one (real data) γ is set to 0.1 and σ is the 0.7-quantile of the pairwise distances of the training data {x i } 1≤i≤n . Quantile levels of interest are τ = (0.25, 0.5, 0.75).

Training time

We aim at comparing the execution time of four approaches for solving (P5): i) Algorithm 1 with active set (PDCD with AS ); ii) Algorithm 1 without active set (PDCD without AS ); iii) an off-the-shelf solver based on an interior-point method for quadratic cone programming (CVXOPT (CONEQP)) and iv) quadratic programming (CVXOPT (QP)) [START_REF] Anderson | CVXOPT: A Python package for convex optimization[END_REF]. For the last approach, we leverage a variational formulation of the sparsity constraint (Bach et al., 2012):

n i=1 α i 2 = inf ∀i∈[n], µ i ∈R + 1 2 n i=1 1 µ i α i 2 2 + µ i ,
and alternate minimization with respect to (α i ) 1≤i≤n and (µ i ) 1≤i≤n . By convexity and differentiability of the objective, alternate minimization converges to an optimal solution [START_REF] Rockafellar | Convex Analysis[END_REF].

We fix 1/(λn) = 10 2 and we use a synthetic dataset for which X ∈ [0, 1.5]. The target Y is computed as a sine curve at 1 Hz modulated by a sine envelope at 1/3 Hz and mean 1. Moreover, this pattern is distorted with a random Gaussian noise with mean 0 and a linearly decreasing standard deviation from 1.2 at X = 0 to 0.2 at X = 1.5.

Figure 3 depicts the dual gap achieved by the competitors with respect to the CPU time in different configurations ( ∈ {0.1, 1, 5} and sample size n ∈ {50, 500}). These curves are obtained by increasing the maximal number of iterations of each solver and recording primal and dual objective values attained along with the CPU time. Given a configuration ( and n fixed), all dual gaps are computed with respect to the lowest primal obtained. This comparison procedure is motivated by the fact that for multiple QR, we cannot compute the exact duality gap (we cannot compute the primal loss ), but only over-estimate it (see Remark 3 in Supplementary A). This explains why plots in Figure 3 do not necessarily converge to 0 (we conjecture that when a method faces no progress, the exact dual gap is almost null). Overall, the curves portray the behavior of the dual objective function rather than the exact duality gap. For small-scale data and small values of , CVXOPT is more efficient than our approach, but its supremacy caves in when increasing and/or the sample size n. The impact of the problem size n on the comparison between interior-point methods (such as CVXOPT) and first-order ones (such as PDCD) is known for a long time, so this first finding is comforting. Also, increasing the sparsity regularizer requires more and more calls to the alternate procedure CVXOPT (QP), so the behavior with respect to was expected.

Analyzing the active-set strategy, we remark that it is never a drawback but becomes a real asset for hard sparsity constraints. This comes from many αi to be null and thus put aside. Yet, it appears that, contrarily to the unidimensional case, only a few α i are on the border of the box constraint. Therefore, the active-set strategy provides only a limited (but real) gain compared to the vanilla version of PDCD for small values of .

Data sparsity

This section investigates the actual data sparsity introduced by the proposed -insensitive loss, along with its impact on the generalization error. For this purpose, we consider the datasets used in (Takeuchi et al., 2006;[START_REF] Zhang | On Quantile Regression in Reproducing Kernel Hilbert Spaces with the Data Sparsity Constraint[END_REF], coming from the UCI repository and three R packages: quantreg, alr3 and MASS. The sample sizes n vary from 38 (CobarOre) to 1375 (heights) and the numbers of explanatory variables vary from 1 (6 sets) to 12 (BostonHousing). The datasets are standardized coordinate-wise to have zero mean and unit variance. In addition, the generalization error is estimated by the mean over 10 trials of the empirical loss

1 m m i=1 (y test i -f (x test i )
) computed on a test set. For each trial, the whole dataset is randomly split in a train and a test set with ratio 0.7-0.3. The parameter 1/(λn) is chosen by a 3-fold cross-validation (minimizing the pinball loss) on a logarithmic scale between 10 -3 and 10 3 (10 values). At last, a training point x i is considered a support vector if α i 2 /p > 10 -3 .

Table 2 reports the average empirical test loss (scaled by 100) along with the percentage of support vectors (standard deviations appear in Supplementary C). For each dataset, the bold-face numbers are the two lowest losses, to be compared to the loss for = 0. We observe first that increasing does promote data sparsity. Second, the empirical loss tends to increase with , as expected, but we do not lose much by requiring sparsity (we can even gain a little).

Conclusion

This paper introduces a novel and unifying framework concerning -insensitive losses, which offers a variety of opportunities for traditional, quantile and expectile regression, for uni and multivariate outcomes. Estimators are based on vector-valued RKHSs [START_REF] Micchelli | On Learning Vector-Valued Functions[END_REF] and benefit from theoretical guarantees concerning generalization as well as an efficient learning algorithm. Future directions of research include improving further the learning procedure and deriving sharper generalization bounds.

where T i : R p → R p are any transformations [START_REF] Alvarez | Kernels for Vector-Valued Functions: a Review[END_REF].

Remark 3 As it is standard for coordinate descent methods, our implementation uses efficient updates for the computation of both n j=1 K(x i , x j )α j and θ l . In addition, convergence of Algorithm 1 can be assessed by duality gap (objective of (P2) minus objective of (P3) in the body of the text). Yet, even though we do not have closed-form expressions for primal loss , the duality gap can be over-estimated by upper-bounding in the following manner:

∀ξ ∈ R p : (ξ) ≤ ξ 1 - min( , ξ 2 ) ξ 2 .
This is true since

min( , ξ 2 ) ξ 2 ξ 2 ≤ .
Remark 4 Contrarily to QR, expectile regression involves a differentiable mapping . Consequently, it can be easily incorporated to the quadratic contribution of (P7) (see body of the text). Nevertheless, it can also be considered jointly with • 2 , in the same manner as for QR. In this case, the differentiable part remains the same for expectile and quantile regression, only the non-differentiable part changes. The proximal operator needed is given in the following proposition.

Proposition 5 Let ψ : y ∈ R p → 1 2 p i=1 τ j -I y j <0 -1 y 2 j . Then ∀y ∈ R p , ∀j ∈ [p] prox λ • 2 +ψ (y) j = 1 + λ µ + λ τ j -I y j <0 -1 -1 y j , if y 2 > λ,
where µ > 0 is solution to:

p j=1 y 2 j µ 1 + λ τ j -I y j <0 -1 + λ 2 = 1, (1)
(such a solution exists) and

prox λ • Since is convex, (ξ, u) → (ξ -u) + χ u 2
≤ is jointly convex with respect to ξ and u. Therefore, ˜ is convex as the coordinate infimum of a jointly convex function [START_REF] Boyd | Convex Optimization[END_REF].

Let us now show that ˜ = . First, since for any ξ, Slater's constraint qualification are satisfied for (2), strong duality holds, that is,

∀ξ ∈ R p , ∃λ ≥ 0 : ˜ (ξ) = inf u∈R p (ξ -u) + λ u 2 -λ ,
and thanks to the lower semi-continuity of the objective, the infimum is attained at, let us say, û. Then, when ξ 2 ≤ , we can chose û = ξ and we get ˜ (ξ) = 0, which is the infimum of . On the other hand, when ξ 2 > , let us consider the Karush-Kuhn-Tucker (KKT) conditions. By complementary slackness, either λ = 0 and û 2 ≤ , or û 2 = . In the first situation (λ = 0 and û 2 ≤ ), ˜ (ξ) = inf u∈R p (ξ -u) = (ξ -û) = (0) = 0 and û = ξ (by uniqueness of the minimizer of ). Thus, ξ 2 ≤ , which is contradictory. Consequently, we have necessarily, û 2 = . To summarize:

∀ξ ∈ R p : ˜ (ξ) =    0 if ξ 2 ≤ inf u∈R p : u 2 = (ξ -u) otherwise,
which is exactly the definition of .

B.2. Dual and representer theorem

Since is convex and can be replaced by (2), Problem (P2) from the body of the text can be reformulated in (Lagrange multipliers are indicated on the right): Bach et al., 2012), a dual to Problem (P1) is:

minimize h∈H, b∈R p , ∀i∈[n], ξ i ∈R p , r i ∈R p λ 2 h 2 H + 1 n n i=1 (ξ i ) s. t.          ∀i ∈ [n], y i -(h(x i ) + b) n = r i + ξ i n : α i ∈ R p r i 2 2 2 n ≤ 2n : µ i ∈ R + . (P1) n i=1 α i 2 (
minimize

∀i∈[n],α i ∈R p 1 n n i=1 (α i ) + 1 2λn 2 n i,j=1 α i , E x i E * x j α j 2 - 1 n n i=1 α i , y i 2 + n n i=1 α i 2 s. t. n i=1 α i = 0. (P2) B.3. Generalization Let P : f ∈ F → E [ (Y -f (X))] and P n : f ∈ F → 1 n n i=1
(Y -f (X)), as well as respective twins P and P n, obtained by substituting to . Let us decompose P f -P f † :

P f -P f † = P f -P n f + P n f -P n f † + P n f † -P f † .
First, by concentration inequalities (Bartlett and Mendelson, 2002;[START_REF] Maurer | A vector-contraction inequality for Rademacher complexities[END_REF]Sangnier et al., 2016), we have, with probability greater that 1 -δ:

P f -P n f ≤ sup f ∈F (P f -P n f ) ≤ 2 √ 2LR n (F) + LM log(1/δ) 2n .
Second, let us decompose P n f -P n f † :

P n f -P n f † = P n f -P n, f + P n, f -P n, f † + P n, f † -P n f † .
By Lipschitz continuity, we have:

∀ξ, u ∈ R p , u 2 ≤ : (ξ) -(ξ -u) ≤ L ξ -(ξ -u) 2 ≤ L .
Consequently, (ξ) -(ξ) ≤ L and P n f -P n, f ≤ L . In addition P n, f -P n, f † ≤ 0 since f is a minimizer of P n, over F, and f † ∈ F. Finally, P n, f † -P n f † ≤ 0 since upper bounds . To summarize the second point, P n f -P n f † ≤ L . Third and last, by Hoeffding's inequality [START_REF] Boucheron | Concentration Inequalities: A Nonasymptotic Theory of Independence[END_REF], with probability at least 1 -δ:

P n f † -P f † ≤ LM log(1/δ) 2n .
Gathering these three points with a union bound concludes the proof.

B.4. Algorithms

Proof (Lemma 3, body of the text) Let φ Proof (Proposition 4, body of the text) The proof is in two part. First, we write optimality conditions for the proximal operator of interest, then we show that [µy] b -a satisfies these optimality conditions when µ is appropriately defined. From now on, let y ∈ R p .

: µ ∈ [0, 1] → 1 + λ [µy] b -a 2 µ. First, φ(1) = 1 + λ [y] b

Optimality conditions

Let x = prox λ • 2 +χ -a • b (y) = arg min -a x b λ x 2 + 1 2 y -x 2 2 .
1. Assume that x = 0. Then, λ • 2 + 1 2 y -• 2 2 is differentiable at x and for each coordinate j ∈ [p], either:

(a) -a j < x j < b j and 1 + λ Gathering Conditions 1a-1c gives x 2 + λ ≤ y 2 . Since x = 0, we get λ < y 2 . Conversely, if y 2 ≤ λ, then x = 0.

2. If x = 0, then ∀δ > 0 such that -a δy b, λ δy 2 + 1 2 y -δy 2 2 ≥ 1 2 y 2 2 , that is λ y 2 ≥ 1 -δ 2 y 2 2 . Thus, by continuity when δ ↓ 0, we have λ ≥ y 2 . To sum up, x = 0 if and only if y 2 ≤ λ.

Proximal solution Let x = [µy] b -a , where µ is defined in Proposition 4 from the body of the text. Assume that y 2 ≤ λ, then µ = 0 and x = 0 satisfies the optimality conditions. On the other hand, if y 2 > λ, then 1 + λ

x 2 µ = 1. As a result, either:

1. -a j < x j < b j , so necessarily x j = µy j (otherwise it would be clipped to b j or -a j ).

Therefore 1 + λ x j ≤ y j ;

3. or x j = -a j , meaning that µy j ≤ -a j . So 1 + λ x 2 (µy j ) ≤ 1 + λ

x 2 (-a j ), that is 1 + λ x 2

x j ≥ y j .

Thus, when y 2 > λ, x satisfies the optimality conditions. This concludes the proof.

Corollary 6 Let two n-tuples A = (a 1 , . . . , a n ) and B = (b 1 , . . . , b n ) of vectors from R p with positive entries. For any n-tuple Y = (y 1 , . . . , y n ) of vectors from R p , let:

(x 1 , . . . , x n ) = prox λ • 1 / 2 +χ -A • B (Y ),
where Y 1 / 2 = n i=1 y i 2 . Then, ∀i ∈ [n]:

x i = prox λ • 2 +χ -a i • b i (y i ).
Proof This is a direct consequence of the separability of λ • 1 / 2 + χ -A • B .

Proof (Proposition 5) The proof is similar to the one for Proposition 4 (see body of the text). Let y ∈ R p .

Optimality conditions Let x = prox λ • 2 +ψ (y) = arg min x∈R p λ x 2 + λψ(x) + 1 2 yx 2 2 . 1. Assume that x = 0. Then, λ • 2 + λψ + 1 2 y -• 2 2 is differentiable at x and for each coordinate j ∈ [p]:

y j = λ x 2 + λ τ j -I x j <0 -1 + 1 x j .
It appears that x j and y j have same sign. Therefore, I x j <0 = I y j <0 and

x j = λ x 2 + λ τ j -I y j <0 -1 + 1 -1 y j .
Now, the previous relation implies: that is λ < y 2 . Conversely, if y 2 ≤ λ, then x = 0.

2. If x = 0, then ∀δ > 0, λ δy 2 + λψ(δy) + 1 2 y -δy 2 2 ≥ 1 2 y 2 2 , that is λ( y 2 + δψ(y)) ≥ 1 -δ 2 y 2 2 . Thus, by continuity when δ ↓ 0, we have λ ≥ y 2 . To sum up, x = 0 if and only if y 2 ≤ λ. Proximal solution If y 2 ≤ λ, then x = 0 is satisfies trivially the optimality conditions. On the other hand, if y 2 > λ, then p Let µ be such a solution and let x ∈ R p such that for each coordinate j ∈ [p],

x j = 1 + λ µ + λ τ j -I y j <0

-1 -1 y j .

Then:

x 2 and x satisfies the optimality conditions. This concludes the proof.

Table 1: Empirical pinball loss ×100 along with percentage of support vectors (the less, the better).

Data set = 0 = 0.05 = 0.1 = 0.2 = 0.5 = 0.75 = 1 = 1.5 = 2 = 3

Figure 1 :

 1 Figure 1: Examples of unidimensional -losses ( = 0.2, τ = 0.25).

Figure 2 :

 2 Figure 2: Multiple -pinball loss and slices of it ( = 0.2, τ = (0.25, 0.6)).

Figure 3 :

 3 Figure 3: Training time (synthetic dataset).

  Second, for µ ≥ 0 sufficiently close to 0, [µy] b -a = µy (since entries of a and b are positive). Therefore φ(0) = lim µ↓0 µ + φ is a continuous mapping on [0, 1] and 1 ∈ [φ(0), φ(1)], then the equation φ(µ) = 1 has a solution in [0, 1].

  or x j = b j and 1 + λ x 2 x j ≤ y j ;(c) or x j = -a j and 1 + λ x 2 x j ≥ y j .

(

  µy j ) = y j ; 2. or x j = b j , meaning that µy j ≥ b j . So 1 + λ x 2 (µy j ) ≥ 1 + λ x 2 b j , that is 1 + λ x 2

  , by continuity, Equation 1 has a solution µ > 0.

  • H is the norm associated to H. ), where the residues ξ i = y i -(h(x i )+ b) are driven close to zero. Here, the prediction function (or regressor) is f :x ∈ X → h(x) + b. Itis made of a functional part h ∈ H and of an intercept b ∈ R p (such as what we can find in SVMs). The last component of Problem (P1), λ 2 • 2 H , is a regularizer, which penalizes functions with high complexities, allowing for generalizing on unknown observations

		Problem (P1)
	consists in minimizing a data-fitting term 1 n	n i=1 (ξ i

Table 2 :

 2 Empirical (test) pinball loss ×100 (percentage of support vectors in parentheses).

	Data set	= 0	= 0.05	= 0.1	= 0.2	= 0.5	= 0.75	= 1	= 1.5	= 2	= 3
	caution	67.50 (100) 67.40 (99) 67.17 (95) 67.54 (87)	69.93 (65) 73.24 (49) 76.80 (39) 83.42 (24) 100.22 (17) 142.19 (6)
	ftcollinssnow	109.07 (100) 109.12 (100) 109.14 (100) 109.15 (100) 109.11 (100) 110.39 (98) 109.05 (82) 110.90 (53) 109.81 (34) 113.50 (8)
	highway	79.29 (100) 78.10 (97)	76.75 (97)	76.66 (91) 75.09 (64) 70.94 (51) 75.10 (37) 97.67 (31) 112.30 (0) 112.09 (1)
	heights	91.05 (100) 91.00 (100) 90.98 (100) 90.98 (100) 91.18 (99) 91.21 (86) 90.98 (67) 91.09 (40) 91.51 (23) 93.34 (5)
	sniffer	32.34 (100) 31.40 (99) 32.31 (97) 31.40 (84) 34.64 (43) 39.84 (24) 41.82 (12) 52.06 (6)	62.21 (4) 103.76 (3)
	snowgeese	49.62 (100) 50.51 (98) 51.25 (87)	51.08 (70)	52.88 (37) 53.81 (23) 62.81 (17) 90.15 (15) 107.53 (14) 94.25 (1)
	ufc	57.87 (100) 57.90 (100) 57.78 (100) 57.84 (100) 57.67 (78) 57.84 (52) 58.19 (35) 61.04 (15) 66.81 (5)	86.23 (2)
	birthwt	99.93 (100) 99.95 (100) 99.93 (100) 99.70 (100) 99.25 (100) 100.50 (87) 99.80 (67) 98.71 (43) 99.56 (24) 103.39 (6)
	crabs	8.59 (100)	8.52 (91)	8.49 (55)	9.44 (16)	19.94 (6)	23.08 (2)	31.44 (3)	44.08 (2)	53.45 (1)	86.91 (2)
	GAGurine	44.30 (100) 44.26 (99) 44.25 (99) 44.86 (87)	46.20 (50) 49.87 (33) 52.88 (22) 57.06 (12) 65.89 (6) 103.32 (2)
	geyser	77.81 (100) 78.15 (100) 78.12 (100) 78.45 (100) 78.40 (92) 78.28 (78) 78.54 (59) 80.55 (32) 85.15 (16) 99.92 (1)
	gilgais	32.96 (100) 33.12 (99) 33.27 (96)	33.42 (81)	35.08 (43) 36.62 (25) 37.94 (14) 48.17 (5)	94.65 (7) 104.12 (0)
	topo	47.49 (100) 48.93 (100) 48.74 (98)	48.17 (94) 41.65 (57) 45.24 (38) 51.19 (26) 53.68 (16) 58.21 (12) 80.57 (6)
	BostonHousing 34.54 (100) 34.68 (99)	34.70 (97) 34.09 (80) 35.27 (35) 37.65 (20) 41.31 (13) 55.04 (7)	73.39 (7) 112.22 (12)
	CobarOre	0.50 (100)	5.05 (38)	8.75 (36)	12.47 (26)	23.84 (18) 35.82 (17) 47.35 (15) 66.15 (14) 84.51 (12) 106.89 (6)
	engel	43.57 (100) 43.50 (100) 43.47 (99) 43.44 (89) 57.36 (39) 43.98 (19) 46.31 (11) 53.15 (4)	69.43 (5) 100.48 (0)
	mcycle	63.95 (100) 63.88 (99) 64.26 (99)	64.90 (98)	65.89 (88) 67.29 (70) 70.11 (51) 74.78 (26) 86.49 (14) 109.79 (2)
	BigMac2003	49.94 (100) 49.97 (98) 50.00 (96)	50.27 (85)	51.16 (56) 51.44 (36) 53.63 (28) 77.40 (18) 106.38 (14) 136.76 (4)
	UN3	71.27 (100) 70.94 (100) 71.03 (100) 71.49 (99)	71.37 (87) 71.53 (65) 72.68 (50) 76.72 (27) 84.50 (13) 109.59 (0)
	cpus	11.31 (100) 13.32 (28) 15.57 (21)	20.16 (15)	25.88 (8)	35.66 (6)	55.27 (6)	65.05 (0)	65.05 (0)	65.02 (0)

c 2017 M. Sangnier, O. Fercoq & F. d'Alché-Buc.

d 1≤i,j≤p (transformable), c 2017 M. Sangnier, O. Fercoq & F. d'Alché-Buc.

+ψ (y) = 0 if y 2 ≤ λ,.Similarly to Equation 1 in the body of the text, the scaling factor µ in Proposition 5 can be easily obtained by a bisection of a Newton-Raphson method.

+ µ i =

Data sparse nonparametric regression with -insensitive losses: supplementary material

Appendix A. Remarks

Remark 1 It is easy to see that, for a unidimensional loss , is obtained in the following manner:

Consequently, when a multivariate loss is separable, that is (ξ) = p j=1 (j) (ξ j ) (for some unidimensional losses (j) ), it is tempting to consider each component separately and to define = p j=1 (j) . Basically, this boils down to replacing • 2 by • ∞ in the general -loss introduced in this paper.

However, this is not a good idea since this definition would result in adding an 1 -norm n i=1 α i 1 instead of an 1 / 2 -norm in the dual. As a consequence, we would obtain sparse vectors α i , which is not the data sparsity we pursue since α i could have null components but could be different from 0, forcing us to keep the points x i for prediction.

Remark 2 In the body of the text, omitting the intercept b in Problem (P2) comes down to removing the linear constraint in Problem (P3). This practice is common for support vector regression (SVR) with a Gaussian kernel, but is excluded for quantile regression (QR) (Takeuchi et al., 2006;Sangnier et al., 2016).

Example 1 Examples of scalar and matrix kernels are:

where d > 0 is the degree (Mohri et al., 2012), and

Appendix B. Technical details

B.1. Convexity and redefinition of

For the sake of simplicity, let us first define:

Let us compute a dual to Problem (P1). The Lagrangian reads:

The objective function of the dual problem to (P1) is obtained by minimizing the Lagrangian with respect to the primal variables h, b, (ξ i ) 1≤i≤n and (r i ) 1≤i≤n . For this purpose, let us remark that minimizing on ξ i boils down to introducing the Fenchel-Legendre transform of : :

. Thus, it remains to compute:

Since H is unbounded in all directions, the minimum of L with respect to h, b and (r i ) n i=1 is obtained by setting the gradients to 0, which leads to h = 1

. Thus, the dual objective reads:

Then, the dual optimization problem consists in maximizing L D subject to the constraints

Appendix C. Numerical experiments

Table 1 reports the average empirical loss (scaled by 100) along with the standard deviations. It completes Talbe 2 from the body of the text. For each dataset, the bold-face numbers are the two lowest losses. These values should be compared to the loss for = 0.


Pascal Ballet
email: pascal.ballet@univ-brest.fr

SimCells, an advanced software for multicellular modeling Application to tumoral and blood vessel co-development

New biomedical advances at cellular level give the possibility to develop more and more accurate computational models of cells. Moreover, the increasing power of graphical processor units allows the simulation of millions of interacting virtual cells. This paper summarizes efforts made to create multicellular simulators, their attended benefits and their inevitable drawbacks. In particular, it presents the new software SimCells that simulate the dynamics of multicellular systems using a graphically programmable multiagent system. A fully functional example of the tumoral and blood vessel codevelopment is also detailed.

Introduction

The biomedical field at cellular scale is a very active field of research, dealing with data from simple molecules up to entire organisms. It works both on the fundamental comprehension of the living and on the design of treatments against complex diseases. Although numerous entities are in interaction at very different scales and inside heterogeneous environment in cellular systems, many data coming from biological experiments can be integrated in computational systems to explain how the real system works. The creation of such systems is made through a long and difficult process entailing interdisciplinary knowledge, computer science and mathematical skills, rigor and creativity. It also demands high computational power which can be found in massively parallel processors. Creativity is a hard and fragile dynamic process, especially during the development of models on complex living systems which usually are non-linear and counter-intuitive. Advanced software can contribute to this creativity by giving intuitive and computational tools. Graphical designing tools enable researchers to focus on their models instead of how they implement it. Computational tools resolve the models and verify their validity. Moreover, they can organize huge amount of data inside formal paradigms which can be easily understood and shared. The modelling and simulation of cellular systems, both at molecular level and at multicellular scale has already been initiated in the early ages of computers with the works of Alan Turing on the chemical basis of morphogenesis and those of John Von Neumann on cellular automata. Moreover, computer aided design has been used since the 1950's in numerous fields of science and technology like mathematics, molecular modelling, car or plane design. Although computer scientists have been interested in modelling multicellular systems, this field has encountered many difficulties both theoretical, due to partial cell comprehension, and practical, with the limitation of computational power. Despite the development of many fundamental algorithms and programs for the modelling of multicellular systems, most of them have been limited to prototypes and dedicated to specific problems. This work presents the recent efforts to develop fully integrated software, that help biologists to create models of multicellular systems. An attention is given to systems with many interacting entities capable of self-assembly. This paper begins with some objectives, expected benefits and the inevitable drawbacks of such software. Then the contemporary contexts of cellular biology, computer science and computational biology are presented. Next, different algorithmic approaches are shown. Afterwards, two advanced software are described. Finally, a SimCells simulation of the angiogenesis process is detailed.

Objectives

The comprehension of living cell is still incomplete and computational approaches can participate in their study. This help can be divided into four short term and two long term objectives.

Short-term objectives

1. simulation of complex cellular systems to improve the understanding and control of those systems 2. formalization of knowledge which can reduce ambiguities inside biological models 3. integration of diverse scientific knowledge to improve the realism of models 4. dissemination of knowledge which could accelerate the rhythm of discoveries by sharing models and simulations 5. the training of students in biology, enabling them to deal with complex systems.

Long-term objectives

1. computer aided design of tissues and organs for medical applications to contribute on their replacement or restoration 2. computer aided design of tissues and organisms for industrial applications to create materials, systems or computers based on multicellular systems 1. These objectives lead to benefits that irrigate science and technology even before their completion. They are briefly described in the following section.

Benefits

The study of multicellular systems using computers is interesting both scientifically and technically and it can also have important implications in society.

Scientifically, it serves the comprehension of fundamental mechanisms governing cells, from single-cell behaviors to embryogenesis of a whole organism, including tissue formation and organs morphogenesis.

Technically, it requires enhancement of computational models, from algorithms to graphical user interface. In return, it can improve the field of computer science.

As regards the implications in society, computation of multicellular systems could help medical researchers to discover new medical treatments. It could also be used, in some cases, as an alternative to animal testing both for research and teaching.

Moreover, the creation of computational models imitating multicellular systems involves many areas of science such as biology and medicine, computer science and computational biology, physics and complex systems, mathematics and pedagogy and may lead to think how to reorganize scientific fields by encouraging interdisciplinary that induces discoveries.

Despite these benefits, it is necessary to clarify the limits of this approach in order not to over-estimate its potential.

Drawbacks, limitations and difficulties

Modelling and simulation of multicellular systems have to deal with many problems. Some are shown below.

Modelling problems

Any model contains approximations, bias and gaps that inevitably denatures the real system. For example, a model contains what is known on a real system but not what is unknown2: in fine, invivo or in-vitro experiments could be necessary. Even what is supposed to be known must be carefully questioned because data coming from experiments are of various kinds, have different experimental conditions and possess margin of error. A model is only valid under certain conditions. It generally has limited predictive capabilities, especially with non-linear systems. Furthermore, a biological system can be described by two or more different models. Living systems are non-optimal, so their models cannot be discriminated by simply using the principle of parsimony. Models can contain false or obsolete knowledge and can slow down the discoveries by propagating nonrelevant data, mechanisms and behaviors. Creating a computational model of multicellular systems requires either a modeler to have an interdisciplinary expertise or specialists to work closely together. In both cases, the time needed for the construction of the model is greatly increased. A model can be difficult to use because it has no simple algebraic solution, it does not fit any existing computational approaches or because it demands too much computation time. More problems related to simulation are described in the next section.

Simulation problems

Multicellular systems are intrinsically multiscale. Spatial scales start from the size of water molecules, 3 x 10 -10 meter up to the size of Sequoia trees that can nearly reach 10 2 meters [START_REF] Flint | To Find the Biggest Tree[END_REF]. Between a molecule of water and a giant sequoia, there are about 12 orders of magnitude which cannot be simulated yet and probably not for a while. In terms of time, the same problem occurs. The time of a water molecule to reorganize with its neighborhood molecules is about 10 -12 second and the oldest living tree is close to 5 x 10 3 years old [START_REF] Brown | OLDLIST: A Database of Maximum Tree Ages[END_REF], which corresponds to 15 orders of magnitude. Again, this cannot be simulated using today's computers. During a simulation, calculation approximations are possible, especially when working with floatingpoint numbers because of truncation. A simulation is based on an implementation of a model and it does not always reflect accurately the model itself. Topology and geometry of computational objects and environments strongly affect simulation results. For example, a cell can be modelled as a solid sphere in a continuous environment or as a cube with integer coordinates and size inside a 3D matrix. An underestimated problem concerns the use of advanced software that may reduce the critical mind of those who use them. It should be remembered that models and simulations are only approximations of the real system, especially for students. Nevertheless, benefits are more important than drawbacks and computational models and simulations contribute to the construction of knowledge in a broad sense.

Properties and control of cells

Living cells are usually seen as the structural and functional basic unit of life. Encyclopaedia Britannica defines the cell as the basic membrane-bound unit that contains the fundamental molecules of life and of which all living things are composed [START_REF] Alberts | Cell[END_REF]. All known living organisms, in their huge diversity 1, are made of cells, from unicellular organisms to humans 2. Since their discovery by Robert Hooke in the middle of the 17 th century 3, they have been one of the most studied scientific objects. Cells have numerous forms 4, functions 5, behaviors 6, interactions 7 and organizations 8 but are based on very similar structures and mechanisms 9. In terms of evolution, the first cell is dated 3.5 billion years old [START_REF] Schopf | Laser-Raman imagery of Earth's earliest fossils[END_REF]. Multicellular systems are more recent and the oldest one known is 2.1 billion-year-old [START_REF] Albani | Large colonial organisms with coordinated growth in oxygenated environments 2.1 Gyear ago[END_REF].

A multicellular system is made of interacting cells that can express numerous shapes and functions. Furthermore, many scales can be observed, starting from tissues that involve several cells, to organs which realize essential functionalities, up-to whole organisms that constitute societies and ecosystems. Every known life kingdom possesses multicellular structures: animal cells (for instance caenorhabditis elegans), plant cells (oak), bacteria (Nostoc for example), fungal cells (Penicillium), protista (without specialized tissues like algae Seaweeds) or archaea (see Methanosarcina acetivorans for instance). Most of them also can perform cellular recognition, a key mechanism that enables the distinction between self and non-self, even on primitive organisms like Porifera [START_REF] Dammer | Binding Strength Between Cell Adhesion Proteoglycans Measured by Atomic Force Microscopy[END_REF] [START_REF] Misevic | Molecular Recognition between Glyconectins as an Adhesion Selfassembly Pathway to Multicellularity[END_REF]. When studying such living structures, many properties can be observed that are biological, mechanical, chemical, thermodynamic, morphological, etc. Therefore, multicellular systems can be modelled with many different theories or approaches. Each of them corresponds to a specific point of view on real multicellular systems which are, by now, too complex to be fully understandable with a single approach. Numerous and fundamental advances in the understanding and control of cells have been made recently, both at molecular and cellular levels. In 2002, Sydney Brenner, H. Robert Horvitz and John E. Sulston obtained the Nobel Prize for their work on genetic regulation of organ development and programmed cell death, describing seminal mechanisms occurring during organism development. Another important work, made in 2006, involves Shinya Yamanaka's team who has artificially generated a pluripotent stem cell from an adult somatic cell. Qualitative and quantitative improvements have been made to this work since this date. About two years later, in 2008, Doris Taylor and her research team built the first beating bioartificial heart using a tissue scaffold from a rat heart and heart cells from newborn rats. More recently, in 2010, Thomas Vierbuchen et al. were able to directly convert fibroblasts into functional neurons, without passing through the stage of stem cell. Also, in 2010, Craig Venter and his team created a bacterial cell controlled by a chemically synthesized genome, that although controversial, is paving the way for the construction of artificial organisms. In 2011, Mototsugu Eiraku et al. showed that embryonic stem cells derived into retinal epithelium self-organize to spontaneously form hemispherical epithelial vesicles corresponding to the optic cup. Besides the fact that such researches lead to a better understanding of cells, they also enable the creation of tissues and organs able to replace those of patients from their own cells, avoiding immune rejection and replacing cells that normally do not easily regenerate by themselves like neurons, retina or cardiac cells. The engineering of tissues and organs is still in its infancy 10, but the examples we have at our disposal show that they are about to revolutionize medicine and probably many other areas.

A key point we can extract from works in biology is that some mechanisms governing multicellular structures are based on rules that can be written in the form of algorithm. Moreover, these mechanisms, though often based on microscopic stochastic phenomena, are reproducible, giving them a strong macroscopic determinism. Given these data, we see that algorithms can be developed to reproduce and therefore to explore some of these phenomena using computers. This is facilitated by recent advances in multicore processors whose computing power overcomes the limitations due to stagnating clock frequencies of processors. Recent advances in quantum computing suggest that it could be possible to drastically improve the speed of some algorithms [START_REF] Mosca | Quantum computer algorithms[END_REF]. Finally, let us note the importance of philosophy and epistemology at the border of biology and computation. These includes the works of Maturana and Varela [START_REF] Maturana | Autopoiesis and cognition: the realization of the living[END_REF] on the autopoietic principle, or those of Jean-Jacques Kupiec with its ontophylogenetics theory [START_REF]The Origin of Individuals[END_REF] which gave useful formalisms to deal with multicellular systems.

Computer science context

Since the beginning of computers in the 1940s, scientists have tried to reproduce living cells in computers. One of the first try has been made by John von Neumann with the concept of cellular automata 11. This has been completed by the reaction-diffusion approach introduced by Alan Turing in 1952. Although computer power was limited at that time, the concepts and their results launched new scientific fields like artificial life and computational biology. This led to many developments, more and more precise, following the fast evolution of computing power. Two of the most interesting research areas for our work is the multiagent approach [START_REF] Minsky | The society of mind[END_REF] where multiple interacting agents evolve in virtual environments and where different models can be executed inside a same simulation. The second interesting approach is called artificial life [START_REF] Langton | Artificial Life: a review[END_REF]. It aims to imitate living things, and furthermore to explore life-as-it-mightbe. Advances in software architecture, algorithms and graphical user interface are also important in a context of interdisciplinary researches. A software should be flexible enough to integrate new ideas and new data that are not initially foreseen, and it should be usable by a large spectrum of persons, from end-user to expert in software development. Another important point is the modelling of multiscale problems which are legion in biology. Despite the difficulties of simulating these systems, several works have been made. For example, Jean-Louis Deneubourg in 1989 focused on examples where microscopic pheromones guide the collective behaviors of agents [START_REF] Deneubourg | Collective patterns and decision-making[END_REF]. The multiagent architecture of Jacques Ferber showed agents made of agents [START_REF] Ferber | Les systèmes multi-agents. Vers une intelligence collective[END_REF]. Generally, multiscale simulations necessitate powerful processors when the number of entities is high. Today, the power of computers continues to grow exponentially even though clock speed of microprocessors has stabilized below 5 GHz 12. The means chosen to overcome this limitation is the use of multicore microprocessors. The first multicore processor called POWER4, was made in 2001 by IBM. [START_REF] Deneubourg | Collective patterns and decision-making[END_REF] The first general public multicore processor was the Intel Pentium D in 2005 just followed by the AMD Athlon 64 x2 in 2006 14. Today, the graphics cards have the highest computing power through their massively parallel architecture. They are programmable thanks to general purpose language close the C language like Cuda [START_REF] Nickolls | Scalable Parallel Programming with CUDA[END_REF] or OpenCL [START_REF] Munshi | The OpenCL Specification by the Khronos OpenCL Working Group[END_REF]. The third point is the interactions between human and computer. The maturity of multi-touch, threedimensional screens, virtual or modified reality and the recent development of screens with tactile sensations, such as friction or roughness, improve the immersion in virtual environments and can help user to create multicellular simulations.

These advances in computer sciences and technologies, are largely used to develop models in the field of computational biology.

Modeling and computational biology context

Computational biology is a vast interdisciplinary field of science where mathematical and algorithmic approaches are applied to biology. If we focus around the cellular scales, we can see that many scientific approaches can be applied to model biological properties, interactions or mechanisms. These approaches can be chemistry, thermodynamics, classic or quantum mechanics, statistic or individual based modelling, genetic or epigenetic, etc. There also exists numerous ways to represent cellular systems. Any biological object can be modelled as a single and simple entity or as a complex dynamical multi-entity and self-assembled system. Furthermore, their environment can be omitted, viewed as static or modelled as a dynamic system. To model cellular systems, there are two principal methods. The first focuses on individual biological objects like molecules or organels and isolate them. This method is called reductionism, for which a complex system is divided into simpler sub-systems or objects. The second method does the opposite by putting together individual biological objects. This method is called integrative biology, where simple sub-systems or objects are combined to explain the functioning of a bigger system. An important difficulty is the development of theories able to integrate these scattered data. This development has already begun inside the integrative computational biology field, also called system biology. For Denis Noble, "Systems biology is about putting together rather than taking apart, integration rather than reduction. It requires that we develop ways of thinking about integration that are as rigorous as our reductionist programs, but different. It means changing our philosophy, in the full sense of the term". The best tool we got to work on the complex is, by far, the human brain. However, when the quantity of data is important, when the diversity of behaviors, properties, structures, parameters and scales become overwhelming, the help of computer can be decisive. This computational assistance can apply both when the data are scattered, to gather them in a comprehensive way, or when the theories are well established, to explore their potential and make predictions. The computer has become what Jol de Rosnay called in 1975 a macroscop, a tool to investigate the infinitely complex, by analogy with the microscope that explores the infinitely small and the telescope which observes the infinitely large [START_REF] De Rosnay | The Macroscope, Principia Cybernetica Web[END_REF]. We briefly see in the next section that each scale of modelling has its own methods, theories, tools and objectives, making the development of multiscale computational models difficult.

-Molecular scale (from 10 -9 to 10 -8 m) A molecule can be seen either as a whole or made of atoms. It can be modelled by using quantum mechanics or Newtonian theories. At this scale, a classical method is the computation of molecular binding, also called molecular docking where the best-fit location and orientation of a ligand that binds a protein is calculated. This method principally aims at the rational design of drugs [START_REF] Kitchen | Docking and scoring in virtual screening for drug discovery: methods and applications[END_REF]. Overviews of docking methods can be seen in [START_REF] Taylor | A review of protein-small molecule docking methods[END_REF] and [START_REF] Lengauer | Computational methods for biomolecular docking[END_REF]. Software to compute docking have also been developed like Arguslab [START_REF] Thompson | ArgusLab 4.0.1[END_REF] or AutoDock [START_REF] Huey | A Semiempirical Free Energy Force Field with Charge-Based Desolvation[END_REF].

-Macromolecular scale (from 10 -8 to 10 -6 m) A macromolecule can be modelled either using interacting molecules or as a single entity. A common objective consists in the computation of protein folding to determine the tertiary structure (3-dimensional structure) of a protein using its primary sequence of amino acids. An overview of simulations can be seen in [START_REF] Scheraga | Protein-Folding Dynamics: Overview of Molecular Simulation Techniques[END_REF]. We can quote some software like Folding at home [START_REF] Pande | Folding@home[END_REF] made to understand protein folding, misfolding and related diseases, using personal computers across the world. Another software called ProteinShop was developed for an interactive protein manipulation [START_REF] Crivelli | ProteinShop: A Tool for Interactive Protein Manipulation[END_REF]. Another example at molecular scale deals with the simulation of lipid membrane. Some methods can be seen in [START_REF] Dopico | Methods in Membrane Lipids[END_REF] and, as for software, we have GROMACS for example, which is able to treat hundreds to millions of particles [START_REF] Hess | GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation[END_REF]. A more macroscopic model has been made to reproduce the behavior of the inner mitochondrial membrane [START_REF] Demongeot | An open problem: the inner mitochondrial membrane (IMM) as a free boundary problem[END_REF].

-Cell scale (from 10 -7 to 10 -5 m) At the cellular scale, like the previous levels studied above, a cell can be seen as a single irreducible unit or as a system made of interacting parts. At this scale level, an important research concerns the cell metabolism. It can be simulated by using very different methods from the ones seen just before. One of the first work was published in 1964 on metabolic control mechanisms [START_REF] Garfinkel | Metabolic control mechanisms. VII. A detailed computer model of the glycolytic pathway in ascites cells[END_REF]. At this scale, biological regulatory network like genetic regulatory network can be applied to model gene interactions and to help biologists to understand and predict emergent cellular behaviors [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF] [START_REF] Bernot | Behaviour Preservation of a Biological Regulatory Network when Embedded into a Larger Network[END_REF]. Software with advanced graphical user interface have been produced like GEPASI [START_REF] Mendes | GEPASI: A software package for modelling the dynamics, steady states and control of biochemical and other systems[END_REF] which works on biochemical systems, or like Virtual Cell [START_REF] Schaff | A general computational framework for modeling cellular structure and function[END_REF] which treats two or three dimensional biochemical problems or finally, like E-Cell [START_REF] Tomita | E-CELL: software environment for whole-cell simulation[END_REF] that enables precise whole cell simulations with activation or knock-out of genes.

-Multicellular scale (from 10 -6 to 10 -3 m) We move up now to the scale of several interacting cells which is our main subject of interest. A multicellular tissue, organ or organism can be abstracted thanks to interacting single cell models or thanks to more global tissue models. Several computational approaches can be used to model and simulate multicellular systems. The more common are briefly presented below.

Algorithmic approaches

Reaction-diffusion

A seminal research started with the work of Alan Turing on the Chemical Basis of Morphogenesis in 1952 which gave birth to the reaction-diffusion field of science. This field was improved by many researchers like Hanz Meinhardt and J. D. Murray. Although this approach focuses on molecular mechanisms, its scale of observation is close to 106 m. It consists in a macroscopic model of interacting and moving molecules and is particularly relevant and useful for the multicellular scale. The next section describes one of the first integrated software developed to simulate multicellular system. This software, called CompuCell3, is based on the Cellular Potts Model for modelling cells and uses a reaction-diffusion system to model molecules.

Ising system

One of the first system that can reproduce selforganized system and which can be easily simulated with a computer is the Ising model, solved in 1925 by Ernst Ising [START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF]. It consists in a mathematical model of ferromagnetism with discrete variables, called spins, which are located in a lattice and interact with their nearest neighbors. It reproduces phase-transition occurring in real substances.

Cellular Potts Model

An evolution of the Ising and Pott systems is the Cellular Pott's Model (CPM). It was created in 1992 by James Glazier and Francois Graner [START_REF] Graner | Simulation of Biological Cell Sorting Using a Two-Dimensional Extended Potts Model[END_REF]. It enables the simulation of multicellular systems by taking into account surface and volume energies for each cell and contact energy between cells. It is used in many ways today, for example to simulate embryogenesis [START_REF] Guillou | A 3D mechanical model of the early mammalian embryo[END_REF], morphogenesis [START_REF] Marée | From Pattern Formation to Morphogenesis Multicellular Coordination in Dictyostelium discoideum[END_REF] or angiogenesis [START_REF] Merks | Modeling Morphogenesis in silico and in vitro: Towards Quantitative, Predictive, Cell-based Modeling[END_REF].

Cellular automata

A classic approach is the cellular automaton, which was designed in the 1940s by John von Neumann and Stanislaw Ulam [START_REF] Beyer | Stanislaw M. Ulam's Contributions to Theoretical Theory[END_REF]. It was initially used to study self-reproducible systems. In this approach, a cell is represented by a value located in a matrix. It was used to create a simple but expressive system called the Game of Life in the middle of the 1960's by John Conway where simple local rules, leading the birth or death of cells, enables complex emerging behaviors [START_REF] Gardner | Mathematical Games: The fantastic combinations of John Conway's new solitaire game life[END_REF]. This approach is also interesting to reproduce many different physical systems [START_REF] Chopard | Cellular Automata Modeling of Physical Systems[END_REF]. However, programming complex biological systems with this paradigm is not easy because the design of rules are not intuitive. The multiagent approach is more relevant in an interdisciplinary context.

Multiagent systems (MAS) and Individual Based Modelling (IBM)

An Individual Based Modelling 19 is required when each entity has to be distinguished from others. It is also useful when there are too few entities so that they cannot be summarized by a continuous value. Moreover, it is interesting when local behaviors are known whereas the global behavior is not. furthermore, this modelling is appropriate when the bottom-up emergent properties are preferred to more directive top-down approaches. Many examples have been developed in MAS. For instance, they were used in 1990 to simulate immune responses thanks to the work of Philip E. Seiden and Franco Celada [START_REF] Seiden | A model for simulating cognate recognition and response in the immune system[END_REF]. This individual based approach can also be combined with global approaches like in the simulation of the earliest stages of embryogenesis in a multicellular organism or brain [START_REF] Fleischer | A Multiple-Mechanism Developmental Model for Defining Self-Organizing Geometric Structures[END_REF] or for the simulation of the angiogenesis in cancer development made in 2005 by A. Stéphanou, S.R. McDougall, A.R.A. Anderson and M.A.J. Chaplain [START_REF] Stéphanou | Mathematical modelling of flow in 2D and 3D vascular network: application to antiangiogenic and chemotherapeutic drug strategies[END_REF].

Mass-spring systems

A cell membrane can be seen as a mass-spring system [START_REF] Rosello | Model Driven Quantification of Individual and Collective Cell Migration[END_REF], eventually combined with a multiagent systems to improve the expressivity of simulations [START_REF] Ballet | Migration de Cellules Virtuelles Déformablesmodélisation biomécanique multiagent de la migration cellulaire[END_REF]. Such approach has also been developed by Szymon Stoma et al. to simulate a plant morphogenesis thanks to a mass-spring system [START_REF] Stoma | Using mechanics in the modelling of meristem morphogenesis[END_REF].

Voronoi diagram

When simulating numerous interacting cells using a mass-spring system, the computation of interacting cell is essential. It can be done by calculating a Delaunay-Voronoi tessellation, as used in the embryomorphic model developed by R. Doursat [START_REF] Doursat | From the Simulation of Complex Biological Systems to the Design of Artificial Morphogenetic Systems, and Back[END_REF].

Other approaches

Fleisher made a model and a simulator where cells can bind and are able to secrete chemoattractants and chemorepellents to generate patterns of virtual proto-organisms [START_REF] Fleischer | A Simulation Testbed for the Study of Multicellular Development: The Multiple Mechanisms of Morphogenesis[END_REF]. Laforge et al. developed an approach based on an equilibrium between the auto-stabilization of stochastic gene expression and the interdependence of cells for proliferation to model embryogenesis and cancer [START_REF] Laforge | Modeling embryogenesis and cancer: an approach based on an equilibrium between the autostabilization of stochastic gene expression and the interdependence of cells for proliferation[END_REF]. Another model of morphogenesis in cellular systems is related to the study of the cell reorganization during in-vitro wound healing [START_REF] Dugnolle | A mechanical model to simulate cell reorganisation during in-vitro wound healing[END_REF]. Other approaches are currently developed and can be seen in the book [START_REF] Anderson | Single-cell-based models in biology and medicine[END_REF].

Software examples

To complete this paper, this section deals with two fully integrated software designed to reproduce multicellular systems.

CompuCell3D (see figure 1)

-Description CompuCell3D is a free open-source software based on the Cellular Potts Model where cells are modelled thanks to several matrix sites that represent their surface and volume. In this approach, cells also interact with each other thanks to adhering properties and chemical signals. The software enables the simulation of cellular and multicellular systems using parallel devices.

-Key features This software uses an original approach, invented by their authors in 1992, which is relevant for the simulation of many different multicellular systems from cell sorting to organism morphogenesis. This system can also be connected with other matrixbased systems like fluid mechanics based on Navier-Stoke to achieve simulations involving, for example, blood coagulation or blood circulation.

-Computational model In this approach, a molecule is a simple value in a matrix. A cell, which is more detailed than a molecule, is a dynamic set of matrix sites. A membrane is simply defined as cell site borders in contact with other cell sites belonging to other cells. We can notice that many cellular behaviors are simulated like migration, division, differentiation and apoptosis.

-Biological scale The usual scales go from 1 micrometer square (or cube when in 3D) for a cell site to 1 millimeter square (or cube) for a whole simulation.

-Computational method The shape and the size of a cell depend on its target surface and volume. Annex developments have been done, like one enabling the design of a target shape for each cell type [START_REF] Tripodi | Selforganization of a virtual multicellular organism by adding a shape model in the Cellular Potts Model[END_REF]. We can note that about 30 plugins can be added to the simulator like Chemotaxis, Mitosis or Secretion.

-Graphical user interface The graphical user interface enables the user to manage the execution of simulations, view the evolution of the system and see the cellular types used. The XML file describing a simulation must be made outside the interface.

-Programming language Two approaches are possible. The first is the use of an XML file to describe a simulation. The second one, which is more powerful, is the use of the Python language to explicit models and simulations. Developers can create their own plugins using the C language.

-Documentation The web site of Compucell3D offers a useful quick start guide, including tutorials and exercises. Two detailed documents are also dedicated to bio-modelers. The first explains the general objectives, the theoretical and algorithmic point of view and contains advanced examples made with Compucell3D. The second describes how to use Python to develop models and simulations. Moreover, a paper explaining how to create plugins in C language is downloadable for developers.

SimCells (see Figure 2)

-Description SimCells is a free software based on a multiagent system where deformable and interacting cells are modelled inside matrices. The user can create different kind of multicellular systems by using a simple graphical language.

-Key features This software uses an original approach, developed for multicore devices. All steps of modeling are graphically made thanks to an advanced graphical user interface, letting nondevelopers to create their own simulations.

-Computational model Like for CompuCell3D, many cellular behaviors can be simulated like migration, chemotaxis, division, differentiation, apoptosis, membrane deformation and adhesion.

-Biological scale SimCells uses a multiscale approach with cells, bricks, molecules and fields.

-Computational method It uses recent Graphical Processor Unit (GPU). The number of simulated cells depends on the power of the graphical processor unit (GPU) used. It can be about one million with a premium GPU in 2018.

-Graphical user interface The graphical user interface enables the user to manage the execution of simulations, view the evolution of the system, create cells, design the initial state, view and export simulation results.

-Programming language The creation of cellular behaviors is made by a graphical language under the form of Conditions Then Actions. Each entity (a cell, brick or molecule) carries its own behaviors. The number of entities is not limited, and they can be explicitly named (like Lymphocyte, Antigen, etc).

-Documentation The web site of SimCells gives a user guide explaining the different entities specificities, how to program them and how to export results and video.

Timeline of multicellular simulation

This section presents a brief history of software built to simulate multicellular systems. A time line of different models and approaches is shown on the table 1. It covers the period from 1965 up to now. A third time line focuses on advanced software made to simulate multicellular systems. It can be seen on the table 2.

Table 1.1: brief history of computational systems for multi-cellular simulation a b c . Table 1.2: brief history of computational systems for multi-cellular simulation a . Table 2: advanced software for multi-cellular simulation a .

SimCells example

The growth of a cancerous tumor and its blood vessel network is a multicellular self-regulated mechanism. Its comprehension is complex because it involves a feedback loop: the tumor is guiding the angiogenesis, and at the same time, the blood vessels supply the cancerous cells with nutriments and oxygen.

A way of studying the dynamics of this system consists in building a computational multiagent simulation that can i) exhibit fundamental processes occurring over time, ii) make forecasts on how the tumor could evolve in the future and iii) give clues about how to control the tumor development.

A multiagent simulation, based on the software SimCells, is described hereafter.

Virtual agents' description

The simulation is made of three principal agents: cancerous cells, blood vessel and tip cells (see figure 3). Moreover, three diffusing fields are required: the supply field, (field 3) generated by the blood vessel, the vascular growth field (field 1) generated by cancerous cells and the matter consumption field (field 2) generated by the growing blood vessels.

Two other agents are used to reproduce inert tissue that can be found close to tumors: bones and healthy tissue. -Cancerous cells They have two states: starvation, when they are too far from a blood vessel and dividing when they are supplied by blood vessels.

-Blood vessels A blood vessel agent has two states: stabilized, when mature and not evolving anymore and growing when developing towards the tumor.

Stabilized vessels have 2 behaviors:

Supply: at each simulation time step, each blood vessel agent produces a diffusing field 3 (supply field).

Angiogenesis: at each time step, on a probability of 0.01 and when the field 1 (vascular growth field) is greater than a threshold and when the field 2 (matter consumption field) is lower than another threshold, then the stabilized blood vessel become a growing blood vessel. -Tip cells Their main behavior consists in orienting towards the positive gradient of the vascular growth field:

Simulation results

A 2D simulation is performed in a 1024x512 matrix. The initial state contains a single vertical blood vessel and a tiny round tumor with about 50 cancerous cells. The figure 4 shows the early stages of angiogenesis. At about 7000 simulation steps, the tumor reaches a bone (see figure 5).

The simulation ended after 27000 simulation steps (see figure 6). At this time, all the available space is occupied by the tumor and its blood vessel network, showing the robustness of this co-developing system. The simulation can also be executed in a 3D environment of size 126x126x126, but the size of the environment must be reduced (see Figure 5) due to memory size of the graphical card. The early stage of the simulation is shown on the figure 7.

Figure 7: same simulation in a 3D environment.

Conclusion

The biomedical field at cellular scale and computer science are independently and together in deep change, opening new ways in the control of the living at the cellular scale. In this context, computer aided design of large and dynamic multicellular tissues becomes possible.

Computer scientists and biologists can cooperate in a complex interdisciplinary field thanks to advanced software, sharing common formal computing languages.

The software, like SimCells, use specific graphical interfaces to exploit computational models in a simple way, keeping creativity of the user at a high level.

The simulation of numerous interacting cells can use multicore devices to drastically decrease the computational time.

Interdisciplinary in science, although hard to implement, is a requirement for the development of multicellular simulators. More scientists and engineers at the interface of different sciences should be formed to be the spearhead of these new advances.

More generally, it will help to support biology and medicine towards more predictive, quantitative and individualized levels.

Figure 1 :

 1 Figure 1: Graphical user interface of the integrated software CompuCell3D.

Figure 2 :

 2 Figure 2: Graphical user interface of the advanced software SimCells.

Figure 3 :

 3 Figure 3: main agents of the simulation.

Figure 4 :

 4 Figure 4: simulation after 1950 simulation steps.The main blood vessel (on the right), has new developing branches (from left to right).

Figure 5 :

 5 Figure 5: simulation at 6930 steps. The tumor continues to grow, inducing the co-development of the blood vessel network. Round gray agents, on the left, represent bones. Grey cells, at the upper-left and bottom-right, represent inert tissues compressed by cancerous cells.

Figure 6 :

 6 Figure 6: final simulation state at 27030 steps. All available space is occupied by dividing cancerous cells and stabilized blood vessels.