Olga Klopp

Sparse high-dimensional varying coefficient model: non-asymptotic minimax study

Keywords: varying coefficient model, sparse model, minimax optimality AMS 2000 subject classification: 62H12, 62J05, 62C20

The objective of the present paper is to develop a minimax theory for the varying coefficient model in a non-asymptotic setting. We consider a high-dimensional sparse varying coefficient model where only few of the covariates are present and only some of those covariates are time dependent. Our analysis allows the time dependent covariates to have different degrees of smoothness and to be spatially inhomogeneous. We develop the minimax lower bounds for the quadratic risk and construct an adaptive estimator which attains those lower bounds within a constant (if all time-dependent covariates are spatially homogeneous) or logarithmic factor of the number of observations.

Introduction

One of the fundamental tasks in statistics is to characterize the relationship between a set of covariates and a response variable. In the present paper we study the varying coefficient model which is commonly used for describing time-varying covariate effects. It provides a more flexible approach than the classical linear regression model and is often used to analyze the data measured repeatedly over time.

Since its introduction by Cleveland, Grosse and Shyu [START_REF] Cleveland | Local regression models[END_REF] and Hastie and Tibshirani [START_REF] Hastie | Varying-coefficient models[END_REF] many methods for estimation and inference in the varying coefficient model have been developed (see, e.g., [START_REF] Wu | Asymptotic confidence regions for kernel smoothing of a varying-coefficient model with longitudinal data[END_REF][START_REF] Hoover | Non-parametric smoothing estimates of time-varying coefficient models with longitudinal data[END_REF][START_REF] Fan | Statistical estimation in varying coefficient models[END_REF][START_REF] Kauermann | On model diagnostics using varying coefficient models[END_REF] for the kernel-local polynomial smoothing approach, [START_REF] Huang | Varying-coefficient models and basis function approximations for the analysis of repeated measurements[END_REF][START_REF] Huang | Polynomial spline estimation and inference for varying coefficient models with longitudinal data[END_REF][START_REF] Huang | Functional coefficient regression models for nonlinear time series: A polynomial spline approach[END_REF] for the polynomial spline approach, [START_REF] Hastie | Varying-coefficient models[END_REF][START_REF] Hoover | Non-parametric smoothing estimates of time-varying coefficient models with longitudinal data[END_REF][START_REF] Chiang | Smoothing spline estimation for varying coefficient models with repeatedly measured dependent variables[END_REF] for the smoothing spline approach and [START_REF] Fan | Statistical methods with varying coefficient models[END_REF] for a detailed discussion of the existing methods and possible applications). In the last five years, the varying coefficient model received a great deal of attention. For example, Wang et al. [START_REF] Wang | Local Rank Inference for Varying Coefficient Models[END_REF] proposed a new procedure based on a local rank estimator; Kai et al. [START_REF] Kai | New efficient estimation and variable selection methods for semiparametric varying-coefficient partially linear models[END_REF] introduced a semi-parametric quantile regression procedure and studied an effective variable selection procedure. Lee et al. [START_REF] Lee | Flexible generalized varying coefficient regression models[END_REF] extended the model to the case when the response is related to the covariates via a link function while Zhu et al. [START_REF] Zhu | Multivariate varying coefficient model for functional responses[END_REF] studied the multivariate version of the model. Existing methods typically provide asymptotic evaluation of the precision of the estimation procedure under the assumption that the number of observations tends to infinity and is larger than the dimension of the problem.

Recently few authors consider still asymptotic but high-dimensional approach to the problem. Wei et al. [START_REF] Wei | Variable Selection and Estimation in High-Dimensional Varying-Coefficient Models[END_REF] applied group Lasso for variable selection, while Lian [START_REF] Lian | Variable selection for high-dimensional generalized varying-coefficient models[END_REF] used extended Bayesian information criterion. Fan et al. [START_REF] Fan | Nonparametric Independence Screening in Sparse Ultra-High Dimensional Varying Coefficient Models[END_REF] applied nonparametric independence screening. Their results were extended by Lian and Ma [START_REF] Lian | Reduced-rank Regression in Sparse Multivariate Varying-Coefficient Models with High-dimensional Covariates[END_REF] to include rank selection in addition to variable selection.

One important aspect that has not been well studied in the existing literature is the non-asymptotic approach to the estimation, prediction and variables selection in the varying coefficient model. Here, we refer to the situation where both the number of unknown parameters and the number of observations are large and the former may be of much higher dimension than latter. The only reference that we are aware of in this setting, is the recent paper by Klopp and Pensky [START_REF] Klopp | Non-asymptotic approach to varying coefficient model[END_REF]. Their method is based on some recent developments in the matrix estimation problem. Some interesting questions arise in this non-asymptotic setting. One of them is the fundamental question of the minimax optimal rates of convergence. The minimax risk characterizes the essential statistical difficulty of the problem. It also captures the interplay between different parameters in the model. To the best of our knowledge, our paper presents the first non-asymptotic minimax study of the sparse heterogeneous varying coefficient model.

Modern technologies produce very high dimensional data sets and, hence, stimulate an enormous interest in variable selection and estimation under a sparse scenario. In such scenarios, penalization-based methods are particularly attractive. Significant progress has been made in understanding the statistical properties of these methods. For example, many authors have studied the variable selection, estimation and prediction properties of the LASSO in high-dimensional setting (see, e.g., [START_REF] Bickel | Simultaneous analysis of Lasso and Dantzig selector[END_REF], [START_REF] Bunea | Aggregation for Gaussian regression[END_REF], [START_REF] Bunea | Sparsity oracle inequalities for the Lasso[END_REF], [START_REF] Van De Geer | High-dimensional generalized linear models and the Lasso[END_REF]). A related LASSO-type procedure is the group-LASSO, where the covariates are assumed to be clustered in groups (see, for example, [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF][START_REF] Bach | Consistency of the group lasso and multiple kernel learning[END_REF][START_REF] Chesneau | Some theoretical results on the grouped variables Lasso[END_REF][START_REF] Meier | The group Lasso for logistic regression[END_REF][START_REF] Meier | High-dimensional additive modeling[END_REF][START_REF] Lounici | Oracle inequalities and optimal inference under group sparsity[END_REF], and references therein).

In the present paper, we also consider the case when the solution is sparse, in particular, only few of the covariates are present and only some of them are time dependent. This setup is close to the one studied in a recent paper of Liang [START_REF] Lian | Variable selection for high-dimensional generalized varying-coefficient models[END_REF]. One important difference, however, is that in [START_REF] Lian | Variable selection for high-dimensional generalized varying-coefficient models[END_REF], the estimator is not adaptive to the smoothness of the time dependent covariates. In addition, Liang [START_REF] Lian | Variable selection for high-dimensional generalized varying-coefficient models[END_REF] assumes that all time dependent covariates have the same degree of smoothness and are spatially homogeneous. On the contrary, we consider a much more flexible and realistic scenario where the time dependent covariates possibly have different degrees of smoothness and may be spatially inhomogeneous.

In order to construct a minimax optimal estimator, we introduce the block Lasso which can be viewed as a version of the group LASSO. However, unlike in group LASSO, where the groups occur naturally, the blocks in block LASSO are driven by the need to reduce the variance as it is done, for example, in block thresholding. Note that our estimator does not require the knowledge which of the covariates are indeed present and which are time dependent. It adapts to sparsity, to heterogeneity of the time dependent covariates and to their possibly spatial inhomogeneous nature. In order to ensure the optimality, we derive minimax lower bounds for the risk and show that our estimator attains those bounds within a constant (if all time-dependent covariates are spatially homogeneous) or logarithmic factor of the number of observations. The analysis is carried out under the flexible assumption that the noise variables are sub-Gaussian. In addition, it does not require that the elements of the dictionary are uniformly bounded.

The rest of the paper is organized as follows. Section 1.1 provides formulation of the problem while Section 1.2 lays down a tensor approach to estimation. Section 2 introduces notations and assumptions on the model and provides a discussion of the assumptions. Section 3 describes the block thresholding LASSO estimator, evaluates the non-asymptotic lower and upper bounds for the risk, both in probability and in the mean squared risk sense, and ensures optimality of the constructed estimator. Section 4 presents examples of estimation when assumptions of the apper are satisfied. Section 5 contains proofs of the statements formulated in the paper.

Formulation of the problem

Let (W i , t i , Y i ), i = 1, . . . , n be sampled independently from the varying coefficient model

Y = W T f (t) + σξ. (1.1)
Here, W ∈ R p are random vectors of predictors,

f (•) = (f 1 (•), . . . , f p (•))
T is an unknown vectorvalued function of regression coefficients and t ∈ [0, 1] is a random variable with the unknown density function g. We assume that W and t are independent. The noise variable ξ is independent of W and t and is such that E(ξ) = 0 and E(ξ 2 ) = 1, σ > 0 denotes the known noise level.

The goal is to estimate vector function f (•) on the basis of observations (W i , t i , Y i ), i = 1, . . . , n.

In order to estimate f , following Klopp and Pensky (2013), we expand it over a basis (φ l (•)), l = 0, 1, . . . , ∞, in L 2 ([0, 1]) with φ 0 (t) = 1. Expansion of the functions f j (•) over the basis, for any t ∈ [0, 1], yields

