Isabelle Charlier
email: ischarli@ulb.ac.be

Davy Paindaveine

Jérôme Saracco
email: jerome.saracco@math.u-bordeaux1.fr

Conditional Quantile Estimation based on Optimal Quantization: From Theory to Practice 1

Keywords: Conditional quantiles, Optimal quantization, Nonparametric regression

Small-sample properties of a nonparametric estimator of conditional quantiles based on optimal quantization, that was recently introduced (J. Statist. Plann. Inference, 156, 14-30, 2015), are investigated. More precisely, (i) the practical implementation of this estimator is discussed (by proposing in particular a method to properly select the corresponding smoothing parameter, namely the number of quantizers) and (ii) its finitesample performances are compared to those of classical competitors. Monte Carlo studies reveal that the quantization-based estimator competes well in all cases and sometimes dominates its competitors, particularly when the regression function is quite complex. A real data set is also treated. While the main focus is on the case of a univariate covariate, simulations are also conducted in the bivariate case.

Introduction

Quantile regression is used to quantify the relationship between a univariate response variable Y and a d-dimensional covariate X. Since the conditional mean only models the average relationship between these, standard (mean) regression is way less informative than quantile regression, that provides a more complete picture of conditional distributions: quantile regression may indeed capture heteroscedasticity, conditional asymmetry, etc. Conditional quantiles allow to construct reference hypersurfaces (curves if the covariate is univariate) and conditional prediction intervals, hence are widely used in many different areas, including medicine, economics or lifetime analysis.

Quantile regression was first introduced by [START_REF] Koenker | Regression quantiles[END_REF] in a linear framework. Throughout this paper, we consider quantile regression in the more general nonparametric regression setting. Many methods have been proposed to estimate conditional quantiles in this context. In particular, nearest-neighbor estimators were investigated in [START_REF] Bhattacharya | Kernel and nearest-neighbor estimation of a conditional quantile[END_REF], and local linear and double-kernel approaches were considered in [START_REF] Yu | Local linear quantile regression[END_REF]. [START_REF] Koenker | Quantile smoothing splines[END_REF] and [START_REF] He | Bivariate quantile smoothing splines[END_REF] developed spline-based methods, while [START_REF] Lee | Efficient semiparametric estimation of a partially linear quantile regression model[END_REF], [START_REF] Horowitz | Nonparametric estimation of an additive quantile regression model[END_REF] or [START_REF] Belloni | Conditional quantile processes based on series or many regressors[END_REF] focused on global nonparametric series methods. Many other estimators have also been investigated; we refer to, e.g., [START_REF] Chaudhuri | Nonparametric estimates of regression quantiles and their local Bahadur representation[END_REF], [START_REF] Fan | Robust nonparametric function estimation[END_REF], [START_REF] Hallin | Local linear spatial quantile regression[END_REF], [START_REF] Gannoun | Reference curves based on non-parametric quantile regression[END_REF] or [START_REF] Yu | Quantile regression: applications and current research areas[END_REF], and to [START_REF] Koenker | Quantile regression[END_REF] for a review.

In this work, we consider the problem of estimating, in a nonparametric way, the conditional quantiles of a scalar response variable Y with respect to a d-dimensional covariate X. The theoretical properties of the method we will focus on were actually derived in a more stringent setup under which

Y = m(X, ε), (1.1)
where ε is a random error term that is independent of X. It is further assumed that both X and ε are absolutely continuous and that the function (x, z) → m(x, z) is of the form m 1 (x) + m 2 (x)z, where m 1 (•) : R d → R and m 2 (•) : R d → R + 0 are Lipschitz functions. In practice, both this link function m and the distribution of ε are unspecified, so that the conditional distribution of Y given X = x is unknown. Therefore it needs to be estimated by using a random sample (X 1 , Y 1 ), . . . , (X n , Y n ) from (1.1). All consistent estimators of conditional quantiles in this context require some localization or penalization technique; four of the most classical estimators, that are used as competitors to our own approach in the sequel, are described in Section 4.1.

The method we will focus was introduced in [START_REF] Charlier | Conditional quantile estimators through optimal quantization[END_REF] -hereafter, CPS15 -and performs localization through optimal quantization, a concept that projects the covariates X i , i = 1, . . . , n onto a finite grid of N quantizers. This projection of the X i 's proceeds according to their belongings to appropriate Voronoi cells that are partitioning the covariate space R d (see [START_REF] Pagès | A space quantization method for numerical integration[END_REF] for more details on this Voronoi partition). The CPS15 estimator can therefore be viewed as a nonparametric partitioning regression estimator. Partitioning ideas were already used in the construction of least square estimators; see, e.g., [START_REF] Györfi | A distribution-free theory of nonparametric regression[END_REF] for a complete overview and [START_REF] Kohler | Rates of convergence for partitioning and nearest neighbor regression estimates with unbounded data[END_REF] or [START_REF] Cattaneo | Optimal convergence rates, Bahadur representation, and asymptotic normality of partitioning estimators[END_REF] for more recent results. Now, the emphasis in CPS15 was almost exclusively on theoretical aspects; the main results there (i) are quantifying how well quantization allows to approximate conditional quantiles as N grows and (ii) are providing the convergence in probability of a sample quantization-based conditional quantile to its population version (see (2.3) and (2.4) below, respectively). Practical implementation, however, was barely touched in CPS15, and finite-sample performances were not investigated. The goal of the present work is to extensively cover these crucial aspects.

The outline of the paper is as follows. For the sake of completeness, we describe in Section 2 the concept of optimal quantization and the quantization-based estimators from CPS15. Then we discuss in Section 3 the choice of the number N of quantizers, that, in the nonparametric setup considered, plays the role of a smoothing parameter. In Section 4, we compare our estimator with well-known nonparametric estimators of conditional quantiles, namely (see Section 4.1) quantile regression splines, k nearest-neighbor (kNN) estimators, and (kernel) local linear/constant estimators. The comparison is based on visual inspection of estimated conditional quantile curves (Section 4.2) and on empirical integrated squared errors (Section 4.3). We provide a real data example in Section 5 and investigate, in Section 6, the case where the number d of covariates is larger than one. We conclude with a summary and some final comments in Section 7.

Conditional quantiles through optimal quantization

Quantization-based approximations of conditional quantiles

Let X be a random d-vector defined on a probability space (Ω, F, P ), such that the L 2 -norm X = (E[|X| 2 ]) 1/2 is finite (here, |•| denotes the Euclidean norm). For any fixed positive integer N , quantization replaces X with the discrete vector

X γ N = Proj γ N (X) obtained by projecting X onto an N -grid γ N of R d (so that γ N ∈ (R d ) N ). Optimal
quantization is achieved if the grid minimizes the quantization error X -X γ N . In the sequel, X N will denote the approximation obtained by projecting X onto an arbitrary optimal N -grid (unlike existence, unicity does not always hold; see [START_REF] Pagès | A space quantization method for numerical integration[END_REF]). Of course, it is expected that this approximation becomes arbitrarily precise as N increases to infinity, which is confirmed by the fact that

X -X N = O(N -1/d ) as N → ∞;
(2.1) see, e.g., [START_REF] Graf | Foundations of quantization for probability distributions[END_REF].

Denoting by I A the indicator function of the set A, consider the well-known check

function z → ρ α (z) = z(α -I [z<0] ) = -(1 -α)zI [z<0] + αzI [z≥0]
; see [START_REF] Koenker | Regression quantiles[END_REF]. In a regression context involving the scalar response Y and the (continuous) covariate d-vector X, this function allows to characterize the conditional α-quantile q α (x)

of Y given X = x as q α (x) = arg min a∈R E ρ α (Y -a)|X = x . In this context, CPS15
proposed approximating q α (x) by

q N α (x) = arg min a∈R E ρ α (Y -a)| X N = x , (2.2)
where x denotes the projection of x onto the optimal N -grid γ N defining X N ; here, N is a fixed integer whose choice in practice will be discussed below. CPS15 proved that, under mild regularity assumptions, there exists a constant C such that q N α (X) -q α (X) ≤ C X N -X 1/2 for N sufficiently large, which, in view of (2.1), implies that

q N α (X) -q α (X) = O(N -1/2d ) as N → ∞.
(2.3)

Since pointwise results are also appealing in quantile regression, they further showed that

q N α (x) → q α (x) as N → ∞, uniformly in x (in the support S X of X).

Quantization-based estimators of conditional quantiles

In real data applications, independent copies (X 1 , Y 1 ), . . . , (X n , Y n ) of (X, Y ) are available, and we have to define a sample version, q N,n α (x) say, of the quantization-based population quantile q N α (x). As described in CPS15, this can be achieved in the following two steps : 4 (S1) First, an optimal grid is obtained by making use of the CLVQ algorithm described below, with n iterations (based on ξ t = X t , t = 1, . . . , n) and an initial N -grid X 0 that results from sampling randomly among the X i 's without replacement under the constraint that an x-value cannot be picked more than once; ) is the optimal grid obtained in (S1), the proposed empirical conditional quantile q N,n α (x) is the sample α-quantile of the Y i 's whose corresponding X N,n i is equal to xN,n . Throughout, we will use the sample quantile concept that is associated with the default value 7 of the type argument in the R function quantile.

(S2) letting X N,n i = Proj γN,n (X i ), i = 1, . . . ,
It was shown in CPS15 that, for any fixed N and x, q N,n α (x)

P → q N α (x), as n → ∞. (2.4)
In other words, the empirical conditional quantile q N,n α (x) is weakly consistent for the -as we have seen, arbitrarily good -approximation q N α (x) of the true conditional quantile q α (x). In principle, the convergence q N α (x) -q α (x) → 0 as N → ∞ could be combined with (2.4) to provide an asymptotic result stating that, if N = N n goes to infinity in an appropriate way, then q Nn,n α (x) -q α (x) P → 0 as n → ∞. Proving formally such a result, however, is extremely difficult, since all convergence results for the CLVQ algorithm below are as n → ∞ with N fixed. In any case, this is of academic interest only; in practice, indeed, the sample size n is of course fixed and the user needs to choose appropriately the size N of the quantization grid. We later on propose a method that allows to choose this smoothing parameter N by achieving a bias-variance trade-off (see the discussion of Figure 1

below).

The so-called Competitive Learning Vector Quantization (CLVQ) algorithm, that provides an optimal grid in (S1) above, is a stochastic gradient algorithm that proceeds as follows. Let (ξ t ) t∈N0 , with N 0 = {1, 2, . . .}, be a sequence of mutually independent random d-vectors that are all distributed as the vector X subject to quantization, and let (δ t ) t∈N0 be a deterministic sequence in (0, 1) such that ∞ t=1 δ t = +∞ and ∞ t=1 δ 2 t < +∞. Starting from an initial N -grid X 0 = x 0 , with pairwise distinct components, the CLVQ algorithm recursively defines the grid X t , t ≥ 1, as

X t = X t-1 - δ t 2 ∇ x d 2 N (X t-1 , ξ t ), (2.5)
where ∇ x d 2 N (x, ξ) stands for the gradient with respect to the x-argument of the so-called local quantization error

