Matthieu Solnon
email: matthieu.solnon@ens.fr

Sylvain Arlot
email: sylvain.arlot@ens.fr

Francis Bach
email: francis.bach@ens.fr

Tong Zhang

Multi-task Regression using Minimal Penalties

Keywords: multi-task, oracle inequality, learning theory c 2012 Solnon, Arlot and Bach

In this paper we study the kernel multiple ridge regression framework, which we refer to as multi-task regression, using penalization techniques. The theoretical analysis of this problem shows that the key element appearing for an optimal calibration is the covariance matrix of the noise between the different tasks. We present a new algorithm to estimate this covariance matrix, based on the concept of minimal penalty, which was previously used in the single-task regression framework to estimate the variance of the noise. We show, in a non-asymptotic setting and under mild assumptions on the target function, that this estimator converges towards the covariance matrix. Then plugging this estimator into the corresponding ideal penalty leads to an oracle inequality. We illustrate the behavior of our algorithm on synthetic examples.

Introduction

A classical paradigm in statistics is that increasing the sample size (that is, the number of observations) improves the performance of the estimators. However, in some cases it may be impossible to increase the sample size, for instance because of experimental limitations. Hopefully, in many situations practicioners can find many related and similar problems, and might use these problems as if more observations were available for the initial problem. The techniques using this heuristic are called "multi-task" techniques. In this paper we study the kernel ridge regression procedure in a multi-task framework.

One-dimensional kernel ridge regression, which we refer to as "single-task" regression, has been widely studied. As we briefly review in Section 3 one has, given n data points (X i , Y i ) n i=1 , to estimate a function f , often the conditional expectation f