f j (t) = L l=0 a jl φ l (t) + ρ j (t) with ρ j (t) = ∞ l=L+1 a jl φ l (t).
(

If φ(•) = (φ 0 (•), . . . , φ L (•)) T and A 0 denotes a matrix of coefficients with elements A (l,j) 0 = a jl , then relation (1.2) can be re-written as f (t) = A T 0 φ(t) + ρ(t), where ρ(t) = (ρ 1 (t), • • • , ρ p (t)) T . Combining formulae (1.1) and (1.2), we obtain the following model for observations (W i , t i , Y i ), i = 1, . . . , n:

Y i = Tr(A T 0 φ(t i )W T i ) + W T i ρ(t i ) + σξ i , i = 1, . . . , n. (1.3)
Below, we reduce the problem of estimating vector function f to estimating matrix A 0 of coefficients of f .

Tensor approach to estimation

Denote a = Vec(A 0 ) and

B i = Vec(φ(t i )W T i ). Note that B i is the p(L + 1)-dimensional vector with components φ l (t i )W (j) i , l = 0, • • • , L, j = 1, • • • , p, where W (j) i is the j-th component of vector W i . Consider matrix B ∈ R n×p(L+1) with rows B T i , i = 1, • • • , n, vector ξ = (ξ 1 , • • • , ξ n ) T and vector b with components b i = W T i ρ(t i ), i = 1, • • • , n. Taking into account that Tr(A T φ(t i )W T i ) = B T i Vec(A)
we rewrite the varying coefficient model (1.3) in a matrix form

Y = Ba + b + σξ. (1.4)
In what follows, we denote

Ω i = W i W T i , Φ i = φ(t i )(φ(t i )) T , Σ i = Ω i ⊗ Φ i , (1.5)
where Ω i ⊗ Φ i is the Kronecker product of Ω i and Φ i . Note that Ω i , Φ i and Σ i are i.i.d. for i = 1, • • • , n, and that Ω i 1 and Φ i 2 are independent for any i 1 and i 2 . By simple calculations, we derive

a T BB T a = n i=1 (B T i a) 2 = n i=1 Tr(A T φ(t i )W T i ) 2 = n i=1 W T i A T φ(t i )φ T (t i )AW i = n i=1 a T (Ω i ⊗ Φ i )a,
which implies

B T B = n i=1 Ω i ⊗ Φ i . (1.6) Let Σ = n -1 B T B = n -1 n i=1 Σ i . (1.7)
Then, due to the i.i.d. structure of the observations, one has

Σ = EΣ 1 = Ω ⊗ Φ with Ω = E(W 1 W T 1 ) and Φ = E(φ(t 1 )φ T (t 1 )). (1.8)
2 Assumptions and notations

Notations

In what follows, we use bold script for matrices and vectors, e.g., A or a, and superscripts to denote elements of those matrices and vectors, e.g., A (i,j) or a (j) . Below, we provide a brief summary of the notation used throughout this paper.

• For any vector a ∈ R p , denote the standard l 1 and l 2 vector norms by a 1 and a 2 , respectively. For vectors a, c ∈ R p , denote their scalar product by a, c .

• For any function q(t), t ∈ [0, 1], q 2 and • , • 2 are, respectively, the norm and the scalar product in the space L 2 ([0, 1]). Also, q ∞ = sup t∈[0,1] |q(t)|.

• For any vector function q(t) = (q 1 (t), • • • , q p (t)) T , denote

q(t) 2 =   p j=1 q j 2 2   1/2
.

• For any matrix A, denote its spectral and Frobenius norms by A and A 2 , respectively.

• Denote the k × k identity matrix by I k .

• For any numbers, a and b, denote a ∨ b = max(a, b) and a ∧ b = min(a, b).

• In what follows, we use the symbol C for a generic positive constant, which is independent of n, p, s and l, and may take different values at different places.

• If r = (r 1 , • • • , r p ) T and r j = r j + 1/2 -1/ν j for some 1 ≤ ν j < ∞, denote r * j = r j ∧ r j and r * min = min j r * j .

• Denote t = (t 1 , • • • , t n ), W = (W 1 , • • • , W n ), (2.1)
i.e., W is the p × n matrix with columns

W i , i = 1, • • • , n.

Assumptions

We impose the following assumptions on the varying coefficient model (1.3).

(A0). Only s out of p functions f j are non-constant and depend on the time variable t, s 0 functions are constant and independent of t and (p -s -s 0 ) functions are identically equal to zero. We denote by J the set of indices corresponding to the non-constant functions f j .

(A1).

Functions (φ k (•)) k=0,...,∞ form an orthonormal basis of L 2 ([0, 1]), and are such that φ 0 (t) = 1 and, for any t ∈ [0, 1], any l ≥ 0 and some

C φ < ∞ l k=0 φ 2 k (t) ≤ C 2 φ (l + 1). (2.2)
(A2). The probability density function g(t) is bounded above and below 0 < g 1 ≤ g(t) ≤ g 2 < ∞. Moreover, the eigenvalues of E(φφ T ) = Φ are bounded from above and below

0 < φ min = λ min (Φ) ≤ λ max (Φ) = φ max < ∞.
Here, φ min and φ max are absolute constants independent of L.

(A3).

Functions f j (t) have efficient representations in basis φ l , in particular, for any j

= 1, • • • , p, one has ∞ k=0 |a jk | ν j (k + 1) ν j r j ≤ (C a ) ν j , r j = r j + 1/2 -1/ν j , (2.3)
for some C a > 0, 1 ≤ ν j < ∞ and r j > min(1/2, 1/ν j ). In particular, if function f j (t) is constant or vanishes, then r j = ∞. We denote vectors with elements ν j and r j , j = 1, • • • , p, by ν and r, respectively, and the set of indices of finite elements r j by J :

J = {j : r j < ∞}. (2.4) (A4). Variables ξ i , i = 1, • • • , n, are i.i.d. sub-Gaussian such that Eξ i = 0, Eξ 2 i = 1 and ξ i ψ 2 ≤ K
where • ψ 2 denotes the sub-Gaussian norm.

(A5) "Restricted Isometry in Expectation" condition. Let W Λ , Λ ∈ {1, . . . , p} be the subvector obtained by extracting the entries of W corresponding to indices in Λ and let

Ω Λ = E W Λ W T Λ .
We assume that there exist two positive constants ω max (ℵ) and ω min (ℵ) such that for all subsets Λ with cardinality |Λ| ≤ ℵ and all v ∈ R ℵ one has

ω min (ℵ) v 2 2 ≤ Ω Λ v 2 2 ≤ ω max (ℵ) v 2 2 .
(2.5)

Moreover, we suppose that E(W (j) ) 4 ≤ V for any j = 1, • • • , p, and that, for any µ ≥ 1 and for all subsets Λ with |Λ| ≤ ℵ, there exist positive constants U µ and C µ and a set W µ such that

W Λ ∈ W µ =⇒ ( W Λ 2 ≤ U µ ) ∩ max j=1,••• ,p |W (j) Λ | ≤ C µ , P (W Λ ∈ W µ ) ≥ 1 -2 p -2µ . (2.6)
Here,

U µ = U µ (ℵ), C µ = C µ (ℵ).

(A6).

We assume that (s + s 0 )(1 + log n) ≤ p and that there exists a numerical constant

C ω > 1 such that log(n) ≥ C ω φ max ω max ((s + s 0 ) log n) φ min ω min ((s + s 0 )(1 + log n)) .
We denote

ω * max = ω max ((s + s 0 )(1 + log n)) , ω * min = ω min ((s + s 0 )(1 + log n)) .
(2.7)

Discussion of assumptions

• Assumptions (A0) corresponds to the case when s of the covariates f j (t) are indeed functions of time, s 0 of them are time independent and (p -s -s 0 ) are irrelevant.

• Assumption (A1) deals with the basis of L 2 ([0, 1]). There are many types of orthonormal bases satisfying those conditions. a) Fourier basis.

Choose φ 0 (t) = 1, φ k (t) = 2 sin(2πkt) if k > 1 is odd, φ k (t) = 2 cos(2πkt) if k > 1 is even. The basis functions are bounded and C φ = 2. b) Wavelet basis. Consider a periodic wavelet basis on [0, 1]: ψ h,i (t) = 2 h/2 ψ(2 h t -i) with h = 0, 1, • • • , i = 0, • • • , 2 h -1. Set φ 0 (t) = 1 and φ j (t) = ψ h,i (t) with j = 2 h + i + 1. If l = 2 J , then the condition (2.2) is satisfied with C φ = ψ ∞ . Observing that, for 2 J < l < 2 J+1 , we have (l + 1) ≥ (2 J+1 + 1)/2 one can take C φ = 2 ψ ∞ .
• Assumption (A2) that φ min and φ max are absolute constants independent of L is guaranteed by the fact that the sampling density g is bounded above and below. For example, if g(t) = 1, one has φ min = φ max = 1.

• Assumption (A3) describes sparsity of the vectors of coefficients of functions f j (t) in basis φ l , j = 1, • • • , p and its smoothness. For example, if ν j < 2, the vector of coefficients a jl of f j is sparse. In the case when basis φ l is formed by wavelets, condition (2.3) implies that f j belongs to a Besov ball of radius C a . If we chose Fourier bases and ν j = 2, then f j belongs to a Sobolev ball of smoothness r j and radius C a . Note that Assumption (A3) allows each non-constant function f j to have its own sparsity and smoothness patterns.

• Assumption (A4) that ξ i are sub-Gaussian random variables means that their distribution is dominated by the distribution of a centered Gaussian random variable. This is a convenient and reasonably wide class. Classical examples of sub-gaussian random variables are Gaussian, Bernoulli and all bounded random variables. Note that Eξ 2 i = 1 implies that K ≤ 1.

• Assumption (A5) is closely related to the Restricted Isometry (RI) conditions usually considered in the papers that employ LASSO technique or its versions, see, e.g., [START_REF] Bickel | Simultaneous analysis of Lasso and Dantzig selector[END_REF]. However, usually the RI condition is imposed on the matrix of scalar products of the elements of a deterministic dictionary while we deal with a random dictionary and require this condition to hold only for the expectation of this matrix.

Note that the upper bound in the condition (2.5) is automatically satisfied with

ω max = Ω 2
where Ω is the spectral norm of the matrix E(WW T ) = Ω. If the smallest eigenvalue of Ω, λ min (Ω), is non-zero, then the lower bound in (2.5) is satisfied with ω min = λ 2 min (Ω). However, since the λ-restricted maximal eigenvalue ω max (λ) may be much smaller than the spectral norm of Ω and ω min (λ) may be much larger then λ min (Ω), using those values will result in sharper bounds for the error. Note that in the case when W has i.i.d. zero-mean entries W j with E W (j) 2 = ν 2 , we have ω max = ω min = ν 2 .

• Assumption (A6) is usual in the literature on the hight-dimensional linear regression model, see, e.g., [START_REF] Bickel | Simultaneous analysis of Lasso and Dantzig selector[END_REF]. For instance, if W has i.i.d. zero-mean entries W j and g(t) = 1, this condition is satisfied for any 1 < C ω ≤ log(n).

3 Estimation strategy and non-asymptotic error bounds

Estimation strategy

Formulation (1.4) implies that the varying coefficient model can be reduced to the linear regression model and one can apply one of the multitude of penalized optimization techniques which have been developed for the linear regression. In what follows, we apply a block LASSO penalties for the coefficients in order to account for both the constant and the vanishing functions f j and also to take advantage of the sparsity of the functional coefficients in the chosen basis.

In particular, for each function f j , j = 1, • • • , p, we divide its coefficients into M + 1 different groups where group zero contains only coefficient a j0 for the constant function φ 0 (t) = 1 and M groups of size d ≈ log n where M = L/d. We denote a j0 = a j0 and a jl = (a j,d(l-1)+1 , • • • , a j,dl ) T the sub-vector of coefficients of function f j in block l, l = 1, • • • , M . Let K l be the subset of indices associated with a jl . We impose block norm on matrix A as follows

A block = p j=1 M l=0 a jl 2 . (3.1)
Observe that A block indeed satisfies the definition of a norm and is a sum of absolute values of coefficients a j0 of functions f j and l 2 norms for each of the block vectors of coefficients a jl ,

j = 1, • • • , p, l = 1, • • • , M .
The penalty which we impose is related to both the ordinary and the group LASSO penalties which have been used by many authors. The difference, however, lies in the fact that the structure of the blocks is not motivated by naturally occurring groups (like, e.g., rows of the matrix A) but rather our desire to exploit sparsity of functional coefficients a jl . In particular, we construct an estimator A of A 0 as a solution of the following convex optimization problem

A = arg min A n -1 n i=1 Y i -Tr(A T φ(t i )W T i ) 2 + δ A block , (3.2)
where the value of δ will be defined later.

Note that with the tensor approach which we used in Section 1.2, optimization problem (3.2) can be re-written in terms of vector α = Vec(A) as

â = arg min α n -1 Y -Bα 2 2 + δ α block , (3.3)
where α block = A block is defined by the right-hand side (3.1) with vectors a jl being subvectors of vector α. Subsequently, we construct an estimator f (t) = ( f1 (t),

• • • , fp (t)) T of the vector function f (t) using fj (t) = L k=0 âjk φ k (t), j = 1, • • • , p. (3.4)
In what follows, we derive the upper bounds for the risk of the estimator â (or A) and suggest a value of parameter δ which allows to attain those bounds. However, in order to obtain a benchmark of how well the procedure is performing, we determine the lower bounds for the risk of any estimator A under assumptions (A0)-(A4). Remark 1 Assumption that K = L/d is an integer is not essential. Indeed, we can replace the number of groups K by the largest integer below or equal to L/d and then adjust group sizes to be d or d + 1 where d = [log n], the largest integer not exceeding log n.