d 2 N (x, ξ) = min 1≤i≤N |x i -ξ| 2 , with x = (x 1 , . . . , x N ) ∈ (R d ) N and ξ ∈ R d . For any ξ, the ith component of this gradient is ∇ x d 2 N (x, ξ) i = 2 (x i -ξ) I [xi=Proj x (ξ)] , i = 1, . . . , N.
The grids X t-1 and X t thus only differ by one point, namely the one corresponding to the non-zero component of the gradient above; see, e.g., [START_REF] Pagès | A space quantization method for numerical integration[END_REF] or [START_REF] Pagès | Optimal quadratic quantization for numerics: the Gaussian case[END_REF] for details.

In the sequel, we will write throughout that the grids provided by the CLVQ algorithm are optimal. Strictly speaking, though, the grid γN,n resulting from this algorithm, for any fixed n, does not minimize the quantization error, hence is not optimal. Yet it is shown in [START_REF] Pagès | A space quantization method for numerical integration[END_REF] that, under mild assumptions, γN,n converges to an optimal grid as the number of iterations (here, the sample size n) tends to infinity (with N fixed).

Clearly, only moderate-to-large values of n are expected to provide good approximations of genuine optimal N -grids.

We end this section by quickly evaluating the performances of the CLVQ algorithm, both for d = 1 and d = 2. To do so, we generated, for the various combinations of n and N considered in

Bootstrapped quantization-based estimators of conditional quantiles

As explained in CPS15, it is more appropriate in practice to use a bootstrap version of the estimator q N,n α (x), particularly so for small to moderate values of n. Obtaining bootstrapped quantization-based estimators can be achieved in the following two steps : The resulting bootstrap estimator qN,n

α,B (x) is qN,n α,B (x) = 1 B B b=1 q(b) α (x).
(2.7)

In the sequel, the notation B = 1 means that no bootstrap is performed.

Figure 1 provides an illustration of the original estimator q N,n α (x) and its bootstrap version qN,n α,B (x) (with B = 50). Based on a random sample of size n = 500 from the model

Y = X 2 + ε, (2.8)
where X = 4Z -2 (with Z ∼ Beta(0.3, 0.3)) and ε ∼ N (0, 1) are independent, we evaluated both q N,n α (x) and qN,n α,B (x) for N = 10, 25, and 50 (actually, these were evaluated only at 300 equispaced points between the minimum and maximum values of X 1 , . . . , X n ).

Clearly, the bootstrap significantly smooths all curves. This happens provided that B is -2 5 (green), 0.75 (orange), and0.95 (red); see (2.8) for the data generating model. not too small. But of course, B should not be taken very large, as it would increase the computational burden (in the sequel, we use B = 50 throughout).

-1 0 1 2 -2 0 2 4 6 X Y (a) N = 10, B = 1 -2 -1 0 1 2 -2 0 2 4 6 X Y (b) N = 10, B = 50 -2 -1 0 1 2 -2 0 2 4 6 X Y (c) N = 25, B = 1 -2 -1 0 1 2 -2 0 2 4 6 X Y (d) N = 25, B = 50 -2 -1 0 1 2 -2 0 2 4 6 X Y (e) N = 50, B = 1 -2 -1 0 1 2 -2 0 2 4 6 X Y (f) N = 50, B = 50
Figure 1 also shows that the number N of quantizers has an important impact on the resulting quantile curves. For small N (N = 10), the curves show a large bias and a small variability. For larger N (N = 50), the bias is smaller but the variability is larger (with or without bootstrap). This bias-variance trade-off is standard in the selection of a smoothing parameter. The intermediate value N = 25, that achieves a nice balance between those two extreme cases, was actually obtained from an original N -selection method we describe in Section 3 below. This method to choose N , along with the computation of the (bootstrap or non-bootstrap) quantization-based estimators of conditional quantiles under investigation, are entirely implemented in a package R called QuantifQuantile, that is available on the CRAN.

Data-driven selection of N

In this section, we first define (Section 3.1) a natural integrated square error (ISE) quantity that is essentially convex in N , hence allows to identify an optimal value N opt for N . ISEs involve the true unknown conditional quantiles, so that N opt cannot be obtained from the data. We therefore propose (Section 3.2) a data-driven selection method for N obtained by replacing the true ISEs with bootstrap versions. We then investigate how close the resulting Nopt is to N opt .

Infeasible selection of N

For any fixed N , we can consider the ISEs (recall that S X denotes the support of X)

S X q N,n α (x) -q α (x) 2 dx and S X qN,n α,B (x) -q α (x) 2 dx,
associated with the non-bootstrap and bootstrap estimators, respectively. These can be approximated by

ISE α,J (N ) = 1 J J j=1 q N,n α (x j )-q α (x j ) 2 and ISE ᾱ,B,J (N ) = 1 J J j=1 qN,n α,B (x j )-q α (x j ) 2 ,
where x 1 , . . . , x J are equispaced between the minimum and maximum values of X 1 , . . . , X n .

It is of course natural to consider optimal a value of N that minimizes these ISEs, which leads to considering

N α,J;opt = arg min N ∈N0
ISE α,J (N ) and N ᾱ,B,J;opt = arg min

N ∈N0 ISE ᾱ,B,J (N ).
These optimal N -values may depend on α, which explains the notation. The dependence on J of these ISEs and of optimal N -values is way less crucial than their dependence on α and B; accordingly, we simply write ISE α(N ), ISE ᾱ,B (N ), N α;opt and N ᾱ,B;opt in the sequel.

To illustrate these definitions, we generated mutually independent random samples of size n = 300 according to the models

(M1) Y = 1 5 X 3 1 + ε, (M2) Y = f (X 2 ) + ε , where X 1 = 6Z 1 -3 (with Z 1 ∼ Beta(0.3, 0.3)), X 2 = 3Z 2 -1.5 (with Z 2 ∼ Beta(2, 2)),
ε ∼ N (0, 1), and ε ∼ χ 2 2 are mutually independent. Denoting the standard normal density as ϕ, the link function f is defined as f

(x) = 1 2 ϕ(x) + 10 3 =1 (0.2) ϕ(10(x - 2 + 1))
, a choice that is inspired by [START_REF] Marron | Exact mean integrated squared error[END_REF]. Obviously, q α (x) =

x 3 5 + Φ -1 (α) for (M1) and q α (x) = f (x) + Ψ -1 2 (α) for (M2), where Φ and Ψ 2 denote the cumulative distribution functions of the N (0, 1) distribution and of the χ 2 2 distribution, respectively. We evaluated the ISEs above with J = 300 points. Figure 2 (with B = 50) for a single sample, while the bottom panels plot the corresponding graphs averaged (pointwise in x) over m = 500 mutually independent samples. For the sake of brevity, Figure 3, that is related to Model (M2), only reports the corresponding plots for bootstrapped estimators (still with B = 50).

Figures 23show that ISE curves are indeed essentially convex in N : while random variability prevents perfect convexity for a single sample, the smoothest case based on 500 replications with qN,n α,B leads to strictly convex I(M)SE curves. Besides, these curves are quite flat in a neighborhood of N α;opt or N ᾱ,B;opt , so that picking a value of N that is close to (but not exactly equal to) the optimal N -value will not have an important impact on the resulting ISE values. It seems that, as expected, N 1-α;opt = N α;opt and N 1-α,B;opt = N ᾱ,B;opt when the error ε has a symmetric distribution. In contrast, asymmetric errors (Figure 3) lead to optimal values of N that are not symmetric in α and depend more strongly of α.

To conclude this section, we stress that the existence of such an optimal value of N does not contradict the theoretical results from Section 2. One might indeed have guessed from (2.3) that the ISE functions above would be monotone decreasing (rather than convex) in N . The result in (2.3), however, (i) involves the population quantile approximation q N α (x) and not its sample counterpart q N,n α (x), and (ii) requires that the projection, for any N , is performed on an optimal quantization grid. As pointed out earlier, the CLVQ algorithm provides a good approximation of an optimal grid only if the number of iterations, that is equal to the sample size n, is large compared to N . Consequently, in the present setup where the sample size n is fixed, increasingly large values of N will result into CLVQ grids that are less and less optimal, which explains the increasing ISE values for such values of N .

Data-driven selection of N

In practice, the population conditional quantile function q α (•) is unknown, so that it is impossible to obtain optimal N -values by minimizing ISE α(N ) and ISE ᾱ,B (N ) as above. We then propose the following approach that consists in replacing the unknown quantile q α (x) by a bootstrap estimate (see [START_REF] Efron | The jackknife, the bootstrap and other resampling plans[END_REF] for more details on such methods).