(X i ) = E[Y i |X i ]
, by minimizing the quadratic risk of the estimator regularized by a certain norm. A practically important task is to calibrate a regularization parameter, that is, to estimate the regularization parameter directly from data. For kernel ridge regression (a.k.a. smoothing splines), many methods have been proposed based on different principles, for example, Bayesian criteria through a Gaussian process interpretation (see, e.g., [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] or generalized cross-validation (see, e.g., [START_REF] Wahba | Spline Models for Observational Data[END_REF]. In this paper, we focus on the concept of minimal penalty, which was first introduced by [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] and [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF] for model selection, then extended to linear estimators such as kernel ridge regression by [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF].

In this article we consider p ≥ 2 different (but related) regression tasks, a framework we refer to as "multi-task" regression. This setting has already been studied in different papers. Some empirically show that it can lead to performance improvement [START_REF] Thrun | Discovering structure in multiple learning tasks: The TC algorithm[END_REF][START_REF] Caruana | Multitask learning[END_REF][START_REF] Bakker | Task clustering and gating for bayesian multitask learning[END_REF]. [START_REF] Liang | Asymptotically optimal regularization in smooth parametric models[END_REF] also obtained a theoretical criterion (unfortunately non observable) which tells when this phenomenon asymptotically occurs. Several different paths have been followed to deal with this setting. Some consider a setting where p ≫ n, and formulate a sparsity assumption which enables to use the group Lasso, assuming all the different functions have a small set of common active covariates (see for instance [START_REF] Obozinski | Support union recovery in high-dimensional multivariate regression[END_REF][START_REF] Lounici | Oracle inequalities and optimal inference under group sparsity[END_REF]. We exclude this setting from our analysis, because of the Hilbertian nature of our problem, and thus will not consider the similarity between the tasks in terms of sparsity, but rather in terms of an Euclidean similarity. Another theoretical approach has also been taken (see for example, [START_REF] Brown | Adaptive multivariate ridge regression[END_REF], [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF] or Ando and [START_REF] Zhang | Learning bounds for kernel regression using effective data dimensionality[END_REF] on semi-supervised learning), the authors often defining a theoretical framework where the multi-task problem can easily be expressed, and where sometimes solutions can be computed. The main remaining theoretical problem is the calibration of a matricial parameter M (typically of size p), which characterizes the relationship between the tasks and extends the regularization parameter from single-task regression. Because of the high dimensional nature of the problem (i.e., the small number of training observations) usual techniques, like cross-validation, are not likely to succeed. [START_REF] Argyriou | Convex multi-task feature learning[END_REF] have a similar approach to ours, but solve this problem by adding a convex constraint to the matrix, which will be discussed at the end of Section 5.

Through a penalization technique we show in Section 2 that the only element we have to estimate is the correlation matrix Σ of the noise between the tasks. We give here a new algorithm to estimate Σ, and show that the estimation is sharp enough to derive an oracle inequality for the estimation of the task similarity matrix M , both with high probability and in expectation. Finally we give some simulation experiment results and show that our technique correctly deals with the multi-task settings with a low sample-size.

Notations

We now introduce some notations, which will be used throughout the article.

• The integer n is the sample size, the integer p is the number of tasks.

• For any n × p matrix Y , we define y = vec(Y ) := (Y 1,1 , . . . , Y n,1 , Y 1,2 , . . . , Y n,2 , . . . , Y 1,p , . . . , Y n,p ) ∈ R np , that is, the vector in which the columns Y j := (Y i,j ) 1≤i≤n are stacked.

• M n (R) is the set of all matrices of size n.

• S p (R) is the set of symmetric matrices of size p.

• S +

p (R) is the set of symmetric positive-semidefinite matrices of size p.

• S ++ p (R) is the set of symmetric positive-definite matrices of size p.

• denotes the partial ordering on S p (R) defined by: A B if and only if B -A ∈ S + p (R).

• 1 is the vector of size p whose components are all equal to 1.

• • 2 is the usual Euclidean norm on R k for any k ∈ N: ∀u ∈ R k , u 2 2 := k i=1 u 2 i .

Multi-task Regression: Problem Set-up

We consider p kernel ridge regression tasks. Treating them simultaneously and sharing their common structure (e.g., being close in some metric space) will help in reducing the overall prediction error.

Multi-task with a Fixed Kernel

Let X be some set and F a set of real-valued functions over X . We suppose F has a reproducing kernel Hilbert space (RKHS) structure [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], with kernel k and feature map Φ : X → F. We observe

D n = (X i , Y 1 i , . . . , Y p i ) n i=1 ∈ (X × R p ) n , which gives us the positive semidefinite kernel matrix K = (k(X i , X ℓ )) 1≤i,ℓ≤n ∈ S + n (R). For each task j ∈ {1, . . . , p}, D j n = (X i , y j i ) n
i=1 is a sample with distribution P j , for which a simple regression problem has to be solved. In this paper we consider for simplicity that the different tasks have the same design (X i ) n i=1 . When the designs of the different tasks are different the analysis is carried out similarly by defining X i = (X 1 i , . . . , X p i ), but the notations would be more complicated.

We now define the model. We assume (f 1 , . . . , f p ) ∈ F p , Σ is a symmetric positivedefinite matrix of size p such that the vectors (ε j i ) p j=1 are i.i.d. with normal distribution N (0, Σ), with mean zero and covariance matrix Σ, and ∀i ∈ {1, . . . , n}, ∀j ∈ {1, . . . , p}, y j i = f j (X i ) + ε j i .

(1)

This means that, while the observations are independent, the outputs of the different tasks can be correlated, with correlation matrix Σ between the tasks. We now place ourselves in the fixed-design setting, that is, (X i ) n i=1 is deterministic and the goal is to estimate f 1 (X i ), . . . , f p (X i ) n i=1 . Let us introduce some notation:

• µ min = µ min (Σ) (resp. µ max ) denotes the smallest (resp. largest) eigenvalue of Σ.

• c(Σ) := µ max /µ min is the condition number of Σ.

To obtain compact equations, we will use the following definition:

Definition 1 We denote by F the n × p matrix (f j (X i )) 1≤i≤n , 1≤j≤p and introduce the vector

f := vec(F ) = (f 1 (X 1 ), . . . , f 1 (X n ), . . . , f p (X 1 ), . . . , f p (X n )) ∈ R np , obtained by stacking the columns of F . Similarly we define Y := (y j i ) ∈ M n×p (R), y := vec(Y ), E := (ε j i ) ∈ M n×p (R)
and ε := vec(E). In order to estimate f , we use a regularization procedure, which extends the classical ridge regression of the single-task setting. Let M be a p × p matrix, symmetric and positivedefinite. Generalizing the work of [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF], we estimate (f 1 , . . . , f p ) ∈ F p by

f M ∈ argmin g∈F p    1 np n i=1 p j=1 (y j i -g j (X i )) 2 + p j=1 p ℓ=1 M j,l g j , g ℓ F    . (2)
Although M could have a general unconstrained form we may restrict M to certain forms, for either computational or statistical reasons.

Remark 2 Requiring that M 0 implies that Equation (2) is a convex optimization problem, which can be solved through the resolution of a linear system, as explained later. Moreover it allows an RKHS interpretation, which will also be explained later.

Example 3 The case where the p tasks are treated independently can be considered in this setting: taking M = M ind (λ) := 1 p Diag(λ 1 , . . . , λ p ) for any λ ∈ R p leads to the criterion

1 p p j=1 1 n n i=1 (y j i -g j (X i )) 2 + λ j g j 2 F , (3)
that is, the sum of the single-task criteria described in Section 3. Hence, minimizing Equation (3) over λ ∈ R p amounts to solve independently p single task problems.

Example 4 As done by [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF], for every λ, µ ∈ (0, +∞) 2 , define

M similar (λ, µ) := (λ + pµ)I p -µ11 ⊤ =    λ + (p -1)µ -µ . . . -µ λ + (p -1)µ    .
Taking M = M similar (λ, µ) in Equation (2) leads to the criterion

1 np n i=1 p j=1 (y j i -g j (X i )) 2 + λ p j=1 g j 2 F + µ 2 p j=1 p k=1 g j -g k 2 F . ( 4
)
Minimizing Equation (4) enforces a regularization on both the norms of the functions g j and the norms of the differences g j -g k . Thus, matrices of the form M similar (λ, µ) are useful when the functions g j are assumed to be similar in F. One of the main contributions of the paper is to go beyond this case and learn from data a more general similarity matrix M between tasks.

Example 5 We extend Example 4 to the case where the p tasks consist of two groups of close tasks. Let I be a subset of {1, . . . , p}, of cardinality 1 ≤ k ≤ p -1. Let us denote by I c the complementary of I in {1, . . . , p}, 1 I the vector v with components v i = 1 i∈I , and Diag(I) the diagonal matrix d with components d i,i = 1 i∈I . We then define

M I (λ, µ, ν) := λI p + µ Diag(I) + ν Diag(I c ) - µ k 1 I 1 ⊤ I - ν p -k 1 I c 1 ⊤ I c .
This matrix leads to the following criterion, which enforces a regularization on both the norms of the functions g j and the norms of the differences g jg k inside the groups I and I c :

1 np n i=1 p j=1 (y j i -g j (X i )) 2 +λ p j=1 g j 2 F + µ 2k j∈I k∈I g j -g k 2 F + ν 2(p -k) j∈I c k∈I c g j -g k 2 F .
(5) As shown in Section 6, we can estimate the set I from data (see [START_REF] Jacob | Clustered multi-task learning: A convex formulation[END_REF] for a more general formulation).

Remark 6 Since I p and 11 ⊤ can be diagonalized simultaneously, minimizing Equation (4) and Equation (5) is quite easy: it only demands optimization over two independent parameters, which can be done with the procedure of [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF].

Remark 7 As stated below (Proposition 8), M acts as a scalar product between the tasks. Selecting a general matrix M is thus a way to express a similarity between tasks.

Following [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF], we define the vector-space G of real-valued functions over X × {1, . . . , p} by G := {g : X × {1, . . . , p} → R / ∀j ∈ {1, . . . , p} , g(•, j) ∈ F} .

We now define a bilinear symmetric form over G,

∀g, h ∈ G , g, h G := p j=1 p l=1 M j,l g(•, j), h(•, l) F ,
which is a scalar product as soon as M is positive semi-definite (see proof in Appendix A) and leads to a RKHS (see proof in Appendix B):

Proposition 8 With the preceding notations •, • G is a scalar product on G. Corollary 9 (G, •, • G ) is a RKHS.
In order to write down the kernel matrix in compact form, we introduce the following notations.

Definition 10 (Kronecker Product) Let A ∈ M m,n (R), B ∈ M p,q (R). We define the Kronecker product A ⊗ B as being the (mp) × (nq) matrix built with p × q blocks, the block of index (i, j) being A i,j • B:

A ⊗ B =    A 1,1 B . . . A 1,n B . . . . . . . . . A m,1 B . . . A m,n B    .
The Kronecker product is a widely used tool to deal with matrices and tensor products. Some of its classical properties are given in Section E; see also [START_REF] Horn | Topics in Matrix Analysis[END_REF].

Proposition 11 The kernel matrix associated with the design X := (X i , j) i,j ∈ X × {1, . . . , p} and the RKHS (G,

•, • G ) is K M := M -1 ⊗ K.
Proposition 11 is proved in Appendix C. We can then apply the representer's theorem [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF] to the minimization problem (2) and deduce that f M = A M y with

A M = A M,K := K M ( K M + npI np ) -1 = (M -1 ⊗ K) (M -1 ⊗ K) + npI np -1 .

Optimal Choice of the Kernel

Now when working in multi-task regression, a set M ⊂ S ++ p (R) of matrices M is given, and the goal is to select the "best" one, that is, minimizing over M the quadratic risk n -1 f Mf 2 2 . For instance, the single-task framework corresponds to p = 1 and M = (0, +∞). The multi-task case is far richer. The oracle risk is defined as inf

M ∈M f M -f 2 2 . ( 6
)
The ideal choice, called the oracle, is any matrix

M ⋆ ∈ argmin M ∈M f M -f 2 2 .
Nothing here ensures the oracle exists. However in some special cases (see for instance Example 12) the infimum of f Mf 2 over the set { f M , M ∈ M} may be attained by a function f * ∈ F p -which we will call "oracle" by a slight abuse of notation-while the former problem does not have a solution.

From now on we always suppose that the infimum of { f Mf 2 } over M is attained by some function f ⋆ ∈ F p . However the oracle M ⋆ is not an estimator, since it depends on f .

Example 12 (Partial computation of the oracle in a simple setting) It is possible in certain simple settings to exactly compute the oracle (or, at least, some part of it).

Consider for instance the set-up where the p functions are taken to be equal (that is,

f 1 = • • • = f p ).
In this setting it is natural to use the set

M similar := M similar (λ, µ) = (λ + pµ)I p - µ p 11 ⊤ / (λ, µ) ∈ (0, +∞) 2 .
Using the estimator f M = A M y we can then compute the quadratic risk using the biasvariance decomposition given in Equation (36):

E f M -f 2 2 = (A M -I np )f 2 2 + tr(A ⊤ M A M • (Σ ⊗ I n )) .
Computations (reported in Appendix D) show that, with the change of variables µ = λ + pµ, the bias does not depend on µ and the variance is a decreasing function of µ. Thus the oracle is obtained when µ = +∞, leading to a situation where the oracle functions

f 1,⋆ , . . . , f p,⋆ verify f 1,⋆ = • • • = f p,⋆ .
It is also noticeable that, if one assumes the maximal eigenvalue of Σ stays bounded with respect to p, the variance is of order O(p -1 ) while the bias is bounded with respect to p.

As explained by [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF], we choose

M ∈ argmin M ∈M {crit(M )} with crit(M ) = 1 np y -f M 2 2 + pen(M ) ,
where the penalty term pen(M ) has to be chosen appropriately.

Remark 13 Our model (1) does not constrain the functions f 1 , . . . , f p . Our way to express the similarities between the tasks (that is, between the f j ) is via the set M, which represents the a priori knowledge the statistician has about the problem. Our goal is to build an estimator whose risk is the closest possible to the oracle risk. Of course using an inappropriate set M (with respect to the target functions f 1 , . . . , f p ) may lead to bad overall performances. Explicit multi-task settings are given in Examples 3, 4 and 5 and through simulations in Section 6.

The unbiased risk estimation principle (introduced by [START_REF] Akaike | Statistical predictor identification[END_REF] requires

E [crit(M )] ≈ E 1 np f M -f 2 2
, which leads to the (deterministic) ideal penalty

pen id (M ) := E 1 np f M -f 2 2 -E 1 np y -f M 2 2 . Since f M = A M y and y = f + ε, we can write f M -y 2 2 = f M -f 2 2 + ε 2 2 -2 ε, A M ε + 2 ε, (I np -A M )f .
Since ε is centered and M is deterministic, we get, up to an additive factor independent of M ,

pen id (M ) = 2E [ ε, A M ε ] np , that is, as the covariance matrix of ε is Σ ⊗ I n , pen id (M ) = 2 tr A M • (Σ ⊗ I n ) np . (7)
In order to approach this penalty as precisely as possible, we have to sharply estimate Σ.

In the single-task case, such a problem reduces to estimating the variance σ 2 of the noise and was tackled by [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF]. Since our approach for estimating Σ heavily relies on these results, they are summarized in the next section.

Note that estimating Σ is a mean towards estimating M . The technique we develop later for this purpose is not purely a multi-task technique, and may also be used in a different context.

Single Task Framework: Estimating a Single Variance

This section recalls some of the main results from [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF] which can be considered as solving a special case of Section 2, with p = 1, Σ = σ 2 > 0 and

M = [0, +∞]. Writing M = λ with λ ∈ [0, +∞], the regularization matrix is ∀λ ∈ (0, +∞) , A λ = A λ,K = K(K + nλI n ) -1 , A 0 = I n and A +∞ = 0; the ideal penalty becomes pen id (λ) = 2σ 2 tr(A λ ) n .
By analogy with the case where A λ is an orthogonal projection matrix, df(λ) := tr(A λ ) is called the effective degree of freedom, first introduced by [START_REF] Mallows | Some comments on C P . Technometrics[END_REF]; see also the work by [START_REF] Zhang | Learning bounds for kernel regression using effective data dimensionality[END_REF]. The ideal penalty however depends on σ 2 ; in order to have a fully data-driven penalty we have to replace σ 2 by an estimator σ 2 inside pen id (λ). For every λ ∈ [0, +∞], define

pen min (λ) = pen min (λ, K) := (2 tr(A λ,K ) -tr(A ⊤ λ,K A λ,K )) n .
We shall see now that it is a minimal penalty in the following sense. If for every C > 0

λ 0 (C) ∈ argmin λ∈[0,+∞] 1 n A λ,K Y -Y 2 2 + C pen min (λ, K) ,
then-up to concentration inequalities-λ 0 (C) acts as a mimimizer of

g C (λ) = E 1 n A λ Y -Y 2 2 + C pen min (λ) -σ 2 = 1 n (A λ -I n )f 2 2 + (C -σ 2 ) pen min (λ) .
The former theoretical arguments show that

• if C < σ 2 , g C (λ)
decreases with df(λ) so that df( λ 0 (C)) is huge: the procedure overfits;

• if C > σ 2 , g C (λ) increases with df(λ) when df(λ) is large enough so that df( λ 0 (C)) is much smaller than when C < σ 2 .

The following algorithm was introduced by [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF] and uses this fact to estimate σ 2 .

Algorithm 14

Input: Y ∈ R n , K ∈ S ++ n (R) 1. For every C > 0, compute λ 0 (C) ∈ argmin λ∈[0,+∞] 1 n A λ,K Y -Y 2 2 + C pen min (λ, K) . 2. Output: C such that df( λ 0 ( C)) ∈ [n/10, n/3].
An efficient algorithm for the first step of Algorithm 14 is detailed by [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF], and we discuss the way we implemented Algorithm 14 in Section 6. The output C of Algorithm 14 is a provably consistent estimator of σ 2 , as stated in the following theorem.

Theorem 15 (Corollary of Theorem 1 of Arlot and Bach, 2011

) Let β = 150. Sup- pose ε ∼ N (0, σ 2 I n ) with σ 2 > 0, and that λ 0 ∈ (0, +∞) and d n ≥ 1 exist such that df(λ 0 ) ≤ √ n and 1 n (A λ 0 -I n )F 2 2 ≤ d n σ 2 ln n n . (8)
Then for every δ ≥ 2, some constant n 0 (δ) and an event Ω exist such that

P(Ω) ≥ 1 -n -δ and if n ≥ n 0 (δ), on Ω, 1 -β(2 + δ) ln n n σ 2 ≤ C ≤ 1 + β(2 + δ)d n ln(n) n σ 2 . ( 9
)
Remark 16 The values n/10 and n/3 in Algorithm 14 have no particular meaning and can be replaced by n/k, n/k ′ , with k > k ′ > 2. Only β depends on k and k ′ . Also the bounds required in Assumption (8) only impact the right hand side of Equation (9) and are chosen to match the left hand side. See Proposition 10 of [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF] for more details.

Estimation of the Noise Covariance Matrix Σ

Thanks to the results developped by [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF] (recapitulated in Section 3), we know how to estimate a variance for any one-dimensional problem. In order to estimate Σ, which has p(p + 1)/2 parameters, we can use several one-dimensional problems. Projecting Y onto some direction z ∈ R p yields

Y z := Y • z = F • z + E • z = F z + ε z , (10)
with

ε z ∼ N (0, σ 2 z I n ) and σ 2 z := Var[ε • z] = z ⊤ Σz.
Therefore, we will estimate σ 2 z for z ∈ Z a well chosen set, and use these estimators to build back an estimation of Σ.

We now explain how to estimate Σ using those one-dimensional projections.

Definition 17 Let a(z) be the output C of Algorithm 14 applied to problem (10), that is, with inputs Y z ∈ R n and K ∈ S ++ n (R).

The idea is to apply Algorithm 14 to the elements z of a carefully chosen set Z. Noting e i the i-th vector of the canonical basis of R p , we introduce Z = {e i , i ∈ {1, . . . , p}} ∪ {e i + e j , 1 ≤ i < j ≤ p}. We can see that a(e i ) estimates Σ i,i , while a(e i + e j ) estimates Σ i,i + Σ j,j + 2Σ i,j . Henceforth, Σ i,j can be estimated by (a(e i + e j )a(e i )a(e j ))/2. This leads to the definition of the following map J, which builds a symmetric matrix using the latter construction.

Definition 18 Let J : R p(p+1) 2

→ S p (R) be defined by

J(a 1 , . . . , a p , a 1,2 , . . . , a 1,p , . . . , a p-1,p ) i,i = a i if 1 ≤ i ≤ p , J(a 1 , . . . , a p , a 1,2 , . . . , a 1,p , . . . , a p-1,p ) i,j = a i,j -a i -a j 2 if 1 ≤ i < j ≤ p .
This map is bijective, and for all B ∈ S p (R)

J -1 (B) = (B 1,1 , . . . , B p,p , B 1,1 + B 2,2 + 2B 1,2 , . . . , B p-1,p-1 + B p,p + 2B p-1,p ) .
This leads us to defining the following estimator of Σ:

Σ := J (a(e 1 ), . . . , a(e p ), a(e 1 + e 2 ), . . . , a(e 1 + e p ), . . . , a(e p-1 + e p )) .

Remark 19 If a diagonalization basis (e ′ 1 , . . . , e ′ p ) (whose basis matrix is P ) of Σ is known, or if Σ is diagonal, then a simplified version of the algorithm defined by Equation (11) is

Σ simplified = P ⊤ Diag(a(e ′ 1 ), . . . , a(e ′ p ))P . ( 12
)
This algorithm has a smaller computational cost and leads to better theoretical bounds (see Remark 24 and Section 5.2).

Let us recall that ∀λ ∈ (0, +∞), [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF] we make the following assumption from now on: ∀j ∈ {1, . . . , p} , ∃λ 0,j ∈ (0, +∞) ,

A λ = A λ,K = K(K + nλI n ) -1 . Following
df(λ 0,j ) ≤ √ n and 1 n (A λ 0,j -I n )F e j 2 2 ≤ Σ j,j ln n n      (13)
We can now state the first main result of the paper.

Theorem 20 Let Σ be defined by Equation (11), α = 2 and assume (13) holds. For every δ ≥ 2, a constant n 0 (δ), an absolute constant L 1 > 0 and an event Ω exist such that

P(Ω) ≥ 1 -p(p + 1)/2 × n -δ and if n ≥ n 0 (δ), on Ω, (1 -η)Σ Σ (1 + η)Σ (14)
where

η := L 1 (2 + δ)p ln(n) n c(Σ) 2 .
Theorem 20 is proved in Section E. It shows Σ estimates Σ with a "multiplicative" error controlled with large probability, in a non-asymptotic setting. The multiplicative nature of the error is crucial for deriving the oracle inequality stated in Section 5, since it allows to show the ideal penalty defined in Equation ( 7) is precisely estimated when Σ is replaced by Σ.

An important feature of Theorem 20 is that it holds under very mild assumptions on the mean f of the data (see Remark 22). Therefore, it shows Σ is able to estimate a covariance matrix without prior knowledge on the regression function, which, to the best of our knowledge, has never been obtained in multi-task regression.

Remark 21 (Scaling of (n, p) for consistency) A sufficient condition for ensuring Σ is a consistent estimator of Σ is

pc(Σ) 2 ln(n) n -→ 0 ,
which enforces a scaling between n, p and c(Σ). Nevertheless, this condition is probably not necessary since the simulation experiments of Section 6 show that Σ can be well estimated (at least for estimator selection purposes) in a setting where η ≫ 1.

Remark 22 (On assumption (13)) Assumption ( 13) is a single-task assumption (made independently for each task ). The upper bound ln(n)/n can be multiplied by any factor 1 ≤ d n ≪ n/ ln(n) (as in Theorem 15), at the price of multiplying η by d n in the upper bound of Equation ( 14). More generally the bounds on the degree of freedom and the bias in (13) only influence the upper bound of Equation ( 14). The rates are chosen here to match the lower bound, see Proposition 10 of [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF] for more details. Assumption (13) is rather classical in model selection, see [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF] for instance. In particular, (a weakened version of ) (13) holds if the bias n

-1 (A λ -I n )F e i 2 2 is bounded by C 1 tr(A λ ) -C 2 , for some C 1 , C 2 > 0.
Remark 23 (Choice of the set Z) Other choices could have been made for Z, however ours seems easier in terms of computation, since |Z| = p(p + 1)/2. Choosing a larger set Z leads to theoretical difficulties in the reconstruction of Σ, while taking other basis vectors leads to more complex computations. We can also note that increasing |Z| decreases the probability in Theorem 20, since it comes from an union bound over the one-dimensional estimations.

Remark 24 When Σ = Σ simplified as defined by Equation (12), that is, when a diagonalization basis of Σ is known, Theorem 20 still holds on a set of larger probability 1κpn -δ with a reduced error η = L 1 (α + δ) ln(n)/n. Then, a consistent estimation of Σ is possible whenever p = O(n δ ) for some δ ≥ 0.

Oracle Inequality

This section aims at proving "oracle inequalities", as usually done in a model selection setting: given a set of models or of estimators, the goal is to upper bound the risk of the selected estimator by the oracle risk (defined by Equation ( 6)), up to an additive term and a multiplicative factor. We show two oracle inequalities (Theorems 26 and 29) that correspond to two possible definitions of Σ.

Note that "oracle inequality" sometimes has a different meaning in the literature (see for instance [START_REF] Lounici | Oracle inequalities and optimal inference under group sparsity[END_REF] when the risk of the proposed estimator is controlled by the risk of an estimator using information coming from the true parameter (that is, available only if provided by an oracle).

A General Result for Discrete Matrix Sets M

We first show that the estimator introduced in Equation ( 11) is precise enough to derive an oracle inequality when plugged in the penalty defined in Equation ( 7) in the case where M is finite.

Definition 25 Let Σ be the estimator of Σ defined by Equation (11). We define

M ∈ argmin M ∈M f M -y 2 2 + 2 tr A M • ( Σ ⊗ I n )
.

We assume now the following holds true:

∃(C, α M ) ∈ (0, +∞) 2 , card(M) < Cn α M . ( 15
)
Theorem 26 Let α = max(α M , 2), δ ≥ 2 and assume (13) and (15) hold true. Absolute constants L 2 , κ ′ > 0, a constant n 1 (δ) and an event Ω exist such that P( Ω) ≥ 1κ ′ p(p + C)n -δ and the following holds as soon as n ≥ n 1 (δ). First, on Ω,

1 np f M -f 2 2 ≤ 1 + 1 ln(n) 2 inf M ∈M 1 np f M -f 2 2 + L 2 c(Σ) 4 tr(Σ)(α + δ) 2 p 3 ln(n) 3 np . ( 16
)
Second, an absolute constant L 3 exists such that

E 1 np f M -f 2 2 ≤ 1 + 1 ln(n) 2 E inf M ∈M 1 np f M -f 2 2 +L 2 c(Σ) 4 tr(Σ)(α + δ) 2 p 3 ln(n) 3 np + L 3 p(p + C) n δ/2 |||Σ||| + f 2 2 np . ( 17
)
Theorem 26 is proved in Section F.

Remark 27 If Σ = Σ simplified is defined by Equation (12) the result still holds on a set of larger probability 1κ ′ p(1 + C)n -δ with a reduced error, similar to the one in Theorem 29.

A Result for a Continuous Set of Jointly Diagonalizable Matrices

We now show a similar result when matrices in M can be jointly diagonalized. It turns out a faster algorithm can be used instead of Equation ( 11) with a reduced error and a larger probability event in the oracle inequality. Note that we no longer assume M is finite, so it can be parametrized by continuous parameters.

Suppose now the following holds, which means the matrices of M are jointly diagonalizable:

∃P ∈ O p (R) , M ⊆ P ⊤ Diag(d 1 , . . . , d p )P , (d i ) p i=1 ∈ (0, +∞) p . ( 18
)
Let P be the matrix defined in Assumption ( 18), Σ = P ΣP ⊤ and recall that A λ = K(K + nλI n ) -1 . Computations detailed in Appendix D show that the ideal penalty introduced in Equation ( 7) can be written as

∀M = P ⊤ Diag(d 1 , . . . , d p )P ∈ M, pen id (M ) = 2 tr A M • (Σ ⊗ I n ) np = 2 np   p j=1 tr(A pd j ) Σ j,j   . (19)
Equation ( 19) shows that under Assumption (18), we do not need to estimate the entire matrix Σ in order to have a good penalization procedure, but only to estimate the variance of the noise in p directions.

Definition 28 Let (e 1 , . . . , e p ) be the canonical basis of R p , (u 1 , . . . , u p ) be the orthogonal basis defined by ∀j ∈ {1, . . . , p}, u j = P ⊤ e j . We then define Σ HM = P Diag(a(u 1 ), . . . , a(u p ))P ⊤ , where for every j ∈ {1, . . . , p}, a(u j ) denotes the output of Algorithm 14 applied to Problem (Pu j ), and

M HM ∈ argmin M ∈M f M -y 2 2 + 2 tr A M • ( Σ HM ⊗ I n ) . ( 20
)
Theorem 29 Let α = 2, δ ≥ 2 and assume (13) and (18) hold true. Absolute constants L 2 > 0, and κ ′′ , a constant n 1 (δ) and an event Ω exist such that P( Ω) ≥ 1κ ′′ pn -δ and the following holds as soon as n ≥ n 1 (δ). First, on Ω,

1 np f M HM -f 2 2 ≤ 1 + 1 ln(n) 2 inf M ∈M 1 np f M -f 2 2 + L 2 tr(Σ)(2 + δ) 2 ln(n) 3 n . (21)
Second, an absolute constant L 4 exists such that

E 1 np f M HM -f 2 2 ≤ 1 + 1 ln(n) 2 E inf M ∈M 1 np f M -f 2 2 +L 4 tr(Σ)(2 + δ) 2 ln(n) 3 n + p n δ/2 f 2 2 np . ( 22
)
Theorem 29 is proved in Section F.

Comments on Theorems 26 and 29

Remark 30 Taking p = 1 (hence c(Σ) = 1 and tr(Σ) = σ 2 ), we recover Theorem 3 of [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF] as a corollary of Theorem 26.

Remark 31 (Scaling of (n, p)) When assumption (15) holds, Equation (16) implies the asymptotic optimality of the estimator f M when

c(Σ) 4 tr Σ p × p 3 (ln(n)) 3 n ≪ inf M ∈M 1 np f M -f 2 2 .
In particular, only (n, p) such that p 3 ≪ n/(ln(n)) 3 are admissible. When assumption (18) holds, the scalings required to ensure optimality in Equation ( 21) are more favorable:

tr Σ × (ln(n)) 3 n ≪ inf M ∈M 1 np f M -f 2 2
.

It is to be noted that p still influences the left hand side via tr Σ.

Remark 32 Theorems 26 and 29 are non asymptotic oracle inequalities, with a multiplicative term of the form 1 + o(1). This allows us to claim that our selection procedure is nearly optimal, since our estimator is close (with regard to the empirical quadratic norm) to the oracle one. Furthermore the term 1 + (ln(n)) -1 in front of the infima in Equations ( 16), ( 21), ( 17) and ( 22) can be further diminished, but this yields a greater remainder term as a consequence.

Remark 33 (On assumption (18)) Assumption (18) actually means all matrices in M can be diagonalized in a unique orthogonal basis, and thus can be parametrized by their eigenvalues as in Examples 3, 4 and 5.

In that case the optimization problem is quite easy to solve, as detailed in Remark 36. If not, solving (20) may turn out to be a hard problem, and our theoretical results do not cover this setting. However, it is always possible to discretize the set M or, in practice, to use gradient descent.

Compared to the setting of Theorem 26, assumption (18) allows a simpler estimator for the penalty (19), with an increased probability and a reduced error in the oracle inequality.

The main theoretical limitation comes from the fact that the probabilistic concentration tools used apply to discrete sets M (through union bounds). The structure of kernel ridge regression allows us to have a uniform control over a continuous set for the single-task estimators at the "cost" of n pointwise controls, which can then be extended to the multitask setting via (18). We conjecture Theorem 29 still holds without (18) as long as M is not "too large", which could be proved similarly up to some uniform concentration inequalities.

Note also that if M 1 , . . . , M K all satisfy (18) (with different matrices P k ), then Theorem 29 still holds for M = K k=1 M k with the penalty defined by Equation (20) with P = P k when M ∈ M k , and P( Ω) ≥ 1 -9Kp 2 n -δ , by applying the union bound in the proof.

Remark 34 (Relationship with the trace norm) Our approach relies on the minimization of Equation (2) with respect to f . [START_REF] Argyriou | Convex multi-task feature learning[END_REF] has shown that if we also minimize Equation (2) with respect to the matrix M subject to the constraint tr M -1 = 1, then we obtain an equivalent regularization by the nuclear norm (a.k.a. trace norm), which implies the prior knowledge that our p prediction functions may be obtained as the linear combination of r ≪ p basis functions. This situation corresponds to cases where the matrix M -1 is singular.

Note that the link between our framework and trace norm (i.e., nuclear norm) regularization is the same than between multiple kernel learning and the single task framework of [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF]. In the multi-task case, the trace-norm regularization, though efficient computationally, does not lead to an oracle inequality, while our criterion is an unbiased estimate of the generalization error, which turns out to be non-convex in the matrix M . While DC programming techniques (see, e.g., [START_REF] Gasso | Recovering sparse signals with non-convex penalties and dc programming[END_REF], and references therein) could be brought to bear to find local optima, the goal of the present work is to study the theoretical properties of our estimators, assuming we can minimize the cost function (e.g., in special cases, where we consider spectral variants, or by brute force enumeration).

Simulation Experiments

In all the experiments presented in this section, we consider the framework of Section 2 with X = R d , d = 4, and the kernel defined by ∀x, y ∈ X , k(x, y) = d j=1 e -|x j -y j | . The design points X 1 , . . . , X n ∈ R d are drawn (repeatedly and independently for each sample) independently from the multivariate standard Gaussian distribution. For every j ∈ {1, . . . , p}, f j (•) = m i=1 α j i k(•, z i ) where m = 4 and z 1 , . . . , z m ∈ R d are drawn (once for all experiments except in Experiment D) independently from the multivariate standard Gaussian distribution, independently from the design (X i ) 1≤i≤n . Thus, the expectations that will be considered are taken conditionally to the z i . The coefficients (α j i ) 1≤i≤m , 1≤j≤p differ according to the setting. Matlab code is available online. 1

Experiments

Five experimental settings are considered: A⌋ Various numbers of tasks: n = 10 and ∀i, j, α j i = 1, that is, ∀j,

f j = f A := m i=1 k(•, z i ).
= • • • = f 5 = f D , f 6 = • • • = f 10 = -f D .
E⌋ Comparison to cross-validation parameter selection: p = 5, Σ = 10 • I 5 , ∀j, f j = f A . The sample size is taken in {10, 50, 100, 250}.

Collections of Matrices

Two different sets of matrices M are considered in the Experiments A-C, following Examples 3 and 4:

M similar := M similar (λ, µ) = (λ + pµ)I p - µ p 11 ⊤ / (λ, µ) ∈ (0, +∞) 2
and M ind := {M ind (λ) = Diag(λ 1 , . . . , λ p ) / λ ∈ (0, +∞) p } .

In Experiment D, we also use two different sets of matrices, following Example 5:

M clus := I⊂{1,...,p},I / ∈{{1,...,p},∅} M I (λ, µ, µ) / (λ, µ) ∈ (0, +∞) 2 ∪ M similar and M interval := 1≤k≤p-1 M I (λ, µ, µ) / (λ, µ) ∈ (0, +∞) 2 , I = {1, . . . , k} ∪ M similar .
Remark 35 The set M clus contains 2 p -1 models, a case we will denote by "clustering". The other set, M interval , only has p models, and is adapted to the structure of the Experiment D. We call this setting "segmentation into intervals".

Estimators

In Experiments A-C, we consider four estimators obtained by combining two collections M of matrices with two formulas for Σ which are plugged into the penalty (7) (that is, either Σ known or estimated by Σ):

∀α ∈ {similar, ind} , ∀S ∈ Σ, Σ HM , f α,S := f M α,S = A M α,S y
where M α,S ∈ argmin

M ∈Mα 1 np y -f M 2 2 + 2 np tr (A M • (S ⊗ I n ))
and Σ HM is defined in Section 5.2. As detailed in Examples 3-4, f ind, Σ HM and f ind,Σ are concatenations of single-task estimators, whereas f similar, Σ HM and f similar,Σ should take advantage of a setting where the functions f j are close in F thanks to the regularization term

j,k f j -f k 2 F .
In Experiment D we consider the following three estimators, that depend on the choice of the collection M:

∀β ∈ {clus, interval, ind} , f β := f M β = A M β y where M β ∈ argmin M ∈M β 1 np y -f M 2 2 + 2 np tr A M • ( Σ ⊗ I n )
and Σ is defined by Equation ( 11).

In Experiment E we consider the estimator f similar, Σ HM . As explained in the following remark the parameters of the former estimator are chosen by optimizing (20), in practice by choosing a grid. We also consider the estimator f similar,CV where the parameters are selected by performing 5-fold cross-validation on the mentionned grid.

Remark 36 (Optimization of (20)) Thanks to Assumption (18) the optimization problem (20) can be solved easily. It suffices to diagonalize in a common basis the elements of M and the problem splits into several multi-task problems, each with one real parameter. The optimization was then done by using a grid on the real parameters, chosen such that the degree of freedom takes all integer values from 0 to n.

Remark 37 (Finding the jump in Algorithm 14) Algorithm 14 raises the question of how to detect the jump of df(λ), which happens around C = σ 2 . We chose to select an estimator C of σ 2 corresponding to the smallest index such that df( λ 0 ( C)) < n/2. Another approach is to choose the index corresponding to the largest instantaneous jump of df( λ 0 (C)) (which is piece-wise constant and non-increasing). This approach has a major drawback, because it sometimes selects a jump far away from the "real" jump around σ 2 , when the real jump consists of several small jumps. Both approaches gave similar results in terms of prediction error, and we chose the first one because of its direct link to the theoretical criterion given in Theorem 15.

Results

In each experiment, N = 1000 independent samples y ∈ R np have been generated. Expectations are estimated thanks to empirical means over the N samples. Error bars correspond to the classical Gaussian 95% confidence interval (that is, empirical standard-deviation over the N samples multiplied by 1.96/ √ N ). The results of Experiments A-C are reported in Figures 2345678. The results of Experiments C-E are reported in Tables 123. The p-values correspond to the classical Gaussian difference test, where the hypotheses tested are of the shape H 0 = {q > 1} against the hypotheses H 1 = {q ≤ 1}, where the different quantities q are detailed in Tables 23.

t 0.01 100 E[ f similar, Σ -f 2 / f ind, Σ -f 2 ]
1.80 ± 0.02 0.300 ± 0.003

E[ f similar, Σ -f 2 ] (2.27 ± 0.38) × 10 -2 0.357 ± 0.048 E[ f similar,Σ -f 2 ] (1.20 ± 0.28) × 10 -2 0.823 ± 0.080 E[ f ind, Σ -f 2 ] (1.26 ± 0.26) × 10 -2 1.51 ± 0.07 E[ f ind,Σ -f 2 ] (1.20 ± 0.24) × 10 -2 4.47 ± 0.13
Table 1: Results of Experiment C for the extreme values of t.

Comments

As expected, multi-task learning significantly helps when all f j are equal, as soon as p is large enough (Figure 1), especially for small n (Figure 6) and large noise-levels (Figure 8 and

E[ f similar, Σ -f 2 / f ind, Σ -f 2 ].
E[ f similar, Σ -f 2 / f ind, Σ -f 2 ]. q E [q] Std[q] p-value for H 0 = {q > 1} f clus -f 2 / f ind -f 2 0.668 0.294 < 10 -15 f interval -f 2 / f ind -f 2 0.660 0.270 < 10 -15 f interval -f 2 / f clus -f 2 1.00 0.165 0.50
Table 2: Clustering and segmentation (Experiment D).

q n E [q] Std[q] p-value for H 0 = {q > 1} f similar, Σ HM -f 2 / f similar,CV -f 2 10 0.35 0.46 < 10 -15 f similar, Σ HM -f 2 / f similar,CV -f 2 50 0.56 0.42 < 10 -15 f similar, Σ HM -f 2 / f similar,CV -f 2 100 0.71 0.34 < 10 -15 f similar, Σ HM -f 2 / f similar,CV -f 2 250 0.87 0.19 < 10 -15
Table 3: Comparison of our method to 5-fold cross-validation (Experiment E).

Table 1). Increasing the number of tasks rapidly reduces the quadratic error with multi-task estimators (Figure 2) contrary to what happens with single-task estimators (Figure 3). A noticeable phenomenon also occurs in Figure 2 and even more in Figure 3: the estimator f ind,Σ (that is, obtained knowing the true covariance matrix Σ) is less efficient than f ind, Σ where the covariance matrix is estimated. It corresponds to the combination of two facts: (i) multiplying the ideal penalty by a small factor 1 < C n < 1+ o(1) is known to often improve performances in practice when the sample size is small (see Section 6.3.2 of [START_REF] Arlot | Model selection by resampling penalization[END_REF], and (ii) minimal penalty algorithms like Algorithm 14 are conjectured to overpenalize slightly when n is small or the noise-level is large [START_REF] Lerasle | Optimal model selection in density estimation[END_REF] (as confirmed by Figure 7). Interestingly, this phenomenon is stronger for single-task estimators (differences are smaller in Figure 2) and disappears when n is large enough (Figure 5), which is consistent with the heuristic motivating multi-task learning: "increasing the number of tasks p amounts to increase the sample size".

Figures 4 and5 show that our procedure works well with small n, and that increasing n does not seem to significantly improve the performance of our estimators, except in the single-task setting with Σ known, where the over-penalization phenomenon discussed above disappears.

Table 2 shows that using the multitask procedure improves the estimation accuracy, both in the clustering setting and in the segmentation setting. The last line of Table 2 does not show that the clustering setting improves over the "segmentation into intervals" one, which was awaited if a model close to the oracle is selected in both cases.

Table 3 finally shows that our parameter tuning procedure outperforms 5-fold crossvalidation.

Conclusion and Future Work

This paper shows that taking into account the unknown similarity between p regression tasks can be done optimally (Theorem 26). The crucial point is to estimate the p × p covariance matrix Σ of the noise (covariance between tasks), in order to learn the task similarity matrix M . Our main contributions are twofold. First, an estimator of Σ is defined in Section 4, where non-asymptotic bounds on its error are provided under mild assumptions on the mean of the sample (Theorem 20). Second, we show an oracle inequality (Theorem 26), more particularly with a simplified estimation of Σ and increased performances when the matrices of M are jointly diagonalizable (which often corresponds to cases where we have a prior knowledge of what the relations between the tasks would be). We do plan to expand our results to larger sets M, which may require new concentration inequalities and new optimization algorithms.

Simulation experiments show that our algorithm works with reasonable sample sizes, and that our multi-task estimator often performs much better than its single-task counterpart. Up to the best of our knowledge, a theoretical proof of this point remains an open problem that we intend to investigate in a future work. n k=1 T i,k f (•, k). We then have:

g, g G = p i=1 p j=1 M i,j g(•, i), g(•, j) F = p i=1 p j=1 p k=1 p l=1 M i,j T i,k T j,l f (•, k), f (•, l) F = p j=1 p k=1 p l=1 T l,j f (•, k), f (•, l) F p i=1 M j,i T i,k = p j=1 p k=1 p l=1 T l,j f (•, k), f (•, l) F (M • T ) j,k = p k=1 p l=1 T l,j f (•, k), f (•, l) F p j=1 T l,j (M • T ) j,k = p k=1 p l=1 f (•, k), f (•, l) F (T • M • T ) k,l = p k=1 f (•, k) 2 F .
This shows that g, g G ≥ 0 and that g, g

G = 0 ⇒ f = 0 ⇒ g = 0. g, Φ(x, l) G = p j=1 p i=1 M j,i g(•, j), Φ(x, l) i F = p j=1 p i=1 p m=1 M j,i M -1 i,m δ m,l g(•, j), Φ(x) F = p j=1 p m=1 (M • M -1 ) j,m δ m,l g(x, j) = p j=1 δ j,l g(x, j) = g(x, l) .
Thus we can write:

k((x, i), (y, j)) = Φ(x, i), Φ(y, j) G = p h=1 p h ′ =1 M h,h ′ M -1 h,i Φ(x), M -1 h ′ ,j Φ(y) F = p h=1 p h ′ =1 M h,h ′ M -1 h,i M -1 h ′ ,j K(x, y) = p h=1 M -1 h,i (M • M -1 ) h,j K(x, y) = p h=1 M -1 h,i δ h,j K(x, y) = M -1 i,j K(x, y) .
• span {e 1 + • • • + e p } corresponding to eigenvalue p,

• span {e 2e 1 , . . . , e pe 1 } corresponding to eigenvalue 0.

Thus, with µ = λ + pµ we can diagonalize in an orthonormal basis any matrix M λ,µ ∈ M as M = P ⊤ D λ, µ P , with D = D λ, µ = Diag{λ, µ, . . . , µ}. Les us also diagonalise in an orthonormal basis K: K = Q ⊤ ∆Q, ∆ = Diag{µ 1 , . . . , µ n }. Thus we can write (see Properties 38 and 39 for basic properties of the Kronecker product):

A M = A M λ,µ = (P ⊤ ⊗ Q ⊤ ) (D -1 ⊗ ∆) (D -1 ⊗ ∆) + npI np -1 (P ⊗ Q) .
We can then note that (D -1 ⊗ ∆) (D -1 ⊗ ∆) + npI np -1 is a diagonal matrix, whose diagonal entry of index (j -1)n + i (i ∈ {1, . . . , n}, j ∈ {1, . . . , p}) is

µ i µ i +npλ if j = 1 , µ i µ i +np µ if j > 1 .
We can now compute both bias and variance.

Bias: We can first remark that (P ⊤ ⊗ Q ⊤ ) = (P ⊗ Q) ⊤ is an orthogonal matrix and that P × 1 = (1, 0, . . . , 0) ⊤ . Thus, as in this setting f 1 = • • • = f p , we have f = 1⊗(f 1 (X 1 ), . . . , f 1 (X n )) ⊤ and (P ⊤ ⊗Q ⊤ )f = (1, 0, . . . , 0) ⊤ ⊗Q(f 1 (X 1 ), . . . , f 1 (X n )) ⊤ .

To keep notations simple we note Q(f 1 (X 1 ), . . . , f 1 (X n )) ⊤ := (g 1 , . . . , g n ) ⊤ . Thus

(A M -I np )f 2 2 = (P ⊗ Q) ⊤ (D -1 ⊗ K) (D -1 ⊗ K) + npI np -1 -I np (P ⊗ Q)f 2 2 = (D -1 ⊗ ∆) (D -1 ⊗ ∆) + npI np -1 -I np × (1, 0, . . . , 0) ⊤ ⊗ (g 1 , . . . , g n ) ⊤ 2 2 .
As only the first n terms of (P ⊗ Q)f are non-zero we can finally write

(A M -I np )f 2 2 = n i=1 npλ µ i + npλ 2 g 2 i .
Variance: First note that

(P ⊗ Q)(Σ ⊗ I n )(P ⊗ Q) ⊤ = (P ΣP ⊤ ⊗ I n ) .
This norm has several interesting properties, some of which we will use are stated below.

Property 41 The norm ||| • ||| is a matricial norm: ∀(A, B) ∈ S p (R) 2 , |||AB||| ≤ |||A||||||B|||.

We will use the following result, which is a consequence of the preceding Property.

∀S ∈ S p (R), ∀T ∈ S ++ p (R), |||T -1 2 ST -1 2 ||| ≤ |||S||||||T -1 ||| .
We also have:

Proposition 42 ∀Σ ∈ S p (R), |||Σ ⊗ I n ||| = |||Σ||| .
Proof We can diagonalize Σ in an orthonormal basis: ∃U ∈ O n (R), ∃D = Diag(µ 1 , . . . , µ p ), Σ = U ⊤ DU . We then have, using the properties of the Kronecker product:

Σ ⊗ I n = (U ⊤ ⊗ I n )(D ⊗ I n )(U ⊗ I n ) = (U ⊗ I n ) ⊤ (D ⊗ I n )(U ⊗ I n ) .
We just have to notice that U ⊗ I n ∈ O np (R) and that: ) .

D ⊗ I n = Diag(µ 1 , . . . ,
This norm can also be written in other forms:

Property 43 If M ∈ M n (R),

E.2 The Proof

We now give a proof of Theorem 20, using Lemmas 46, 48 and 49, which are stated and proved in Section E.3. The outline of the proof is the following:

1. Apply Theorem 15 to problem (10) for every z ∈ Z in order to 2. control sζ ∞ with a large probability, where s, ζ ∈ R p(p+1)/2 are defined by s := (Σ 1,1 , . . . , Σ p,p , Σ 1,1 + Σ 2,2 + 2Σ 1,2 , . . . , Σ i,i + Σ j,j + 2Σ i,j , . . .) and ζ := (a(e 1 ), . . . , a(e p ), a(e 1 + e 2 ), . . . , a(e 1 + e p ), a(e 2 + e 3 ), . . . , a(e p-1 + e p )) .

3. Deduce that Σ = J(ζ) is close to Σ = J(s) by controlling the Lipschitz norm of J.

Proof 1. Apply Theorem 15: We start by noticing that Assumption (13) actually holds true with all λ 0,j equal. Indeed, let (λ 0,j ) 1≤j≤p be given by Assumption (13) and define λ 0 := min j=1,...,p λ 0,j . Then, λ 0 ∈ (0, +∞) and df(λ 0 ) since all λ 0,j satisfy these two conditions. For the last condition, remark that for every j ∈ {1, . . . , p}, λ 0 ≤ λ 0,j and λ → (A λ -I)F e j 2 2 is a nonincreasing function (as noticed in Arlot and Bach, 2011 for instance), so that

1 n (A λ 0 -I n )F e j 2 2 ≤ 1 n (A λ 0,j -I n )F e j 2 2 ≤ Σ j,j ln(n) n . (23)
In particular, Equation ( 8) holds with d n = 1 for problem (10) whatever z ∈ {e 1 , . . . , e p }.

Let us now consider the case z = e i + e j with i = j ∈ {1, . . . , p}. Using Equation ( 23) and that F e i +e j = F e i + F e j , we have

(B λ 0 -I n )F e i +e j 2 2 ≤ (B λ 0 -I n )F e i 2 2 + (B λ 0 -I n )F e j 2 2 +2 (B λ 0 -I n )F e i , (B λ 0 -I n )F e j .
The last term is bounded as follows:

2 (B λ 0 -I n )F e i , (B λ 0 -I n )F e j ≤ 2 (B λ 0 -I n )F e i • (B λ 0 -I n )F e j ≤ 2 n ln(n) Σ i,i Σ j,j ≤ n ln(n)(Σ i,i + Σ j,j ) ≤ (1 + c(Σ)) n ln(n)(Σ i,i + Σ j,j + 2Σ i,j ) = (1 + c(Σ)) n ln(n)σ 2 e i +e j ,
because Lemma 46 shows 2(Σ i,i + Σ j,j ) ≤ (1 + c(Σ))(Σ i,i + Σ j,j + 2Σ i,j ) .

Therefore, Equation (8) holds with d n = 1 + c(Σ) for problem (10) whatever z ∈ Z. 2. Control sζ ∞ : Let us define

η 1 := β(2 + δ)(1 + c(Σ)) ln(n) n .
By Theorem 15, for every z ∈ Z, an event Ω z of probability greater than 1n -δ exists on which, if n ≥ n 0 (δ), (1

-η 1 )σ 2 z ≤ a(z) ≤ (1 + η 1 )σ 2 z . So, on Ω := z∈Z Ω z , ζ -s ∞ ≤ η 1 s ∞ , (24)
and P(Ω) ≥ 1p(p + 1)/2 × n -δ by the union bound. Let

Σ ∞ := sup i,j |Σ i,j | and C 1 (p) := sup Σ∈Sp(R) Σ ∞ |||Σ||| .
Since s ∞ ≤ 4 Σ ∞ and C 1 (p) = 1 by Lemma 48, Equation ( 24) implies that on Ω, By Lemma 49, C 2 (p) ≤ 3 2 p. By Equation ( 25), on Ω,

ζ -s ∞ ≤ 4η 1 Σ ∞ ≤ 4η 1 |||Σ||| . ( 25
||| Σ -Σ||| = |||J(ζ) -J(s)||| ≤ C 2 (p) ζ -s ∞ ≤ 4η 1 C 2 (p)|||Σ||| . ( 26
) Since |||Σ -1 2 ΣΣ -1 2 -I p ||| = |||Σ -1 2 (Σ -Σ)Σ -1 2 ||| ≤ |||Σ -1 ||||||Σ -Σ||| ,
and |||Σ||||||Σ -1 ||| = c(Σ), Equation ( 26) implies that on Ω,

|||Σ -1 2 ΣΣ -1 2 -I p ||| ≤ 4η 1 C 2 (p)|||Σ||||||Σ -1 ||| = 4η 1 C 2 (p)c(Σ) ≤ 6η 1 pc(Σ) .
To conclude, Equation ( 14) holds on Ω with

η = 6pc(Σ)β(2 + δ)(1 + c(Σ)) ln(n) n ≤ L 1 (2 + δ)p ln(n) n c(Σ) 2 (27)
for some numerical constant L 1 .

Remark 44 As stated in [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF], we need n 0 (δ)/ ln(n 0 (δ)) ≥ 504 and n 0 (δ)/ ln(n 0 (δ)) ≥ 24(290 + δ).

Remark 45 To ensure that the estimated matrix Σ is positive-definite we need that η < 1, that is, n ln(n) > 6β(2 + δ)pc(Σ) (1 + c(Σ)) .

E.3 Useful Lemmas

Lemma 46 Let p ≥ 1, Σ ∈ S ++ p (R) and c(Σ) its condition number. Then,

∀1 ≤ i < j ≤ p , Σ i,j ≥ - c(Σ) -1 c(Σ) + 1 Σ i,i + Σ j,j 2 , ( 28
)
Remark 47 The proof of Lemma 46 shows the constant c(Σ)-1 c(Σ)+1 cannot be improved without additional assumptions on Σ.

Proof It suffices to show the result when p = 2. Indeed, (28) only involves 2×2 submatrices Σ(i, j) ∈ S ++ 2 (R) for which

1 ≤ c( Σ) ≤ c (Σ) hence 0 ≤ c( Σ) -1 c( Σ) + 1 ≤ c(Σ) -1 c(Σ) + 1 . So, some θ ∈ R exists such that Σ = |||Σ|||R ⊤ θ DR θ where R θ := cos(θ) sin(θ) -sin(θ) cos(θ) D = 1 0 0 λ and λ := 1 c(Σ) . Therefore, Σ = |||Σ||| cos 2 (θ) + λ sin 2 (θ) 1-λ 2 sin(2θ) 1-λ 2 sin(2θ) λ cos 2 (θ) + sin 2 (θ) .
So, Equation ( 28) is equivalent to

(1 -λ) sin(2θ) 2 ≥ - 1 -λ 1 + λ 1 + λ 2 ,
which holds true for every θ ∈ R, with equality for θ ≡ π/2 (mod. π).

Lemma 48 For every p

≥ 1, C 1 (p) := sup Σ∈Sp(R) Σ ∞ |||Σ||| = 1 . Proof With Σ = I p we have Σ ∞ = |||Σ||| = 1, so C 1 (p) ≥ 1. Let us introduce (i, j) such that |Σ i,j | = Σ ∞ .
We then have, with e k being the k th vector of the canonical basis of R p ,

|Σ i,j | = |e ⊤ i Σe j | ≤ |e ⊤ i Σe i | 1/2 |e ⊤ j Σe j | 1/2 ≤ ( Σ 1/2 2 ) 2 . Lemma 49 For every p ≥ 1, let C 2 (p) := sup ζ∈R p(p+1)/2 |||J(ζ)||| ζ ∞ . Then, p 4 ≤ C 2 (p) ≤ 3 2 p .
Proof For the lower bound, we consider

ζ 1 = (1, . . . , 1 p times , 4, . . . , 4 p(p-1) 2 times ) , then J(ζ 1 ) =    1 . . . 1 . . . . . . . . . 1 . . . 1    so that |||J(ζ)||| = p and ζ ∞ = 4.
For the upper bound, we have for every ζ ∈ R p(p+1)/2 and z ∈ R p such that z 2 = 1

z ⊤ J(ζ)z = 1≤i,j≤p z i z j J(ζ) i,j ≤ 1≤i,j≤p |z i | |z j | |J(ζ)| ≤ J(ζ) ∞ z 2 1 . By definition of J, J(ζ) ∞ ≤ 3/2 ζ ∞ . Remarking that z 2 1 ≤ p z 2 2 yields the result. Lemma 54 Let C A = 2, C B = 1, C C = 2, C D = 1, C E = 306.25, C F = 306.25 .
Suppose that (18) holds. Then P(Ω

x (M, C A , C B , C C , C D , C E , C F )) ≥ 1 -pe 1027+ln(n) e -x .
Suppose that (15) holds. Then P(Ω

x (M, C A , C B , C C , C D , C E , C F )) ≥ 1 -6p card(M)e -x . . Proof First common step. Let M ∈ M, P M ∈ O p (R) such that M = P ⊤ M DP M , with D = Diag(d 1 , . . . , d p )
. We can write:

A M = A d 1 ,...,dp = (P M ⊗ I n ) ⊤ (D -1 ⊗ K) D -1 ⊗ K + npI np -1 (P M ⊗ I n ) = Q ⊤ A d 1 ,...,dp Q , with Q = P M ⊗ I n and A d 1 ,...,dp = (D -1 ⊗ K)(D -1 ⊗ K + npI np ) -1 . Remark that A d 1 ,...,
= Qε = ( ε 1 ⊤ , . . . , ε p ⊤ ) ⊤ and f = Qf = ( f 1 ⊤ , . . . , f p ⊤ ) ⊤ we can write |δ 1 (M )| = ε, A d 1 ,...,dp ε -E ε, A d 1 ,...,dp ε , |δ 2 (M )| = A d 1 ,...,dp ε 2 2 -E A d 1 ,...,dp ε 2 2 , |δ 3 (M )| = 2 A d 1 ,...,dp ε, ( A d 1 ,...,dp -I np ) f , |δ 4 (M )| = 2 ε, (I np -A d 1 ,...,dp ) f .
We can see that the quantities δ i decouple, therefore

|δ 1 (M )| = p i=1 ε i , A pd i ε i -E [ ε i , A pd i ε ] , |δ 2 (M )| = p i=1 A pd i ε i 2 2 -E A pd i ε i 2 2 , |δ 3 (M )| = p i=1 2 A pd i ε i , (A pd i -I n ) f i , |δ 4 (M )| = p i=1 2 ε i , (I n -A pd i ) f i .
Supposing (18). Assumption (18) implies that the matrix P used above is the same for all the matrices M of M. Using Lemma 9 of [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF], where we have p concentration results on the sets Ω i , each of probability at least 1e 1027+ln(n) e -x we can state that, on the set p i=1 Ω i , we have uniformly on

M |δ 1 (M )| ≤ p i=1 θ 1 Var[ ε i ] tr(A ⊤ pd i A pd i ) + (C A + C B θ -1 1 )x Var[ ε i ] , |δ 2 (M )| ≤ p i=1 θ 2 Var[ ε i ] tr(A ⊤ pd i A pd i ) + (C C + C D θ -1 2 )x Var[ ε i ] , |δ 3 (M )| ≤ p i=1 θ 3 (I n -A pd i ) f i 2 2 + C E θ -1 3 x Var[ ε i ] , |δ 4 (M )| ≤ p i=1 θ 4 (I n -A pd i ) f i 2 2 + C F θ -1 4 x Var[ ε i ] .
Supposing (15). We can use Lemma 8 of [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF] where we have p concentration results on the sets Ω j,M , each of probability at least 1 -6e -x we can state that, on the set p j=1 M ∈M Ω i , we have uniformly on M the same inequalities written above.

Final common step. To conclude, it suffices to see that for every i ∈ {1, . . . , p}, Var[ ε i ] ≤ |||Σ|||.

F.2 Intermediate Result

We first prove a general oracle inequality, under the assumption that the penalty we use (with an estimator of Σ) does not underestimate the ideal penalty (involving Σ) too much.

Proposition 55 Let C A , C B , C C , C D , C E ≥ 0 be fixed constants, γ > 0, θ S ∈ [0, 1/4) and K S ≥ 0. On Ω γ ln(n) (M, C A , C B , C C , C D , C E ), for every S ∈ S ++ p (R) such that tr A Mo(S) • ((S -Σ) ⊗ I n ) ≥ -θ S tr A Mo(S) • (Σ ⊗ I n ) inf M ∈M b(M ) + v 2 (M ) + K S ln(n)|||Σ||| v 1 (M ) (34)
and for every θ ∈ (0, (1 -4θ S )/2), we have:

1 np f Mo(S) -f 2 2 ≤ 1 + 2θ 1 -2θ -4θ S inf M ∈M 1 np F M -F 2 2 + 2 tr (A M • ((S -Σ) + ⊗ I n )) np + 1 1 -2θ -4θ S (2C A + 3C C + 6C D + 6C E + 2 θ (C B + C F ))γ + θ S K S 4 ln(n)|||Σ||| np (35)

F.3 The Proof Itself

We now show Theorem 26 as a consequence of Proposition 55. It actually suffices to show that Σ does not underestimate Σ too much, and that the second term in the infimum of Equation ( 35) is negligible in front of the quadratic error (np) -1 f Mf 2 . Proof On the event Ω introduced in Theorem 20, Equation ( 14) holds. Let γ = pc(Σ) (1 + c(Σ)) .

By Lemma 56 below, we have:

inf M ∈M b(M ) + v 2 (M ) + K S ln(n)|||Σ||| v 1 (M ) ≥ 2 K S ln(n)|||Σ||| n tr(Σ) .
We supposed Assumption (15) holds. Using elementary algebra it is easy to show that, for every symmetric positive definite matrices A, M and N of size p, M N implies that tr(AM ) ≥ tr(AN ). In order to have M o ( Σ) satisfying Equation ( 34), Theorem 20 shows that it suffices to have, for every θ S > 0, We now take θ S = θ = (9 ln(n)) -1 . Let Ω be the set given by Theorem 20. Using Equation (35) and requiring that ln(n) ≥ 6 we get, on the set Ω = Ω∩Ω (α+δ) ln(n) (M, C A , C B , C C , C D , C E , C F ) of probability 1 -(p(p + 1)/2 + 6pC)n -δ , using that α ≥ 2:

1 np f M -f 2 ≤ 1 + 1 ln(n) inf M ∈M    1 np f M -f 2 2 + 2 tr A M • (( Σ -Σ) + ⊗ I n ) np    + 1 - 2 3 ln(n) -1 2C A + 3C C + 6C D + 6C E + ln(n) 18C B + 18C F + 729β 2 γ 2 tr(Σ) 2 4|||Σ||| 2 × (α + δ) 2 ln(n) 2 |||Σ||| np .
Using Equation ( 27) and defining

η 2 := 12β(α + δ)γ ln(n) n ,
we get

1 np f M -f 2 ≤ 1 + 1 ln(n) inf M ∈M 1 np f M -f 2 2 + η 2 tr(A M • (Σ ⊗ I n )) np + 1 - 2 3 ln(n) -1 2C A + 3C C + 6C D + 6C E + ln(n) 18C B + 18C F + 729β 2 γ 2 tr(Σ) 2 4|||Σ||| 2 ×(α + δ) 2 ln(n) 2 |||Σ||| np . (42)
Now, to get a classical oracle inequality, we have to show that η 2 v 1 (M ) = η 2 tr(A M • (Σ ⊗ I n )) is negligible in front of f Mf 2 . Lemma 56 ensures that:

∀M ∈ M , ∀x ≥ 0 , 2 x|||Σ||| n tr(Σ) v 1 (M ) ≤ v 2 (M ) + x|||Σ||| .
With 0 < C n < 1, taking x to be equal to 72β 2 ln(n)γ 2 tr(Σ)/(C n |||Σ|||) leads to

η 2 v 1 (M ) ≤ 2C n v 2 (M ) + 72β 2 ln(n)γ 2 tr(Σ) C n . ( 43
)
Then, since v 2 (M ) ≤ v 2 (M ) + b(M ) and using also Equation ( 36), we get

v 2 (M ) ≤ f M -f 2 2 + |δ 2 (m)| + |δ 3 (M )| .
On Ω we have that for every θ ∈ (0, 1), using Equation ( 31) and (32),

η 2 v 1 (M ) ≤ 1 + 4C n θ 1 + 2θ f M -f 2 2 + 2C n C C + (C D + C E )θ -1 1 + 2θ (α + δ) ln(n)|||Σ||| + 72β 2 ln(n)γ 2 tr(Σ) C n .
Taking θ = 1/2 then leads to From this we can deduce Equation ( 16) by noting that γ ≤ 2pc(Σ) 2 . Finally we deduce an oracle inequality in expectation by noting that if n -1 f Mf 2 ≤ R n,δ on Ω, using Cauchy-Schwarz inequality

η 2 v 1 (M ) ≤ (1 + C n ) f M -f
E 1 np f M -f 2 2 = E 1 Ω np f M -f 2 2 + E 1 Ω c np f M -f 2 2 ≤ E [R n,δ ] + 1 np 4p(p + 1) + 6pC n δ E f M -f 4 2 . ( 44
)
We can remark that, since

|||A M ||| ≤ 1, f M -f 2 2 ≤ 2 A M ε 2 2 + 2 (I np -A M )f 2 2 ≤ 2 ε 2 2 + 8 f 2 2 . So E f M -f 4 2 ≤ 12 np|||Σ||| + 4 f 2 2 2 ,
together with Equation ( 42) and Equation ( 44), induces Equation ( 17), using that for some constant L 3 > 0, To get to the oracle inequality in expectation we use the same technique than above, but we note that P( Ω c ) ≤ L 4 × p/n δ/2 . We can finally define the constant L 4 by: L 3 tr(Σ)(2 + δ) 2 p ln(n) 3 np + p n δ/2 |||Σ||| ≤ L 4 γ 2 tr(Σ)(α + δ) 2 p ln(n) 3 np .

  The number of tasks is varying: p ∈ {2k / k = 1, . . . , 25}. The covariance matrix is Σ = 10 • I p . B⌋ Various sample sizes: p = 5, ∀j, f j = f A and Σ = Σ B has been drawn (once for all) from the Whishart W (I 5 , 10, 5) distribution; the condition number of Σ B is c(Σ B ) ≈ 22.05. The only varying parameter is n ∈ {50k / k = 1, . . . , 20}. C⌋ Various noise levels: n = 100, p = 5 and ∀j, f j = f A . The varying parameter is Σ = Σ C,t := 5t • I 5 with t ∈ {0.2k / k = 1, . . . , 50}. We also ran the experiments for t = 0.01 and t = 100. D⌋ Clustering of two groups of functions: p = 10, n = 100, Σ = Σ E has been drawn (once for all) from the Whishart W (I 10 , 20, 10) distribution; the condition number of Σ E is c(Σ E ) ≈ 24.95. We pick the function f D := m i=1 α i k(•, z i ) by drawing (α 1 , . . . , α m ) and (z 1 , . . . , z m ) from standard multivariate normal distribution (independently in each replication) and finally f 1

Figure 1 :

 1 Figure1: Increasing the number of tasks p (Experiment A), improvement of multi-task compared to single-task:E[ f similar, Σf 2 / f ind, Σf 2 ].

Figure 2 :Figure 3 :

 23 Figure 2: Increasing the number of tasks p (Experiment A), quadratic errors of multi-task estimators (np) -1 E[ f similar,Sf 2 ]. Blue: S = Σ. Red: S = Σ.

Figure 4 :Figure 5 :

 45 Figure 4: Increasing the sample size n (Experiment B), quadratic errors of multi-task estimators (np) -1 E[ f similar,Sf 2 ]. Blue: S = Σ. Red: S = Σ.

Figure 6 :Figure 7 :

 67 Figure 6: Increasing the sample size n (Experiment B), improvement of multi-task compared to single-task: E[ f similar, Σf 2 / f ind, Σf 2 ].

Figure 8 :

 8 Figure 8: Increasing the signal-to-noise ratio (Experiment C), improvement of multi-task compared to single-task: E[ f similar, Σf 2 / f ind, Σf 2 ].

2 -

 2 |δ 2 (M )|+|δ 3 (M )| ≤ 2θ f Mf 2 |δ 2 (M )| -|δ 3 (M )| +(C C +(C D +C E )θ -1 )(α+δ) ln(n)|||Σ||| , which leads to |δ 2 (M )| + |δ 3 (M )| ≤ 2θ (C D + C E )θ -1 1 + 2θ (α + δ) ln(n)|||Σ||| .Now, combining this equation with Equation (43), we get

+

  C n (C C + 2(C D + C E ))(α + δ) ln(n)|||Σ||| + 72β 2 ln(n)γ tr(Σ) C n .We now take C n = 1/ ln(n). We now replace the constants C A , C B , C C , C D , C E , C F by their values in Lemma 54 and we get, for some constant L 2 , (n)γ 2 tr(Σ)C n ≤ L 2 ln(n)γ 2 tr(Σ) 2 |||Σ||| 2

.

  Lemma 56 Let n, p ≥ 1 be two integers, x ≥ 0 and Σ ∈ S ++ p (R). Then,inf A∈Mnp(R),|||A|||≤1 tr(A ⊤ A • (Σ ⊗ I n )) + x|||Σ||| tr(A • (Σ ⊗ I n )) ≥ 2 x|||Σ||| n tr(Σ)Proof First note that the bilinear form onM np (R), (A, B) → tr(A ⊤ B • (Σ ⊗ I n )) is a scalar product. By Cauchy-Schwarz inequality, for every A ∈ M np (R), tr(A • (Σ ⊗ I n )) 2 ≤ tr(Σ ⊗ I n ) tr(A ⊤ A • (Σ ⊗ I n )) .Thus, since tr(Σ⊗ I n ) = n tr(Σ), if c = tr(A • (Σ ⊗ I n )) > 0, tr(A ⊤ A • (Σ ⊗ I n )) ≥ c 2 n tr(Σ).

  the operator norm M 2 := sup t∈R n \{0}

		M t 2 t 2	is equal
	to the greatest singular value of M :	ρ(M ⊤ M ). Henceforth, if S is symmetric, we have
	|||S||| = S 2

  dp is block-diagonal, with diagonal blocks being B d 1 , . . . , B dp using the notations of Section 3. With ε

Matlab code can be found at http://www.di.ens.fr/ ~solnon/multitask_minpen_en.html.

Acknowledgments

This paper was supported by grants from the Agence Nationale de la Recherche (Detect project, reference ANR-09-JCJC-0027-01) and from the European Research Council (SIERRA Project ERC-239993).

We give in Appendix the proofs of the different results stated in Sections 2, 4 and 5. The proofs of our main results are contained in Sections E and F.

Appendix A. Proof of Proposition 8

Proof It is sufficient to show that •, • G is positive-definite on G. Take g ∈ G and S = (S i,j ) 1≤i≤j≤p the symmetric postive-definite matrix of size p verifying S 2 = M , and denote T = S -1 = (T i,j ) 1≤i,j≤p . Let f be the element of G defined by ∀i ∈ {1 . . . p}, g(•, i) =

Appendix B. Proof of Corollary 9

Proof If (x, j) ∈ X × {1, . . . , p}, the application (f 1 , . . . , f p ) → f j (x) is clearly continuous. We now show that (G, •, • G ) is complete. If (g n ) n∈N is a Cauchy sequence of G and if we define, as in Section A, the functions f n by ∀n ∈ N, ∀i ∈ {1 . . . p}, g n (•, i) = p k=1 T i,k f n (•, k). The same computations show that (f n (•, i)) n∈N are Cauchy sequences of F, and thus converge. So the sequence (f n ) n∈N converges in G, and (g n ) n∈N does likewise.

Appendix C. Proof of Proposition 11

Proof We define

with δ i,j = 1 i=j being the Kronecker symbol, that is, δ i,j = 1 if i = j and 0 otherwise. We now show that Φ is the feature function of the RKHS. For g ∈ G and (x, l) ∈ X × {1, . . . , p}, we have:

Appendix D. Computation of the Quadratic Risk in Example 12

We consider here that f 1 = • • • = f p . We use the set M similar :

Using the estimator f M = A M y we can then compute the quadratic risk using the biasvariance decomposition given in Equation ( 36):

Les us denote by (e 1 , . . . , e p ) the canonical basis of R p . The eigenspaces of p -1 11 ⊤ are:

We can also note that Σ := P ΣP ⊤ is a symmetric positive definite matrix, with positive diagonal coefficients. Thus we can finally write

As noted at the end of Example 12 this leads to an oracle which has all its p functions equal.

D.1 Proof of Equation

. The computations detailed above also show that the ideal penalty introduced in Equation ( 7) can be written as

Appendix E. Proof of Theorem 20

Theorem 20 is proved in this section, after stating some classical linear algebra results (Section E.1).

E.1 Some Useful Tools

We now give two properties of the Kronecker product, and then introduce a useful norm on S p (R), upon which we give several properties. Those are the tools needed to prove Theorem 20.

Property 38 The Kronecker product is bilinear, associative and for every matrices A, B, C, D such that the dimensions fit,

Definition 40 We now introduce the norm ||| • ||| on S p (R), which is the modulus of the eigenvalue of largest magnitude, and can be defined by

Appendix F. Proof of Theorem 26

The proof of Theorem 26 is similar to the proof of Theorem 3 in [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF]. We give it here for the sake of completeness. We also show how to adapt its proof to demonstrate Theorem 29. The two main mathematical results used here are Theorem 20 and a gaussian concentration inequality from [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF].

F.1 Key Quantities and their Concentration Around their Means Definition 50 We introduce, for S ∈ S ++ p (R),

C F be fixed nonnegative constants. For every x ≥ 0 we define the event

on which, for every M ∈ M and θ 1 , θ 2 , θ 3 , θ 4 ∈ (0, 1]:

Of key interest is the concentration of the empirical processes δ i , uniformly over M ∈ M. The following Lemma introduces such a result, when M contains symmetric matrices parametrized with their eigenvalues (with fixed eigenvectors).

Proof The proof of Proposition 55 is very similar to the one of Proposition 5 in [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF]. First, we have

Combining Equation ( 29) and ( 37), we get:

On the event Ω γ ln(n) , for every θ ∈ (0, 1] and M ∈ M, using Equation ( 30) and ( 33) with

Using Equation ( 31) and ( 32) with θ 2 = θ 3 = 1/2 we get that for every M ∈ M Equation

Combining Equation ( 39) and (40), we get

With Equation (38), and with

Using Equation ( 34) we can state that

which then leads to Equation (35) using Equation ( 40) and (41).

Therefore inf

≥ 2 x|||Σ||| n tr(Σ) .

F.4 Proof of Theorem 29

We now prove Theorem 29, first by proving that Σ HM leads to a sharp enough approximation of the penalty.

Lemma 57 Let Σ HM be defined as in Definition 28, α = 2, κ > 0 be the numerical constant defined in Theorem 15 and assume (13) and (18) hold. For every δ ≥ 2, a constant n 0 (δ), an absolute constant L 1 > 0 and an event Ω exist such that P(Ω HM ) ≥ 1pn -δ and for every n ≥ n 0 (δ), on Ω HM , for every

where

Proof Let P be defined by ( 18). Let M ∈ M, and (d 1 , . . . , d p ) ∈ (0, +∞) p such that M = P ⊤ Diag(d 1 , . . . , d p )P . Thus, as shown in Section D, we have with Σ = P ΣP ⊤ :

tr(A pd j ) Σ j,j .

let σ j be defined as in Definition 28 (and thus Σ HM = P Diag( σ 1 , . . . , σ p )P ⊤ ), we then have by Theorem 15 that for every j ∈ {1, . . . , p} an event Ω j of probability 1κn -δ exists such that on Ω j | Σ j,jσ j | ≤ η Σ j,j . Since


Thierry Dumont
email: thierry.dumont@math.u-psud.fr

Sylvain Le Corff
email: y.le-corff@warwick.ac.uk

Supplement paper to 'Nonparametric estimation in hidden Markov models'

This document is a supplementary material to the article "Nonparametric estimation using partially observed Markov chains". It provides additional proofs of some technical results given in the original paper. Section 1 recalls the model, the definitions and the assumptions used in the paper. Section 2 provides proofs of some results stated in the paper and Section 3 gives details on the algorithm used to perform the Expectation-Maximization based estimation.

Model and definitions

In this section, we recall the model and the assumptions given in [START_REF] Dumont | Nonparametric estimation in hidden Markov models[END_REF]. Comments on these assumptions can be found in [Dumont and Le Corff, 2012, Section 2]. Let and m be positive integers and K be a subset of R m . The main statistical problem considered in this paper is the estimation of an unknown target function f : K → R when observing a process {Y k } k∈N such that for any k ≥ 0, Y k belongs to R and satisfies

Y k def = f (X k ) + k .
{ k } k∈N is assumed to be an i.i.d Gaussian process with common distribution N (0, σ 2 I ), I being the identity matrix of size and σ 2 a fixed positive parameter. Denote by ϕ the probability distribution of 0 , i.e.

∀z ∈ R , ϕ(z) def = 2πσ 2 -/2 exp - z 2 2σ 2 ,
where • is the euclidean norm on R m . {X k } k∈N is assumed to be a non observed Markov chain, taking its values in K and independent of { k } k∈N . In

= C a (x)q x -x a , (1)
where q is a known, positive, continuous and strictly monotone function on R + and where

C a (x) def = K q x -x a dx -1 , ( 2
)
where dx is a shorthand notation for µ(dx ). The Markov transition kernel associated with q a is denoted by Q a . Assume the existence of an unknown parameter a > 0 such that H1 {X k } k∈Z is a stationary Markov chain with transition kernel Q a .

Assume the following statement on the set K:

H2 (i) K is compact. (ii) K is homeomorphic to a convex subset of R m .
(iii) K has a local Lipschitz boundary.

As an immediate consequence of the compactness of K and of the continuity of q, there exists 0 < σ -(a) < σ + (a) < +∞ such that, for all x, x ∈ K, σ -(a) ≤ q a (x, x ) ≤ σ + (a) .

(3)

For any a > 0, Q a is a ψ-irreducible and recurrent Markov kernel and then, it has a unique invariant probability distribution, see [Meyn and Tweedie, 2009, Theorem 10.0.1]. By the symmetry of the kernel (x, x ) → q x-x a , the finite measure on K with density function x → C -1 a (x) is Q a -invariant. Therefore, the unique invariant probability of Q a has a density given by

∀x ∈ K, ν a (x) def = K q x -x a dx K 2 q x -x a dx dx . ( 4
)
Let s ∈ N and p ≥ 1.

Remark 1.1. i) By H2(iii) and the Stein Theorem [Adams and Fournier, 2003, Theorem 5.24], there exists a positive constant C such that any bounded [Adams and Fournier, 2003, Theorem 6.3], assuming that K satisfies H2(i) and H2(iii) and s > m/p+k, W s,p (K, R) is compactly embedded into the subspace of bounded functions in

function f in C 1 ( • K) can be extended by a function f in C 1 (R m ), with f C 1 (R m ) ≤ C f 0, by
C k ( • K), • C k ( • K)
. Provided that s > m/p + 1, and arguing component by component, W s,p is compactly embedded into the subspace of bounded functions C 1 ( • K, R ). Moreover, the identity function id : W s,p → C 1 ( • K, R ) being linear and continuous, there exists a positive coefficient κ such that, for any f ∈ W s,p ,

f C 1 ( • K,R ) ≤ κ f W s,p , ( 5
) thus f is a bounded function in C 1 ( • K, R
) and, by i), can be extended by a function in C 1 (K, R ) shortly denoted by C 1 , and

f C 1 ≤ κ f W s,p . (6)
H3 s > m/p + 1.

For any f ∈ C 1 and any x ∈ K, the Jacobian of f at x, is defined by

J 2 f (x) def = Det D f (x) T D f (x) ,
where D f (x) is the ×m gradient matrix of f at x defined, for any j ∈ {1, . . . , } and any i ∈ {1, . . . , m}, by

D f (x) j,i def = ∂f j ∂x i (x) . H4 (i) f ∈ W s,p . (ii) f : K → Im(f ) is a diffeomorphism.
Consider the following assumption on υ.

H5 υ > 2 .

Additional proofs

This section is devoted to the proof of [Dumont and Le Corff, 2012, Propositon 5.1]. First of all, for any M ≥ 1 the Sobolev ball of radius M centred in 0 is denoted by W s,p M . Recall that, for any probability density function p on R 2 ,

g p def = 1 2 ln p + p f ,a 2p f ,a .
Define the following collections of functions on R 2 :

G M def = g p f,a ; f ∈ W s,p M , a > 0 and G M def = {g -E [g(Y 0 , Y 1 )] ; g ∈ G M } ,
where E is the expectation under the distribution P . Under H1, H2(i), H2(iii), H4(i) and H3, [Dumont and Le Corff, 2012, Propositon 5.1] states that there exist some positive constants K 1 , K 2 , C and c, depending on f and a such that, for any M ≥ 1, any n ≥ 1 and any t ≥ Cn -1/2 , P sup

g∈G M |ν n (g)| ≥ cE sup g∈G M |ν n (g)| + M t ≤ K 1 e -K2t 2 + e -K2t . (7) Let Z def = {Z k } k≥0 be the Markov chain, defined, for any k ≥ 0, by Z k def = (X 2k , Y 2k , X 2k+1 , Y 2k+1 ). For any z = (x, y) in K 2 × (R ) 2
, we denote by P z the conditional version of P where the starting distribution of the Markov chain Z is the Dirac distribution in z.

The proof of ( 7) is obtained by integration of the following result with respect to the invariant distribution of the Markov chain Z.

Proposition 2.1. Assume that H2(i), H4(i) and H3 hold. There exist some positive constants K 1 , K 2 , C and c, depending on f and a such that, for any Lemma 2.2 is proved in [Dumont and Le Corff, 2012, Lemma A.1].

M > 1, any z ∈ (K × R ) 2 , any n ≥ 1 and any t ≥ Cn -1/2 , P z sup g∈G M |ν n (g)| ≥ cE sup g∈G M |ν n (g)| + M t ≤ K 1 e -K2t 2 + e -K2t . (8)
Lemma 2.2. Assume H2(i), H2(iii) H3 and H4(i).There exists a constant C > 0 such that, for any

y ∈ R 2 , G M (y) ≤ C G (1 + M y ) and G M (y) ≤ C (1 + M y ) .
Proposition 2.1 is then an application of [Adamczak and Bednorz, 2012, Theorem 7] to the class {ḡ/M ; ḡ ∈ G M }. Indeed, Lemma 2.2 gives an upper bound for G M /M which is independent from M . This allows us to apply [Adamczak and Bednorz, 2012, Theorem 7] where all the constants in the upper bound of

P z sup g∈G M |ν n (g)| ≥ cE sup g∈G M |ν n (g)| + M t not depend on M .
By [Adamczak and Bednorz, 2012, Section 3.2], it is sufficient to prove that there exists a small set D, see [Meyn and Tweedie, 2009, Section 5.2], such that i) there exists κ > 1 satisfying sup

z∈D E z [κ τ D ] < +∞, with τ D def = min{k ≥ 1; Z k ∈ D}.
ii) The extended chain satisfies a drift condition: there exists a function V :

(K × R ) 2 → R + and b > 0 such that Q a V (z) -V (z) ≤ -exp G M (y)/M + b1 D (z) ,
where Q a is the Markov transition kernel of the extended chain Z. By [Meyn and Tweedie, 2009, Theorem 14.2.3 and Theorem 14.2.4], ii) is satisfied if

sup z∈D E z τ D -1 k=0 exp G M M (Y k ) < +∞ , (9)
where

Y k def = (Y 2k , Y 2k+1 ).
In this case, we can choose

V (z) def = E z σ D k=0 exp G M M (Y k ) ,
where σ D def = min{k ≥ 0; Z k ∈ D}. By Lemma 2.2 there exists K > 0 (independent from M ) such that the function V is upper bounded by K on D. Therefore, [Adamczak and Bednorz, 2012, Theorem 7] states the existence of constants K 1 , K 2 , c and C such that, for any t ≥ Cn -1/2 , any z ∈ (K × R ) 2 and any n > 1,

P z sup g∈G M |ν n (g)| ≥ cE sup g∈G M |ν n (g)| + M t ≤ K 1 e -K2t 2 + e -K2t √ n/ log(n) + e -K2t √ n + e -K2t ,
which concludes the proof of Proposition 2.1.

We now turn to the proof of i) and (9). By (3), it can be proved that the transition kernel Q a of the extended chain Z also satisfies a strong mixing condition. Therefore any subset of (K × R ) 2 is a small set for this extended chain. i) and ( 9) can be established by a proper choice of D. By H4(i), there exists M < +∞ such that f W s,p = M . Furthermore, by Remark 1.1, we have, for all x ∈ K, f (x) ≤ √ κM . Consider the set