Lower bounds for the risk

In this section we will obtain the lower bounds on the estimation risk. We consider a class F = F s 0 ,s,ν,r (C a ) of vector functions f (t) such that s of their components are non-constant with coefficients satisfying condition (2.3) in (A3), s 0 of the components are constant and (p -s -s 0 ) components are identically equal to zero. We construct the lower bound for the minimax quadratic risk of any estimator f of the vector function f ∈ F s 0 ,s,ν,r (C a ). Let ω * max be given by formula (2.7). Denote r max = max{r j : j ∈ J } and

n low = 2 σ 2 κ C 2 a ω * max φ max max 1, 6 s 2rmax+1 , (3.5)
∆ lower (s 0 , s, n, r) = max

   κ σ 2 s 0 4n ω * max φ max , 1 8 j∈J C 2 2r j +1 a σ 2 κ n ω * max φ max 2r j 2r j +1    . (3.6)
Then, the following statement holds.

Theorem 1 Let s ≥ 1 and s 0 ≥ 3. Consider observations Y i in model (1.3) with W i , i = 1, • • •
, n and t satisfying assumptions (A5) and (A2), respectively. Assume that, conditionally on W i and t i , variables ξ i are Gaussian N (0, 1) and that n ≥ n low . Then, for any κ < 1/8 and any estimator f of f , one has

inf f sup f ∈F P f -f 2 2 ≥ ∆ lower (s 0 , s, n, r) ≥ √ 2 1 + √ 2 1 -2κ - 2κ log 2 . ( 3.7)
Note that condition s 0 ≥ 3 is not essential since, for s 0 < 3, the first term in (3.6) is of parametric order. Condition n ≥ n low is a purely technical condition which is satisfied for the collection of n's for which upper bounds are derived. Observe also that inequality (3.7) immediately implies that

inf f sup f ∈F E f -f 2 2 ≥ ∆ lower (s 0 , s, n, r) √ 2 1 + √ 2 1 -2κ - 2κ log 2
.

(3.8)

Adaptive estimation and upper bounds for the risk

In this section we derive an upper bound for the risk of the estimator (3.2). For this purpose, first, we shall show that, with high probability, the ratio between the restricted eigenvalues of matrices Σ defined in (1.7) and Σ = E Σ is bounded above and below. This is accomplished by the following lemma.

For any Λ ∈ {1, . . . , p}, we denote by

(Ω i ) Λ = (W i ) Λ ((W i ) Λ ) T , Σ Λ = n -1 n i=1 (Ω i ) Λ ⊗Φ i and Σ Λ = E Σ Λ .
For some 0 < h < 1 and 1 ≤ ℵ ≤ p, we define

N (ℵ) = 64 µ ℵ (L + 1) log(p + L) C 2 φ U 2 µ (ℵ) φ max ω max (ℵ) h 2 φ 2 min ω 2 min (ℵ) . (3.9)
Lemma 1 Let n ≥ N (ℵ) and µ in (2.6) be large enough, so that

p µ ≥ max √ 2V n 8µ √ ℵ U 2 µ (ℵ) log(p + L) , 2n , (3.10)
where V is defined in Assumption (A5). Then, for any Λ ∈ {1, . . . , p}

inf Λ: |Λ|≤ℵ P Σ Λ -Σ Λ < h φ min ω min (ℵ) ∩ W ⊗n µ ≥ 1 -2p -µ , (3.11)
where W µ is the set of points in R p such that condition (2.6) holds and W ⊗n µ is the direct product of n sets W µ .

Moreover, on the set W ⊗n µ , with probability at least 1 -2p -µ , one has simultaneously

inf Λ: |Λ|≤ℵ λ min ( Σ Λ ) ≥ (1-h)φ min ω min (ℵ), sup Λ: |Λ|≤ℵ λ max ( Σ Λ ) ≤ (1+h)φ max ω max (ℵ). (3.12)
Lemma 1 ensures that the restricted lowest eigenvalue of the regression matrix Σ is within a constant factor of the respective eigenvalue of matrix Σ. Since p may be large, this is not guaranteed by a large value of n (as it happens in the asymptotic setup) and leads to additional conditions on the relationship between parameters L, p and n.

By applying a combination of LASSO and group LASSO arguments, we obtain the following theorem that gives an upper bound for the quadratic risk of the estimator (3.2). We set

U µ = U µ (s + s 0 ). Define N = max 64 µ (s + s 0 )C 2 φ U 2 µ (L + 1)φ max ω * max log(p + L) h 2 φ 2 min (ω * min ) 2 , U 2 µ C 2 φ (L + 1)µ log p g 2 ω max (s) , 3 C 2 a g 2 s ω max (s) (3.13) and δ = 2 (σC ω K √ µ + 1)) (1 + h)φ max ω max (1) log p n . (3.14) Theorem 2 Let min k (r k ∧ r k ) ≥ 2, L + 1 ≥ n 1/2
and n ≥ N. Let µ in (2.6) be large enough, so that

p µ ≥ max √ 2V n 8µ √ s + s 0 U 2 µ log(p + L) , 2L log n , 2n . (3.15)
If â is an estimator of a obtained as a solution of optimization problem (3.3) with δ = δ, and the vector function f is recovered using (3.4), then, one has

P f -f 2 2 ≤ ∆(s 0 , s, n, r) ≥ 1 -8 p -µ (3.16)
where

∆(s 0 , s, n, r) = C 2 a s n 2 + C B (1 + h)ω * max φ max (1 -h)ω * min φ min C ω σ 2 K 2 µ + 1 (s 0 + s) log p n(1 -h)ω * min φ min + j∈J C 2/(2r j +1) a C ω σ 2 K 2 µ + 1 n(1 -h)ω * min φ min 2r j 2r j +1 (log n) (2-ν j ) + -2ν j r j ν j (2r j +1) (log p) 2r j 2r j +1   .
Note that construction (3.3) of the estimator â does not involve knowledge of unknown parameters r and ν or matrix Σ, therefore, estimator â is fully adaptive. Moreover, conclusions of Theorem 2 are derived without any asymptotic assumptions on n, p and L.

In order to assess the optimality of estimator â, we consider the case of the Gaussian noise, i.e. K = 1. Observe that, under Assumption (A2), the values of φ min and φ max are independent of n and p, so that the only quantities in (3.16) where C is an absolute constant independent of n, p and σ 2 and we use (2 -ν j ) + -2ν j r j ≤ 0. Inequality (3.17) implies that, for any values of the parameters, the ratio between the upper and the lower bound for the risk (3.6) is bounded by C log(p) ω * max 2 /ω * min 2 . Note that ω * max /ω * min is the condition number of matrix Ω Λ with |Λ| = (s + s 0 )(1 + log n). Hence, if matrix Ω Λ is well conditioned, so that ω * max /ω * min is bounded by a constant, the estimator f attains optimal convergence rates up to a log p factor.

Consider the case when all functions f j (t) are spatially homogeneous, i.e., min j ν j ≥ 2 and n is large enough, i. e. there exists a positive β such that n β ≥ p. Then, the estimator f attains optimal convergence rates up to a constant factor, if s/s 0 is bounded or nω * min /σ 2 is relatively large. In particular, if all functions in assumption (A3) belong to the same space, then the following corollary is valid.

Corollary 1 Let conditions of Theorem 2 hold with r j = r and ν j = ν, j -1, • • • , p, and matrix Ω be well conditioned, i.e. ω * max /ω * min is bounded by for some absolute constant independent of n, p and σ. Then,

∆ upper (s 0 , s, n, r) ∆ lower (s 0 , s, n, r) ≤    C log p, if σ 2 (s/s 0 ) 2r+1 ≥ nω * min , C(log p) 2r 2r+1 , if σ 2 (s/s 0 ) 2r+1 < nω * min , 1 ≤ ν < 2, C, if σ 2 (s/s 0 ) 2r+1 < nω * min , ν ≥ 2.
(3.18)

Adaptive estimation with respect to the mean squared risk

Theorem 2 derives upper bounds for the risk with high probability. Suppose that an upper bound on the norms of functions f j is available due to physical or other considerations:

max 1≤j≤p f j 2 2 ≤ C 2 f . (3.19)
Then, a 2 ≤ p C 2 f and â given by (3.3) can be replaced by the solution of the convex problem â = arg min

a n -1 Y -Ba 2 2 + δ a block s.t. a 2 ≤ p C 2 f (3.20)
with δ = δ where δ is defined in (3.14), and estimators fj of f j , j = 1, • • • , p, are constructed using formula (3.4). Choose µ in (2.6) large enough, so that

16nC 2 f ≤ p µ-1 . (3.21)
Then, the following statement is valid.

Theorem 3 Under the assumptions of the Theorem 2, and for any µ satisfying condition (3.21), one has

E f -f 2 2 ≤ C∆ upper (s 0 , s, n, r) (3.22)
where C is an absolute constant independent of n, p and σ.

Examples and discussion

Examples

In this section we provide several examples when assumptions of the paper are satisfied. For simplicity, we assume that g(t) = 1, so that φ min = φ max = 1.

Example 1 Normally distributed dictionary Let W i , i = 1, • • • , n, be i.i.d. standard Gaussian vectors N (0, I p ). Then, Ω = I p , so that ω min = ω max = 1. Moreover, W (j) are independent standard Gaussian variables and (W i ) T Λ (W i ) Λ are independent chi-squared variables with |Λ| = ℵ degrees of freedom. Using inequality (see, e.g., [START_REF] Birgè | Minimal Penalties for Gaussian Model Selection[END_REF], page 67)

P χ 2 ℵ ≤ ℵ + 2 √ ℵx + 2x ≥ 1 -e -x , x > 0,
for any µ 1 ≥ 0, derive

P (W 1 ) T Λ (W 1 ) Λ ≤ ( √ ℵ + √ 2µ 1 ) 2 ≥ 1 -exp(-µ 2 1 ).
Choose any µ > 0 and set µ 2 1 = 2µ log(p). Then, using a standard bound on the maximum of p Gaussian variables one obtains that Assumption (A5) holds with

C µ = 2 log p, U 2 µ = ( √ ℵ + 2 µ log p) 2 .
Example 2 Symmetric Bernoulli dictionary Let W (j)

i , i = 1, • • • , n, j = 1, • • • , p, be independent symmetric Bernoulli variables P(W (j) i = 1) = P(W (j) i = -1) = 1/2.
Then, Ω = I p , ω min = ω max = 1 and, for any µ,

C µ = 1, U 2 µ (ℵ) = ℵ.
In both cases, N in (3.13) is of the form N = C(s + s 0 ) 2 (L + 1) log(p). Under conditions of Theorem 2, the upper bounds for the risk are of the form

∆ upper (s 0 , s, n, r) = C   σ 2 s 0 log n n + j∈J σ 2 n 2r j 2r j +1 C 2/(2r j +1) a (log n) (2-ν j ) + -2ν j r j ν j (2r j +1) (log p) 2r j 2r j +1  
where C is a numerical constant, so that it follows from Corollary 1 that the block LASSO estimator is minimax optimal up to, at most logarithmic factor of p.