Consider a grid {x 1 = X (1) , x 2 , . . . , x J-1 , x J = X (n) } of equispaced points between the minimal and maximal observed X-values. We generate B bootstrap samples of size n from X 1 , . . . , X n and use them as stimuli (that is, as the ξ t 's) in the CLVQ algorithm, which provides B resulting quantization grids. We then consider where the cardinality of N (⊂ N 0 ) is finite (and may be chosen as a function of n). As for their infeasible counterparts, we will not stress in the sequel the dependence on J in these sample ISEs and optimal N -values, nor the dependence on B, that we choose equal to 30 throughout. to over-estimate the optimal N -values, with a bias of about 10. Note that this bias is not alarming since Figure 2 shows that virtually the same performances in terms of ISE are obtained with or without bias (particularly so for intermediate values of α). Finally, we point out that, when considering a single sample (instead of 500 replications), the estimated optimal N -value stays in the same interval (the plots are not provided here).

ISE α, B,J (N ) = 1 J J j=1 1 B B b=1 q N,n α (x j ) -q( b) α (x j ) 2 and ISE ᾱ,B, B,J (N ) = 1 J J j=1 1 B B b=1 qN,n α,B (x j ) -q( b) α (x j ) 2 , ( 3
To sum up, the above simulation results suggest that the sample ISE functions allow us to obtain a data-driven value of N that, for all practical purposes, is as satisfactory as the infeasible optimal N -value based on the original ISE functions. Note that, in a nonparametric framework, the smoothing parameter is often rather selected by minimizing the mean square error of the estimator at hand. In the context of this paper, however, deriving explicit expressions of the bias and variance (or standard errors) of quantization-based quantile estimators turns out to be extremely challenging. The main reason is that the randomness of observations appears twice in the definition of these estimators : (i) when obtaining the quantization grid (through the CLVQ algorithm, on which few theoretical results are available) and (ii) when calculating the empirical quantile of the selected observations. This explains why we need to rely on the selection procedure described in this section.

Comparison with some classical competitors

The aim of this section is to investigate how the proposed estimator of conditional quantile behaves compared to its main competitors, that are introduced in Section 4.1.

We both compare the estimated quantile curves themselves (Section 4.2) and the respective ISEs, averaged over m = 500 mutually independent random samples (Section 4.3).

The competitors considered

The first competing estimator we consider is associated with quantile regression spline methods; see [START_REF] Koenker | Quantile smoothing splines[END_REF]. More precisely, this estimator, q sp α , say, is defined as

q sp α = q sp α,λ = arg min g∈G n i=1 ρ α (Y i -g(X i )) + λP (g) ,
where λ is a nonnegative real number, G is an appropriately chosen space of functions, and the penalty P (g) = 1 0 |g (x)| dx is the total variation of the first derivative of g.

It is well-known that q sp α is always piecewise linear. As usual in penalization methods, λ governs the trade-off between fidelity to the sample and smoothness of the resulting estimator. Several methods to select an appropriate value of λ are available in the literature. Below, we will use the AIC criterion, which selects the value λ = λα;opt

minimizing AIC(λ) = log 1 n i ρ α Y i -ĝλ (X i ) + p λ n ,
where p λ is the number of knots defining the spline. The reader can refer to Koenker et al.

(1994) and [START_REF] Koenker | Penalized triograms: total variation regularization for bivariate smoothing[END_REF] for more details on how to choose λ (or G). When implementing this method, we used the rqss function from the R package quantreg to compute q sp α and performed the λ-selection via the R function AIC (the package quantreg does not propose an automatic λ-selection procedure).

Another competitor is the kNN estimator, which is defined as follows; see Bhattacharya and [START_REF] Bhattacharya | Kernel and nearest-neighbor estimation of a conditional quantile[END_REF]. For any integer k ∈ {1, . . . , n}, this estimator, qk α;kNN (x) say, is defined as the α-quantile of Y 1 , . . . , Y k after permuting the observations (X i , Y i ), i = 1, . . . , n in such a way that |X i -x| is increasing in i. In other words, qk α;kNN (x) is the α-quantile of the Y i 's associated with the k X i 's that are closest to x. There does not seem to exist an efficient and well-accepted data-driven method to choose k in the literature. In the sequel, we throughout considered the kNN estimator qα;kNN (x) = qkα;opt α;kNN (x), with

k α;opt = arg min k∈K 1 J J j=1 q k α;kNN (x j ) -q α (x j ) 2 ,
where {x 1 , . . . , x J } is chosen as in Section 3.2 and K denotes some set of possible values for k. Clearly, k cannot be chosen this way in practice since the theoretical quantiles q α (x j ) are unknown. To the best of our knowledge, no R package allows to compute q α;kNN (x) and we therefore wrote our own implementation to conduct the simulations below.

The last estimators we consider are the kernel (local linear or local constant) estimators introduced in Yu and Jones (1998). The local linear estimator is of the form q α;YJ (x) = â, with

(â, b) = arg min (a,b)∈R×R n i=1 ρ α Y i -a -b(X i -x) K X i -x h ,
where K is a kernel function and h is the bandwidth. In the sequel, K will be the standard normal density, and we choose

h = ĥα;opt = α(1 -α) ϕ(Φ -1 (α)) 2 h mean ,
where ϕ and Φ are respectively the standard normal density and distribution functions, and where h mean is the optimal choice of h for mean regression, selected through crossvalidation; see [START_REF] Yu | Local linear quantile regression[END_REF]. The local constant version of this estimator is defined as q α;YJc (x) = â, with

â = arg min a∈R n i=1 ρ α Y i -a K X i -x h ,
where K and h will throughout be chosen as for the local linear estimator. As for q α;kNN (x), the simulations below are based on our own R implementation of q YJ α (x) and q YJc α (x).

Comparison of estimated quantile curves

We now compare our quantization-based estimators with their competitors described above. Since we saw in Section 2.3 that the bootstrap estimators qN,n α,B are to be favored over their original versions q N,n α , we restrict to qN,n α,B below, with the corresponding datadriven value of N , namely Nᾱ ,B;opt , that was proposed in Section 3.2. In this section, we also do not consider the local constant estimator q YJc α (x) since the results in Section 4.3 below show that it is usually outperformed by its local linear version in terms of ISEs.

We start the comparison by investigating estimated quantile curves computed from n independent observations generated according to the models

(M1) Y = 1 5 X 3 1 + ε, (M2) Y = f (X 2 ) + ε , (M3) Y = sin(X 3 ) + 0.5 + 1.5 sin 2 ( π 2 X 3 ) ε,
where

X 1 = 6Z 1 -3 (with Z 1 ∼ Beta(0.3, 0.3)), X 2 = 3Z 2 -1.5 (with Z 2 ∼ Beta(2, 2)), X 3 = 6Z 2 -3, ε ∼ N (0, 1)
, and ε ∼ χ 2 2 are mutually independent. The link function f is the same as in Section 3.1. For n = 300, the resulting estimated quantile curves are plotted in Figures 567for Models (M1)-(M3), respectively. For the sake of comparison, Figure 8 provides the results for Model (M1) and n = 1, 000.

The quantization-based quantile curves are smooth and adapt well to the polynomial or more complex nature of the link function at hand. In contrast, while the piecewise linear curves obtained from the spline estimator are beneficial for polynomial link functions (see Figure 5), they hurt for more complex link functions (Figure 6, e.g., shows that the spline-based curves miss some of the underlying bumps). The curves resulting from kNN estimation are relatively close to the theoretical ones in each model, but show some peaks or constant pieces, hence are less pleasant from a visual point of view. Eventually, the curves associated to q YJ α (x) lack smoothness for medium values of the covariate in Figure 5 (further numerical experiments, not reported here, reveal that this follows from the non-uniform covariate distribution). Also, the local linear estimator q YJ α (x) does not catch correctly the link function in Figure 7.

The values of the various tuning parameters are given in Table 2. For the sake of brevity, we only comment on the results for Model (M2). We see that the optimal value for N decreases with α, which is due to the chi-square error term, that implies that estimation of conditional α-quantiles should be based on increasingly many data points as α grows. In accordance with this, the optimal number of neighbours k α;opt in kNN estimation increases with α. Such a monotonic pattern is not observed for the bandwidth ĥα;opt used in kernel estimators since it is chosen such a way that ĥα;opt = ĥ1-α;opt . Finally, the parameter λα;opt , that is a monotone decreasing function of the number of knots of the splines, is related to α in a quite unclear fashion.

Comparison of the ISEs

Obtaining well-behaved curves is of course desirable, particularly so in applied statistics, but this should not be achieved at the expense of efficiency. That is why we now compare the various estimators in terms of ISEs. To do so, we generated 500 independent samples from Models (M1)-(M3) with sample sizes n = 300 and n = 1, 000. In each case, we evaluated the ISEs corresponding to the quantization-based estimators and to their four competitors defined in Section 4.1. For each model, sample size, and quantile order α considered, this provides a series of 500 observed ISEs for each estimator. -3 The lighter curves correspond to population conditional quantiles.

-2 -1 0 1 2 3 -5 0 5 X Y (a) qN,n α,B -3 -2 -1 0 1 2 3 -5 0 5 X Y (b) q sp α -3 -2 -1 0 1 2 3 -5 0 5 X Y (c) q α;kNN -3 -2 -1 0 1 2 3 -5 0 5 X Y (d) q YJ α
-1.0

-0.5 0.0 0.5 1.0 1.5 0 2 4 6 8 10 X Y (a) qN,n α,B -1.0 -0.5 0.0 0.5 1.0 1.5 0 2 4 6 8 10 X Y (b) q sp α -1.0 -0.5 0.0 0.5 1.0 1.5 0 2 4 6 8 10 X Y (c) q α;kNN -1.0 -0.5 0.0 0.5 1.0 1.5 0 2 4 6 8 10 X Y (d) q YJ α Figure 6:
The same plots as in Figure 5, but for Model (M2) and n = 300.

20

-3

-2 -1 0 1 2 3 -4 -2 0 2 4 X Y (a) qN,n α,B -3 -2 -1 0 1 2 3 -4 -2 0 2 4 X Y (b) q sp α -3 -2 -1 0 1 2 3 -4 -2 0 2 4 X Y (c) q α;kNN -3 -2 -1 0 1 2 3 -4 -2 0 2 4 X Y (d) q YJ α Figure 7:
The same plots as in Figures 56, but for Model (M3) and n = 300.

21

-3 ISEs are, as expected, smaller for n = 1, 000 than for n = 300, but otherwise are quite similar.

-2 -1 0 1 2 3 -5 0 5 X Y (a) qN,n α,B -3 -2 -1 0 1 2 3 -5 0 5 X Y (b) q sp α -3 -2 -1 0 1 2 3 -5 0 5 X Y (c) q α;kNN -3 -2 -1 0 1 2 3 -5 0 5 X Y (d) q YJ α
For most models and quantile orders considered, the quantization-based estimators, and n = 1, 000 (bottom). The estimators considered are the quantization-based estimator qN,n α,B (in blue), the spline estimator q sp α (in orange), the kNN estimator estimator q α;kNN (in green), the local linear estimator q YJ α (in purple), and the local constant estimator q YJc α (in red).

Section 4.2 (note, however, that an exception also appears in Figure 11(b)).

Since the computational burden is also an important issue, we gather in Table 3 the computation times (in seconds) used by each estimator to produce Figure 9(a). To study computation times in the light of efficiency, we further report there, for each estimator, a global measure of efficiency (Eff), defined as the sum, over the five α-values considered in that figure, of the medians of the 500 fixed-α observed ISEs. Since the choice of the grid N plays an important role in the selection of the optimal N -value (at least from a computational point of view), we considered two grids, namely N 1 = {5, 6, . . . , 29, 30}

and N 5 = {5, 10, . . . , 25, 30} (the latter one, that may seem too coarse, is actually the one that led to the good ISE performances in Figures 9). Table 3 confirms that there is no free lunch, as it shows that the gain in terms of efficiency (for Model (M1) and n = 300) has a price in terms of computation time. This price, however, is quite reasonable. Moreover, and N 5 actually differ across models but not their size.

A real data example

In this section, we illustrate the behavior of qN,n α,B on a real data set, that was kindly sent to us by Dr. Keming Yu. This data set, of size n = 298, corresponds to the serum concentration, in grams per liter, of immunoglobulin-G in children aged from 6 months to 6 years. This data set was already investigated in [START_REF] Yu | Local linear quantile regression[END_REF], where a kernel local linear estimator and a double kernel estimator were considered; we therefore refer to that paper for (smoothed versions of) the corresponding quantile curves. Also, we do not plot the quantile curves resulting from the kNN estimator q α;kNN since the selection of the parameter k is not data-driven, hence cannot be achieved on a real data set. Consequently, we only compare the proposed quantization-based quantile curves with those obtained from their main competitor, namely the spline estimator q sp α . We now describe how we obtained quantization-based quantile curves in the present context, for each α = 0.05, 0.1, 0.25, 0.5, 0.75, 0.9 and 0.95.

(i) The first step consists in choosing the optimal number N of quantizers. To do so, we adopted the method proposed in Section 3. The resulting quantile curves are plotted in Figure 13(a). As announced, we compare these curves with those associated with the spline estimator q sp α . For each α considered, the parameter λ was selected according to the AIC procedure described in Section 4.1,

Summary and final comments

In this paper, we investigated the empirical performances of the quantization-based estimators of conditional quantile that were recently introduced in CPS15. This led us to design an efficient data-driven method to select the number N of quantizers to be used, that plays the role of a smoothing parameter. From extensive simulations, we concluded that quantization-based estimators compete well with alternative conditional quantile estimators and sometimes dominate its main competitor, spline estimators, particularly when the link function is complex. We treated a real data example, in which we showed that the proposed methodology provides very satisfactory conditional quantile curves, and we saw that the good properties of quantization-based estimators extend to the bivariate covariate case. This should make conditional quantile estimation based on quantization of interest to practitioners; in this spirit, we wrote an R package, named

QuantifQuantile, that allows to compute in a straightforward way the proposed estimators (the data-driven selection of N is included) and to plot the resulting quantile curves/surfaces. This package is already available on the CRAN.

We conclude with two final comments.

(i) Since conditional quantiles q α (x) are monotone increasing in α for any fixed x, population conditional quantile hypersurfaces cannot cross. Yet, like many of their competitors, the proposed quantization-based quantile hypersurfaces might cross. No such crossings occurred in the illustrations of this paper, nor in all other numerical experiments we have conducted. Practitioners, however, usually do not like such crossings, that violate their interpretation of quantiles. An easy way to exclude crossings is to select an optimal N -value that is independent of α. A natural way to do so is to minimize

N → ISE B (N ) = AVE α ISE ᾱ,B (N ) , (7.1)
where the average extends over all α-values considered. Figure 16 plots the resulting ISE curve for the real data example from Section 5, which leads to the optimal value N = 5.

Another procedure that allows to avoid crossings is the rearrangement procedure from [START_REF] Chernozhukov | Quantile and probability curves without crossings[END_REF].

(ii) Both in CPS15 and in this work, quantization was applied to the covariate only.

One may wonder whether or not it may be of interest to move entirely to a discrete ) is the optimal value for N X at the given N Y (i.e., the one that minimizes the ISE function N → ISE ᾱ,B (N, N Y ). It is seen that, irrespective of N Y and α, the optimal number of quantizers for X stays close to 25. However, there is no optimal value for N Y : the larger N Y , the smaller the averaged ISE. In conclusion, an optimal choice for (N X , N Y ) appears to be N X close to 25 and N Y → ∞. Since letting the size of the quantization grid go to infinity is equivalent to not quantizing at all, this leads to favouring the solution adopted in CPS15 and in the present work over double quantization. Quantizing X only also makes perfect sense in the problem of estimating nonparametrically a given conditional quantile q α (x). Most estimation methods indeed consist in selecting the observations whose X-part is the closest to x (where the meaning of "close" depends on the method at hand) and to take the sample quantile of the corresponding Y 's. Quantizing Y thus does not seem natural, as it adds an unnecessary approximation error.

(

  S1) We first generate B(∈ N 0 ) samples of size n from the original sample with replacement. Performing (S1) on each of them leads to B optimal grids; ( S2) each grid obtained in ( S1) then yields, through (S2), an estimator q. . . , B, obtained from the projection of the original sample {(X i , Y i ), i = 1, . . . , n} on the corresponding grid.

Figure 1 :

 1 Figure1: Estimated conditional quantile curves x → q N,n α (x) (left) and their bootstrap counterparts x → qN,n α,B (x) for B = 50 (right), based on N = 10 (top), N = 25 (middle), and N = 50 (bottom). In all cases, the sample size is n = 500, and the quantile orders considered are α = 0.05 (red), 0.25 (orange), 0.5 (green), 0.75 (orange), and 0.95 (red); see (2.8) for the data generating model.

  focuses on Model (M1). Its top panels plot the graphs of N → ISE α(N ) and N → ISE ᾱ,B (N )

Figure 2 :Figure 3 :

 23 Figure 2: Plots of the mappings N → ISE α(N ) (top left) and N → ISE ᾱ,B (N ) with B = 50 (top right)for a random sample of size n = 300 from Model (M1). The bottom panels report the corresponding plots obtained by averaging these mappings over m = 500 mutually independent random samples (all mappings are actually only evaluated atN = 5, 6, 7, . . . , 29, 30, 35, 40, . . . , 150).

  .1) where q( b) α (x j ), b = 1, . . . , B, stands for the estimator from (S2) in Page 5 computed from the original sample (X 1 , Y 1 ), . . . , (X n , Y n ) and the bth grid obtained above. Evaluating ISE ᾱ,B, B,J (N ) thus requires generating B + B bootstrap samples of size n : B for the construction of qN,n α,B (x j ), and B to obtain q( b) α (x j ), b = 1, . . . , B. These sample ISEs are to be minimized in N . Since not all values of N can be considered in practice, we rather consider Nα , B,J;opt = arg min N ∈N ISE α, B,J (N ) and Nᾱ ,B, B,J;opt = arg min N ∈N ISE ᾱ,B, B,J (N ),

Figure 4 Figure 4 :

 44 Figure4plots the mappings N → ISE α(N ) and N → ISE ᾱ,B (N ) and, for the sake of comparison, the (infeasible) mappings N → ISE α(N ) and N → ISE ᾱ,B (N ), in the setup of Model (M1) with sample size n = 300 (more precisely, the average of the corresponding plots, over 500 mutually independent random samples, are showed there). It is seen that, for large N , the (averaged) sample ISE functions quite poorly estimate their theoretical versions. Our primary interest, however, rather is in the agreement of the corresponding argmins, which are associated with the feasible or infeasible optimal N -values. In that respect, the results show that the values of N minimizing the sample ISE functions tend 13

  Figures 9-11 draw the boxplots of those 500 ISE values for Models (M1)-(M3), respectively.

Figure 5 :

 5 Figure 5: For a random sample of size n = 300 generated according to Model (M1), conditional quantile curves obtained from quantization-based estimation (upper left), spline methods (upper right),nearest-neighbour estimation (lower left), and local linear kernel methods (lower right). In all cases, the quantile levels considered are α = 0.05 (red), 0.25 (orange), 0.5 (green), 0.75 (orange), and 0.95 (red).

Figure 8 :

 8 Figure 8: The same plots as in Figures 5-7, but for Model (M1) and n = 1, 000.

Figure 9 :

 9 Figure 9: Boxplots, for α = 0.05, 0.25, 0.50, 0.75, and 0.95, of the ISEs of various conditional α-quantile estimators obtained from 500 independent random samples from Model (M1), with size n = 300 (top)

Figure 10 :

 10 Figure10: The same boxplots as in Figure9, but for Model (M2).

Figure 11 :

 11 Figure 11: The same boxplots as in Figures 9-10, but for Model (M3).

Figure 12 :

 12 Figure 12: In the context of the real data example of Section 5, this figure plots, for various α-values, N → ISE ᾱ,B (N ) over N = {5, 10, . . . , 35, 40} (left) and over N = {4, 5, . . . , 9, 10} (right).

  2, based on the minimization of ISE ᾱ,B (N ) (see (3.1)), for which we considered throughout B = 50, B = 30, and a grid of J = 300 equispaced points between the minimum and maximum values of the X-part of the sample. We first evaluated ISE ᾱ,B (N ) for all N ∈ N = {5, 10, . . . , 35, 40}. Figure 12(a) plots the resulting sample ISE curves for each α. Irrespective of α, these curves are monotone increasing, which meansthat the tested values of N are too large. We therefore did the same exercise for N = {4, 5, . . . , 9, 10}, which led to Figure12(b). For each α, this provides an optimal N (that is equal to 5 or 6).(ii) The second step of course consists in obtaining the estimated quantile curves themselves, based on the selected values of N above. These are the plots of the mappings x → qN,n α,B (x) in (2.7), where we chose B = 50.

Figure 13 :Figure 14 :Figure 15 :

 131415 Figure 13: Conditional α-quantile curves for the real-data example of Section 5 and α = 0.05, 0.1, 0.25, 0.75, 0.9 and 0.95. The estimators considered are the quantization-based estimator qN,n α,B with B = 50 (left) and the spline estimator q sp α (right).

Figure 16 :=

 16 Figure 16: Plot of the mapping N → ISE B (N ) in (7.1) for the real data example of Section 5.

Figure 17 :

 17 Figure 17: Plots of the mapping (N X , N Y ) → ISE ᾱ,B (N X , N Y ) (for the estimator qY,N X ,N Y ,n α,B

  n, and xN,n = Proj γN,n (x), where γN,n =

	(x N,n 1 , . . . , xN,n N

Table 1

 1

	, 50 mutually independent random samples of size n from the
	uniform distribution over [-2, 2] d and performed in each case the CLVQ algorithm with
	the corresponding value of N . For each (n, N ), this leads to 50 CPU times to completion
	and 50 sample quantization errors
	e =	1 n	n i=1	X i -Proj γN,n (X i )	2 .	(2.6)
	Table 1 reports the corresponding averaged CPU times, along with the averages and
	standard deviations of the various collections of 50 quantization errors. Since the number
	of iterations of the algorithm is equal to n, it is no surprise that the error decreases with n
	while the computation time increases with n.

Table 2 :

 2 Values of the various tuning parameters involved in the conditional quantile estimators consid-

		,B;opt	20	23	23	18	18
	(M1)	λα;opt	0.595	0.971	1.075	0.638	1.033
	n = 300	k α;opt	25	32	30	28	27
		ĥα;opt	0.202	0.170	0.164	0.170	0.202
		Nᾱ ,B;opt	23	23	19	16	15
	(M2)	λα;opt	0.556	0.672	0.632	0.504	0.716
	n = 300	k α;opt	16	38	42	43	55
		ĥα;opt	0.135	0.113	0.109	0.113	0.135
		Nᾱ ,B;opt	13	10	9	10	9
	(M3)	λα;opt	0.518	1.127	1.490	0.840	0.516
	n = 300	k α;opt	70	61	38	28	34
		ĥα;opt	0594	0.498	0.482	0.498	0.594
		Nᾱ ,B;opt	32	43	43	41	33
	(M1)	λα;opt	1.122	1.117	0.927	1.120	1.208
	n = 1, 000	k α;opt	72	54	84	74	94
		ĥα;opt	0.148	0.124	0.120	0.124	0.148
	ered in Figures 5-8.

Table 3 :

 3 (First line:) Computation times (CPU, in seconds) used by each estimator to obtain Fig-Mac and Linux versions of our R package QuantifQuantile offer an option for parallel computing, which divides the computing times shown in Table3by a factor of 4, which clearly makes our estimators very competitive in this respect. In addition, the computation times for the nearest-neighbour estimator should be considered with care, since kNN procedures here are not based on a (typically always quite computationally intensive) data-driven selection of smoothing parameters. This table also gathers the efficiency measures Eff associated with Figures10(a) and 11(a) (we do not provide the corresponding CPU times, that barely depend on the model considered). The grids N 1

	ure 9(a). (Subsequent lines:) ISE-based global efficiency measures (Eff) associated with Figures 9(a),
	10(a), and 11(a); see Section 4.3 for details.

the

For the interpretation of the references to color in some of the figures, the electronic version of the paper should be consulted.

Acknowledgments

The authors are grateful to an Associate Editor and two anonymous referees for their careful reading and insightful comments that led to a significant improvement of the original manuscript. The first author's research is supported by a Bourse FRIA of the Fonds National de la Recherche Scientifique, Communauté française de Belgique.

The second author's research is supported by an A.R.C. contract from the Communauté Française de Belgique and by the IAP research network grant P7/06 of the Belgian government (Belgian Science Policy).

regression function are avoided, d should be small enough to avoid the usual curse of dimensionality.

For illustration purposes, we focus on the bivariate case d = 2 to be able to plot the resulting conditional quantile hypersurfaces. We first generated a random sample of size n = 3, 000 from the model

where Z 1 , Z 2 and ε are mutually independent. Conditional quantile hypersurfaces of this model for α =0.05, 0.25, 0.5, 0.75, and 0.95 are plotted in Figure 14(a). Figure 14(b) provides the corresponding sample quantile surfaces obtained from our estimator qN,n α,B , with B = 50 and the optimal N -values that are obtained through our data-driven Nselection procedure for the various values of α. The quantile surfaces were readily obtained from the R package QuantifQuantile, that can therefore also be used for d > 1.

Clearly, qN,n

α,B provides very smooth surfaces that properly catch the link function, even though the amplitude, for extreme values of α, is a bit under-estimated.

The results of Section 4 suggest restricting the comparison to spline-based estimators.

For this purpose, we consider the triogram-based bivariate splines from [START_REF] Koenker | Penalized triograms: total variation regularization for bivariate smoothing[END_REF], that is implemented in the rqss function of the R package quantreg. The resulting estimated quantile surfaces, based on the same λ-selection method as for d = 1 in Section 4, are plotted in Figure 14(c). Clearly, these quantile surfaces exhibit more variability than the quantization-based ones. Incidentally, note that, unlike the proposed estimator, this spline-based estimator does not extend easily to d > 2.

Parallel to our investigation of the case d = 1, we complement the results above with some efficiency results in terms of ISE. To do so, we generated 50 independent random samples from Model (M4) above with sample size n = 3, 000 and computed, for each sample and various values of α, the quantization-based and spline-based conditional quantile estimators. The corresponding 50 resulting ISEs for each estimator and every α considered, that are defined in a similar way as in Section 4.3, are represented in Figure 15.

In line with the d = 1 results in Figure 11, this shows that, for the complex link function considered, qN,n α,B dominates q sp α in terms of efficiency. Note that this is particularly true for extreme values of α.


G Brèthes

O Allain

A Dervieux

A Mesh-Adaptive Metric-Based Full-Multigrid for the Poisson problem

Keywords: anisotropic mesh adaption, full multi-grid, finite element, stopping criterion, Poisson problem

This paper studies the combination of the Full-Multi-Grid (FMG) algorithm with an anisotropic metric-based mesh adaption algorithm. For the sake of simplicity, the case of an elliptic two-dimentional Partial Differential Equation (PDE) is studied. Meshes are unstructured and non-embedded, defined through the metricbased parametrisation. A rather classical MG preconditionner is applied, in combination with a quasi-Newton fixed point. An anisotropic metric-based mesh adaptation loop is introduced inside the FMG algorithm. FMG convergence stopping test is re-visited. Applications to a few 2D continuous and discontinuouscoefficient elliptic model problems show the efficiency of this combination.

INTRODUCTION

Multi-grid methods (MG) can produce fast and robust solution algorithms. They apply to a large variety of models and approximations in Computational Mechanics.

MG uses intensively approximation properties. A simple local iteration is applied on the given grid. The iterative convergence is then accelerated by means of a set of coarser-grid corrections, ranging typically from a just twice coarser level, to a coarsest level with just a few dozens degrees of freedom. The approach generally produces an iterative convergence which, when expressed in terms of (logarithm of) residual norm decreasing with iterations, is more or less of constant slope. Further this slope does not depend on mesh size. In particular, an accurate enough discrete solution with N degrees of freedom is obtained with a number of operations bounded by K.N.Log(N), a complexity nearly optimal. Lastly, MG can be combined with a nested iteration producing the Full-Multi-grid (FMG) algorithm: FMG involves n phases, working from the 1-st coarsest mesh to the nth finest mesh. The j-th phase of FMG solves the appproximate PDE on the j-th mesh. This j-th phase starts from an interpolation of the result of j -1-th mesh and applies a certain number k j of MG cycles with the available j coarser meshes, from 1-st to j-th. In an ideal case, the number k j of cycles in each phase is the same. Then FMG has then an optimal complexity of K.N, predicted by theory and observed on many practical examples, see the reference book [START_REF] Trottenberg | Multigrid[END_REF].

But this rosy picture needs some rectifications. First, in many cases, the single local iteration of MG is frequently not sufficient to deal with singular or stiff configurations like discontinuities, or boundary layers. It becomes necessary to use more sophisticated less local iterations, and/or more sophisticated coarse mesh definition, as proposed by Algebraic MG [START_REF] Brandt | Algebraic Multigrid for automatic multigrid solutions with application to geodetic computations[END_REF] [START_REF] Stüben | Appendix A: An introduction to Algebraic Multigrid[END_REF] or by anisotropic mesh coarsening [START_REF] Mulder | A new multigrid approach to convection problems[END_REF], [START_REF] Mavriplis | Adaptive mesh generation for viscous flows using Delaunay triangulation[END_REF], [START_REF] Francescatto | Méthodes multigrilles par agglomération directionnelle pour le calcul d'écoulements turbulents[END_REF], and/or more adapted inter-grid transfers. Indeed, in some case, the directly-coarser grid correction is not able to complement the fine grid iteration, or may even work in a defavourable way. As a consequence, the best set of grids to apply for MG acceleration is not necessarily the best set of grid to apply for the FMG process. In this paper, we apply an isotropic refinement for the FMG nested iteration and an anisotropic mesh coarsening for MG acceleration.

Second, for a lot of complex applications, it has been remarked that FMG does not work, in the sense that the usual stopping criterion produces a discrete final solution with an accuracy deteriorated by an insufficient iterative resolution. See for example [START_REF] Carré | On the application of FMG to variational approximation of flow problems[END_REF]. Let us examine a possible FMG failure scenario. The two assumtions in FMG theory are (1) a MG convergence which does not depend of level fineness, (2) the asymptotic high-order convergence of the discrete solution to continuous one on the different meshes of FMG, including coarse ones. Assuming the MG cycling convergence is good, FMG failure can then be explained by the lack of asymptotic convergence to continuous, either because meshes are still too coarse, or because solution involves small details or singularities. As a consequence, when the basic FMG algorithm involving a fixed number of cycles per phase is applied, the solution produced at end of FMG may be inaccurate. Adjusting the number of cycles to the necessary convergence, if not done accurately, may result in a computational cost much larger than the one which the theory of FMG would let expect and may result in loosing the K.N complexity. Stopping criteria for iterative solvers have been the topic of many published works. In the case of quasi-Newton iteration, a typical work is [START_REF] Wathen | Stopping criteria for iterations in finite element methods[END_REF]. It is commonly admitted that the best criterion is to stop the iteration when the iteration error between iterated approximated solution and converged approximated solution is smaller than the approximation error between converged approximated solution and exact solution [START_REF] Trottenberg | Multigrid[END_REF], [START_REF] Wathen | Stopping criteria for iterations in finite element methods[END_REF]. But computing the approximation error is computationally costly. In [START_REF] Wathen | Stopping criteria for iterations in finite element methods[END_REF] its evaluation is replaced by a assumption concerning O(h 2 ) mesh convergence. In the present paper, we propose to use the a posteriori residual as reference for stopping MG cycles.

MG and FMG have been combined very frequently with mesh adaption. Let us cite a pioneering work of R. Bank, [START_REF] Bank | PLTMG: A software package for solving Elliptic Partial Differential Equations[END_REF], and a few more recent ones, such as [START_REF] Rüde | Fully adaptive multigrid methods[END_REF], [START_REF] Braaten | Three-dimensional unstructured adaptive multigrid scheme for the Navier-Stokes equations[END_REF], [START_REF] Nastase | A parallel hp-multigrid solver for threedimensional Discontinuous Galerkin discretizations of the Euler equations[END_REF], [START_REF] Mitchel | The hp-multigrid method applied to hp-adaptive refinement of triangular grids[END_REF]. Adaptive works are most frequently based on mesh refinement by local division, producing embedded meshes. More generally, unstructured nonembedded MG and FMG have been penalized during years by the difficulty in building and managing multiple coarse and fine unstructured meshes in particular for industrial applications. This difficulty is more easy to address today, with the recent progress of mesh generation and adaptation, see e.g. [START_REF] Frey | Mesh generation. Application to finite elements[END_REF]. Due to this progress, novel anisotropic strongly mesh-adaptive algorithms are now available. By strongly mesh-adaptive we mean that an anisotropic mesh adaption is strongly coupled with the solver thanks to a nonlinear fixed point iteration. Anisotropic mesh adaptors have been observed as carrying two important advantages. First not only many computations are performed in much better conditions than with traditional methods, but also they allow computations which were simply not feasible without anisotropic adaptation, like the propagation of a sonic boom from aircraft to ground [START_REF] Loseille | Adaptive anisotropic simulations in aerodynamics[END_REF]. Second, anisotropic mesh adaptors provide mesh convergence at high-order for singular problems [START_REF] Loseille | Adaptation de maillage anisotrope 3D multi-échelles et ciblée à une fonctionnelle pour la mécanique des fluides. Application à la prédiction haute-fidélité du bang sonique[END_REF]. For non-singular problems but rather heterogeneous problems, non-adaptive methods will produce higher order convergence only with very fine meshes. Anisotropic adaption will give a high order numerical mesh convergence with a much smaller number of nodes.

The plan is as follows. The next section introduces Riemannian metrics for defining what we call a continuous mesh model and the fixed-point mesh adaptation algorithm. Section 3 combines MG and mesh adaption. Section 4 presents FMG and proposes a stopping criterion for it. Section 5 defines the complet proposed algorithm, combining FMG and anisotropic mesh adaption. The paper is completed by several test cases and a discussion.

Mesh parametrization and mesh adaptation loop

Continuous mesh model

We recall shortly the continuous mesh framework, introduced in [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework. Part II: validations and applications[END_REF]. This framework lies in the class of metric-based methods. A continuous mesh M of the computational domain Ω is identified to a Riemannian metric field [START_REF] Berger | A panoramic view of Riemannian geometry[END_REF] M = (M (x)) x∈Ω . For all x of Ω, M (x) is a symmetric 2 × 2 matrix. Its diagonalisation writes:

M (x) = d(x) R(x) r -1 1 (x) r -1 2 (x) t R(x), (1)
The total number of vertices C is defined as:

C (M ) = Ω d(x) dx = Ω det(M (x)) dx.
A discrete mesh H of the same domain Ω is a unit mesh with respect to M , if, to simplify, each edge e = ab verifies 1 0

t ab M (a + t ab) ab dt ∈ 1 √ 2 , √ 2 .
Given a smooth function u, to each unit mesh H with respect to M corresponds a local interpolation error |u -Π H u|. In [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework. Part II: validations and applications[END_REF], it is shown that this interpolation error is well represented by the so-called continuous interpolation error related to M , which is locally expressed in terms of the Hessian H u of u as follows:

|u -π M u|(x) = 1 10 trace(M -1 2 (x) |H u (x)| M -1 2 (x)) (2)
where |H u | is deduced from H u by taking the absolute values of its eigenvalues. We define as optimal metric the one which minimizes the right-hand side under the constraint of a total number of vertices C equal to a parameter N. After solving analytically this optimization problem, this defines the unique optimal (M L p opt (x)) x∈Ω as:

M L p opt = D L p (det |H u |) -1 2p+2 |H u | with D L p = N 2 2 Ω (det |H u |) p 2p+2 -2 2 , (3)
where D L p is a global normalization term set to obtain a continuous mesh with complexity N and (det

|H u |) -1
2p+2 is a local normalization term accounting for the sensitivity of the L p norm. In the sequel we choose p = 2.

Fixed-point mesh adaptation

In the case where the function u is the solution of a Partial Differential Equation, the Hessian-based method extends heuristically as follows. Given a discrete solution u h to the PDE, a Hessian H(u h ) is defined from it. The so-called optimal mesh M pde is defined by:

M pde = M opt (H(u h (M pde )))
Where for any metric M , u h (M ) is the discrete PDE solution computed on a unit mesh for M . We solve the non-linear problem giving the optimal mesh M pde by applying the following loop:

Fixed point for adaptive PDE approximation 1-compute the PDE approximate solution u h on current mesh M 2-compute an approximate Hessian H(u h ) 3-adapt with N nodes according to this Hessian, obtain M = M opt (H(u h )) 4-go to 1.

For remeshing phases, we used indifferently MeshGems-Adapt of Distene and in-house versions of Yams ( [START_REF] Frey | About surface remeshing[END_REF][START_REF] Frey | Yams, a fully automatic adaptive isotropic surface remeshing procedure[END_REF]). Due to the discrete and noisy character of remeshing, a strictly-fixed point cannot be obtained, but instead the iteration needs be stopped when further work would be useless. In contrast to the adaptation to an analytic function, the deviation to a target is not available. However, it is possible to rely on the approximation of the interpolation error given by integral [START_REF] Bank | PLTMG: A software package for solving Elliptic Partial Differential Equations[END_REF]. In [START_REF] Alauzet | High order sonic boom modeling by adaptive methods[END_REF], it is proposed to stop the iteration when the difference between to approximate solution field u h is smaller than a positive quantity to choose cleverly since this difference does not converge to zero. The simple option of a fixed number of adaption iterations is also a rather secure one.

The MG Anisotropic fixed-point

Let us assume now that we want to solve our mesh-adaptive discrete PDE by means of a MG algorithm. This means that the PDE to solve, used for finding the approximate solution u h , is replaced by the problem of finding the couple (M h , u h ) such that: M h is adapted to u h and u h is computed on mesh M h .

As noted in previous section this is the solution of a non-linear coupling system, but the dependancy of M h with respect to u h is explicit, in the sense that the cost of systems to solve is neglectible with respect to the cost of computation of the solution u h on a given mesh. Therefore, we propose to apply the adaptive loop as an external one, the MG resolution by MG being an internal loop. The resulting algorithm is depicted in Figure 2. The external loop of mesh adaption is iterated five times for convergence of the coupling bertween PDE solution and adapted mesh. We now define in more details the ingredients of this loop.

Anisotropic coarsening

Let us examine how to build coarser meshes in order to apply MG. The adopted stanpoint is to use the metric based mesh parametrization. Firstly, we specify the number of nodes N k of the adapted fine mesh of the current phase, phase number k. We have to choose an initial metric M :

M (x, y) = R t (x, y) 1 ∆ξ 2 (x,y) 0 0 1 ∆η 2 (x,y) R(x, y)
with: ∆ξ (x, y) = mesh size in the first caracteristic direction ∆η(x, y) = mesh size in the second caracteristic direction R(x, y) = matrix of eigenvectors.

The specification of the number of nodes of this fine mesh writes:

(∆ξ ∆η) -1 dxdy = N k (4)
where the integral is taken over the computational domain. Then coarser metrics are build using the metric-based embedding:

M coarser 1 (x, y) = R t (x, y) 1 4∆ξ 2 (x,y) 0 0 1 4∆η 2 (x,y) R(x, y)
and even coarser:

M coarser 2 (x, y) = R t (x, y) 1 16∆ξ 2 (x,y) 0 0 1 16∆η 2 (x,y) R(x, y)
etc. In particular the first coarser mesh has about N k /4 nodes, the next coarser has about N k /16 nodes, etc.

Anisotropic coarsening can also be applied. For example, assuming that the ordering of eigenvalues satisfies ∆ξ ≤ ∆η, we can coarsen in an equivalent way to [START_REF] Francescatto | Méthodes multigrilles par agglomération directionnelle pour le calcul d'écoulements turbulents[END_REF]:

M coarser 1 (x, y) = R t (x, y) 1 4∆ξ 2 (x,y) 0 0 1 Max(∆η 2 (x,y),4∆ξ 2 (x,y))

R(x, y).

This option has been tested but did not improve the results for the test cases we present, which involve meshes which are not enough stretched.

MG

The above metric coarsening produces a sequence coarse meshes H 1 , ...H max to be used together with the fine initial mesh as levels + 1, ..., max for a MG cycle. Those are kept during the MG cycles and regenerated during the adaptation phase when the fine mesh is adapted. For applying the MG cycle, transfers are defined as follows: correction transfers from coarse to fine are P 1 interpolated in triangles, and residual transfers, from fine to coarse, are accumulated on coarse nodes with barycentric weighting. A saw-tooth V-cycle with 10 Jacobi sweeps as pre-smoothing and without post-smoothing is applied.

Global fixed point

In the mesh-adapted MG, the adapted solution with a prescribed number of nodes N k is obtained by encapsulating the MG cycle into the adaption loop. For the adaptation convergence, we have chosen to uniformly apply 5 adaptations.

The FMG algorithm

The adaptive FMG is the succession of adaptive MG phases with transfer of the solution between each phase. At phase k φ , the number of nodes is prescribed to be equal to N k φ . We have chosen the usual option of a new mesh size two times smaller in next FMG phase:

N k φ +1 = 4 N k φ . (5)
The resulting approximation error will be presumably 4 times smaller. With N k φ +1 = 2 N k φ the error would be 2 times smaller, and so for the necessary cycles but two phases would be necessary for the same final accuracy, with a similar global cpu effort. Therefore the choice of N k φ +1 /N k φ is not a sensible one.

Global FMG under O2 convergence assumption

FMG can be defined as the combination of a MG loop with a nested iteration. A first coarse mesh is used for a first evaluation of the solution. On the coarse mesh, in principle, a coarser level for acceleration is not necessary since the convergence of a standard iterative solution algorithm is rather fast. A finer mesh is built, generally by uniformly refining the first mesh. The previous solution is transfered to the new mesh, typically by interpolation. The two meshes are available for playing the role of two levels in order to solve fastly the problem with a two-grid iteration on the new mesh, starting with a good initial condition. This process is reiterated with a 3-grid solution on next mesh etc. We call FMG k φth phase the k φ -th nested iteration phase, using with k φ -grid cycling. Due to the initialization by the previous phase, a sufficient convergence at each phase k φ can be obtained by a small number k c of k φ -grid cycles. Due to the ability of MG in exhibiting a convergence rate quasi-independant of number of level and grid size, an important gain is reachable. The prescription in stopping iterative solutions is generally based on the following criterion:

Stopping criterion: Iteration error should be smaller than the approximation error. In the FMG theory presented in [START_REF] Hackbusch | Multigrid methods and applications[END_REF], it is assumed that there exists a constant K such that for any phase k φ , we have

||u -u k φ || ≤ Kh 2 k φ . (6)
The error for u k φ is four times larger than for u k φ -1 , which suggests a residual reduction by 4 at each phase. A more accurate analysis ( [START_REF] Hackbusch | Multigrid methods and applications[END_REF]) shows that it is enough to converge each FMG cycling phase by dividing the residual by 10, for example:

||A h u k c h -f h || ≤ 1 10 ||A h u 0 h -f h ||. (7)
The weakness of this theory lies in the central assumption [START_REF] Carré | On the application of FMG to variational approximation of flow problems[END_REF]. An evident first remark is that ( 6) is a mesh convergence assumption, which can be established only for a fine enough mesh, that is for grids/levels with large enough. For close to unity, the number of nodes is small and the above estimate is clearly not true. But the main disaster arises when [START_REF] Carré | On the application of FMG to variational approximation of flow problems[END_REF] does not holds for the last fine grids.

In that latter case, the dilemma of FMG is: Either the convergence test [START_REF] Francescatto | Méthodes multigrilles par agglomération directionnelle pour le calcul d'écoulements turbulents[END_REF] based on a fixed decrease of residual at each phase, being not sufficient, will produce an inaccurate solution on the finer mesh, or, assuming we have found a sufficiently severe way to stop convergence in each phase, we get an accurate solution, but we have increased the cost in a rather unpredictible way and have presumably lost the optimal K •N complexity of FMG.

We are also motivated by a second issue: if a mesh adaption loop is applied, the MG cycling between mesh updates will be initialised by a candidate solution obtained with a just slightly less adapted mesh of same fineness. This candidate solution can be already very close to the converged solution. Then, in order to avoid a rather large amount of unnecessary computing time, we need to recognize it and stop early the MG cycling.

Next two paragraphs deal with introducing a measure of approximation error and a control of iteration error in FMG, in order to improve the cycling stopping test for FMG, consolidating the accuracy of FMG. The second issue, efficiency, will be addressed by introducing mesh adaption, in Section 7.

Approximation error estimate

Let u k h (x) = ∑ u k h,i N i (x) be the approximate function at iteration k of a given iterative solver (N i holds for the finite-element basis function related to node i). For a second-order accurate approximation, Arioli and co-workers propose in [START_REF] Wathen | Stopping criteria for iterations in finite element methods[END_REF] a stopping criterion based on:

||A h u k h -f h || H -1 ≤ h 2 || f h || H -1
where H -1 holds for the dual of the Sobolev space H 1 0 . Using the h 2 factor assumes that the numerical mesh-convergence is close to scheme asymptotic convergence.

We propose here a method which does not rely on mesh convergence. Let us introduce some notations: Let V = H 1 0 (Ω), Ω being the computational domain. The continuous PDE system is written in short:

Au = f or u ∈ V ∀ φ ∈ V a(u, φ ) = ( f , φ )
Let V h be the usual P 1 -continuous finite-element approximation space. The finiteelement discretisation is written:

u h ∈ V h ∀ φ h ∈ V h a(u h , φ h ) = ( f , φ h )
We introduce the components of u h for the finite-element basis (N i ) as follows, with T : R n → V :

u h = ∑ u h,i N i ⇔ u h = T u h
Let us denote T * the adjoint of T :

T * : V → R n [T * f ] i = ( f , N i ).
The variational discretization:

a( ∑ u h, j N j , N i ) = ( f , N i ) ∀i ⇔ ∑ a(N j , N i )u h, j = ( f , N i ) ∀i.
transforms into an algebraic one:

A h u h = f h , where [A h ] i j = a(N j , N i ) and f h = T * f . (8)
The exact a posteriori estimate:

u -u h = A -1 ( f -Au h )
can approximated as:

u -u h ≈ T A -1 h T * ( f -Âu h ) where Âu h denotes a smoothed approximation of Au h ( Âu h , φ ) = ∑ i j |D i j | |i j| D i j φ [∇u h ] i j • n i j dv
where :

φ is an arbitrary function of V , -the sum ∑ is taken for all internal edge i j of the mesh (2D), -[] holds for the jump of quantity inside bracket across the edge i j, -the integral D i j is taken over surface (2D) of the diamond quadrilateral D i j:

D i j = iG i j jG ji
where G i j and G ji are the centroids of the two elements (triangles) having i j as common edge.

4.3. Iterative-Error-controlled FMG Now the algebraic system ( 8) is solved by a number k c of cycles:

u 0 h = 0 ; u k c +1 h = Iterate(u k c h ) ; u ∞ h = u h
The iterative error can be evaluated by solving the system with a right-hand side equal to the local iteration residual:

u h -u k c h = A -1 h (f h -A h u k c h ).
It remains to compare the iterative error with the above approximation error. Remember first that the approximation error above derivation did not use the assumption that u h is the solution of the discrete system. In particular, the same estimate holds for the result u k c h of the incomplete iterative resolution. Let us introduce the element of V h :

u k c h = T u k c h = ∑ [u k c
h ] i N i , Then, using the a posteriori error estimate:

u -u k c h ≈ T A -1 h T * ( f -Âu k c h ). (9)
It is now useful to transform the algebraic iterative residual A h u k c h -f in similar terms. We start from:

A h (u k c h -u h ) = A h u k c h -f h -(A h u h -f h ) = A h u k c h -f h , thus u k c h -u h = A -1 h (A h u k c h -f h ) ⇔ u k c h -u h = T A -1 h (A h u k c h -f h ). (10)
Heuristics: Assuming that, in some norm to specify later, ||uu h || ≤ ε, ε small and positive, and that after k c solver iterations, we have

||u h -u k c h || ≤ 0.1 ||u -u k c h ||, (11)
then ||u -u k c h || ≤ 1.1 ||u -u h ||. ( 12
)
As already mentioned, the last statement ( 12) is for us an acceptable iterative convergence stopping criterion. We observe that stopping criterion ( 11) is realizable, since as iteration number increases, the norm ||uu k c h || is supposed to converge to ||uu h ||, assumed to be not zero, while ||u hu k c h || can be driven to machine zero by iterating over k. Unfortunately, evaluating the two terms of [START_REF] Hackbusch | Multigrid methods and applications[END_REF] involve solving two discrete systems with matrix A, a computation which is more or less computationally as costly as the original system to solve, and is therefore too costly. A possible solution is to solve approximatively the discrete error system with a coarse grid. Here, we propose to decrease the cost of the stopping test, with some risk of decreasing its accuracy, by taking the following l 1 norm of the right-hand side:

||f h -A h u k c h || l 1 = ∑ i |[f h -A h u k c h ] i |
We rely on the observation that MG-cycles decrease many different norms of the residual with about the same slope. Therefore, we do not claim that the proposed stopping criterion is able to work adequately when associated with another iteration than MG.

Stopping test 1: Assume that after k solver iterations

||f h -A h u k c h || l 1 ≤ ε ||f h -A h u 0 h || l 1 , ε , small and positive, ( 13
)
and that we have

||f h -A h u k c h || l 1 ≤ 0.1 ||T * ( f -Âu k c h )|| l 1 , (14)
then stop the iteration.

According to [START_REF] Francescatto | Méthodes multigrilles par agglomération directionnelle pour le calcul d'écoulements turbulents[END_REF], ε is chosen to be 1 10 . In practice, since the computation of the RHS of ( 14) may need more cpu than a cycle, the test [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF] means that several iterations are performed in order to decrease the iterative residual to satisfy (13) before the second test ( 14) is evaluated. If test ( 14) is negative, several iterations are again performed before a second test of ( 14) is again evaluated, etc. This splitting allows for a lower CPU cost. This device is inspired by an analog one proposed by Arioli and co-workers [START_REF] Wathen | Stopping criteria for iterations in finite element methods[END_REF]. Our final formulation is as follows:

Stopping test 2: 1. Iterate cycling until ||f h -A h u k c h || l 1 ≤ 1 10 ||f h -A h u 0 h || l 1 2. If ||f h -A h u k c h || l 1 > 0.1 ||T * ( f -Âu k c h )|| l 1 then u 0 h = u k c h go to 1. 3.Stop the iteration.

Application to the proposed FMG

The MG algorithm which we use is built from: -a sequence of unstructured meshes which are not necessarily embedded. the basic inter-grid transfers are classically the P 1 interpolation (for a transfer from a coarse mesh to a finer mesh) and an accumulation weighted with barycentric ceofficients (for a tranfer from a fine mesh to a coarse mesh).

-a saw-tooth V-cycle with 10 damped-Jacobi relaxation a a smoother.

-an encapsulation of this MG cycle as a preconditioner of a GMRES loop. When we shall talk about "a cycle", we shall mean the combination of the MG V-cycle with the GMRES updating. The rather high number of sweeps, together with the use of GMRES is the price we pay in order to get a robust convergence for high-density ratio case.

Three contexts are now examined for showing how works the combination of FMG with the stopping criterion.

-Fig. 4: the function to compute is not important, but to fix the ideas, it is the circular test case described in the sequel. The initial solution is uniform. We apply 90 GMRES-MG cycles. The preconditioned (by MG) residual l 1 norm (started at iteration 1) indeed decreases momotonely and fastly from 1 to 10 -8 (multiply-shaped marked curve). Marked by plus-covered-by-multiply, the equation residual norm ||f h -A h u k c h || l 1 starts from less than 1, shows an increasing phase and then decreases to about 2.5 10 -3 . Marked by plus-symbols, the approximation residual norm ||T * ( f -Âu k c h )|| l 1 also first increases in a similar way and then goes down to a non-zero limiting value. The two above above curves intersect at about 70 GMRES-MG cycles. This is probably too many cycles, i.e. rather conservative. In the chosen example we know the exact solution and can also depict the approximation error norm ||uu k c h || L 1 , which starts from a number close to 1 and decreases to a limiting level ||uu h || L 1 of about 5 10 -5 . We get confirmation that with 70 cycles, this level is tightly approached.

-Fig. 5: the main change is that the initial solution is provided by interpolating the discrete solution computed on previous coarser mesh. Cycles are stopped at 34. The preconditioned (by MG) residual l 1 norm starts with a small level, and decreases fastly to 5 10 -8 (multiply-shaped marked curve). Marked by pluscovered-by-multiply, the equation residual norm ||f h -A h u k c h || l 1 starts from about 0.5, and monotonely (this time) decreases to about 2 10 -3 . Marked by plussymbols, the approximation residual norm 0.1||T * ( f -Âu k c h )|| l 1 starts from a lower value then previous case and monotonely decreases in a similar way to equation residual, going down to a non-zero limiting value. The two above above curves intersect at about 22 GMRES-MG cycles. In contrast to the previous case, the approximation error norm ||uu k c h || L 1 , starts from a low level of about 10 -4 and decreases to about 5 10 -5 , indicating (1) that the numerical convergence is not so good (first-order) between the two meshes and (2) that again the iterative error at iteration 22 is much smaller than the approximation one.

Fig. 6: we give an example of FMG sequence for a slightly easier problem (Laplace equation with uniform meshes).The convergence on four sucessive meshes is shown. The approximation error is numerically converging at second order. With mesh 2 and mesh 4, the second test is negative and a second MG convergence is applied. The four phases are complete with a total of 16 cycles, that is a reasonable mean number of 4 cycles per FMG phase.

The FMG Anisotropicly adaptive algorithm

The synthesis of the above sections is the FMG anisotropicly adaptive algorithm. We insert the adaptation loop as an intermediate loop between FMG phases and MG cycles. Concerning the process of going to a larger number of nodes, we keep the previous meshes and define a finer one by a simple division of each element into four elements of same area. The global algorithm is sketched in Fig. 7. Let us re-visit the ways these loops are stopped. The external loop, FMG phase will increase the number of mesh nodes. Theoretically it should stop when some norm of the approximation error |uu h | is smaller than a number prescribed by user. This option is yet rather far from practice and is not studied in this work. The intermediate loop, mesh adaption is stopped after 5 iterations. As concerns inner loop, the cycling loop, it is controlled by the stopping criterion defined in Section 6. The stopping criterion is used when (1) changing from a first mesh of N nodes to a finer mesh of 4N nodes between two FMG phases, as well as (2) changing from a first mesh of N nodes to a second mesh of same number N of nodes, but more adapted.

Examples

The proposed Adaptive Anisotropic FMG is applied to three test cases and compared with a pure FMG method applied with a sequence of embedded uniform meshes. Due to our restriction to a Poisson-like model, test cases will be toy problems. However, we choose them in order to represent the three following typical difficulties of multiphase incompressible flows: -boundary layers, -discontinuous phase changes, -Dirac layer source terme from capillarity. Our three simplified representations of these difficulties, in combination with the specification of mesh sizes (number of vertices) constitute a small benchmark for the performance of mesh adaptive methods, which could be used for the evaluation of various sensors.

In cases where the test case has an analytic solution, we shall call total approximation error or simply approximation error the error between the analytic solution from one side, and, from the other side, the discrete solution produced the algorithm: since the GMRES-MG algorithm is not converged to machinezero, our approximation error combines numerical scheme approximation error and iterative error.

A smooth boundary layer test case

For modelling the stiffnes of a boundary layer, we consider a Poisson problem with a smooth solution presenting some anisotropic local variation. Let: rhs(x, y) = α 2 (exp(1/α) -1) -1 exp(x/α) with α = 0.03.

We solve

-∆u = rhs with ∂ u ∂ y (x, 0) = ∂ u ∂ y (x, 1) = 0 and u(0, y) = u(1, y) = 0. Then u(x, y) = [exp(1/α) -1] -1 exp(x/α) + x + [exp(1/α) -1] -1 .
An example of approximated solution u h is shown in Figure 8. That allows us to compute directly We observe a convergence of order two for the non-adaptive case and a similar convergence in the adaptive case. Since the criterion for adaption which we adopted postulates that the interpolation error is a good representation of approximation error, it is interesting to examine the convergence of the interpolation error of the exact solution on the meshes we used, see Figure 10. We observe that indeed both convergence are similar, but that the interpolation error decreases to values which are smaller by a factor larger than two orders of magnitude. This observation will be also done for the other test cases and deliver a strong message saying that the interpolation error and its companion the Hessian criterion, while providing rather good adaptations, are not faithful representations of the actual approximation error. The difference between Figure 9 and Figure 10 measures the relative inadequacy of the Hessian-based option. Comparing the non-adaptive case and the adaptive one (Figure 9) for the same number of vertices, we observe that the error ||uu h || L 1 is notably smaller in the adaptive case. We also draw this error in function of the CPU time in Figure 11. We distinguish mesh division phase with steep slopes from the mesh adaption ones with less steep slopes.

Thanks to the stopping criterion, the iterations 2 to 5 of adaption consume less cycles and therefore less CPU. For about 1000 seconds of a workstation, the obtained accuracy is 610 -3 for the non-adaptive case, and 10 -5 for the adaptive one. the accuracy fo the non-adaptive calculation with 750 seconds is obtained by adaption with 10 seconds.

A non-smooth internal layer case

The second test case exemplifies the singularity which is met in the simulation of multiphase non-mixed flows with a large deviation between the physical properties of each phase. Let us consider the equation of Poisson -div( 1ρ ∇u) = rhs with discontinuous coefficient ρ and a right hand side rhs which are strongly discontinuous on the domain, as it is shown in Figure 12. The solution u of the homogeneous Dirichlet prolem has discontinuous gradients along the coefficient discontinuity. A mesh-adaptive approximation u h is depicted in Fig. 13. Since this time an analytic solution is not available, we compute the L 1 -norm of the solution u h : ||u h || L 1 = Ω |u h |dxdy and compare with an interpolated evaluation on uniform mesh. Figure 14 shows this norm in function of the number of points. It is also shown in Tables 1 and2. The uniform-mesh approach is definitively penalized by the singularity of the solution. We can expect first-order convergence and indeed the observed numerical convergence order of the non-adaptive L1-norms is 0.96 for the finest computations. Since the L 1 norm is an integral, we can try an extrapolation of it, which gives ||u|| L 1 ≈ 0.82, rather close to our mesh adaptive results. But the L 1 norm obtained with more than a 100, 000-node uniform mesh and a CPU time of 1811 seconds still show an error of more than 8%. With adaptation, numerical convergence is rather noisy. However, a L 1 norm at less than 1% from fine-mesh one is already obtained with 552 nodes and a CPU time of 57 seconds. Figure 15 shows the same norm in function of the CPU-time.

Circular test case

Capillary models exhibit, along the interfaces, Dirac layer source terms for the pressure equation. These terms imply discontinuous pressures. For example the pressure could be is equal to 1 on a disk at center and equal to 0 in the rest of the domain. Instead of considering a strictly discontinuous solution, we approach it by defining a thickness ε of the layer between the two uniform phases as shown in Figure 16 given by: u c (x, y)

= 1 2 [1 + ψ ε + 1 π sin( πψ ε )] with ψ = R -(x C -x) 2 + (y C -y) 2 .
The value of ε controls the thickness of the transition between u c = 1 and u c = 0 and is chosen equal to 0.02. Let rhs = ∆u c . We consider the Dirichlet problem ∆u = rhs in Ω, u = 0 on ∂ Ω. The right-hand side rhs is close to a Dirac distribution concentrated along the circle limiting the disk. In practical nonlinear situations as capillary models, that kind of feature is not a priori known. Then we choose in our variational formulation to integrate the discrete RHS (rhs, N i ) (N i : finite element basis function) on the given mesh without particular care of the quasi singularity of rhs, which means that a coarse mesh may produce a very inaccurate solution. Indeed, in our computations we observe in Fig. 22 that very large errors are produced by coarse uniform meshes. With 10, 000 nodes, a 100 % L 1 error is still produced. Full second-order asymptotic convergence seems be reached only after 30000 nodes are used. This behavior can be a strong handicap for 3D calculations where the number of nodes cannot be much increased. In contrast, the mesh adaptive computation produces much smaller errors with coarse meshes and always perform as well or better. An mesh-adapted approximate solution u h is shown in Figure 17. A second remark is that most gain of adaption is obtained at about 1000 -2000 nodes while adaption for finer meshes seem to bring no further acceleration (slope remains second-order). This phenomenon is observed also for the interpolation error of the exact solution onto the different meshes which are used. Our interpretation is that once the layer around the circle is captured, no further adaption in needed, since a fine enough, but uniform mesh is quasi-optimal in the vicinity of the layer described by a sinus function. Then the subsequent efforts in adaption are useless, which explains that after an interesting performance for 2000 nodes, the adaption option looses its CPU advantage.

Thinner circular test case

Now, the value of ε controlling the thickness of the transition between u = 1 and u = 0 and is chosen equal to 0.001. An approximate (adapted) solution u h is shown in Figure 21 (right). Now, due to the very thin definition of the Diractype right-hand-side, the brut force use of an uniform mesh of 100, 000 nodes does not allow the computation of a good solution, see Fig. 21 (left). Probably, a good solution is obtainable when a sufficiently fine uniform mesh is considered in a subsequent nested-iteration phase, but in that latter case, MG convergence on that mesh would need be iterated during many iterations and the FMG ideal complexity is lost. In the same figure is depicted the approximate mesh-adaptive solution with the same number of nodes. The adaptive L 1 approximation error norm is 0.3 with 10, 000 nodes, and 0.04 with 100, 000 nodes.

Concluding remarks

We have proposed a combination of the well-established FMG method with an anisotropic mesh adaption method. The mesh adaption fixed point loop is introduced in the FMG process. In order to master the extra computational complexity, an improved stopping criterion for MG cycling inside FMG is proposed. We believe that it is important for robustness and efficiency. We emphasize that this stopping criterion assumes that the iterative MG convergence rate is more or less norm-independant, a property surely not enjoyed by many other (non-MG) solution algorithms, for which our stopping criterion is definitively not recommanded.

The number of adaption iterations is fixed once for all. The overall anisotropic adaptive FMG is of rather high programming/algorithmic complexity, due to the higher number of tests and embedded loops. The central question is therefore: does it enjoy a robust computational efficiency, that is, is it computationally efficient for difficult problems.

The few numerical experiments tend to promote a positive answer. The four test cases are run with a unique set of parameters, i.e. without parameter tuning. In contrast to the non-adaptive case, with the use of anisotropic adaptation, we observe the early capturing of many different scales. A significative comparison between AFMG and FMG relies on the total approximation error as a function of CPU time. In most test cases, for a same CPU time, the mesh adaptive computation produces a lower or much lower approximation error. For a same CPU time, the number of nodes is much lower. A second measure concerns the asymptotic behavior in terms of number of unknowns and CPU. When it works, our FMG indeed shows O(N) complexity, namely the considered norm of the total approximation error is proportional to the logarithm of the number of nodes N. This also approximatively holds for the AFMG version. For FMG, with some variations, the norm of total approximation error is also proportional to the logarithm of CPU time. Because of our choice of limiting adaption iteration to 5, we expect that this also holds for AFMG, and we observe it, again approximatively. The discontinuous coefficient case desserves a particular mention since even for the L 1 norm of solution, the convergence on uniform meshes is first-order. The improvement in that case is of two orders of magnitude.

We have not proposed a smart stopping criterion for the adaption loop. We plan to discuss this issue in a forthcoming paper in combination with a different mesh adaption criterion.

Indeed, the present study relies on the Hessian-based mesh adaptive criterion. This simple and robust option has some limits, which we have measured by comparing the convergence of the interpolation error and the convergence of the actual approximation error. The Hessian-based criterion is designed exactly for the interpolation error, and, by the way, the interpolation error converges fastly to small values. This shows that the different approximations of u by u h and of the Hessian of u h do not introduce a too important penalty. The -expected-bad news is that the approximation error does not decrease as fastly and as low. We interpret this as an effect of the lack of consistency between the interpolation error and the approximation error. Introducing adjoint-based adaption criteria may improve this issue. We plan to discuss this in a forthcoming paper using the benchmark proposed in the present work.

Of course, we are not sure that the extension to 3D will enjoy the same qualities, but the present results are encouraging. Then, such an evaluation is now in progress. Also, in our opinion, the extension of the proposed methods to other models of Continuum Mechanics can be envisaged as far as the application of MG works satisfactorily.

Figure 1 :

 1 Figure 1: Mesh adaption loop

Figure 2 :

 2 Figure 2: Mesh adaptation loop with Multigrid

Figure 3 :

 3 Figure 3: Each phase k φ of the Full Multi-Grid algorithm is made of several k φ -grid cycles

Figure 4 :

 4 Figure 4: Convergence of the GMRES-MG-iterative l 1 residual norm ||A h u k hf h || l 1 (+ ×), the approximation l 1 residual norm || Âu k hf || l 1 (+), the preconditioned residual norm (×), the norm ||uu h || L 1 of deviation to exact ( ), for a Poisson problem, starting from a uniform field u 0 h = 1 at iteration 0.

Figure 5 :

 5 Figure 5: Convergence of the GMRES-MG-iterative l 1 residual norm||A h u k hf h || l 1 (+ ×), the approximation l 1 residual norm || Âu k hf || l 1 (+), the preconditioned residual norm (×), the norm ||uu h || L 1 of deviation to exact ( ), for a Poisson problem, starting from a coarser-grid interpolated solution at iteration 0.

Figure 6 :

 6 Figure 6: Convergence of the FMG-iterative l 1residual norm ||A h u k hf h || l 1 (+ ×), the approximation l 1 residual norm || Âu k hf || l 1 (+), the norm ||uu h || L 1 of deviation to exact ( ), for a Poisson problem, starting from the coarsest-grid (grid 1) solution, and performing four FMG phases from grid 2 to grid 5.

Figure 7 :

 7 Figure 7: Mesh adaptive Full Multi-Grid

Figure 8 :

 8 Figure 8: Boundary layer test case solution and adapted mesh

Figure 9 :

 9 Figure 9: Boundary layer test case. Approximation error ||uu h || L 1 as a function of the number of mesh nodes. (+) non-adaptive FMG, (×) adaptive FMG. The straight line shows the second-order slope.

Figure 10 :

 10 Figure 10: Boundary layer test case. Corresponding behavior of the interpolation error of exact solution ||u -Π h u|| L 1 on the same meshes as in Fig.9. The straight line shows the second-order slope.

Figure 11 :

 11 Figure 11: Boundary layer test case. Approximation error ||uu h || L 1 as a function of the CPU time. (+) non-adaptive FMG, (×) adaptive FMG.

Figure 12 :Figure 13 :

 1213 Figure 12: Stiff layer test case domain

Figure 14 :

 14 Figure 14: Stiff layer test case results: L1-norm ||u h || L 1 of the approximate solution as a function of the number of points.(×) non-adaptive FMG, (+) adaptive FMG.

Figure 15 :

 15 Figure 15: Stiff layer test case results: L1-norm ||u h || L 1 of the approximate solution as a function of the CPU time.(×) non-adaptive FMG, (+) adaptive FMG.

Figure 16 :

 16 Figure 16: Circular-test-case-domain

Figure 17 :Figure 18 :

 1718 Figure 17: Circular layer test case solution and adapted mesh

Figure 19 :

 19 Figure 19: Circular test case. Behavior of the interpolation error of exact solution ||u -Π h u|| L 1 as a function of the number of mesh nodes, on the same meshes. The straight line shows the second-order slope.

Figure 20 :

 20 Figure 20: Circular test case: errors as a function of the CPU time. Approximation error ||u-u h || L 1 as a function of the CPU time. (×) non-adaptive FMG, (+) adaptive FMG. The straight line shows the second-order slope.

Figure 21 :

 21 Figure 21: Thin-circular layer test case uniform-mesh solution and mesh-adaptive solution.

Figure 22 :

 22 Figure 22: Thin circular test case: errors as functions of the number of mesh nodes. Approximation error ||uu h || L 1 as a function of the number of mesh nodes. (×) non-adaptive FMG, (+) adaptive FMG. The straight line shows the second-order slope.

Figure 23 :

 23 Figure 23: Thin circular test case: errors as a function of the CPU time. Approximation error ||uu h || L 1 as a function of the CPU time. (×) non-adaptive FMG, (+) adaptive FMG. The straight line shows the second-order slope.

Table 1 :

 1 . If (x, y) is located inside the thickness of the layer, u c (x, y) is Stiff internal layer: convergence of L 1 -norm of the approximate solution for a series of embedded uniform meshes.

	Numb. nodes	121	441	1681	6561	25921 103041
	L 1 norm	0.1354 0.01929 0.03806 0.05679 0.06869 0.07488
	Numb. nodes	142	552	2089	9243	36126
	L 1 norm	0.07512 0.08211 0.08292 0.0831 0.08376

Table 2 :

 2 Stiff internal layer: convergence of L 1 -norm of the approximate solution for a series of adapted meshes.

Acknowledgements

We thank Adrien Loseille and Frédéric Alauzet for the fruitful discussions we had together and help in using mesh software. Thanks also for Franc ¸oise Loriot from Distene. This work has been supported by French National Research Agency (ANR) through COSINUS program (projects ECINADS n o ANR-09-COSI-003 and MAIDESC n o ANR-13-MONU-0010). HPC resources from GENCI-[CINES] (Grant 2010-x2010026386 and 2010-c2009025067) are also gratefully acknowledged.