D def = K × K × B(0, √ κM + ρ) × B(0, √ κM + ρ) ,
where B(0, √ κM + ρ) def = {y ∈ R ; y ≤ √ κM + ρ} and where ρ > 0 is a constant to be chosen later.

Lemma 2.3. For all k ≥ 0 and all z ∈ D,

P z {τ C > k} ≤ exp {-λ(ρ)k} , where λ(ρ) def = t 2 - 2 -ln 2 - 2 ln t 2 and t 2 def = ρ 2 σ 2 . ( 10
)
Proof. For all i ≥ 1,

P X2i Y 2i / ∈ B(0, √ κM + ρ) = P X2i f (X 2i ) + 2i ≥ √ κM + ρ ≤ P X2i { 2i ≥ ρ} . Since 2i is N (0, σ 2 I ), if t 2 > , P X2i { 2i ≥ ρ} ≤ exp 2 ln t 2 - t 2 - 2 .
This concludes the proof.

Then, for κ > 1 and z ∈ D,

E z [κ τ D ] = k≥1 P z {τ D = k} κ k ≤ e λ(ρ) k≥1 e -(λ(ρ)-ln κ)k .
The right hand side of the last equation is finite if ρ is chosen sufficiently large. This concludes the proof of i).

Proof of (9). Let z ∈ D. By Lemma 2.2, there exists a constant C such that