The two examples above illustrate the situation when estimator (3.4) attains nearly optimal convergence rates when p > n. This, however, is not always possible. Note that our analysis of the performance of the estimator (3.4) relies on the fact that the eigenvalues of any sub-matrix Σ Λ are close to those of matrix Σ Λ (Lemma 1). The latter requires n ≥ N where N depends on the nature of vectors W i . The next example shows that sometimes n < p prevents Lemma 1 from being valid.

Example 3 Orthonormal dictionary Let W i , i = 1, • • • , n, be uniformly distributed on a set of canonical vectors e k , k = 1, • • • , p. Then, Ω = I p /p, so that ω min = ω max = 1/p. Moreover, W 1 2 2 = 1 and |W (j) | ≤ 1. Therefore, for any µ > 0,

C µ = 1, U 2 µ (ℵ) = 1.
In the case of the orthonormal dictionary, N in (3.13) is of the form N = C(s + s 0 )(L + 1)p log(p). Under conditions of Theorem 2, the upper bound for the risk is of the form

∆ upper (s 0 , s, n, r) = C   σ 2 ps 0 log n n + j∈J pσ 2 n 2r j 2r j +1 C 2/(2r j +1) a (log n) (2-ν j ) + -2ν j r j ν j (2r j +1) (log p) 2r j 2r j +1   ,
so, n ≥ N implies n > Cσ 2 p(s 0 + s) log n which also guarantees that the risk of the estimator is small. This, indeed, coincides with one's intuition since one would need to sample more than p vectors in order to ensure that each component of the vector has been sampled at least once.

Discussion

In the present paper, we provided a non-asymptotic minimax study of the sparse high-dimensional varying coefficient model. To the best of our knowledge, this has never been accomplished before. An important feature of our analysis is its flexibility: it distinguishes between vanishing, constant and time-varying covariates and, in addition, it allows the latter to be heterogeneous (i.e., to have different degrees of smoothness) and spatially inhomogeneous. In this sense, our setup is more flexible than the one usually used in the context of additive or compound functional models (see, e.g., [START_REF] Dalalyan | Statistical inference in compound functional models[END_REF] or [START_REF] Raskutti | Minimax-optimal rates for sparse additive models over kernel classes via convex programming[END_REF]).

The adaptive estimator is obtained using block LASSO approach which can be viewed as a version of group LASSO where groups do not occur naturally but are rather driven by the need to reduce the variance, as it is done, for example, in block thresholding. Since we used tensor approach for derivation of the estimator, we believe that the results of the paper can be generalized to the case of the multivariate varying coefficient model studied in [START_REF] Zhu | Multivariate varying coefficient model for functional responses[END_REF].

An important feature of our estimator is that it is fully adaptive. Indeed, application of the proposed block LASSO technique does not require the knowledge of the number of the non-zero components of f . It only depends on the highest diagonal element of matrix Ω which can be estimated with high precision even when n is quite small due to Lemma 1.

Note that, even when p is larger than n, the vector function f is completely identifiable due to Assumption (A5), as long as the number of non-zero components of f does not exceed ℵ in condition (2.5) and the number of observations n is large enough. The examples of the paper deal with the dictionaries such that (2.5) holds for all ℵ = 1, • • • , p. The latter ensures identifiability of f provided n ≥ N where N is specified for each type of the random dictionary and depends on the sparsity level of f . On the other hand, large values of p ensure great flexibility of the choice of f , so one can hope to represent the data using only few components of it.

Finally, we want to comment on the situation when the requirement n ≥ N is not met due to lack of sparsity or insufficient number of observations. In this case, f is not identifiable and one cannot guarantee that f is close to the true function f . However, this kind of situations occur in all types of high-dimensional problems.

Proofs

Proofs of the lower bounds for the risk

In order to prove Theorem 1, we consider a set of test vector functions f ω

(t) = (f 1,ω , • • • , f p,ω ) T indexed by binary sequences ω with components f k,ω (t) = ω k0 u k + 2l 0k -1 l=l 0k ω kl v k φ l (t), (5.1)
where

ω kl ∈ {0, 1} for l = l 0k , l 0k + 1, • • • , 2l 0k -1, k = 1, • • • , p. Let K 0 and K 1 , respectively, be the sets of indices such that K 0 ∩ K 1 = ∅ and u k = u if k ∈ K 1 and u k = 0 otherwise, v k = v if k ∈ K 0 and v k = 0 otherwise, so that u k v k = 0.
In order assumption (2.3) holds, one needs u ≤ C a and

v ν j 2l 0k -1 l=l 0k (l + 1) ν j r j ≤ (C a ) ν j , j ∈ Υ. (5.2)
By simple calculations, it is easy to verify that condition (5.2) is satisfied if we set

u ≤ C a , v = C a (2l 0k ) -(r k +1/2) , (5.3)
where the constancy of v implies that l 0k in (5.3) are different for different values of k.

Consider two binary sequences ω and ω and the corresponding test functions f (t) = f ω (t) and f (t) = f ω (t) indexed by those sequences. Then, the total squared distance in L 2 ([0, 1]) between f ω (t) and f ω (t) is equal to

D 2 = u 2 k∈K 1 |ω k0 -ωk0 | + v 2 k∈K 0 2l 0k -1 l=l 0k |ω kl -ωkl |.
(5.4)

Let P f and P f be probability measures corresponding to test functions f and f , respectively. We shall consider two cases. In Case 1, the first s 0 functions are constant and the rest of the functions are equal to identical zero. In Case 2, the first s functions are time-dependent and the rest of the functions are equal to identical zero. In both cases, f (t) and f (t) contain at most s + s 0 non-zero coordinates. Using that conditionally on W i and t i , the variables ξ i are Gaussian N (0, 1), we obtain that the Kullback-Leibler divergence K(P f , P f ) between P f and P f satisfies

K(P f , P f ) = (2σ 2 ) -1 E n i=1 Q i (f ) -Q i ( f ) 2 = (2σ 2 ) -1 n E Q 1 (f ) -Q 1 ( f ) 2 where Q i (f ) = W T i f (t i
) and, due to conditions (A2) and (A5),

E Q 1 (f ) -Q 1 ( f ) 2 = E (f -f ) T (t 1 )W 1 W T 1 (f -f )(t 1 ) = E (f -f ) T (t 1 ) Ω (f -f )(t 1 ) ≤ ω * max E (f -f )(t 1 ) 2 2 ≤ ω * max φ max D 2 ,
where ω * max and D 2 are defined in (2.7) and (5.4), respectively. In order to derive the lower bounds for the risk, we use Theorem 2.5 of [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF] which implies that, if a set Θ of cardinality M + 1 contains sequences ω 0 , • • • , ω M with M ≥ 2 such that, for any j = 1, • • • , M , one has f ω 0 -f ω j | ≥ D > 0, P ω j << P ω 0 and K(P

f j , P f 0 ) ≤ κ log M with 0 < κ < 1/8, then inf ω sup f ω ,ω∈Θ P ( f ω -f ω 2 ≥ D/2) ≥ √ M 1 + √ M 1 -2κ - 2κ log M . (5.5)
Now, we consider two separate cases.

Case 1. Let the first s 0 functions be constant and the rest of the functions be equal to identical zero. Then, v = 0 and K 1 = {1, • • • , s 0 }. Use the Varshamov-Gilbert Lemma (Lemma 2.9 of [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]) to choose a set Θ of ω with card(Θ) ≥ 2 s 0 /8 and D 2 ≥ u 2 s 0 /8. Inequality K(P f j , P f 0 ) ≤ (2σ 2 ) -1 n ω * max φ max u 2 s 0 /8 ≤ κ log (card(Θ)) holds if u 2 = 2 σ 2 κ/(n ω * max φ max ). Then, D 2 = (s 0 σ 2 κ)/(4n ω * max φ max ) and u ≤ C a provided n ≥ 2 σ 2 κ/(C 2 a ω * max φ max ).

Case 2. Let the first s functions be time-dependent and the rest of the functions be equal to identical zero. Then u = 0, v is given by formula (5.3) Use Varshamov-Gilbert Lemma to choose a set Θ of ω with card(Θ) ≥ 2 L/8 and D 2 ≥ v 2 L/8.

, K 0 = {1, • • • , s}. Let r k , k ∈ K 0 ,
Inequality K(P f j , P f 0 ) ≤ κL/8 holds if v 2 ≤ (σ 2 κ)/(4n ω * max φ max ),
which, together with (5.3) and (5.4) imply that

l 0k = 1 2 4 C 2 a n ω * max φ max σ 2 κ 1 2r k +1 + 1, D 2 ≥ C 2 a 16 s k=1 4 C 2 a n ω * max φ max σ 2 κ - 2r k 2r k +1
where x denotes the integer part of x. Condition L ≥ 3 is satisfied for any s ≥ 1 provided n ≥ n low .

Proofs of the upper bounds for the risk

Proof of Theorem 2.

For any α ∈ R p(L+1) , one has

n -1 Bâ -Y 2 2 + δ â block ≤ n -1 Bα -Y 2 2 + δ α block .
Consider a set F 1 ⊆ W ⊗n µ such that (3.12) holds for any t and W ∈ F 1 , and a set F 2 ⊆ W ⊗n µ such that, (5.36) hold. Let F = F 1 ∩ F 2 . Lemma 2 implies that on the event F

2 â -α, B T b n ≤ δ â -α block .
Consider a set Ξ of values of the vector ξ such that 2σn -1 âα, B T ξ ≤ δ âα block for ξ ∈ Ξ.

(5.6)

Using (1.4), for ξ ∈ Ξ, any t and W ∈ F, one obtains

B(â -a 0 ) 2 2 n ≤ B(α -a 0 ) 2 2 n -2 δ â block + 2 δ α block + δ â -α block . (5.7)
Since α is an arbitrary vector, setting α = a in (5.7) yields δ â block -δ a block ≤ 0.5 δ âa block .

Let the set J 0 contain the indices of nonzero blocks of a 0 :

J 0 = {(j, l) : a jl 2 = 0} .

(5.8)

Then, the last inequality implies

(i,j)∈J C 0 (a -â) ij ≤ 3 (i,j)∈J 0 (a -â) ij (5.9)
From Lemma 1 it follows that

λ min ( Σ) ≥ (1 -h)φ min ω * min and λ max ( Σ) ≤ (1 + h)φ max ω * max
where ω * min and ω * max are defined in (2.7). For 1 ≤ j ≤ p consider sets

G 00 = {j : 1 ≤ j ≤ p, α j0 = 0}, G 01 = {j : 1 ≤ j ≤ p, α j0 = 0}, G j0 = {l : 1 ≤ l ≤ M, α jl 2 = 0}, G j1 = {l : 1 ≤ l ≤ M, α jl 2 = 0}
We choose α j0 = a j0 if a j0 = 0 and α j0 = 0 otherwise. Let sets G j1 be so that l ∈ G j1 iff

a jl 2 2 > ε = 8 2 σ 2 C ω K 2 µ + 1 log p nλ min ( Σ) .
We set α jl = a jl if j ∈ J and l ∈ G j1 and α jl = 0 otherwise where J is the set of indices corresponding to non-constant functions f j .

With δ = 2 δ, inequality (5.9) and Lemma 3 guarantee that

B(â -a) 2 2 n ≥ C B λ min ( Σ) â -a 2 2 .
(5.10)

On the other hand, Lemma 1 and the definition of α imply that

B(α -a) 2 2 n ≤ λ max ( Σ) α -a 2 2 .
(5.11)

Then, using (5.10) and (5.11), we rewrite inequality (5.7) as

C B λ min ( Σ) â -a 0 2 2 ≤ λ max ( Σ) (α -a 0 )
C B λ min ( Σ) ,
and similar inequality applies to the first sum in (5.12). By subtracting 0.5

C B λ min ( Σ) â -a 2 2
from both sides of (5.12) and plugging in the values of δ and α, derive:

(â -a) 2 2 ≤ C B λ max ( Σ) λ min ( Σ)   p j=1 l∈G j0 a jl 2 2 + 8 2 σ 2 C ω K 2 µ + 1 (s 0 + s) log p nλ min ( Σ) (5.13)
+ j∈J 8 2 σ 2 C ω K 2 µ + 1 card(G j1 ) log p nλ min ( Σ)   .
Observe that for any 1 ≤ j ≤ p one has

l∈G j0 a jl 2 2 + 8 2 σ 2 C ω K 2 µ + 1 card(G 1 j) log p n λ min ( Σ) ≤ M l=1 min a jl 2 2 , 8 2 σ 2 C ω K 2 µ + 1 log p n λ min ( Σ) . Application of inequality (5.31) with ε = 8 2 σ 2 C ω K 2 µ + 1 log p nλ min ( Σ) log n (see Lemma 4) yields â -a 2 2 ≤ C B λ max ( Σ) λ min ( Σ) C ω σ 2 K 2 µ + 1 (s 0 + s) log p nλ min ( Σ) (5.14) + j∈J C 2/(2r j +1) a C ω σ 2 K 2 µ + 1 nλ min ( Σ) 2r j 2r j +1 (log n) (2-ν j ) + -2ν j r j ν j (2r j +1) (log p) 2r j 2r j +1   
Denote r * j = r j ∧ r j and r * min = min j r * j . Choose L + 1 ≥ n 1/2 and note that r * min ≥ 2. Using (5.33), we obtain f -

f 2 2 ≤ â -a 2 2 + C 2 a s(L + 1) -2r * min , which implies f -f 2 2 ≤ â -a 2 2 + C 2 a s n 2 .
Now, (5.14) together with Lemmas 1,2 and 4 imply

P f -f 2 2 ≤ ∆(s 0 , s, n, r) ≥ 1 -8 p -µ .
Proof of Theorem 3. Let sets F be such that (5.14) holds. Then, P ( F) ≥ 1 -8p -µ and (5.14) yields

E â -a 2 2 ≤ E â -a 2 2 I( F) + E â -a 2 2 I( FC ) ≤ C ω * max ω * min   σ 2 s 0 log n nω * min + j∈J C 2/(2r j +1) a σ 2 nω * min 2r j 2r j +1 (log n) (2-ν j ) + -2ν j r j ν j (2r j +1) (log p) 2r j 2r j +1   + 16 C 2 f p 1-µ ≤ C ∆ upper (s 0 , s, n, r),
due to (3.21).

Proof of Lemma 1.

In order to simplify the notations, we set ω max (ℵ) = ω max and U µ = U µ (ℵ). Let W µ be the set of points described in condition (2.6) of Assumption (A5). Denote the direct product of n sets W µ by W ⊗n µ . Then, P(W ⊗n µ ) ≥ 1 -2 n p -2µ . Consider random matrices

Z i = (Σ i ) Λ -Σ Λ = (Ω i ) Λ ⊗ Φ i -Ω Λ ⊗ Φ, ζ i = (Σ i ) Λ I(W µ ) -E((Σ i ) Λ I(W µ )).
Then, ζ i are i.i.d. with Eζ i = 0. We apply the matrix version of Bernstein's inequality, given in Tropp [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF]:

Proposition 1 (Theorem 1.6, Tropp (2011)) Let ζ 1 , . . . , ζ n be independent random matri- ces in R m 1 ×m 2 such that E(ζ i ) = 0. Define σ ζ = max    1 n n i=1 E ζ i ζ T i 1/2 , 1 n n i=1 E ζ T i ζ i 1/2    .
and suppose that ζ i ≤ T for some T > 0. Then, for all t > 0, with probability at least 1 -e -t one has

1 n n i=1 ζ i ≤ 2 max σ ζ t + log(d) n , T t + log(d) n , (5.15
)

where d = m 1 + m 2 .
In order to find σ ζ , note that

1 n n i=1 E(ζ i ζ T i ) = E(ζ 1 ζ T 1 ) ≤ E((Σ 1 ) Λ (Σ 1 ) T Λ I(W µ ) + E((Σ 1 ) Λ I(W µ ) E((Σ 1 ) T Λ I(W µ ) = E [((Ω 1 ) Λ ⊗ Φ 1 )((Ω 1 ) Λ ⊗ Φ 1 )I(W µ )] + E [((Ω 1 ) Λ ⊗ Φ 1 )I(W µ )] 2 = E [((Ω 1 ) Λ (Ω 1 ) Λ ) ⊗ (Φ 1 Φ 1 )I(W µ )] + E((Ω 1 ) Λ I(W µ )) 2 E(Φ 1 ) 2 ≤ E(Φ 1 Φ 1 ) E [((Ω 1 ) Λ (Ω 1 ) Λ )I(W µ )] + E((Ω 1 ) Λ I(W µ )) 2 E(Φ 1 ) 2 ,
and, similarly,

1 n n i=1 E(ζ T i ζ i ) ≤ E(Φ 1 Φ 1 ) E [((Ω 1 ) Λ (Ω 1 ) Λ )I(W µ )] + E((Ω 1 ) Λ I(W µ )) 2 E(Φ 1 ) 2 .
Here,

E(Φ 1 Φ 1 ) ≤ C 2 φ (L + 1)Φ = C 2 φ (L + 1)φ max and E [((Ω 1 ) Λ (Ω 1 ) Λ )I(W µ )] = E W Λ W T Λ W Λ W T Λ I(W µ ) ≤ U 2 µ E(W Λ W T Λ ) = U 2 µ Ω Λ = U 2 µ ω max , so that σ 2 ζ ≤ 2 C 2 φ U 2 µ (L + 1)φ max ω max . (5.16)
Now, observe that, since matrix E((Σ i ) Λ I(W µ )) is non-negative definite and matrices Φ i and (Ω i ) Λ have rank one for any i, one has

T = sup ζ 1 ≤ 2 sup (Σ 1 ) Λ I(W µ )) = 2 sup (Ω 1 ) Λ I(W µ ) Φ 1 ≤ 2 C 2 φ U 2 µ (L + 1
). (5.17)

Apply Bernstein inequality (5.15) with σ 2 ζ and T given by formulae (5.16) and (5.17), respectively. Then, we obtain for any t > 0, with probability at least 1 -e -t (5.18) In order to apply inequality (5.18) to Z i , observe that Z i -

ζ i = (Σ i ) Λ I(W c µ ) -E((Σ i ) Λ I(W c µ ) and E((Σ i ) Λ I(W c µ )) 2 = E (Ω i ) Λ I(W c µ ) ⊗ E(Φ i ) 2 ≤ E (Ω i ) Λ I(W c µ ) 2 E(Φ i ) 2 ≤ E (Ω i ) Λ I(W c µ ) 2 2 E Φ i 2 2 ≤ 2 C 4 φ V ℵ p -2µ (L + 1) 2 , (5.19)
due to the fact that E(W (j) ) 4 ≤ V for any j = 1, • • • , p.

Proof. First, let us show that for J ≥ 1 |a k | ν k r ν ≤ C 2 a (J + 1) -2r , so that (5.33) is valid. Now, using (5.33), we prove (5.31). Again, we consider cases ν ≥ 2 and 1 ≤ ν < 2, separately. If ν ≥ 2, then partitioning the sum into the portion for l ≤ J and l > J (which corresponds to k > Jd), we derive Minimization of the last expression with respect to J yields (5.31).

In order to prove (5.32), observe that for any J > 1 and r * > 3/2, one has

f -f J ∞ ≤ sup t∈[0,1] ∞ l=1 J(l+1)-1 k=Jl a 2 k J(l+1)-1 k=Jl φ 2 k (t) ≤ C a C φ ∞ l=1 √ Jl(Jl + 1) -r * ≤ C a C φ J -(r * -1/2)
which completes the proof.

Lemma 5 Let r * j = r j ∧ r j , r * = min j r * j ≥ 2 in assumption (A3). Let µ in (2.6) be large enough, so that p µ ≥ 2 n (5.34)

and n be such that n ≥ U 2 µ C 2 φ (L + 1)µ log p g 2 ω max (s) (5.35) where U µ = U µ (s + s 0 ). Then, one has

P n -1 b 2 2 ≤ 3 g 2 ω max (s) s C 2 a (L + 1) -4 ∩ W ∈ W ⊗n µ ≥ 1 -2 p -µ .
(5.36)

Here b is the vector with components b i = W T i ρ(t i ), i = 1, • • • , n, where ρ(t) = (ρ 1 (t), • • • , ρ p (t)) T and ρ i (t) are defined in (1.2).

Proof. In order to estimate b 2 2 we apply Bernstein inequality to the centered random variables β i = b 2 i I(W µ ) -E b 2 i I(W µ ) . We start with establishing an upper bound for Eb 2 i :

Eb 2 i = E n -1 b 2 2 = E ρ(t) Ω 1 ρ(t) ≤ ω max (s)E ρ(t) 2 2 .

(5.37)

Using (5.33) and the orthonormality of the basis, we derive E ρ(t) It is also easy to see that, for any i and W ∈ W ⊗n µ , by (5.32), one has,

U b = max(b 2 i ) ≤ W J

  )

  coincide with the values of finite components of vector r. Denote L = s k=1 l 0k .

1 (J + 1 )

 11 -2r . Since (ν -2)/(2rν) < 1 for ν ≥ 2, inequality (5.33) holds. If 1 ≤ ν < 2, then

2 2min( a l 2 2min( a l 2 2

 222 , εd) ≤ Jεd + C 2 a (Jd) -2r .Minimizing the last expression with respect to J, we obtain (5.31) without the log-factor. If1 ≤ ν < 2, then ∞ l=1 , εd) ≤ Jεd + (dε) 1-ν/2 , εd) ≤ Jεd + C ν a (dε) 1-ν/2 (Jd) -r ν .

σ 2 b

 2 and U b given by (5.39) and (5.40), respectively, and t = µ log(p), one obtains, that with probability at least 1 -p -µ1 n n i=1 β i ≤ 2 max σ b µ log p n , U b µ log p n .(5.41)

  2 2 ≤ s g 2 C 2 a (L + 1) -2r * .(5.38)Plugging (5.38) into (5.37), we obtainEb 2 i ≤ g 2 ω max (s) s C 2 a (L + 1) -2r * .Now, we use the upper bound (5.32) to establish an upper bound for the variance σ 2 b of b 2 i :σ 2 b ≤ E b 4 i I(W µ ) ≤ U 4 µ E ρ(t) 4 2 ≤ g 2 U 4 µ s 2 C 4 a C 4 φ (L + 1) -4r * +2 .(5.39)

Acknowledgments

Marianna Pensky was partially supported by National Science Foundation (NSF), grant DMS-1106564. The authors want to thank Alexandre Tsybakov for extremely valuable suggestions and discussions.

Now, we use the union bound over all Λ such that |Λ| ≤ ℵ, the inequality p ℵ ≤ ep ℵ ℵ and choose t = 2 µ ℵ log( ep ℵ ). Combining (5.18) and (5.19), for all Λ such that |Λ| ≤ ℵ, we derive inf

for any z such that

Note that, under condition (3.10), one has

It is easy to check that, whenever n ≥ N (λ) where N (λ) is defined in (3.9), condition (5.20) is satisfied with

which, together with condition p µ ≥ 2n, implies that inf

In order to complete the proof, observe that

Proofs of supplementary lemmas

Lemma 2 Let vector α ∈ R p(L+1) be partitioned into subgroups the way it has been done for vector a = Vec(A). Let L + 1 ≥ n 1/2 and n ≥ N where N is defined in (3.13). Suppose that Assumptions (A5), (A4) are valid and µ in (2.6) is large enough, so that condition (3.15) holds. Then, (i)

where δ is defined in (3.14).

Proof. In order to prove (i), note that

Fix vectors W 1 , • • • , W n and t. Using Hoeffding-type inequality for sub-gaussian random variables (see, e.g., Proposition 5.2 in [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF]) we obtain that, for any µ ≥ 1,

where C > 0 is an absolute constant and |Λ| ≤ 1.

For the second maximum in (5.24), we use a corollary of the Hanson-Wright inequality for sub-gaussian random vectors (see Theorem 2.1 in [START_REF] Rudelson | Hanson-Wright inequality and sub-gaussian concentration[END_REF]), which states that, for any matrix Q ∈ R n×d and any µ ≥ 1, one has

We apply this inequality with matrix Q = n -1/2 B j,l , where B j,l is a n × d sub-matrix of matrix

where we used M = L/ log(n). Now, consider a set F 1 ⊆ W ⊗n µ such that (3.12) holds for any t and any W ∈ F 1 . By Lemma 1, obtain P(F 1 ) ≥ 1 -2 p -µ . Then,

with probability at least 1 -5 p -µ , due to condition (3.15). Inequality (ii) is proved in a similar manner. Indeed, one has

Consider a set F 1 ⊆ W ⊗n µ such that (3.12) holds for any t and any W ∈ F 1 . By Lemma 1, obtain that P(F 1 ) ≥ 1 -2p -µ . Then, on the event F 1 , the following inequalities are valid

(5.28)

Now consider a set F 2 ⊆ W ⊗n µ such that (5.36) holds. Let F = F 1 ∩ F 2 . Lemmas 1 and 5 imply that on the event F

Now, the statement of the Lemma is valid due to assumption (L + 1) 2 ≥ n.

Lemma 3 Let S ∈ {1, . . . , p} be a subset of indices with |S| ≤ (s+s 0 ) and J = S ×{0, . . . , L} ⊂ {1, . . . , p} × {0, . . . , L}. We define A to be the following set of vectors

Suppose that the assumptions of Lemma 1 and Assumption (A6) hold. Then, there exists a numerical constant C B such that, where ω * min is defined in (2.7). Proof. The proof follows the lines of the proof of Lemma 4.1 in [START_REF] Bickel | Simultaneous analysis of Lasso and Dantzig selector[END_REF].

for some C a > 0, 1 ≤ ν < ∞ and r > min(1/2, 1/ν). Let a l be blocks of coefficients a k of length d, so that a l = (a (l-1)d+1 , • • • , a ld ). Then, for any ε > 0 one has

(5.31)

Moreover, if r ≥ 2 and basis {φ k } satisfies assumption (A1), then

.32)