E z τ D -1 k=0 exp G M M (Y k ) ≤ e C E z τ D -1 k=0 exp {C Y k } ≤ e C E z τ D exp C τ D -1 k=0 Y k ≤ e C E z τ 2 D 1/2 E z exp 2C τ D -1 k=0 Y k 1/2
. By Lemma 2.3, sup z∈D E z τ 2 D < +∞. For the second term we write

E z exp 2C τ D -1 k=0 Y k = p≥1 E z 1 τ D =p exp 2C p-1 k=0 Y k ≤ p≥1 P z {τ D = p} 1/2 E z exp 4C p-1 k=0 Y k 1/2 , (11)
where

E z exp 4C p-1 k=0 Y k ≤ exp {4C y } × 2πσ 2 - exp - y 2 2σ 2 + y √ 2 κM σ 2 + 4C dy p-1
.

Let H be the Hausdorff measure on R 2 of order 2 -1 restricted to S 2 -1 , where

S 2 -1 def = {x ∈ R 2 ; x = 1}. Then, 2πσ 2 - exp - y 2 2σ 2 + y √ 2 κM σ 2 + 4C dy ≤ 2πσ 2 - R + ×S 2 -1 exp - ru 2 2σ 2 + ru √ 2 κM σ 2 + 4C H(du)r 2 -1 dr = H(S 2 -1 ) 2πσ 2 - R + exp - r 2 2σ 2 + r √ 2 κM σ 2 + 4C r 2 -1 dr ≤ H(S 2 -1 ) exp      √ 2 κM + 4Cσ 2 2 2σ 2      I 2 -1 √ 2 κM + 2 ,
where, for any c ∈ R, the sequence {I k (c)} ∞ k=1 is given by :