Here, r * = min(r, r ), (x) + = x if x > 0 and zero otherwise.


Baptiste Hemery

Jean-Jacques Schwartzmann

Christophe Rosenberger

Study on color spaces for single image enrolment face authentication


Zakaria Guelilia
email: zakaria.guelilia@insa-rennes.fr

Renaud Loison

Raphael Gillard

MacroModel based DG-FDTD for Calculating Local Dosimetry in a Variable and Highly Multiscale Problem

This paper proposes a method to estimate human exposure to electromagnetic field radiation in a variable and highly multiscale problem. The electromagnetic field is computed using a combination of two methods: a rigorous time domain and multiscale method, the DG-FDTD (Dual Grid Finite Difference Time Domain) and a fast substitution model based on the use of transfer functions. The association of these methods is applied to simulate a scenario involving an antenna placed on a vehicle and a human body model located around it. The purpose is to assess the electromagnetic field in the left eye of the human body model. It is shown that this combination permits to analyse many different positions in a fast and accurate way.

INTRODUCTION.

The estimation of human exposure to electromagnetic radiation is of growing interest with the development and the daily presence of communication objects. Technologies must comply with legislations in which exposure limits often come from ICNIRP guidelines [1]. Numerical dosimetry is an essential approach for the assessment of human exposure to electromagnetic radiation. When dealing with actual exposure scenarios involving a fixed transmitter placed in a realistic environment, numerical methods face several difficulties. Firstly, it is necessary to consider a large computational volume while describing precisely the human body or parts of it (for local dosimetry typically). Secondly, the exposure has to be assessed for all possible body's positions with respect to the transmitter in a defined area. Finally, the exposure must be assessed for different body's morphologies. In this paper, we propose a numerical technique to solve the first two difficulties, namely multiscale and position variability aspects.

For numerical dosimetry, the FDTD (Finite Difference Time Domain) method is widely used. Indeed, it has proved its usefulness for treating electromagnetic problems involving complex and heterogeneous structures like the human body [START_REF] Taflove | Computational Electrodynamics: The finite-Difference Time-Domain Method, Third Edition[END_REF]. However, if a high resolution is required for the description of parts of the body such as eyes, the uniform meshing scheme of the FDTD technique leads to large oversampled areas. This dramatically increases the computation time and memory requirements. Moreover, for each new position of the body within the environment, a new entire simulation must be performed.

In order to solve these issues, several techniques based on the FDTD method have been presented.

To address the multiscale problem, several works highlight the use of subgridding FDTD [START_REF] Johansson | Wave-propagation from medical implants-influence of body shape on radiation pattern[END_REF][START_REF] Wiart | The analysis of human head interaction with a hand-held mobile using the non-uniform FDTD[END_REF]. This kind of methods uses finer resolution for the description of specific areas in a single volume. Although these approaches reduce the simulation time and computing resources, they have some disadvantages such as instabilities [START_REF] Berenger | A Huygens subgridding for the FDTD method[END_REF] or parasitic reflections at grids interfaces [START_REF] Chilton | Explicit 3D FDTD subgridding with provable stability and conservative properties[END_REF].

As an alternative to classical subgridding approaches, the Dual Grid FDTD has been proposed in [START_REF] Pascaud | Dual-grid fnite difference timedomains cheme for the fast simulation of surrounded antennas[END_REF]. It proposes to perform sequentially FDTD simulations each with its own resolution. A coarse FDTD simulation of the environment is combined with the fine FDTD simulation of the human body model. Thanks to that, instabilities are avoided and reflections are reduced. This method has shown that it can effectively treat a highly multiscale dosimetric problem [START_REF] Miry | An Efficient Bilateral Dual-Grid-FDTD Approach Applied to On Body Transmission Analysis and Specifc Absorption Rate Computation[END_REF]. Although this method is fast and robust, it is not suitable to address the variability issue since the multiplication of costly fine EM simulations leads to a prohibitive computation time.

To handle the problem of variability, several solutions have been proposed. Most of them are based on the construction of an alternative model for rapidly assessing exposure. Many of these models are based on statistical approaches. First of all, EANN (Evolutionary Artificial Neural Networks) is a tool often used in dosimetric problem to treat variability [START_REF] Colak | Estimation of local SAR level using RBFNN in threelayer cylindrical human model[END_REF][START_REF] Ortiz-Rodriguez | A Computational Tool Design for Evolutionary Artificial Neural Networks in Neutron Spectrometry and Dosimetry[END_REF]. However, it has some drawbacks in the optimal determination of the network topology. This kind of methods need a long process of data acquisition to construct the model. This learning phase requires numerous EM simulations and thus its implementation is time-consuming. Stochastic collocation [START_REF] Kientega | A surrogate model to assess the whole body SAR induced by multiple plane waves at 2.4 GHz[END_REF] or Monte carlo simulations [START_REF] Cheol | Monte Carlo simulation to measur light dosimetry within the biological tissue[END_REF][START_REF] Habachi | Statistical analysis of whole-body absorption depending on anatomical human characteristics at a frequency of 2.1 GHz[END_REF][START_REF] Aiouaz | Uncertainty analysis of the specific absorption rate induced in a phantom using a stochastic spectral collocation method[END_REF] are often proposed to face parameters variability. However, the construction or implementation of these tools requires large number of costly electromagnetic simulations. To finish, a method presented in [START_REF] Hirata | Estimation of Whole-Body Averaged SARs in Human Models for Far-Field Exposures in Whole-Body Resonance and GHz Frequency Regions[END_REF] proposes to replace the human body model by an antenna. This substitution model permits to estimate the Specific Absorption Rate (SAR) for the whole body but is not suitable for local SAR computation.

In this paper, we propose a new strategy wich combines the DG-FDTD and a fast and rigorous substitution model based on the macromodel proposed in [START_REF] Vairavanathan | A Reciprocity-Based Framework for the Efficient Modeling of Antenna-Wireless Channel Interaction[END_REF]. This combination is called MacroModel based DG-FDTD method (MM-DG-FDTD).

This new method is particularly suitable for problems in which the field has only to be known at a few points in the computation volume. This is definitively the case in dosimetry problems where only the local exposure has to be assessed. In such cases, a fast and rigorous substitution model using quite simple transfer functions will this be built to replace costly fine EM simulation steps. This paper is organised as follows. Firstly, the addressed scenario is described in section 2. In section 3, the MM-DG-FDTD approach is detailed by presenting the substitution model, its principle and its implementation. Section 4 presents the validation of this new DG-FDTD approach. To finish, an example of its application is presented in section 5.

PRESENTATION OF THE ADDRESSED SCENARIO

In order to illustrate the new method, we consider the multiscale and variable problem described in figure 1. The purpose is to calculate the electric field in the eye of a whole human body located near an antenna onboard a vehicle.

The human eye diameter is about 2.5 cm whereas the analysed scene has a length and a width of 10.24 m. This very high contrast demonstrates the multiscale feature of the considered problem. The calculation of the field must be carried out for different positions of the person around the vehicle. This aspect reveals the concept of variability.

The vehicle and the person are located in the near field of a monopole antenna matched at the frequency of 60 MHz. We consider this frequency as it is used for military communication systems and, furthermore, the whole body SAR reaches its maximum at this frequency. This shows the need to use a full wave method to analyse the addressed scenario. As the study objective is only to show our method effectiveness the power emitted in the scenario does not correspond to a real system. The amplitude of the transmitter voltage generator is simply fixed to 10 V . The human body model used in the simulations comes from the "Visible Human Project" [START_REF]The Visible Human Project[END_REF]. It corresponds to the Hugo homogeneous model except eyes which have their own specifications (table 1). The vehicle and the ground are modelled as PEC (Perfect Electric Conductor). The characteristics of the scenario are:

• h = 0.416 × λ 60M Hz = 2.08 m, • d = 0.384 × λ 60M Hz = 1.92 m,
• l = 0.896 × λ 60M Hz = 4.48 m, • D = L = 2.048 × λ 60M Hz = 10.24 m.

MACROMODEL BASED DG-FDTD METHOD

The aim of this section is to present the principle and the implementation of the MacroModel based on DG-FDTD method. The building of the substitution model and its insertion in the DG-FDTD method are described.

Principle of the MM-DG-FDTD method

The principle of the MM-DG-FDTD method is to decompose the electromagnetic calculation into two steps performed sequentially (figure 2).

The first step is a coarse FDTD simulation whose objective is to characterize the whole environment including the human body model. The resolution is mainly governed by the smallest considered wavelength. In our case, the body is thus described coarsely. A closed near field surface S is placed around the body in order to store the incident field incoming on the body. This field includes, at a coarse level, all EM interactions between the body and the environment.

In the second step, instead of performing a fine FDTD simulation of the body as in the classic DG-FDTD, the MM-DG-FDTD method proposes to use a rigorous substitution model. The input data are the N field components on the near field surface collected in step 1 { E step1 ; H step1 }, whereas the output is the electric field at one point in the eye ( E step2 ) :

E step2 u ( r out ) = N i=1 [ F E u ( r out , r in i ). E step1 ( r in i ) + F H u ( r out , r in i ). H step1 ( r in i )] (1)
where:

• u = x, y or z,

• r in i belongs to the near field surface S,

• r out locates the point in the eye where the output field has to be computed,

• N is the number of cells over the surface S,

• F E u and F H u are transfer functions.

Computation of transfer functions

To build the substitution model, an equivalent problem is firstly considered. The field { E step1 ; H step1 } on the near field surface S is replaced by surface sources { J step1 S ; M step1 S } using Huygens theory [START_REF] Jin | Theory and computation of electromagnetic fields[END_REF]:

J step1 S = -H step1 ∧ n (2)
M step1 S = E step1 ∧ n (3)
To determine the transfer functions F E u and F H u , a reciprocal problem is considered in which a source J R is located at r out . This source generates a field { E R ; H R }.

In the frequency domain and according to the reciprocity theorem [START_REF] Jin | Theory and computation of electromagnetic fields[END_REF], sources and fields of the inital problem are related to those of the reciprocal one as follows:

S ( E R . J step1 S -H R . M step1 S )dS = V ( E step2 . J R )dV ( 4
)
where V is the volume inside the near field surface S. This relation shows that if the fields of the reciprocal problem { E R ; H R } in response to a given J R are known on S then the transfer functions of equation 1 can be derived easily. To do so, a single fine FDTD simulation of the reciprocal problem is performed. For this simulation, the fine FDTD computation volume is the same as that used in the second step of the classical DG-FDTD (figure 2) but it is excited by the localized electrical source J R ( r out ). During this simulation, { E R ; H R } are recorded on the near field surface S. Expressed in the frequency domain thanks to FFTs, J R , E R and H R are used to determine the transfer functions as it is explained below.

For the sake of simplicity, only one face of the near field surface is considered in the following equations. This face is discretized into M cells located. All these cells have equal surface area ∆S and equal unitary normal vectors n. The currents are constant on each cell, thus the left side of 4 becomes :

M i=1 ∆S[ E R ( r in i ). J step1 S ( r in i ) -H R ( r in i ). M step1 S ( r in i )] (5)
The source J R is oriented along e u and is constant in the FDTD cell of volum ∆V located in r out , thereby, the right side of 4 is given by :

E step2 u ( r out ).J R u ( r out )∆V (6)
Combining ( 2), ( 3), ( 4), ( 5) and ( 6) leads to :

F E u ( r out , r in i ) = - ∆S ∆V . n ∧ H R ( r in i ) J R u ( r out ) (7)
F H u ( r out , r in i ) = - ∆S ∆V . n ∧ E R ( r in i ) J R u ( r out ) (8)
One fine FDTD simulation of the reciprocal problem permits to determine the transfer functions { F E u ( r out , r in i ); F H u ( r out , r in i )}. This single simulation is then enough to determine the field component E step2 u ( r out ) in response to any incident fields { E step1 ; H step1 } using (1).

Method overview

The method starts with the computation of transfer functions. In order to determine the three output field components ( E step2 ( r out )), it is sufficient to perform three simulations of the reciprocal problem with the three source orientations : J R x ( r out ), J R y ( r out ) and J R z ( r out ). The transfer functions can therefore be calculated and be used again for any incident field { E step1 ; H step1 }. Once the transfer functions calculated, the electromagnetic analysis of the addressed scenario for each position of the body is derived as follow

• The first step of the DG-FDTD method is run and the incident fields { E step1 ; H step1 } is recorded over the near field surface S. • After being expressed in the frequency domain, the incident field is combined with transfer functions in order to obtain fastly and precisely the desired output field E step2 ( r out ).

Next section is about the validation of this MM-DG-FDTD method.

VALIDATION OF THE MACROMODEL BASED DG-FDTD METHOD

This section presents the validation of the MM-DG-FDTD method. Firstly, a canonical case is considered for the validation of the substitution model. Next, the new method is applied to the addressed scenario of figure 1 in order to assess its performance.

Validation of the substitution model for a canonical case

To validate the substitution model, a canonical problem is first considered (Figure 3). A volume is illuminated by a z polarized plane wave of amplitude 10 V /m which propagates along the x direction. The time domain excitation covers a 100MHz frequency bandwith ([5;115]MHz). A dielectric block is placed in the simulation box. It has the same specifications as the human body model (table 1). The objective is to estimate the three electric field components at an output point located in the dielectric block. Two simulations are carried out: a standard FDTD one which serves as a reference and a simulation using the substitution model. Table 2 gives the parameters of the reference FDTD simulation and the one used to calculate the transfer functions. The relative difference between the two results is calculated using this relation:

∆E u (f ) = | E subs u (f ) -E F DT D u (f ) E F DT D u (f )
| % with : 4 presents the relative difference between FDTD and substitution model for the three electric field components versus frequency. Even if the substitution model is rigorous, a slight relative error smaller than 2% is observed. This can be explained by the difference in the numerical implementations of the two methods. Anyhow, this agreement is found sufficient for dosimetry applications, as will be shown in the next section.

• u = x,
The analysis of the canonical case has shown that in a simple problem, the substitution model works well. Next subsection assesses its performance in a variable and multiscale problem.

Application and performance of the MM-DG-FDTD method

We consider the addressed scenario (figure 1) where the aim is to assess the three field components at the center of the eye in response to an illumination resulting from the radiation of an antenna located on a vehicle. The study is done at the frequency of 60 MHz and several positions of the body around the vehicle are considered.

To assess the performance of our approach, a comparison between a classical DG-FDTD and the MM-DG-FDTD method is done. The description of the simulation steps are presented in figure 2. The 3. Concerning the new method, the simulation parameters and the computation time associated to the calculation of the transfer functions are given in tables 4 and 5 respectively. The mesh is selected according to the smallest detail in the corresponding computation volume (i.e the eye). Such a grid density is required to describe it with a sufficient accuracy. 6 and7 show the computation time obtained with the two approaches for the study of W different body positions. Compared to a classical DG-FDTD method, the MM-DG-FDTD requires an incompressible 748 min CPU time to calculate the transfer functions. Thus, for a single position study, the use of the substitution model is useless compared to the classical DG-FDTD. However, the study of a new position only requires an additionnal 0.79 min CPU time with the new method whereas the additional cost of classical DG-FDTD is 357 min. Therefore, the MM-DG-FDTD method rapidly outperforms the classical DG-FDTD. As figure 5 shows, the use of substitution model presents an interest in terms of calculation time for the study of more than two positions. As an example, an 8-position study requests a time of 2856 min whereas the new method only needs 754.3 min.

It is important to note that the study is done for a single output point. A study conducted on several output points would require more transfer functions and thus an important increase of the computation time (748 min for each new output point).

Next section presents a study in which the new method is applied for all possible positions of the addressed scenario.

MM-DG-FDTD EXPLOITATION

The MM-DG-FDTD method is applied to analyse the addressed scenario for all possible positions of the body around the vehicle. In the coarse description of the environment, this corresponds to 888 positions and the study is done at the frequency of 60 MHz. In order to verify the validity of the MM-DG-FDTD simulation, table 8 compares the results of the classical DG-FDTD simulation with those of the new method for one particular position in the addressed scenario. The MM-DG-FDTD method provides results with an error smaller than 1%. Figure 6 presents the space orientation of the body and a field cartography. The cartography represents the module of the electric field at the center of the left eye for all 888 positions around the vehicle. Each pixel gives the amplitude of the electric field in the eye for a given position of the body. For readability, the body is not represented on this map. This cartography is normalized to its maximum. The whole study has been carried out in 1.4495 × 10 3 min with the proposed method. For comparison, it would require 3.1705 × 10 5 min with the classical DG-FDTD.

CONCLUSION.

This paper presents a new method based on the DG-FTD approach. It combines the DG-FDTD method with a substitution model to proceed highly multiscale and variable dosimetric problems. It is demonstrated that the proposed MM-DG-FDTD method provides quick results with a very good accuracy. Reciprocity principle, used in the calculation of transfer functions, works for isotropic materials and passive structures. Furthermore, the DG-FDTD has already been used successfully to analyse heterogeneous body models [START_REF] Miry | An Efficient Bilateral Dual-Grid-FDTD Approach Applied to On Body Transmission Analysis and Specifc Absorption Rate Computation[END_REF]. Thus, the MM-DG-FDTD can be also applied to a heterogeneous body models. To conclude, the MM-DG-FDTD is ready to be used to compute real life dosimetric problems.

CITATIONS

Figure 1 :

 1 Figure 1: Addressed scenario

Figure 2 :

 2 Figure 2: Application of classical DG-FDTD and the MM-DG-FDTD on the addressed scenario

  y or z, • E subs u represents the result for the substitution model, • E F DT D u represents the result for the FDTD simulation.

Figure 3 :

 3 Figure 3: Canonical case

Figure 4 :

 4 Figure 4: Canonical case's results: relative difference between FDTD and substitution model

Figure 5 :

 5 Figure 5: Computation time for both methods.

Figure 6 :

 6 Figure 6: Normalized results using the MM-DG-FDTD.

Table 1 :

 1

		Human body parameters
	F=60MHz	Permittivity (F/m) Conductivity (S/m)
	Human body	0.3970	45.55
	Eyes	0.8870	76.68

Table 2 :

 2 Simulation parameters

	Time Step dt	14.626 ps
	Spatial resolution	(8 × 8 × 8)mm
	(∆x × ∆y × ∆z)	≈ ( λ 60M Hz 625 ) 3 mm
	Volume size
	(N x × N y × N z ) cells	(50 × 50 × 30)
	Observation time T obs	200 ns
	Figure

Table 3 :

 3 Simulation parameters of the DG-FDTD method for one position

		First step	Second step
	Time Step dt	292.5 ps	14.626 ps
	Spatial resolution	(160 × 160 × 160) mm	(8 × 8 × 8)mm
	(∆x × ∆y × ∆z)	≈ ( λ 60M Hz 31	) 3 mm	≈ ( λ 60M Hz 625 ) 3 mm
	Volume size
	(N x × N y × N z ) cells	(90 × 90 × 90)	(150 × 150 × 300)
	Observation time T obs	300 ns	300 ns
	parameters of the DG-FDTD simulation are reported in table