I k (c) = 2πσ 2 - R + exp - 1 2σ 2 [r -c] 2 r k dr .
If ξ denotes a Gaussian random variable with mean c and variance σ 2 , we have

I k (c) = (2πσ 2 ) -+1/2 E ξ k 1 ξ>0 ≤ E |ξ| k .
Then,

I k (c) ≤ (2πσ 2 ) -+1/2 E |ξ -c + c| k ≤ (2πσ 2 ) -+1/2 k i=0 k i c i E |ξ -c| k-i . If B(k) def = (2πσ 2 ) -+1/2 max 0≤i≤k E |ξ -c| k-i (which is independent from c), then I k (c) ≤ B(k) k i=0 k i c i ≤ B(k)(1 + c) k .
This yields,

E z exp 4C p-1 i=0 Y i ≤ exp {4C y } B(2 -1)H(S 2 -1 ) p-1 × exp ( √ 2 κM + 4Cσ 2 ) 2 2σ 2 p (1 + √ 2 κM + 4Cσ 2 ) (2 -1)p .
Finally, for all p ≥ 1 and all z def = (x, y) ∈ C,

E z exp 4C p-1 i=0 Y i ≤ exp {4C y } exp {η(4C)(p -1)} ,
where

η(4C) def = ln(κ( )) + √ 2 κM + 4Cσ 2 2 2σ 2 + (2 -1) ln 1 + √ 2 κM + 4Cσ 2
(12) and where κ( ) is a constant depending only on . Therefore, by ( 11) and Lemma 2.3, this concludes the proof for a sufficiently large ρ.

Numerical experiments

Let n be a positive integer, in this section, we denote by f the estimator defined as a maximizer the function T defined by

T : W s,2 → R f → 1 n n-1 k=0 log p f,a (Y 2k , Y 2k+1 ) -λ 2 n ||f || 2 W s,2 .
The HMM framework suggests to use an Expectation-Maximization (EM) type procedure, see [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. This algorithm iteratively produces a sequence of estimates { f p } p≥0 . Assume the current parameter estimate is given by f p . The estimate f p+1 is defined as one of the maximizer of the function Q defined by

f → Q(f, f p ) def = 1 n n-1 k=0 E f p log p f,a (X 2k , Y 2k , X 2k+1 , Y 2k ) Y 2k , Y 2k+1 -λ 2 n ||f || 2 W s,2
, where E f p [•] denotes the expectation under the law of the stationary HMM parameterized by f p and where

p f,a (x, y, x , y ) = ν a (x)q a (x, x )ϕ(y -f (x))ϕ(y -f (x )) . The differential of f → Q(f, f p ) is given, for any f, h ∈ W s,2 , by d f Q(•, f p )(h) = S n,1 ( f p , f, h) + S n,2 ( f p , f, h) -2λ 2 n 0≤|α|≤s D α f, D α h L2 ,
where

S n,1 ( f p , f, h) def = 1 nσ 2 n-1 k=0 E f p [ h(X 2k ), f (X 2k ) -Y 2k |Y 2k:2k+1 ] , S n,2 ( f p , f, h) def = 1 nσ 2 n-1 k=0 E f p [ h(X 2k+1 ), f (X 2k+1 ) -Y 2k+1 |Y 2k:2k+1 ] .
f p+1 is then defined as the function f ∈ W s,2 such that for any h ∈ W s,2 , d f Q( f p , •)(h) = 0. In the sequel, we choose s = 2 and K = [0, 1], therefore, this implies, for any h ∈ W ([0, 1], R),

S n,1 ( f p , f, h) + S n,2 ( f p , f, h) -2λ 2 n 2 α=0 f (α) , h (α) L2 = 0 . ( 13
)
This equation can be applied to any function h in W 2,2 0 def = {h ∈ W ([0, 1], R); h(0) = h(1) = 0}. Using integration by parts, this yields, for any component f j and any x ∈ [0, 1],

1 + 1 2nλ 2 n σ 2 n-1 k=0 φ f p ,a 2k|2k:2k+1 (x) + φ f p ,a 2k+1|2k:2k+1 (x) f j (x) -f (2) j (x)+f (4) j (x) = 1 2nλ 2 n σ 2 n-1 k=0 Y 2k φ f p ,a 2k|2k:2k+1 (x) + Y 2k+1 φ f p ,a 2k+1|2k:2k+1 (x) , (14)
where φ f p ,a 2k|2k:2k+1 and φ f p ,a 2k+1|2k:2k+1 are the filtering distributions defined by

φ f p ,a 2k|2k:2k+1 (x) def = ν a (x)q a (x, x )ϕ(Y 2k -f p (x))ϕ(Y 2k+1 -f p (x ))dx p f p ,a (Y 2k , Y 2k+1 ) , φ f p ,a 2k+1|2k:2k+1 (x ) def = ν a (x)q a (x, x )ϕ(Y 2k -f p (x))ϕ(Y 2k+1 -f p (x ))dx p f p ,a (Y 2k , Y 2k+1 ) .

Numerical approximations

Let N ≥ 1 be an integer. The differential system ( 14) is solved using a discretization of the state space [0, 1] by i N N i=0 . Let q a be the transition probability associated to this discretization: for any i

, j ∈ {0, • • • , N }, q a (i, j) def = q a ( i N , j N ) N j=0 q a ( i N , j N )
and ν a the invariant distribution of q a on {0, . . . , N }. Define ν a as the invariant distribution of q a on {0, . . . , N }. The filtering distributions φ , defined by

φ f p ,a k (x) def = N -1 i=0 1 [ i N , i+1 N [ (x) ϕ f p i,k and φ f p ,a k (x) def = N -1 i=0 1 [ i N , i+1 N [ (x) ϕ f t i,k ,
where, for any i, j ∈ {0, . . . , N -1},

ϕ f p i,k def = N -1 j=0 ν a (i) q a (i, j)ϕ f p i N -Y 2k ϕ f p j N -Y 2k+1 N -1 i ,j =0 ν a (i ) q a (i , j )ϕ f p i N -Y 2k ϕ f p j N -Y 2k+1 , (15)
ϕ f p j,k def = N -1 i=0 ν a (i) q a (i, j)ϕ f p i N -Y 2k ϕ f p j N -Y 2k+1 N -1 i ,j =0 ν a (i ) q a (i , j )ϕ f p i N -Y 2k ϕ f p j N -Y 2k+1 . (16)
The equation ( 14) is solved on each interval

[ i N , i+1 N [, i ∈ {0, • • • , N -1}
, which is straightforward since the coefficients are constant and the equation is linear.

Computation of f p+1

Let i ∈ {0, • • • , N -1} and j ∈ {1, • • • , }. We denote by f i the solution of ( 14) on i N , i+1 N . We have, for each component f j , j ∈ {1, • • • , },

f (4) j,i (x) -f (2) j,i (x) + (1 + α i ) f j,i (x) = β j,i ,
where

α i def = 1 2nλ 2 n σ 2 n-1 k=0 ϕ f n i,k + ϕ f n i,k , β j,i def = 1 2nλ 2 n σ 2 n-1 k=0 Y 2k,j ϕ f n i,k + Y 2k+1,j ϕ f n j,k
.

Therefore, there exist c i,1 , c i,2 , s i,1 and s i,2 such that, for any x ∈ i N , i+1 N , f j,i (x) = e ηix [c 1,i cos (γ i x) + s 1,i sin (γ i x)]

+ e -ηix [c 2,i cos (γ i x) + s 2,i sin (γ i x)] + β j,i 1 + α i , where, if

r i def = √ 1 + α i η i def = √ 1 + 2r i 2 and γ i def = √ 2r i -1 2 .
Therefore, 4N parameters have to be chosen to uniquely determine the solution

f p+1 j = N -1 i=0 1 [ i N , i+1
N [ f j,i . The C 3 -regularity conditions for each boundary provides 4(N -1) equations: for any i ∈ {0, • • • , N -2},

f j,i i + 1 N = f j,i+1 i + 1 N , f j,i i + 1 N = f j,i+1 i + 1 N , f (2) j,i i + 1 N = f (2) j,i+1 i + 1 N , f (3) j,i i + 1 N = f (3) j,i+1 i + 1 N , (17)

  Proposition 2.1 is an application of[Adamczak and Bednorz, 2012, Theorem 7] and relies on an intermediate lemma on the envelope functions of the sets G M and G M defined, for any y ∈ R 2 , by G M (y) def = sup g∈G M g(y) and G M (y)

where, for any x ∈ i N , i+1 N ,

Solving (13) with h(x) = 1, h(x) = x, h(x) = x 2 and h(x) = x 3 leads to four other linear equations which conclude the computation of f p+1 j .