Table 4 :

 4 Simulation parameters of the FDTD simulation used for the calculation of transfer functions

	Time Step dt	14.626 ps
	Spatial resolution	(8 × 8 × 8)mm
	(∆x × ∆y × ∆z)	≈ ( λ 60M Hz 625 ) 3 mm
	Volume size
	(N x × N y × N z ) cells (150 × 150 × 300)
	Observation time T obs	200 ns
	Tables

Table 5 :

 5 Computation time for the calculation of transfer functions Simulation time (220 × 3) min

				Post-processing		88 min
					Total time		748 min
	Table 6: Computation time for the classical DG-FDTD approach
										W positions
				Step one (160 mm)			W × 0.56 min
	Spatial interpolation (time domain)	W × 30 min
				Step two (8 mm)			W × 326.48 min
				Total time				W × 357.04 min
	Table 7: Computation time for the MM-DG-FDTD method
											W positions
	Transfer Functions step (8mm)			748 min
				Step one (160mm)				W × 0.56 min
	Spatial interpolation (frequency domain)		W × 0.067 min
		Step two (linear system)				W × 0.16 min
				Total time				748 + W × 0.79 min
		4000
						Computation time for the DG-FDTD method
		3500				Computation time for the MM-DG-FDTD method
											X: 8
		3000									Y: 2856
	Computation time (min)	1500 2000 2500 1000		X: 1 Y: 748.8							X: 8 Y: 754.3
				X: 1
		500		Y: 357
		0	0	1	2	3	4	5	6	7	8	9	10
							Wth positions

Table 8 :

 8 Results for one position in the addressed scenario Relative difference between DG-FDTD and MM-DG-FDTD |E x | 0.8526 % |E y | 0.48189 % |E z | 0.89392 % |E total | 0.73927 %


Guelilia Zakaria

Loison Renaud

Raphael Gillard

Combination of DG-FDTD with a Substitution Model for Calculating Local Dosimetry in a Variable and Highly Multiscale Problem

Keywords: EM analysis, multiscale, numerical dosimetry

This paper proposes a method to estimate human exposure to electromagnetic field radiation in a variable and highly multiscale problem. The electromagnetic field is estimated using a combination of two methods: a rigorous time domain and multiscale method, the DG-FDTD (Dual Grid Finite Difference Time Domain) and a substitution model based on the superposition principle and the use of transfer functions. The association of these methods is applied to simulate a scenario involving an antenna placed on a vehicle and a human body model located around it. The purpose is to assess the electromagnetic fields in the left eye of the human body model. It is shown that this combination permits to analyse many different positions in a fast and accurate way.

I. INTRODUCTION

Numerical dosimetry is an essential approach for the assessment of human exposure to electromagnetic radiation. In order to be more realistic to the actual exposure conditions, it is necessary to consider different environments involving large structures. In this context, it is common to treat very small elements compared to the overall size of the scenario. Then, we clearly face a contrast of scale which poses difficulties for electromagnetic simulations. In addition to these multiscale aspects, dosimetry problems have naturally a variable feature as the position of the person or its orientation. This variability is a challenge for numerical dosimetry which must be able to assess the exposure for a variety of configurations of the problem.

To solve the multiscale problem, the choice of the DG-FDTD [START_REF] Pascaud | Dual-grid fnite difference time-domainscheme for the fast simulation of surrounded antennas[END_REF], [START_REF] Miry | An Effcient Bilateral Dual-Grid-FDTD Approach Applied to On Body Transmission Analysis and Specifc Absorption Rate Computation[END_REF] turns out to be judicious. Indeed, this method has shown that it can effectively treat a highly multiscale dosimetry problem [START_REF] Guelilia | Utilisation de la DG-FDTD pour un calcul de dosimtrie dans un problme fortement multichelle : dtermination du DAS il pour une personne situe proximit d'une source HF/VHF embarque sur un vhicule[END_REF]. Although this method is fast, this approach is not suitable to solve the second issue because the multiplication of simulations leads to prohibitive computation time.

To handle the problem of variability, several solutions have been proposed. Most of them are based on the construction of an alternative model for rapidly assessing exposure. Many of these models are based on statistical approaches such as neural networks [START_REF] Colak | Estimation of local SAR level using RBFNN in threelayer cylindrical human model[END_REF] or stochastic collocation [START_REF] Silly-Carette | Dtermination de lincertitude sur le DAS en utilisant une mthode stochastique et la FDTD[END_REF]. However, the construction or implementation of these tools requires a large number of costly IT resources including electromagnetic In this paper, we propose a fast substitution model whose construction requires few electromagnetic simulations. It is composed of transfer functions and operates using the superposition theorem. In addition, we combine the DG-FDTD with this model to treat effectively multiscale scenarios with variable parameters. Firstly, the selected scenario is described. Then, The classic DG-FDTD approach is detailed. The next part presents the substitution model and it association with the DG-FDTD method. An application is performed for one position. To finish, the electrectromagnetic fields is assessed for several positions of the human body model using the classic DGFDTD method and the combination of the methods.

II. PRESENTATION OF THE SELECTED SCENARIO

The purpose of this study is to calculate the electromagnetic fields in an eye of a whole human body located near an antenna onboard a vehicle (figure [?]).

The calculation of the fields must be carried out for different positions of the person around the vehicle. The vehicle and the person are located in the near field of an antenna. The eye is very small compared to the rest of the environment. Thus, we face a multiscale problem. The studied frequency point of study is at 60 MHz. With this is scenario, 888 positions are considered.

III. PRINCIPLE OF THE DG-FDTD METHOD AND APPLICATION FOR ONE POSITION

A. Principle of the DG-FDTD method

The principle of the DG-FDTD method is to decompose a FDTD calculation in several steps performed sequentially [START_REF] Miry | An Effcient Bilateral Dual-Grid-FDTD Approach Applied to On Body Transmission Analysis and Specifc Absorption Rate Computation[END_REF], [START_REF] Guelilia | Utilisation de la DG-FDTD pour un calcul de dosimtrie dans un problme fortement multichelle : dtermination du DAS il pour une personne situe proximit d'une source HF/VHF embarque sur un vhicule[END_REF]. Each step involves a FDTD volume with its own resolution. The transfer of electromagnetic datas is realised through a sampling surface and an injection one.

B. Application for one position

The description of the simulation is presented in figure 2 and the parameters are given in table I. "Electromagnetic zoom" tool proposed by the DG-FDTD is used on the study.

The first step consists of a simulation with a large mesh and describing the whole problem. In fact, The mesh is chosen in a suitable way in order to describe all the elements of the scenario: the vehicle, the transmitter and the human body model. A near field surface is inserted around the person. The field collected in the first step is used in the second step to excite the volume describing more finely the human body model. This last step permits to compute in a more precise way the fields in the left eye.

The computation times associated with each step are reported in Table II. It is clear that step 2 is the most consuming in time and numerical resources. With the purpose of processing a new position of the person around the vehicle using DG-FDTD method, it is necessary to restart the two

IV. SUBSTITUTION MODEL

A. Principle

The substitution model intends to rigorously replace the second steps of the DG-FDTD calculation presented in the previous section. The input data of this model are the N field components on the near field surface collected in step 1. The output data are the three components of the field at the center of the human body eye model (step 2). In the frequency domain and using the superposition principle, the last step of the DG-FDTD can be replaced by a linear system compound of N inputs and three outputs. Indeed, the superposition principle permits to cumpute the output field through the sum of contributions of input field components (figure 3 and figure 4). To determine the N transfer functions effectively without performing N electromagnetic simulations, the reciprocity 5). It should be noted that the approach could be generalized to more outputs (several field points, for example). The three electromagnetic simulations necessary to determine the transfer functions are performed with the FDTD method. We limit the construction of the alternative model to replace only the simulation of step 2 of the initial problem (the most expensive step).

The multiplication of these transfer functions with the input field values can estimate in a short time the value of the three component field at the output point. These transfer functions can be used again to study several positions of the human body model around the vehicle.

B. Application and results

The study is done at the frequency of 60 MHz. The transfer functions are firstly determined using three FDTD simulations in accordance to figure 5. Three simulations are executed for the three electromagnetic field components. The parameters of these simulations are given in Table III. The transfer functions are then used to calculate the field components at the center of the eye in response to an illumination resulting from a complete simulation of the scenario (160 mm description).

Table IV compares the results of the classic DG-FDTD simulation (160 to 8 mm) with results of the substitution model. The combination of DG-FDTD and substitution model provides results with an error less than 1%.

The Tables V and VI show the computation time obtained with the two approaches for the study of one body position in the environment and for the study of ten body positions.

These results indicate that the use of substitution model presents an interest in terms of time calculation for the study of multiple positions. The study of a new position requires only 0.54 min.

V. APPLICATION ON THE ADDRESSED SCENARIO

Combination of DG-FDTD method and substitution model is applied on the selected scenario. The study is done at the frequency of 60 MHz. 888 positions are considered around the vehicle.

Figure 6 presents the orientations of the body in the space and a field cartography. The cartography represents the module of the electric fields in the center of the left eye when the body is located in the environment for 888 positions around the vehicle. This cartography is normalized to its maximum.

VI. CONCLUSION

This paper proposes to combine the DG-FDTD method with a substitution model to process a highly multiscale and variable problem for a numerical dosimetry. The combination of these two techniques is applied on a scenario and provides quick results with a very good accuracy.

Fig. 1 .

 1 Fig. 1. Selected scenario

Fig. 2 .

 2 Fig. 2. Application for one position

Fig. 3 .

 3 Fig. 3. Architecture of the substitution model

Fig. 4 .

 4 Fig. 4. Sketch of the substitution model

Fig. 6 .

 6 Fig. 6. Results normalized using the combination of methods.

TABLE I SIMULATION

 I PARAMETERS OF DG-FDTD METHOD

		First step	Second step
	Time Step	292.5 ps	14, 626 ps
	Observation time		300 ns	300 ns
		TABLE II
	RESULTS OF DG-FDTD SIMULATIONS
			Computing time
	First step		0.56 min
	Second step	356.48 min

Spatial resolution (dx × dy × dz)

(160 × 160 × 160) mm (8 × 8 × 8)mm

TABLE III SIMULATION

 III PARAMETERS FOR THE CALCULATION OF TRANSFER Spatial resolution (dx × dy × dz) (8 × 8 × 8) mm

	FUNCTIONS
	Time step	14, 626 ps
	Observation time	200 ns
	Simulation time	(220 × 3) min

TABLE V COMPUTATION

 V TIME WITH THE DG-FDTD APPROACH

		One position	Ten positions
	Step one (160 mm)	0.56 min	5.6 min
	Step two (8 mm)	356.48 min	3564.8 min
	Total	357.04 min	3570.4 min
		TABLE VI
	COMPUTATION TIME WITH THE SUBSTITUTION PATTERN
		One position	Ten positions
	Calcul of transfer Functions (8mm)	748 min	748 min
	Step one (160mm)		0.56 min	5.6 min
	Step two (8mm)		0.2 min	2 min
	Total	748.76 min	755.6 min