Didier Chauveau
email: didier.chauveau@univ-orleans.fr

Pierre Vandekerkhove

The Nearest Neighbor entropy estimate: an adequate tool for adaptive MCMC evaluation

Keywords: Adaptive MCMC algorithms, Bayesian model, entropy, Kullback divergence, Metropolis-Hastings algorithm, nearest neighbor estimation, nonparametric statistic

 but used mostly in univariate situations until recently when entropy estimation has been considered in other fields like neuroscience. Theoretically, we prove that, under certain uniform control conditions, the successive densities of a generic class of Adaptive Metropolis-Hastings algorithms to which most of the strategies proposed in the recent literature belong can be estimated consistently with our method. We then show that in MCMC setup where moderate to large dimensions are common, this estimate seems appealing for both computational and operational considerations, and that the problem inherent to a non neglictible bias arising in high dimension can be overcome. All our algorithms for MCMC simulation and entropy estimation are implemented in an R package taking advantage of recent advances in high performance (parallel) computing.

Introduction

A Markov Chain Monte Carlo (MCMC) method generates an ergodic Markov chain for which the stationary distribution is a given probability density function (pdf) f . For common Bayesian inference, f is a posterior distribution of the model parameter θ over a state space Θ ⊆ R d . This posterior is typically known only up to a multiplicative normalizing constant, and simulation or integration w.r.t. f are approximated by ergodic averages from the chain. The Metropolis-Hastings (MH) algorithm [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF][START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF] is one of the most popular algorithm used in MCMC methods. Another commonly used method is the Gibbs sampler introduced by [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the bayesian restoration of images[END_REF].

Each step of a MH algorithm at a current position θ t is based on the generation of the proposed next move from a general proposal density q(•|θ t ). Historically, two popular MH strategies used to be (i) the Independence Sampler (MHIS), which uses a proposal distribution independent of the current position, and (ii) the Random Walk MH algorithm (RWMH), for which the proposal is a random perturbation of the current position, most often drawn from a Gaussian distribution with a fixed variance matrix that has to be tuned.

To actually implement a MCMC algorithm, many choices for the proposal density are possible, with the goal of improving mixing and convergence properties of the resulting Markov chain. For instance running a RWMH strategy requires the determination of a "good" scaling constant, since the mixing depends dramatically on the variance matrix of the perturbation [START_REF] Roberts | Optimal scaling for various Metropolis-Hastings algorithms[END_REF]. As a consequence, a growing interest in new methods appeared this last decade, which purpose is to optimize in sequence the proposal strategy in MCMC algorithms on the basis of the chain(s) history; see, e.g., [START_REF] Andrieu | A tutorial on adaptive MCMC[END_REF] for a recent survey. These approaches called adaptive Markov Chain Monte Carlo (AMCMC) can be described (not in an entirely general way) as follows: let f be the pdf of interest and suppose that we aim to simulate efficiently from f given a family of Markov kernels {P ϑ , ϑ ∈ E}. This can be done adaptively using a joint process (θ t , ϑ t ) t≥0 such that the conditional distribution of θ t+1 given the information available up to time t is a kernel P ϑ t where ϑ t is an Euclidean parameter tuned over time to fit a supposed relevant strategy. Some general sufficient conditions insuring convergence (essentially ergodicity and the strong law of large numbers) of such algorithms have been established by various authors, see [START_REF] Andrieu | A tutorial on adaptive MCMC[END_REF]. These conditions are informally based on the two following ideas. Containment: for any (θ 0 , ϑ 0 ), and any ε > 0, the stochastic process (M ε (θ t , ϑ t )) t≥0 is bounded in probability, where

M ε (θ, ϑ) = inf t ≥ 1 : P t ϑ (θ, •) -f (•) T V ≤ ε
is the "ε-time to convergence". Diminishing Adaptation: for any (θ 0 , ϑ 0 ), lim t→∞ D t = 0 in P θ 0 ,ϑ 0 -probability, where

D t = sup θ∈Θ P ϑ t+1 (θ, •) -P ϑ t (θ, •) T V ,
represents the amount of adaptation performed between iterations t and t + 1. Note that in [START_REF] Bai | On the containment condition for adaptive markov chain monte carlo algorithms[END_REF] two examples are provided to show that either Diminishing Adaptation or Containment is not necessary for ergodicity of AMCMC, and diminishing Adaptation alone cannot guarantee ergodicity. See also the very simple four-state Markov Chain Example 1 in [START_REF] Rosenthal | Coupling and ergodicity of adaptive mcmc[END_REF], which illustrates the fact that ergodicity is not an automatic heritage when adapting a Markov Chain from its past.

These various and sometimes experimental algorithmic choices are associated in general to unknown rates of convergence because of the complexity of the kernel, and the difficulty in computing, when available, the theoretical bounds of convergence. For instance, [START_REF] Bai | Divide and conquer: A mixture-based approach to regional adaptation for mcmc[END_REF] compare two AMCMC strategies in dimension d ≤ 5, and [START_REF] Vrugt | Accelerating markov chain monte carlo simulation by differential evolution with self-adaptive randomized subspace sampling[END_REF] compare two AMCMC's against some benchmark in dimension d = 10. More recently [START_REF] Fort | A central limit theorem for adaptive and interacting markov chains[END_REF] define the best interacting ratio for a simple equi-energy type sampler, by minimizing the corresponding limiting variance involved in the Central Limit Theorem (see Fig. 1 in [START_REF] Fort | A central limit theorem for adaptive and interacting markov chains[END_REF]). There are also recent works proposing tools or methods for MCMC comparisons, showing that these questions are crucial in nowaday MCMC application and research. [START_REF] Thompson | SamplerCompare: A framework for comparing the performance of MCMC samplers[END_REF] proposes the R package SamplerCompare for comparing several MCMC's differing by a single tuning parameter, using standard evaluation criterion.

In this paper, we propose a methodological approach, and corresponding software tool, only based on Monte Carlo simulation (i.e. not requiring a theoretical study typically MCMC and/or target-specific) with two goals: (i) For MCMC users to easily select a good sampler among possible candidates; (ii) For researchers to better understand which (A)MCMC methods perform best in which circumstances. Let

H(p) := p log p = E p (log p) (1)
be the differential entropy of a probability density p over Θ, and p t be the marginal density of the (A)MCMC algorithm at "time" (iteration) t. Our approach is grounded on a criterion which is the evolution of the Kullback-Leibler divergence between p t and f ,

t → K(p t , f ) := p t log p t f = H(p t ) -p t log f.
This Kullback "distance" is indeed a natural measure of the algorithm's quality and has strong connections with ergodicity of Markov chains and rates of convergence, see [START_REF] Harremoes | Convergence of Markov chains in information divergence[END_REF] for recent results. In MCMC setup, [START_REF] Chauveau | Smoothness of Metropolis-Hastings algorithm and application to entropy estimation[END_REF] showed that if the proposal density of a Metropolis-Hastings algorithm satifies a uniform minorization condition implying its geometric convergence as in [START_REF] Holden | Geometric convergence of the Metropolis-Hastings simulation algorithm[END_REF], then K(p t , f ) also decreases geometrically.

In Section 2, we detail our approach which is methodological but relies on estimation techniques that have been proved to be theoretically consistent in simple cases like Gaussian unidimensional RWMH or independent samplers [START_REF] Chauveau | Smoothness of Metropolis-Hastings algorithm and application to entropy estimation[END_REF]. Our estimation of H(p t ) is grounded on the simulation of N parallel (iid) chains. In Section 3, we prove the consistency of our entropy estimate based on Nearest Neighbor (NN) estimate from [START_REF] Kozachenko | Sample estimate of entropy of a random vector[END_REF], adapted to our adaptive MCMC setup. Section 4 illustrates the good behavior of our criterion on synthetic multi-dimensional examples. These examples also allow us to show the difficulty due to a bias coming from the curse of dimension in nonparametric statistical estimation. In Section 5 we detail our methodological solution for handling that bias problem, in such a way that our approach being still usable even in large dimension, in practice for Bayesian models with dozens of parameters.

Entropy and Kullback estimation in MCMC context

Recent motivations for entropy estimation in other fields like molecular science appeared recently in the literature (see, e.g. [START_REF] Singh | Nearest neighbor estimate of entropy[END_REF], and are concerned by estimation of H(p) for multivariate densities p. Most of the estimation techniques proved to be consistent under various conditions are based on iid samples from p. There exists some results about entropy estimation for dependent sequences, but these heavily rely on the mixing properties of these sequences themselves, that are precisely what we want to capture by our simulationbased approach without theoretical investigations concerning mixing properties of the Markov kernel. More importantly, these approaches could be used to estimate H(f ) but cannot estimate H(p t ) for each t.

Simulation of iid copies of the (A)MCMC algorithm

Our approach is consequently based on the simulation of N parallel (iid) copies of (Adaptive) Markov chains started from a diffuse initial distribution p 0 and using the transition kernel defined by the MCMC strategy under investigation. The N chains, started from θ 0 1 , . . . , θ 0 N iid∼ p 0 , are denoted chain # 1 :

θ 0 1 → θ 1 1 → • • • → θ t 1 ∼ p t → • • • . . . . . . chain # N : θ 0 N → θ 1 N → • • • → θ t N ∼ p t → • • •
where "→" indicates the (eventually non-homogeneous) Markov dependence. At "time" (iteration) t, the locations of the N simulated chains θ t = (θ t 1 , . . . , θ t N ) forms a N -sample iid∼ p t . In an experimental framework where one wants to evaluate a new (A)MCMC algorithm the target f often corresponds to a benchmark example, hence is completely known (as e.g., in [START_REF] Vrugt | Accelerating markov chain monte carlo simulation by differential evolution with self-adaptive randomized subspace sampling[END_REF]. In this case a strongly consistent estimate of p t log f is given by Monte Carlo integration and the Strong Law of Large Numbers,

pt N (log f ) = 1 N N i=1 log f (θ t i ), (2)
so that estimation of K(p t , f ) is in turn accessible provided H(p t ) is. However, if the objective is to evaluate an experimental MCMC method for an actual Bayesian model for which f is a posterior density proportional to the likelihood, say f (•) ∝ φ(•) where the normalization constant is not known, only pt N (log φ) is accessible. We will see that this is not really a flaw since φ itself retains all the specificity (shape, modes, tails, . . . ) of f , and since we are mostly interested in the stabilization in t of K(p t , f ), not necessarily in knowing its limiting value, as will be detailed in Section 5. In addition, the normalization problem can be eliminated by comparing the MCMC under study to a benchmark MCMC algorithm (e.g., a gaussian RWMH) for the same target f . Indeed, considering two MCMC strategies leading to two sequences of marginal densities, say (p t 1 ) t≥0 and (p t 2 ) t≥0 allows the difference of the divergences to be accessible to estimation since

D(p t 1 , p t 2 , f ) = K(p t 1 , f ) -K(p t 2 , f ) = H(p t 1 ) -H(p t 2 ) + E p t 2 [log φ] -E p t 1 [log φ]. (3)
The Kullback criterion is the only usual divergence insuring this property and, in addition to its connection with ergodicity, it motivates our choice. Note also that the Kullback divergence is currently used as a criterion in other simulation approaches, see [START_REF] Douc | Convergence of adaptive mixtures of importance sampling schemes[END_REF]. The choice of this estimate also has the advantage of avoiding numerical integration in moderate or high dimensional spaces (replaced by Monte Carlo integration), in contrary to other criterion such as the L 1 -distance.

For estimating the entropy H(p t ) a classical, plug-in approach, is to build a nonparametric kernel density estimate of p t , and to compute the Monte Carlo integration of this estimate. Techniques based on this approach have been suggested by [START_REF] Ahmad | A nonparametric estimation of the entropy for absolutely continuous distributions[END_REF], and studied by many authors under different assumptions (see, e.g., the survey paper [START_REF] Beirlant | Nonparametric entropy estimation, an overview[END_REF]. Several consistency and asymptotic normality results pertaining to this approach have been proved (see references in [START_REF] Eggermont | Best asymptotic normality of the kernel density entropy estimator for smooth densities[END_REF]. However, most of these are not suitable to estimate H(p t ) even in the simplest MH cases, either because they do not apply to multivariate densities, or because they require smoothness conditions that are far too restrictive to be proved for the sequences of densities p t we have to consider here. Up to our best knowledge, the unique consistency result applicable in this MCMC simulation context is the one proved in [START_REF] Györfi | An entropy estimate based on a kernel density estimation[END_REF], that essentially requires a Lipschitz type smoothness condition. Indeed, for that approach, [START_REF] Chauveau | Smoothness of Metropolis-Hastings algorithm and application to entropy estimation[END_REF] have proved that adequate smoothness and tail conditions on the "input ingredients" of the MH algorithm (namely p 0 , q and f ) propagate a Lipschitz condition to the successive marginals p t , t = 1, . . . , n, so that the sequence of (H(p t )) t=1,...,n can be consistently estimated. These technical conditions have been proved to hold in simple univariate IS and RWMH cases, but are not meant to be verified in general, since it would require tedious (and often unfeasible) calculations.

Estimates based on nearest neighbor distances

The plug-in estimate presented above requires the tuning of several parameters: a certain threshold for truncating the data over the tails of p t , the choice of the kernel and the difficult issue of the appropriate bandwidth matrix, particularly in high dimensions. All these issues motivated us to find an alternative, and study the behavior of the somehow simpler Nearest Neighbor (NN) estimate initiated by [START_REF] Kozachenko | Sample estimate of entropy of a random vector[END_REF] (see also [START_REF] Beirlant | Nonparametric entropy estimation, an overview[END_REF], for a survey on these entropy estimates). In our setup, based on the sample θ t iid∼ p t in dimension d, this NN entropy estimate is

ĤN (p t ) = 1 N N i=1 log(ρ d i ) + log(N -1) + log(C 1 (d)) + C E , (4)
where

C E = - ∞ 0 e -u log u du ≈ 0.5772 . . . is the Euler constant, C 1 (d) = π d/2 Γ(d/2+1
) and where

ρ i = min{ρ(θ t i , θ t j ), j ∈ {1, 2, . . . , N }, j = i}
is the (Euclidean) distance ρ(•, •) from the ith point to its nearest neighbor in the sample θ t . [START_REF] Kozachenko | Sample estimate of entropy of a random vector[END_REF] proved the mean square consistency of the NN entropy estimate ĤN (h) as in (4) for a density h and any dimension d under the following Peak and Tail (P&T) conditions.

Definition 1. A density h over R d satisfies the P&T conditions if there exists α > 0 such that

| log h(x)| 2+α h(x) dx < +∞ (5) | log ρ(x, y)| 2+α h(x)h(y) dx dy < +∞. ( 6
)
This NN estimate seems more appealing than kernel density estimates in our situation, both from an operational point of view (no tuning parameters like the threshold and bandwidth), and from a computational point of view (the nearest distance can be computed faster than a multivariate kernel density estimate in high dimension). Until recently, this nearest neighbor approach has been used and studied mostly in univariate or bivariate (d = 2) situations, like in image processing. One interest of this study is to investigate its behavior in higher dimensions and for MCMC sequences of marginals.

Peak and tail conditions in adaptive MH case

In this section we introduce a generic class of Adaptive Metropolis-Hastings (AMH) algorithms to which most of the AMCMC strategies proposed in the recent literature belong (see, e.g., [START_REF] Andrieu | A tutorial on adaptive MCMC[END_REF]. We prove that, under certain uniform control conditions, the P&T conditions from Definition 1 required to estimate the entropy of successive densities p t by the NN approach from [START_REF] Kozachenko | Sample estimate of entropy of a random vector[END_REF], hold for each fixed iteration t ∈ N.

As recalled in the introduction, an AMCMC algorithm relies on a family of Markov kernels based on a joint process (θ t , ϑ t ) t≥0 such that the conditional distribution of θ t+1 given the information available up to time t is a kernel P ϑ t where ϑ t is a Euclidean parameter depending on the past. In the case of a generic Adaptive MH processes (X t ) t≥0 valued in Ω ⊆ R d , each MH step at time t is based on the generation of the proposed next move y from an adapted proposal density q ϑt (y) ∈ F := {q ϑ |ϑ ∈ Θ}, where ϑ t := ϑ(x t 0 ) is a strategically tuned parameter possibly integrating the whole past trajectory denoted x t 0 = (x 0 , . . . , x t ). For a starting value x 0 ∼ p 0 , the t-th step x t → x t+1 of the AMH algorithm is as follows:

1. generate y ∼ q ϑt (•) 2. compute α ϑt (x t , y) = min 1, f (y)q ϑt (x t ) f (x t )q ϑt (y)
3. take x t+1 = y with probability α ϑt (x t , y)

x t with probability 1 -α ϑt (x t , y).

The proposition below gives the convergence of our NN entropy estimates for the successive AMH marginal densities.

Proposition 1. If there exist nonnegative functions (ϕ 1 , ϕ 2 ) both defined on Ω, a constant a ∈ (0, 1) and α > 0 such that:

i) C 1 = ϕ 1 (x)dx < ∞, C 2 = ϕ 2 1 (x) ϕ 2 (x) dx < ∞, and C 3 = ϕ 2 (x)dx < ∞ ii) ϕ 1 ≤ p 0 ≤ ϕ 2 , and ϕ 1 ≤ f iii) af ≤ q ϑ ≤ ϕ 2 for all ϑ ∈ Θ iv) log C ϕ 2 1 (x) ϕ 2 (x) 2+α ϕ 2 (x)dx < +∞ for any constant C > 0 v) |log (Cϕ 2 (x))| 2+α ϕ 2 (x)dx < +∞ for any constant C > 0 vi) | log ρ(x, y)| 2+α ϕ 2 (x)ϕ 2 (y)dxdy < +∞,
then the successive densities of the Adaptive MH algorithm described above satisfy the P&T conditions (5)-( 6).

Note that the constant a introduced above can be understood as the minorization constant used in, e.g., [START_REF] Holden | Geometric convergence of the Metropolis-Hastings simulation algorithm[END_REF] and [START_REF] Mengersen | Rates of convergence of the Hastings and Metropolis algorithms[END_REF] conditions of geometric ergodicity.

Proof. For all t ≥ 0, we define P ϑt (x t , •), the generic adaptive transition kernel depending on ϑ t = ϑ(x t 0 ): P ϑt (x t , dy) = q ϑt (y)α ϑt (x t , y))dy + 1 -q ϑt (z)α ϑt (x t , z)dz δ x t (dy).

We denote as before by p t the marginal density of the AMH algorithm at iteration t. Define first the two nonnegative functions controlling p t from Lemma 1 in Appendix A.2, Equations ( 15) and ( 16) using conditions (i)-(iii) of Proposition 1:

A t (x) := a 2t C 1 C t-1 2 ϕ 2 1 (x) ϕ 2 (x) B t (x) := 2(C 3 + 1) t-1 ϕ 2 (x) A t (x) ≤ p t (x) ≤ B t (x), t ≥ 1.
To prove that p t , for a fixed step t, satisfies the first P&T condition (5) from Definition 1, set the domain

D t = {x ∈ R d : p t (x) < 1}, Dt = R d \ D t , and let β = (2 + α). We have | log p t (x)| β p t (x) dx = D t | log p t (x)| β p t (x) dx + Dt | log p t (x)| β p t (x) dx. (7) Since for x ∈ D t , 0 ≤ | log p t (x)| = -log p t (x) ≤ -log A t (x) = | log A t (x)|, D t | log p t (x)| β p t (x) dx ≤ D t | log A t (x)| β B t (x) dx.
This last integral is finite from condition (iv) since

D t | log A t (x)| β B t (x) dx ≤ C log C ϕ 2 1 (x) ϕ 2 (x) β ϕ 2 (x) dx < ∞,
where C = a 2t C 1 C t-1 2 and C = 2(C 3 + 1) t-1 are both finite from condition (i). For the rightmost integral in (7), it suffices to note that

Dt | log p t (x)| β p t (x) dx ≤ Dt | log B t (x)| β B t (x) dx ≤ C | log( Cϕ 2 (x))| β ϕ 2 (x) dx,
where the rightmost integral is finite from condition (v). The P&T condition (6) for p t is obtained straightforwardly from condition (vi).

Experiments and simulations

All the estimation techniques and MCMC evaluation criterion presented in the previous Sections are based on intensive simulations and computations for which we provide a software tool implemented in the EntropyMCMC package for the R statistical software (R Core Team, 2013), taking advantage of recent advances in High Performance Computing, that will be publicly available in a near futur. This package includes some predefined target distributions and standard MCMC samplers, easy definition of additional ones, functions for running simulations, estimating entropy and Kullback divergences, results visualization and sampler comparison. For instance, Figs 2 and 3 have been done using simply a default plot() command from this package. The parallel simulations can be done from inside the package, or imported from external files.

Multidimensional Gaussian target density and iid sampler

Several authors, mostly from biology, statistics and information theory have recently shown some evidence that estimation of functionals like the entropy suffers from the curse of dimensionality (see for instance, [START_REF] Stowell | Fast multidimensional entropy estimation by k-d partitioning[END_REF] and [START_REF] Sricharan | Ensemble estimators for multivariate entropy estimation[END_REF]). In these studies, the bias seems much more affected by the dimension than the variance.

To confirm this and also validate our software tool, we ran some experiments in Gaussian situations. Let us denote by N µ,Σ d the pdf of a d-multivariate Gaussian N d (µ, Σ) with mean vector µ and covariance matrix Σ. One reason for choosing Gaussian targets is that the true entropy is known in this case,

H(N µ,Σ d ) = - d(log(2π) + 1) + log(det(Σ)) 2 . ( 8
)
In this section and in Section 5, we use centered Gaussian distributions with spherical covariances matrices, and then simply denote by

N σ d the pdf of N d (01 d , σ 2 I d )
, where 1 d is the unit column vector of size d, and I d is the d × d identity matrix. For these experiments we coded a "fake" MCMC algorithm, i.e. we defined as a MCMC algorithm a simple iid sampler from a completely known target such as N σ d . In other words, p t ≡ N σ d for t ≥ 1, so that running n iterations of N iid "chains" from this algorithm corresponds exactly to simulating n replications of N iid observations ∼ N σ d , since there is no dependence with time. Here, we try our entropy estimate on the simplest possible target: the pdf of the standard Gaussian N 1 d . This allowed us to compare the entropy estimate based on the Monte-Carlo (MC) integration of the known log(N 1 d ) from ( 2) and the NN estimates of H(N 1 d ) from ( 4). In addition to knowing the true entropy from (8), another advantage is that we can reasonably assume that the gaussian density p t ≡ N 1 d satisfies the P&T conditions required by the NN entropy estimates (see a proof in Appendix A.3 for the univariate case). Our results also show numerical evidence of consistency.

We ran this model for dimensions between 3 and 30, and N ranging from 500 to 30, 000 observations. Examples of typical results for d = 20 are provided in Table 1. We can see that the bias of the NN estimation is much more affected by the dimension than the variance. Remember that it is "unfair" to compare the MC estimator with the NN estimator, since the former is a simple application of the Law of Large Numbers for a known log(N 1 d ), whereas in the latter the density itself is estimated from the sample. We also observed that in this case, the bias is negative for large d, leading to the under-estimation of H(N 1 d ). Bias considerations in actual (A)MCMC setups will be discussed more deeply in Section 5.

Multidimensional Gaussian mixtures

We illustrate here our methodology on a synthetic but more complex target density: a 3components mixture of multivariate, d-dimensional Gaussian distributions,

f (x) = 3 j=1 λ j N µ j ,Σ j d (x),
The three weights are set to λ j = 1/3, the mean vectors and (spherical ) covariance matrices are set to (with the notations introduced in Section 4.1)

µ 1 = 01 d , µ 2 = 41 d , µ 3 = -41 d , Σ j = jI d , j = 1, 2, 3.
The advantage of such a synthetic model is that it is defined for any dimension, but the complexity of the target increases with d (the distance between modes increases and the modes get more and more separated and spiked). Note that here the normalizing constant of f is known, so that theoretically the sequence of marginals p t from a proper (converging) MCMC satisfies K(p t , f ) → 0 as t → ∞. We compare several standard, well known MCMC algorithms for recovering this target:

• Two RWMH with Gaussian proposals N d (θ t ; σ 2 I d ) resulting in slow or fast MCMC's depending on the magnitude σ 2 of their (spherical) variance matrix that we set to σ 2 = 1 (algorithm called RW1 in the sequel) and 4 (RW4).

• Three Independence Samplers (IS) with Gaussian proposals N σ d and the choices σ 2 = 2 (IS2 in the sequel), σ 2 = 9 (IS9) and 16 (IS16).

The idea driving the choices above for the tuning parameters of the candidate MCMC's is that we want to compare fast, slow and even MCMC's not converging in practice (i.e., in a feasible amount of iterations). For this Gaussian mixture target in small dimensions (d ≤ 2), it is easy to figure out how to obtain such algorithms. For instance, Fig. 1 displays the target and five proposal densities for the one-dimensional case. For the IS's, the practical support of the proposal density, within which each proposed next move lies, is suppose to include the region of interest of the target. This is definitely not the case for σ 2 = 2, hence IS2 can be viewed in practice as a non converging algorithm. Larger variances should lead to better algorithms, but it is not easy to tell which value of σ 2 is the best choice. Also, the IS converges geometrically if the proposal have heavier tails than the target, which is not the case here.

Our method then shows how these strategies behave in higher dimension, where the three components get more and more separated and spiked. We detail some results of our experiments, associated to Figs 2 and3.

Fig.2:

In small dimensions like d = 2 (top panel), our criterion delivers the right answer straightforwardly, since there is no noticeable bias in the estimates, even with N = 500 chains. All the convergent MCMC's stabilized well before n = 1000 iterations. IS2 is not converging and is almost stabilized to a non-zero value in these n = 1000 iterations. Similar runs for more iterations show clearly its non convergence, and increasing N up to say 1000 reduces the variance, resulting in easier-to-read diagnostics. The bottom panel displays a similar experiment, but now in dimension d = 10, and for N = 5000 chains to illustrate the variance reduction resulting is smooth curves (so many chains are not needed to get a readable diagnostic). The two RW's converge, but note that RW4 is less performant than RW1 in this higher dimension. IS2 is stabilized away from zero, indicating non convergence. It is hard to tell in 1000 iterations wether IS9 and IS16 will converge at some point, but 1000 iterations are sufficient to tell that these are slow, comparing with the RW's fast stabilizations. Actually running the experiment up to n = 10, 000 show that IS16 is slow but converging, whereas IS9 convergence is not clear. This is because the proposal "almost covers" the region of interest of the mixture.

Fig. 3: As expected, the scale of the number of iterations required to detect stabilization increases with the difficulty associated to the dimension, hence in this figure two simulations in the d = 20 dimensional case have been ran up to n = 10, 000 iterations. The difference between top and bottom panels is just the number of iid chains, N = 500 (top) and N = 10, 000 (bottom). This last number has been chosen intentionnally huge to illustrate the variance reduction, bias difficulty, and the fact that running so many chains is feasible but not needed to obtain a meaningful criterion. Note that this example with d = 20, n = 10, 000 and N = 500 only requires about about 25mn of CPU time per algorithm on a 12-cores single workstation. The same example in dimension d = 30 required 32mn of CPU time.

One purpose of Fig. 3 is precisely to illustrate the bias problem. Indeed, if we look at the top panel and the three IS's only, we conclude that IS16 and IS9 are non or very slowly converging, and that IS2 quickly stabilizes near (and above) 0. Hence if we were only comparing algorithms for stabilization near 0, we would falsely conclude that IS2 is the preferable algorithm. But we already know that IS2 is not a convergent MCMC in this synthetic example. This bizarre results is due to a bias in H(p t ) estimation. Looking now at the curves for the RW's, we see that these stabilize on a common negative value, which is theoretically impossible since K(p t , f ) ≥ 0. Hence there is a negative bias, more prominent in the top panel with just N = 500 chains. Looking at the bottom panel, we see that estimating the entropy from samples of size N = 10, 000 has the effect of reducing the variance giving more accurate curves, but only slightly reducing the bias: all the curves are just shifted from a small amount (enough to tell that IS2 is actually not convergent). The difference between these two plots illustrate the fact that the bias reduces dramatically slowly with N , and that since all the stabilization values are biased, we cannot rely on stabilization near 0 to assess convergence or lack of convergence.

There are two reasons allowing us to evaluate the competing algorithms in view of Fig. 3 top: (i) the property of the Kullback divergence says that K(p t , f ) decreases when p t → f , so that if the bias is of the same order in all the algorithms, the faster decay and smaller stabilization value is associated to the best algorithm; (ii) we have a prior knowledge that the RW's are convergent MCMC's, without knowing their rates of convergence. The two RW's can be viewed here as convergent benchmarks, and since they stabilize at the same value, it means that their biases are quite identical. To summarize, we can conclude that RW1 is preferable, RW4 is convergent but slower (probably as a consequence of the more spiked component modes so that its proposal variance is comparatively getting too large), IS2 is non convergent since it stabilizes above the benchmark RW1, and IS9 and IS16 may be convergent but even much slower. This diagnostic comes from the analysis of the top panel only, i.e. N = 500 is enough to get a conclusion. The next section details all these methodological questions.

How to handle the bias in large dimension

The experiments from Section 4 show that some care shoud be taken in the analysis of the plots of the estimates of t → K(p t , f ) delivered by our techniques, particularly in high dimension (say d ≥ 10) where the bias becomes visible. This effect of the dimension on the bias has actually been already noticed in recent literature since nowadays applications of entropy estimation in other fields require moderate to high dimensions. Our results are in accordance with, for instance, [START_REF] Stowell | Fast multidimensional entropy estimation by k-d partitioning[END_REF] and [START_REF] Sricharan | Ensemble estimators for multivariate entropy estimation[END_REF]. These studies show that in H(p) estimation, one can expect the variance to decreases as O(N -1 ) whereas the bias only decreases as O(N -1/(d+1) ), which these authors called a "glacially slow" rate, and this phenomenon occurs both for Kernel density and NN-based estimates. These behaviors has been confirmed in our case using an iid sampler for a Gaussian target of known entropy (Section 4.1). Hence, trying to achieve in practice the asymptotic unbiasedness guaranteed by the theory by "just" increasing N is hopeless when d gets large.

Another difficulty comes from the unknown normalizing constant. In our experiment the target f in the Monte-Carlo estimate (2) was entirely known, but in practical situations like Bayesian inference for a parameter θ, f is a posterior distribution only known up to a multiplicative constant, f (θ) = Cφ(θ) where C = ( φ(θ) dθ) -1 is the (unknown) normalizing constant. Hence what can be actually estimated by our method is K(p t , φ), and we have

K(p t , φ) = H(p t ) -E p t (log φ) = K(p t , f ) + log C.
This is why in (3) we noticed that in the (estimate of the) difference between the Kullback divergences issued from two MCMC strategies with marginal densities p t 1 and p t 2 , the unknown log C cancels out and D(p t 1 , p t 2 , f ) → 0 as t → ∞ if both strategies are converging.

More precisely, consider first the asymptotic in t, i.e. along the iterations of the algorithm. Define by p t b the successive marginals of a converging MCMC algorithm (e.g., a benchmark algorithm), and by p t n the successive marginals of a non converging MCMC. This means that, as t → ∞, p t b → f and p t n → g = f so that

K(p t b , φ) → log C (9) K(p t n , φ) → K(g, f ) + log C > log C. (10)
Consider now the estimation KN (•, φ) based on a N iid copies of each MCMC. For d small (say d < 10), the consistency result applies so that, for any fixed t ≥ 1,

KN (p t , φ) → N →∞ K(p t , φ).
Thus (for N large enough) our criterion provides a graphical tool comparing convergent MCMC's from the decays of their KN (p t , φ)'s. But it can also be used as a convergence assessment tool due to ( 9)-( 10) which imply that a non-converging MCMC will stabilize above a converging one. For large d, each Kullback estimate is biased. As seen previously this bias, which comes from the estimation of the entropy H(p) of any d-dimensional pdf p, depends on p and d and will be denoted bias N (p, d). We can sketch the behavior of the estimates for fixed but large enough N so that the variance becomes neglictible but the bias still remains, which is what is achievable in practice for d large. Informally, for fixed t, the estimate stabilizes at the theoretical value plus this bias,

KN (p t , φ) ≈ K(p t , f ) + log C + bias N (p t , d).
We assume that bias N (p t , d) → bias N (f, d) when p t → f , which seems fairly reasonnable and implies that bias N (p t , d) ≈ bias N (f, d) for p t and f "close enough". This is needed in order to compare decays of t → K(p t , f ) before convergence, and is supported by numerical evidence from the experiments in Section 4 and other experiments detailed below. Intuitively, even a slow MCMC, if converging, should never been too far from its target, leading to similar biases. Finally, for our two strategies above, we have for both N and t large enough:

KN (p t b , φ) ≈ log C + bias N (f, d) (11) KN (p t n , φ) ≈ K(g, f ) + log C + bias N (g, d). (12)
Unfortunately, we cannot assume in this case that K(g, f ) + bias N (g, d) > bias N (f, d) without additionnal assumptions on g and f ; in particular the biases can be of any sign. We thus provide below an additional experiment which purpose is to evaluate, in a Gaussian situation and for several choices of distinc distributions f and g, the typical biases and K(g, f ) we can expect. The idea is to check wether, for f = g, typical biases may be such that bias N (g, d) ≈ bias N (f, d), and both smaller in comparison with K(g, f ).

Let N σ d denote the pdf of the d-dimensional centered multivariate Gaussian with spherical covariance matrix of diagonal variances σ 2 as in Section 4.1. For σ = 1, the target density N 1 d corresponds to the standard multivariate Gaussian for which bias N (N 1 d , d) has been investigated in Table 1 from a Monte-Carlo experiment only using the code provided in the EntropyMCMC package. The same code can be used similarly to estimate any H(N σ d ). Moreover, the Kullback divergence between two multivariate Gaussian distribution is known and available in closed form: For N

µ j ,Σ j d = N d (µ j , Σ j ), j = 1, 2, we have K(N µ 2 ,Σ 2 d , N µ 1 ,Σ 1 d ) = 1 2 Tr Σ -1 1 Σ 2 + (µ 1 -µ 2 ) T Σ -1 1 (µ 1 -µ 2 ) -d -log det(Σ 2 ) det(Σ 1
) .

(13) Our experiments only involve centered Gaussian for simplicity, and compare K(g, f )+bias N (g, d) against bias N (f, d). Since the entropy H(N µ,Σ d ) is not affected by the mean, but the Kullback divergence is, differences between means would result in even larger K(g, f )'s.

Tables 2 and3 show the typical results we obtain in the Gaussian case, for dimensions 20 and 40. In these tables, the case σ = 1 corresponds to the target f ≡ N 1 d , and the other Gaussian pdf's with larger variances to densities g = f . The column rel.bias gives the relative bias of the entropy estimates

( ĤN (N σ d ) -H(N σ d )) H(N σ d )
.

It is clear that for any fixed dimension, the biases are similar and neglictible in comparison with the Kullback distances. For instance, in the case d = 40, g ≡ N 5 d and running our criterion for N = 5000, we obtain (Table 3, row 5) bias N (f, d) = -3.3139 almost equal to bias N (g, d) = -3.3152, whereas K(g, f ) = 23.9. To summarize, in all the experiments we did, bias N (g, d) ≈ bias N (f, d), and both smaller in comparison with K(g, f ) even for small differences between f and g's.

To summarize, in large dimension, our experiments provide numerical evidence that biases in the entropy estimates are of the same order for different densities with roughly the same shape (this means for convergent p t b ≈ f but also for non convergent p t n ≈ g = f ) and neglictible in comparison with K(g, f )'s. Our criterion based ultimately on ( 11)-( 12) can thus compare efficiency of convergent (A)MCMC's, and detect non-convergent algorithms.

Conclusion

We have proposed a methodological approach to evaluate (A)MCMC efficiency and control of convergence on the basis of intensive simulation only. The diagnostic is based on a practical, easy-to-understand graphical criterion. To evacuate the difficulty induced by the bias in high dimensions we have introduced a benchmark convergent MCMC which indicates when stabilization means convergence. Since our method requires intensive simulations that may be computationally demanding, all our algorithms have been implemented in a package named EntropyMCMC for the R statistical software (R Core Team, 2013) that will be publicly available in a near futur. This package takes advantage of recent advances in parallel, High Performance Computing (HPC) using the Rmpi package [START_REF] Yu | Rmpi: Interface (Wrapper) to MPI (Message-Passing Interface[END_REF]. All the examples shown in this paper have been ran with this package on multicore workstations and the regional cluster CCSC1 . These simulations (or part of it) from the best sampler are recyclable after comparisons, or can be re-used in the fly for statistical inference.

Theoretically, we have shown that the peak and tail conditions required for the NNestimates consistency are satisfied by the successive densities of a generic class of MH including Adaptive MH MCMC algorithms, under uniform conditions controlling the adaptation. Like in [START_REF] Chauveau | Smoothness of Metropolis-Hastings algorithm and application to entropy estimation[END_REF] for a Kernel density based entropy estimate, the difficulty often comes from the fact that a MH kernel has a point mass at the chain's current position. We have then assumed that these conditions were satisfied for the MCMC algorithms used in our experiment, and this have been supported by numerical evidence of consistency.

Recent researchs extend the NN idea to a k-th nearest neighbor distance estimate [START_REF] Singh | Nearest neighbor estimate of entropy[END_REF], see also [START_REF] Wang | Divergence estimation for multidimensional densities via k-nearest-neighbor distances[END_REF]. There are some hope that these extensions together with recent computing strategies for computing approximate k-NN estimates reduce the bias in entropy estimation. However, it also brings back a tuning parameter (how to choose k) that plays somehow the role of the bandwidth in the kernel density estimation approach. These considerations are perspective for futur work.

A Appendix: Controlling the successive marginals

We provide in this Appendix some technical results which allow us to control the successive marginals of some generic AMH algorithms, in order to prove the P&T conditions ( 5)-( 6) in Proposition 1.

A.1 The MH independence sampler case

To help intuition, we start here by showing how successive marginals of a simple independent MH sampler can be controlled using the assumptions (i)-(iii) of Proposition 1 (where q does also simply not depends on the past). Recall that "independent" here means that the proposal density q does not depend on the current position of the chain. Let us denote the probability of accepting the move y from x, α(x, y) = min 1, f (y)q(x) f (x)q(y) .

Then

p 1 (y) = p 0 (x)q(y)α(x, y)dx + p 0 (y)q(z)(1 -α(y, z))dz ≥ aϕ 1 (y) ϕ 1 (x)α(x, y)dx .
We have also

α(x, y) ≥ min 1, af (y) q(y) ≥ min 1, aϕ 1 (y) ϕ 2 (y) = a ϕ 1 (y) ϕ 2 (y)
since aϕ 1 ≤ q ≤ ϕ 2 . This leads to

p 1 (y) ≥ a 2 ϕ 2 1 (y) ϕ 2 (y) ϕ 1 (x)dx = a 2 ϕ 2 1 (y) ϕ 2 (y) C 1 .
Iterating, we have

p 2 (y) ≥ q(y) p 1 (x)α(x, y)dx ≥ a 4 C 1 ϕ 2 1 (x) ϕ 2 (x) dx ϕ 2 1 (y) ϕ 2 (y) . (14)
By induction we prove that

p t (y) ≥ q(y) p t-1 (x)α(x, y)dx ≥ a 2t ϕ 1 (x)dx . ϕ 2 1 (x) ϕ 2 (x) dx t-1 ϕ 2 1 (y) ϕ 2 (y) = a 2t C 1 C t-1 2 ϕ 2 1 (y) ϕ 2 (y) .
To majorize p 1 (y) we can simply notice that p 1 (y) ≤ q(y) + p 0 (y) ≤ 2ϕ 2 (y) and iterate to get p t (y) ≤ (t + 1)ϕ 2 (y). However this will not hold in the adaptive case.

A.2 The Adaptive MH (AMH) case

We turn now to the case of the AMH generic algorithm defined in Section 3. For more obvious notations, we will not use the common description of an adaptive MCMC algorithm through a Markov kernel indexed by ϑ t = ϑ(x t 0 ) as we did previously, but directly by the trajectory from all the past x t 0 to indicate dependence. Lemma 1. Let (ϕ 1 , ϕ 2 ) be nonnegative functions satisfying conditions (i)-(iii) of Proposition 1, and q x t 0 (y) be an adaptive proposal density depending on the past such that af ≤ q x t 0 ≤ ϕ 2 for any x t 0 ∈ Ω t+1 . Then, for all y ∈ Ω,

a 2t C 1 C t-1 2 ϕ 2 1 (y) ϕ 2 (y) ≤ p t (y) (15) and p t (y) ≤ 2(C 3 + 1) t-1 ϕ 2 (y), (16)
where the constants a, C 1 , C 2 , C 3 , are defined in Proposition 1.

Proof. For all t ≥ 1, we define the generic AMH transition a kernel depending on the past x t-1 0 = (x 0 , . . . , x t-1 ):

P x t-1 0 (x t-1 , dy) = q x t-1 0 (y)α x t-1 0 (x t-1 , y)dy + q x t-1 0 (z) 1 -α x t-1 0 (x t-1 , z) dz δ x t-1 (y) dy (17)
where

α x t-1 0 (x t-1 , y) = min 1, f (y)q x t-1 0 (x t-1 ) f (x t-1 )q x t-1 0 (y)
is the probability of accepting the move y from x t-1 in the MH step. We handle first the minorization part (15). The technique is similar to the simplest independence sampler case of Appendix A.1, except that here we need to minorize the transition kernel itself as follows:

P x t-1 0 (x t-1 , dy) ≥ q x t-1 0 (y)α x t-1 0 (x t-1 , y) dy.
Similarly to the independent sampler case we have:

α x t-1 0 (x t-1 , y) ≥ min 1, af (y) q x t-1 0 (y)
≥ a ϕ 1 (y) ϕ 2 (y) , which implies P x t-1 0 (x t-1 , dy) ≥ a 2 ϕ 2 1 (y) ϕ 2 (y) dy.

Proceeding in that way we have the following minorization for the densities:

p t (y) dy = p 0 (x 0 ) dx 0 P x 0 (x 0 , dx 1 )P x 1 0 (x 1 , dx 2 ) . . . P x t-1 0 (x t-1 , dy)

≥ a 2t ϕ 1 (x 0 ) ϕ 2 1 (x 1 ) ϕ 2 2 (x 1 ) . . . ϕ 2 1 (x t-1 ) ϕ 2 (x t-1 ) dx t-1 0 ϕ 2 1 (y) ϕ 2 (y) dy = a 2t C 1 C t-1 2 ϕ 2 1 (y) ϕ 2 (y) dy.

To obtain the majorization (16) of the densities, we notice from (17) that: P x t-1 0 (x t-1 , dy) ≤ q x t-1 0 (y)dy + δ x t-1 (y) dy ≤ ϕ 2 (y)dy + δ x t-1 (y) dy = Φ(x t-1 , dy), where Φ(x, dy) := ϕ 2 (y)dy + δ x (y) dy is a non-normalized transition kernel, i.e. Φ(x, dy) = C 3 + 1. This leads to p t (y) dy = p 0 (x 0 ) dx 0 P x 0 (x 0 , dx 1 )P x 1 0 (x 1 , dx 2 ) . . . P x t-1 0 (x t-1 , dy) ≤ p 0 (x 0 ) dx 0 Φ(x 0 , dx 1 )Φ(x 1 , dx 2 ) . . . Φ(x t-1 , dy).

We can now study separately the right hand side term of the above inequality. For the first step we have:

p 1 (x 1 ) dx 1 = p 0 (x 0 ) dx 0 P x 0 (x 0 , dx 1 ) ≤ p 0 (x 0 ) dx 0 Φ(x 0 , dx 1 ) = p 0 (x 0 ) ϕ 2 (x 1 )dx 1 + δ x 0 (x 1 ) dx 1 dx 0 = ϕ 2 (x 1 ) p 0 (x 0 )dx 0 dx 1 + p 0 (x 0 )I {x 1 } (x 0 ) dx 0 dx 1 ≤ ϕ 2 (x 1 ) dx 1 + p 0 (x 1 )dx 1 ≤ 2ϕ 2 (x 1 ) dx 1 .

Similarly for the second step (the integrals being w.r.t. dx 0 and dx 1 ), p 2 (x 2 ) dx 2 = p 0 (x 0 )P x 0 (x 0 , dx 1 ) dx 0 P x 1 0 (x 1 , dx 2 ) ≤ 2ϕ 2 (x 1 ) dx 1 Φ(x 1 , dx 2 ) ≤ ϕ 2 (x 1 ) dx 1 2ϕ 2 (x 2 ) dx 2 + 2ϕ 2 (x 2 ) dx 2 = 2(C 3 + 1)ϕ 2 (x 2 ) dx 2 .

So that, by induction, p t (x t ) dx t ≤ 2(C 3 + 1) t-1 ϕ 2 (x t ) dx t .

This bound degenerates as t → +∞ but it is finite for each fixed iteration t.

A.3 P&T conditions in the Gaussian case

We check here that the P&T conditions are satisfied when f is the pdf of the standard normal N (0, 1). P&T condition (5) is obvious. For P&T condition (6), let

I

Figure 1 :

 1 Figure 1: Target density (bold grey), IS proposal densities (solid lines), and RWHM proposal densities (dashed lines, located on the arbitrary current position * , for d = 1.

  := R×R log(|x -y|) 2+α f (x)f (y) dx dy. (18)We remark thatR×R | log(|x -y|)| 2+α e -1 2 (|x| 2 +|y| 2 ) dxdy = I 1 + I 2 ,whereI 1 := R |x-y|<ε | log(|x -y|)| 2+α e -1 2 (|x| 2 +|y| 2 ) dy dx,andI 2 := R |x-y|≥ε | log(|x -y|)| 2+α e -1 2 (|x| 2 +|y| 2 ) dx.I 2 is trivially convergent. For I 1 , by the change of variable u = y-x and a symmetry argument we have:I 1 = R |u|<ε | log(|u|)| 2+α e -1 2 (|x| 2 +|x-u| 2 ) du dx ≤ R |u|<ε | log(|u|)| 2+α du e -1 2 (|x| 2 ) dx < ∞,where the last integral in convergent because, for ε < 1/e, |u|<ε | log(|u|)| 2+α du = 2

Figure 2 :

 2 Figure 2: Kullback estimates t → KN (p t , f ) for n = 1000 iterations of the 5 MCMC strategies RW1 (black), RW4 (red), IS2 (green), IS9 (blue), IS16 (light blue). Top: d = 2 and N = 500 chains; bottom: d = 10 and N = 5000 chains.

Figure 3 :

 3 Figure 3: Kullback estimates t → KN (p t , f ) for d = 20 and n = 10, 000 iterations of the 5 MCMC strategies RW1 (black), RW4 (red), IS2 (green), IS9 (blue), IS16 (light blue). Top: N = 500 chains; bottom: N = 10, 000 chains.

Table 1 :

 1 100×Bias and standard dev's for estimation of H(N 1 d ) over n = 100 replications of samples of size N , for the d = 20 Gaussian target. Here the true H(N 1 d ) = -28.38.

	N	MC bias MC sd	NN bias	NN sd
	500	1.6296	13.4395 -101.0047 15.2934
	1000	-0.6200 10.4458 -88.8441 12.1195
	5000	-0.0637	4.4483	-60.6133	5.4353
	10,000 -0.4023	3.0028	-52.0266	3.3988
	20,000 -0.0048	2.2053	-44.0156	2.6230
	30,000	0.3701	1.9592	-39.2813	2.2468

Table 2 :

 2 True entropies for some target pdf N σ d and Kullback distances to the reference N 1 d , together with estimated biases and standard deviations (sd) of the NN entropy estimate based on n = 100 replications of N parallel chains, in dimension d = 20.

	N	σ 2 H(N σ d ) K(N σ d , N 1 d ) bias N (N σ d , d) rel.bias(%) NN sd
		1	-28.4	0	-0.8652	3.0489	0.1097
	1000	5	-44.5	23.9	-0.8685	1.9528	0.1160
		10	-51.4	66.97	-0.8913	1.7338	0.1069
		1	-28.4	0	-0.6061	2.1356	0.0518
	5000	5	-44.5	23.9	-0.6099	1.3713	0.0558
		10	-51.4	66.97	-0.6019	1.1709	0.0537
		1	-28.4	0	-0.5084	1.7914	0.0376
	10,000 5	-44.5	23.9	-0.5138	1.1553	0.0397
		10	-51.4	66.97	-0.5207	1.0129	0.0364

Table 3 :

 3 True entropies for some target pdf N σ d and Kullback distances to the reference N 1 d , together with estimated biases and standard deviations (sd) of the NN entropy estimate based on n = 100 replications of N parallel chains, in dimension d = 40.

	N	σ 2 H(N σ d ) K(N σ d , N 1 d ) bias N (N σ d , d) relbias(%) NN sd
		1	-56.8	0	-4.0243	7.09	0.1823
	1000	5	-88.9	47.8	-3.9921	4.49	0.1644
		10	-103	133.9	-4.0222	3.91	0.1636
		1	-56.8	0	-3.3139	5.8387	0.0781
	5000	5	-88.9	47.8	-3.3152	3.7272	0.0765
		10	-103	133.9	-3.3147	3.2242	0.0749
		1	-56.8	0	-3.0610	5.3932	0.0500
	10,000 5	-88.9	47.8	-3.0784	3.4610	0.0448
		10	-103	133.9	-3.0574	2.9739	0.0524

Centre de Calcul Scientifique en région Centre


Clément Crémoux
email: c.cremoux@adveez.com

Reliability of active RFID tag immersed in water for anti

Keywords:

In this paper, a method is developed to the reliability of active RFID tags immersed in water. This situation exists in particular when infants kidnaping tag take bath. This study discusses electromagnetic radiation of an emitting tag water volume and then proposes a statistical method to predict the radio coverage of the immersed tag. Finally it demonstrates that the RFID technology is barely affected by this issue even though other RFID technologies do not support

INTRODUCTION

To prevent abduction for infant, on RFID is a growing application. To control where and when infant is in the maternity are attached to infant's ankle and located by stations installed in the building. technologies or strategies have been and locate the RFID tags. The first strategy consists in measuring the RSSI of the ID signal regularly (e.g. every quarter of minute infant's positions by RSSI triangulation [ strategy uses the Time Of Flight (TOF) triangulation usually based on UWB signal [2]. A third strategy is called polling LF-RF: the tag sends radiofrequency (RF) link only if it is a station placed closed to an access door which covers a limited range by a low frequency (LF) signal Technologies with passive tags and security gate are affected by weaker ranges than technologies active tag.

One of the limits of technologies links is due to the immersion of the tag in example during a bath. The degradation of the RF link can result in false triggering of alarms which discredits the system. The proposed study discusses about impact of immersion of an active RFID tags based on a RF-LF link. A statistical method predict the radio coverage of RFID detection control if the deployment of stations is efficient.

PRESENTATION OF THE ANTI-RAPT

Functional principle of the system

Provide abduction is a hands free access control application. Fig. 1 describes the operation of the system. Stations which are placed near access secured areas emit permanently a low frequency (LF) signal at 125 KHz. The strength of the emitted signal __________________________________________________________________________ Labrousse, Toulouse, France, CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse, France 31400

Reliability of active RFID tag immersed in water for anti applications

Alexandre Boyer The proposed study discusses about the impact of immersion of an active RFID tags based on thod is proposed to of RFID detection and deployment of stations is efficient.

RAPT SYSTEM

Functional principle of the system hands free access control . Fig. 1 describes the operation of the near access doors of emit permanently a low frequency (LF) strength of the emitted signal is adjustable to cover a detection to 3 meters. When a tag enters in a detection area, it sends its ID to the station with band. In our application, 434 and 908 MHz ISM frequencies are used. The access is blocked and alarms are sent if unauthorized persons arrive area. The low power consumption of tag of the detection area and the weak influence of the site configuration on the installation parameters main advantages of this technology.

Compared to usual hands free access contro applications, some changes are provided:

-Tags are attached with they have to be carried all time.

-A permanent RF link between the station is reused. It sends warning if the wristbands is wrenched or if the battery level tag location. Further, RF polling is sent every fifteen minutes to ensure that the tag is working fine.

Limitations

It is unlikely that a bath is situated in detection area of a station. Therefore reliability of out in this study. The immersion of the tag may be harmful for the RF range. Limitations of can be caused by desensitization the throughput is about 20kbits/s, dela rejected. Desensitization of links parameters:

-The radiated power emitted by tags parameter is affected by tag immersion -The sensitivity of station noise, which is supposed constant antennas are also supposed quasi omnidirectional

Reliability of active RFID tag immersed in water for anti-rapt

Karim Ben Dhia 1,4 is adjustable to cover a detection area ranging from 1 When a tag enters in a detection area, it with a RF signal in the UHF . In our application, 434 and 908 MHz ISM . The access is blocked and alarms if unauthorized persons arrive in a forbidden ow power consumption of tag, the exactness the weak influence of the site configuration on the installation parameters are the main advantages of this technology.

Compared to usual hands free access control applications, some changes are provided:

Tags are attached with wristband because they have to be carried all time.

between the tag and the t sends warning if the wristbands is if the battery level is low, whatever the RF polling is sent every fifteen tag is working effectively nctional principle of the anti-rapt system a bath is situated in detection area reliability of LF signal is ruled The immersion of the tag may be Limitations of the RF range can be caused by desensitization and delay spread. As the throughput is about 20kbits/s, delay spread is of links depends on two main adiated power emitted by tags. This parameter is affected by tag immersion.

sensitivity of station and the ambient supposed constant and receiving quasi omnidirectional.

EFFECT OF THE IMMERSION ON AN

The bathtub is considered as a parallelepiped filled with water. The volume of water is equivalent to a dielectric resonator. The tag antenna can excite some of the cavity modes of the dielectric resonator. The resonant frequencies of cavity modes equation ( 1 The volume of water is a dielectric parallelepiped, with a relative permittivity and an electrical conductivity equal to 72 and 650µS/m respectively. It is representative of bath water at 37°C. modeled by an equivalent magnetic dipole with a magnetic moment equal to 7.54e-6 A.m empty volume enwraps this source. The whole is placed within the water volume.

To illustrate the influence of immersion on radiated field, 2D far field patterns various situations are compared in Fig. 2. curve is the reference case without water other cases are the immersion; only tag inside the volume is changed.

Fig. 2. Farfield on level water at 908

Depending on the frequency, the dimensions of the bathtub, the orientation and the position of the tag within the bathtub and the electrical characteristic of the water, the variation of the far-f given direction can exceed 20dB, without relation between the far-field emission and these parameters. This complex behavior confirmed by radiated measurements of an immersed tag in anechoic chamber.

AN RFID TAG

The bathtub is considered as a parallelepiped filled

The volume of water is equivalent to a dielectric resonator. The tag antenna can excite some cavity modes of the dielectric resonator. The resonant frequencies of cavity modes are given by (1) mode number and (a,b,c) the

The relative electric 72 and depends on The high electric permittivity value induces high reflection amount at dielectric boundaries. It enhances the dielectric resonator behavior of the volume of water. Although radiation diagram of the tag is similar to the radiation pattern of a magnetic dipole, the bathtub resonances give rise to a complex field distribution

. Electromagnetic ] are performed in order to

The volume of water is a dielectric parallelepiped, with a relative permittivity and an electrical conductivity equal to 72 and 650µS/m respectively. It is representative of bath water at 37°C. Tag antenna is nt magnetic dipole with a 6 A.m 2 . A cylindrical empty volume enwraps this source. The whole is influence of immersion on the patterns at 908 MHz in in Fig. 2. The first without water and the two ; only the position of the 908 MHz frequency, the dimensions of the bathtub, the orientation and the position of the tag within the bathtub and the electrical characteristic of field radiation in a , without any obvious field emission and these behavior has been confirmed by radiated measurements of an immersed

The bathtubs can have different size and baby's tag can be placed everywhere inside this volume. Since these parameters are unknown influence the radiated emission approach is the best solution to stud tag's immersion.

STOCHASTIC DISTRIBUTION OF FIELD RADIATED BY AN IMMERSED

Methodology of probabilistic approach

The field radiated by the immersed tag cannot be precisely determined since the dimensions of the bathtub, the position and the orientation of the tag remain unknown. However, if these parameters are considered as random variables, the statistical distribution of the field radiated established. This distribution will be integrated to the indoor radio propagation model presented in part V order to predict the effect of immersion on the tag radio coverage. A comparison with of a non-immersed tag will allow evaluate the degradation due to the tag immersion.

The method consists to simulate a large sample with different values of unknown electromagnetic model presented in the previous part. Results are the radiated electric field in all directions from the tag emission frequency (434 and 3645 simulations are performed in about 4 days on a machine equipped with an Intel® Pentium® Dual CPU 1,80GHz processor. Post applied to extract the statistical distribution of the electric field radiated by the immersed tag

Cumulative distribution of the radiated emission of tag

In order to highlight the significant impact of the immersion on the RF link strength fluctuations, the distributions of electric field radiated by a tag immersed in water or placed in free space compared. The statistical distribution of the electric field in this last case can be computed analytically. The electric field produced by a magnetic dipole in far field at a distance R is given by (2) [

sin ⁄ ! " #
With E0 which depends on the excitation and the size of the magnetic dipole the normalized amplitude of the electric field. not depend on the distance R but is influenced by the dipole orientation which is unkn the distribution of the electric field radiated by the magnetic dipole in all directions following cumulative distribution function (CDF

$ %" # & '' ) " # * 1 , 1 , -
The simulation results of immersion for various water volume dimensions, tag position and orientation give the distributions of the The bathtubs can have different size and baby's tag can be placed everywhere inside this volume. Since these parameters are unknown, and as they can the radiated emission strongly, a stochastic best solution to study the impact of the

ON OF THE ELECTRIC AN IMMERSED TAG

Methodology of probabilistic approach

The field radiated by the immersed tag cannot be precisely determined since the dimensions of the nd the orientation of the tag remain unknown. However, if these parameters are considered as random variables, the statistical distribution of the field radiated by the tag can be This distribution will be integrated to the ation model presented in part V in the effect of immersion on the tag comparison with the radio coverage will allow evaluate the degradation due to the tag immersion.

to simulate a large sample unknown parameters from the presented in the previous part. radiated electric field in all directions e tag emission frequency (434 and 908 MHz).

ons are performed in about 4 days on a machine equipped with an Intel® Pentium® Dual Post-processing is then applied to extract the statistical distribution of the electric field radiated by the immersed tag

Cumulative distribution of the radiated

In order to highlight the significant impact of the immersion on the RF link strength fluctuations, the distributions of electric field radiated by a tag immersed in water or placed in free space are compared. The statistical distribution of the electric field in this last case can be computed analytically.

lectric field produced by a magnetic dipole in far is given by (2) [6].

" # ⁄ ! (2)
depends on the amplitude of the magnetic dipole. Efar is litude of the electric field. It does but is influenced by the dipole orientation which is unknown. Theoretically, the distribution of the electric field radiated by the in all directions is given by the following cumulative distribution function (CDF). An offset 7 is added to obtain model. The parameters of the distribution are given in Table I at 434 and908 MHz. Fig. 3. Comparison between CDF of magnetic dipole in free space and immersed in water At 434 MHz, the amplitude of the electric field radiated by the magnetic dipole is larger than the field produced by the same dipole in free s This result shows that the RF coverage will better in these conditions when the However that is not the case at 908 MHz

PROPAGATION MODEL

Radio planning solver

The solver is based on the Multi Frequency Domain ParFlow (MR-FDPF) method [7 The starting point of this algorithm is a approach proposed by Chopard and advantages of the method are its speed accuracy. Attenuation, reflection and diffraction phenomena are modeled. Radio coverage at 434 MHz over a 1500 square meter building floor obtained in only 3 seconds on a machin with an Intel® Core™ I3 clocked at As the coverage is computed only in the study will be limited to a building models are easy to construct: the positions of the walls and their reflectance and absorption wall are requested in input. The i solver are the building map, the emitter positions and the wall proprieties. These parameters from the electric permittivity and loss tangent regular building materials. The solver computes the radiated power density and the propagation loss Lp within the structure. at 434 and 908 MHz, as shown in can be approximated by a Rayleigh ] to model field reverberation chambers. The solver is based on the Multi-Resolution FDPF) method [7]. ng point of this algorithm is a ParFlow and al. [8]. The are its speed and its Attenuation, reflection and diffraction

. Radio coverage at 434 MHz building floor can be on a machine equipped clocked at 2,53GHz. only in a 2D plane, the study will be limited to a building floor. Building the positions of the walls absorption coefficients of its

The input data of the emitter positions and These parameters are determined and loss tangent of The solver computes the radiated power density and the propagation loss Lp

Validity of model

To verify the exactness of the results given by th radio propagation solver, on-site measurements been completed and compared The measurements of the received electric field are performed with a calibrated handheld RF field strength meter in the different rooms of the installation site.

Simulated results are interpolated to grid of measurement points. The d the measured and simulated evaluated. Fig. 4 presents a comparison measured and the simulated propagation losses corridor. This model of loss propagation in indoor environment is validated with 0,947.

STATISTICAL DISTRIBUTION

COVERAGE

The statistical distributions of the tag far field radiation determined in part IV radio planning solver to determine distribution of the radio coverage building. A comparison between immerged and a non-immerged order to evaluate the impact of risk of RF link loss.

Integration of the CDF of the radiated emission of the tag

The aim of the radio planning solver is to evaluate the probability that a RF link can be established between a tag and a receiving signal to noise ratio of the received signal is beyond a given threshold. In other word, a RF link is lost if the propagation loss exceeds a threshold noted L In part IV, the CDF of the electric field radiated by a tag have been established. However, they cannot be exploited from the radio planning solver results directly since they are given in term of loss. It is necessary to convert the CDF of the radiated electric field into the CDF of propagation loss. propagation loss can be computed according to equation ( 5): the exactness of the results given by the site measurements have with simulation results. easurements of the received electric field are performed with a calibrated handheld RF field strength meter in the different rooms of the installation site. The aim of the radio planning solver is to evaluate link can be established nd a receiving station, i.e. that the signal to noise ratio of the received signal is beyond a

In other word, a RF link is lost if the threshold noted Lp max .

CDF of the electric field radiated by been established. However, they cannot be exploited from the radio planning solver results since they are given in term of propagation It is necessary to convert the CDF of the radiated electric field into the CDF of propagation loss. The propagation loss can be computed according to ;< = = > ? > ?# With: -Lp, loss propagation determined by MR -P r is the sensitivity threshold the gain of receiving antenna -P t and G t are the transmitted power and the antenna gain of the tag. In far field, they are the electric field by ( 6) " # @A B C D When, K is equal to 60 for the The probability of RF link loss is given by the complementary CDF of the propagation loss F

$ E ;< 1 , $ % &. F
Other statistical parameters can be used (body influence, attenuation of crowd, different environment noise, jamming …) but they are disregarded

Results of statistical coverage

The radio planning solver provides predictions of the average propagation loss within the building. The CDF of the immersed tag and the tag in free space are used to compute the probability that the propagation loss exceeds the maximum allowed propagation loss on each point of the building. In accordance with previous results, is really better when the tag is immersed as shown in Fig. 6. In this case, the probability that the RF link establishes correctly exceeds 9 of the building surface. When the tag is not immersed, this area falls down to 48 %. At 908 MHz differences are less significant, in accordance with Fig. 6. Percentage of building correctly

Fig. 1 .

 1 Fig. 1. Functional principle of the anti

  ) [3]: √ Where (m,n,p) are the mode number dimensions of the bathtub. The r permittivity ε r of water is about 72 temperature and frequency [4]. The high electric permittivity value induces high reflection amount at dielectric boundaries. It enhances the dielectric resonator behavior of the volume of water. Although the radiation diagram of the tag is similar to the radiation pattern of a magnetic dipole, the bathtub resonances give rise to a complex around the volume of water. Electromagnetic simulation with FEKO© [5] are performed in order to study this field distribution.

  results of the magnetic dipole for various water volume dimensions, tag the distributions of the electric field F Efar at 434 and 908 MHz, as shown in Fig 3. Their CDF can be approximated by a distribution, as suggested in [

  Comparison between CDF of magnetic dipole , the amplitude of the electric field E far larger than the field produced by the same dipole in free space conditions. coverage

Fig. 4 .

 4 Fig. 4. Comparison between the simulation of the propagation loss 908 MHz

  Simulated results are interpolated to fit with the The differences between propagation loss are esents a comparison between the propagation losses along a his model of loss propagation in indoor nment is validated with a correlation factor of Comparison between the measurement and propagation loss along a corridor atDISTRIBUTION OF THE RADIOtatistical distributions of the tag far field IV are introduced in the to determine a statistical coverage of a tag within a between the coverage of an immerged tag is provided in of the immersion on the of the CDF of the radiated

  Fig. coverage map for both tag configurations. station is placed, at the green point.

Fig. 5 .

 5 Fig. 5. Comparison of the RF coverage between a non-immersed and an immersed tag at

In this paper, a method is developed to evaluate immersed in water. This n exists in particular when infants wearing anti- This study discusses the effects on the ing tag plunged into a hen proposes a statistical method to predict . Finally it demonstrates affected by this issue even t support immersion.

		2,3	Karim Ben Dhia
	localization;	radio
	, technology based
	is a growing application. To control where
	maternity ward, RFID tags
	and are identified and
	stations installed in the building. Several
	technologies or strategies have been proposed to detect
	strategy consists in
	signal sent by the tag
	minute) and deducing
	by RSSI triangulation [1]. Another
	(TOF) triangulation
	A third strategy is
	its ID through a
	is in the vicinity of
	placed closed to an access door which covers
	range by a low frequency (LF) signal.
	and security gate are
	technologies based on
	One of the limits of technologies based only on RF
	to the immersion of the tag in water, for
	example during a bath. The degradation of the RF link
	can result in false triggering of alarms which discredits

TABLE I .

 I

	PARAMETERS OF CUMULATIVE DISTRIBUTION FUNCTIO
	Frequencies	E0 (V)	σ (V)
	434 MHz	0,0187	0,0250
	908 MHz	0,0819	0,0522


A Boyer

B Vrignon

J Shepherd

M Cavarroc

Evaluation of the Near-Field Injection Method at Integrated Circuit Level

Keywords: near-field scan injection, electromagnetic attack, susceptibility I

Near-field injection is a promising method in order to induce local faults in integrated circuits. This method can be used for various applications such as electromagnetic attacks on secured circuits or susceptibility investigations. This paper aims at evaluating the ability of near-field scan injection to induce local disturbances in integrated circuits. The study relies on measurements performed by on-chip voltage sensors, which provide an accurate method to characterize the induced voltage fluctuations.

INTRODUCTION

Near-field injection constitutes a promising method for various applications such as electromagnetic attacks on secured circuits [START_REF] Schmidt | Optical and EM Fault -Attacks on CRTbased RSA: Concrete Results[END_REF] [START_REF] Poucheret | Local and Direct EM Injection of Power into CMOS Integrated Circuits[END_REF] or investigations of integrated circuits (ICs) susceptibility to electromagnetic disturbances [START_REF] Boyer | Characterization of the Electromagnetic Susceptibility of Integrated Circuits using a Near Field Scan[END_REF] [START_REF] Dubois | Near-Field Electromagnetic Characterization and Perturbation of Logic Circuits[END_REF]. The method consists in placing a miniature near-field probe in the vicinity of an electronic device and exciting the probe with a disturbance signal in order to induce a local voltage fluctuation and trigger a failure on the device under test. This method has been mainly used for analysis of IC susceptibility with injection coupled at package level. However, the injection at die level offers the advantage to locate precisely the origin of the induced faults. It can improve not only the efficiency of intentional electromagnetic injection on secured ICs such as differential fault injection (DFA) but also the diagnosis of susceptibility of an IC block.

In this paper, the ability of the near-field injection method to induce local disturbances at IC level is investigated. The proposed study aims at evaluating the possibility of the method in terms of resolution and fault generation with basic handmade magnetic and electric field probes. Moreover, it attempts to clarify the nature of the coupling between the injection probe and the circuit under test. The understanding and modeling of the coupling is essential not only to optimize the injection method, but also to understand the origins of triggered failures and the development of efficient countermeasures against intentional electromagnetic attacks.

This study is based on near-field injection measurements performed on a test chip which integrates several on-chip voltage sensors. They ensure an accurate and high frequency characterization of the voltage induced along circuit interconnects. The first part of the paper describes the nearfield injection bench and the near-field probes used for this study. The second part introduces the on-chip sensor and describes the test chip. The third part presents the near-field injection experimental results. Some conclusions about the injection method efficiency and the nature of the coupling between the probe and the circuit under test are provided.

II. NEAR-FIELD INJECTION TEST BENCH

A. Presentation of the test bench

The near-field injection consists in placing a miniature near-field probe in the vicinity of a device and injecting disturbance signals in the probe in order to induce a local voltage fluctuation and trigger a failure. It relies on a near-field scan table in order to control precisely the position of the nearfield probe compared to the device under test. The typical test bench is described in Fig. 1.

RF signal control

Positioning control

Near field probe

IC

IC monitoring

Scan table

RF signal generator

Directional coupler

Power amplifier

IC status

Figure 1. Near-field injection test bench

The near-field probe can be excited by various types of signal (harmonic or pulse). In this study, harmonic signals are considered for frequency analysis purpose. They are produced by a RF synthesizer followed by a 50 W power amplifier. By experience, this range of power is required to induce failures in the IC. However, one drawback of the use of harmonic signals is the major power dissipation that the probe and cables have to withstand. This constraint can be reduced with short pulse or on-off keying modulated sinus waveforms. With harmonic signals, a directional coupler is advised to measure the amplitude of the forward power delivered by the power amplifier. As near-field probe are totally mismatched, the excitation current or voltage is not easily predictable.

The injection at IC level requires a precise control of the position of the near-field probe compared to the IC under test. For this purpose, our test bench is equipped with a laser telemeter, which has a theoretical accuracy of 12.5 µm. First, a preliminary calibration is necessary to measure the coordinates of the near-field probe in the telemeter reference frame. Then, the laser telemeter scans the IC under test and provides a 3D topology of its surface. Finally, the scan table control software is able to locate the position of the near-field probe compared to the IC under test. The scan can be done at constant height above the IC under test. Because of the noise introduced by stepper motors of the scan table, the accuracy of the positioning system is about 50 µm.

B. Near-field probes

Near-field probes are miniature and wideband antennas usually dedicated to the sense of either electric or magnetic fields in the vicinity of electronic components, circuit boards or antennas. They are usually based on small loops and opened tips. If they are excited by an electrical signal, they produce a strong electric or magnetic fields in their near-field region, which can couple on a device under test (DUT) to produce important voltage fluctuations locally. Near-field injection relies on probe with a good resolution and a sufficient efficiency to induce a perceptible effect on the DUT.

The development of optimized near-field probes is outside the scope of this work. In this study, two usual handmade probes made from RG405 semi-rigid coaxial cables are evaluated. Fig. 2 gives the dimensions of the probes used in this study: a tangential magnetic H and a normal electric field E probe with similar dimensions.

C. Near-field probe characterization

Near-field probe for injection are characterized by two main figures of merit: the resolution and the injection efficiency. The resolution quantifies the ability to produce an important field locally. In practice, it means that the coupling between a high resolution probe and two identical neighbor lines separated by a short distance can be distinguished. The resolution is defined as the distance that the injection probe has to travel between a first position where the coupling with the DUT is optimal and a second position where the coupling has been reduced significantly (e.g. by 3 dB for 3 dB resolution). The efficiency is the coupled voltage induced on a given DUT for a given probe excitation. It aims at comparing the effect of different injection probe.

In order to measure the resolution and the efficiency of our injection probes, they are placed at 0.5 mm above a narrow microstrip line with a width of 0.15 mm. The probes are excited by a RF signal synthesizer which delivers 13 dBm forward power. The resolutions of both probes are similar, as shown in Table I. Fig. 3 presents the measured evolution of the voltage coupled on the microstrip line during E and H field injection. The H field probe is 15 dB more efficient than the E field probe up to 1 GHz. Basic models of the near-field radiation of both probes have been constructed using elementary electric and magnetic dipole models. The evolution of the electric and magnetic fields in the vicinity of the probe can be computed. The voltages V E or V H induced on the terminal of a line (e.g. a microstrip line) by an electric or magnetic field probe respectively can be deduced by using the field to transmission line equations [START_REF] Taylor | The response of terminated two-wire transmission line excited by a non uniform electromagnetic field[END_REF].

∫ + = line L O line Z line L L L L E dl h E C j Z Z Z Z V ω 2 1 2 1 (1) ∫ + = line L O line L L L H dl h H j Z Z Z V tan 0 2 1 1 ωµ (2)
From the dimensions of the tested injection probe and microstrip line, the voltage induced by the injection on the microstrip line terminal is simulated. Simulation results are superimposed on measurement curves, as shown in Fig. 3. This basic model is sufficient to obtain a good estimation of probe couplings for relatively short separation distances. The coupling increases with frequency by 20 dB/dec up to 1 or 2 GHz.

III. PRESENTATION OF THE TEST CHIP

Measuring the local voltage fluctuations in an IC induced by a near-field probe is not obvious. First, the measurement of high frequency voltage fluctuations through I/O pads and IC package introduces a non negligible attenuation and distortion. High frequency probes placed directly on IC pads disturbs the measurements and are partially coupled with the injection probe. A convenient method to measure on-chip rapid voltage fluctuation accurately, with a minimum of distortions and coupling with the injection probe relies on the use of on-chip voltage sensors (OCS) [START_REF] Ben Dhia | A new method for measuring signal integrity in CMOS ICs[END_REF] [START_REF] Ben Dhia | On-Chip Noise Sensor for Integrated Circuit Susceptibility Investigations[END_REF]. A dedicated test chip with canonical structures and OCS has been developed for the characterization of near-field injection at IC level.

A. Presentation of the test chip

A test chip has been designed with Freescale® in 0.25 µm SMARTMOS 8 technology with 4 metal layers in order to study the near-field injection on basic interconnect and bus structures. Twenty four OCS are disseminated within the test chip to monitor local voltage fluctuations induced by the nearfield injection.

The presented study is focused on one structure called Struct1, described in Fig. The test chip is mounted in CQFP64 package with a removable metallic lid in order to place the near-field probe as close as possible to the die surface.

B. Presentation of the on-chip sensor

In order to monitor voltage fluctuations induced by an injection near-field probe along IC interconnects, on-chip voltage sensors are disseminated within the test chip. This sensor is able to measure the waveform of voltage bounce across non accessible nodes with a precise time resolution (up to 100 ps for a development with a CMOS 0.25 µm technology). Its analog bandwidth is equal to 2.7 GHz. The acquisition principle is based on a sequential equivalent-time sampling (Fig. 5). Its principle and its implementation in the test chip are explained in [START_REF] Ben Dhia | A new method for measuring signal integrity in CMOS ICs[END_REF]. In order to prevent noise coupling to sensors, they are supplied by an internal voltage regulator connected to a dedicated power supply and are isolated from the bulk substrate by a deep Nwell. Moreover the OCS is routed with only three metal layers. A complete shielding of the OCS is done with the top metal layer, avoiding the OCS to be disturbed by the radiated fields. The sensor is controlled by an external acquisition card which sends sampling clock, delay control command and converts the analog signals delivered by the sensor. The waveform of a signal can be reconstructed only if the signal is reproducible and synchronized with the sampling command of the sensor. During the near-field injection, the RF signal generator is synchronized with the OCS sampling command. Fig. 6 presents an example of reconstructed harmonic signal induced within the test chip during a near-field injection. The frequency of the probe excitation signal is 900 MHz. The step time is set to 150 ps.

B. Injection with the magnetic field probe

The tangential magnetic field probe is oriented along the Y axis defined in Fig. 7 (orientation 0°) in order to optimize the magnetic coupling with Struct1 lines.

1) Evaluation of the localization/resolution

A preliminary experience consists in displacing the H field probe above the bonding wire of V ref (Scan1). The excitation frequency is set to 1.4 GHz. The power amplifier delivers a 45 dBm forward power. Fig. 8 compares the voltage coupled only on four lines of the Struct1 for clarity purpose. The maximum induced voltage reaches up to 220 mV when the probe is placed just above the bonding wire. The same voltage is induced on the eight lines because the coupling is localized on the bonding wire V ref .

The same experiment is repeated but the probe moves along Scan2 area. Fig. 9 presents the evolution of the voltage coupled on lines 0, 0.455, 120 and 320 µm according to the magnetic field probe position. The sensors do not measure exactly the same voltage for a given probe position. The positions of maximum coupling on these four lines are indicated by arrows on Fig. 9. They are shifted along the probe position. The maximum coupled voltage is similar for all the lines except line 0 µm. As line 0 µm is at the border of Struct1, the presence of tiling could modify its impedance and thus the coupling of the magnetic field produced by the injection probe. This result shows that the injection can be distinguished between two lines separated by more than 100 µm. At 300 µm above the IC line, the 1 dB resolution of the magnetic field probe is about 250 µm. In order to verify that the coupling is mainly magnetic, the near-field probe is rotated by 90° and the previous experiment is repeated. Fig. 10 compares the evolution of the voltage coupled on line 120 µm for both probe orientations (0° and 90°). The coupled voltage is divided by two for orientation 90° because the magnetic field direction is not perpendicular to the lines of Struct1. However, as a non negligible voltage is still measured by the sensor, the magnetic field is certainly coupled on the V ref plane placed under the lines.

V ref

2) Effect of the near-field scan altitude

The effect of the near-field probe height above the die surface on the voltage coupled on the lines of Struct1is also measured. Fig. 11 presents the evolution of the voltage measured by the sensor placed on line 120 µm according to the probe height which ranges from 200 to 1000 µm. The probe is placed just above line 120 µm. Obviously, when the scan altitude is reduced, the coupled voltage increases rapidly, especially for the altitude smaller than 500 µm. At such a short distance, the coupling with the line under test is mainly due to the nearest part of the injection probe.

3) Frequency characterization of the magnetic field coupling

The OCS placed on line 120 µm is used to measure the evolution of the coupled voltage on the line depending on the frequency of the excitation signal of the probe. The measurements are done with a constant forward power amplitude (43 dBm measured on a 50 Ω load). Fig. 12 presents the evolution of the voltage amplitude coupled on line 120 µm depending on the frequency. This result shows that the coupling is globally improved at high frequency as expected. Voltage fluctuations of several hundreds of millivolts can be generated on a 1 mm long IC interconnect, which can be just enough to start triggering failures in digital or analog ICs. However, the coupling with the IC interconnect does not increase as continuously and rapidly as the coupling with a PCB microstrip line (Fig. 3). Moreover, the coupling is surprisingly high if the geometry of line 120 µm and the basic coupling model used in part II.C are considered. The structure of the line is certainly more complicated than a basic microstrip line. Moreover, the parasitic coupling between the line and the circuit substrate should be taken into account to model the evolution of the coupled voltage. Further measurements and simulations have to be carried at clarify the physical interaction between the injection probe and a circuit interconnect.

C. Injection with the electric field probe

The normal electric field probe is placed at 300 µm above Struct1 along the Scan2 area. The excitation frequency is set to 1.4 GHz. The power amplifier delivers a 45 dBm forward power. Fig. 13 presents the evolution of the voltage coupled on lines 0, 0.455, 120 and 320 µm according to the electric field probe position. The coupling between the IC lines and the electric field probe is nearly 7 times weaker than the coupling with the magnetic field probe. The measurement of the voltage fluctuation amplitude is partially disturbed by a residual noise of 5 mV. Due to the weak coupling, the four measurement curves are nearly superimposed and no conclusions about the resolutions of electric field probes can be deduced. Fig. 14 presents the evolution of the voltage amplitude coupled on line 120 µm according to the frequency. The measurements are done with a constant forward power amplitude (43 dBm measured on a 50 Ω load). At a given probe height, the coupling between the IC line and the electric field probe is about 7 times weaker than the coupling with the magnetic field probe. It confirms the efficiency difference between both probes measured in part II.C. The coupling increases slightly with frequency, but not as continuously and rapidly as the coupling with a PCB microstrip line (Fig. 3). Once again, this evolution can be explained by a more complicated line structure than a basic microstrip line and the complex parasitic coupling between the line and the circuit substrate. The paper has presented a original experimental work dedicated to the evaluation of the near-field injection to induce local disturbances at IC level. This method is very promising for various applications such as electromagnetic attacks on secured ICs or investigations of susceptibility to electromagnetic disturbances. The measurements rely on onchip voltage sensor to overcome the parasitic effect of IC pads and package interconnects on high frequency signal. Two handmade and basic electric and magnetic field probes have been tested. Their dimensions are in the millimeter range.

Magnetic field rpobes provide a resolution sufficient to distinguish the coupling between two lines separated by 100 µm, for a scan altitude of 300 µm. The improvement of the resolution is necessary. It can be provided by a miniaturization of the probes and a reduction of the scan altitude. It requires a test bench with a precise control of the altitude probe to prevent contact with the die surface. The measurement results have shown that even with similar dimensions the magnetic field probe produces voltage fluctuations about 7 times larger along an IC interconnect than the electric field probe. If the magnetic field probe is connected to a RF power amplifier which can deliver several tens of watts, the amplitude of the induced voltage fluctuation can reach several hundred of millivolts. It can be enough to cause triggering failures in digital or analog circuits. The use of a harmonic RF power source has a major drawback: the heating of the injection probe. In order to increase the amplitude of IC voltage fluctuations without excessive heating of the injection probe, it is recommended to use pulsed signals.

However, the results have shown that the physical coupling between the injection probe and an IC interconnect remains complex. It has to be clarified to propose a valid model. A simple coupling model of the electric or magnetic field produced by the injection probe on a canonical line structure such as a microstrip line underestimates the actual voltage fluctuations measured on IC interconnects. Further works have to be done to improve of field coupling on IC interconnects.

Figure 2 .

 2 Figure 2. Dimensions of the near-field probes: tangential magnetic field (on the left) and normal vertical electric field probes (on the right)

Figure 3 .

 3 Figure 3. Comparison between measurement and simulation of the coupling of the near-field E and H probes on a microstrip line

Figure 4 .

 4 Figure 4. Description of Struct1

Figure 5 .

 5 Figure5. Principle of signal reconstruction with the on-chip sensor[START_REF] Ben Dhia | A new method for measuring signal integrity in CMOS ICs[END_REF]

Figure 6 .

 6 Figure 6. Example of reconstructed harmonic signal by the OCS

Figure 7 .

 7 Figure 7. Scanned area above the test chip

Figure 8 .

 8 Figure 8. Evolution of the voltage coupled on four lines of Struct1 when the magnetic field probe moves along Scan1

Figure 10 .

 10 Figure 10. Evolution of the voltage coupled on line 120 µm of Struct1 for two perpendicular orientations of the magnetic field probe

Figure 11 .

 11 Figure 11. Evolution of the voltage coupled on line 120 µm of Struct1 according to the magnetic field probe altitude

Figure 12 .

 12 Figure 12. Evolution of the voltage coupled on line 120 µm of Struct1 according to the magnetic field probe excitation frequency

Figure 13 .

 13 Figure 13. Evolution of the voltage coupled on the four lines of Struct1 when the electric field probe is positionned at 300 µm above Struct1

Figure 14 .

 14 Figure 14. Evolution of the voltage coupled on line 120 µm of Struct1 according to the electric field probe excitation frequency

TABLE I .

 I

	1 dB	0.3 mm	0.4 mm
	3 dB	0.65 mm	0.8 mm
	6 dB	0.95 mm	0.95 mm

RESOLUTION OF THE NEAR-FIELD PROBES

Resolution H field probe E field probe

This work is supported by the French Agence Nationale de la Recherche under contract ANR-12-INSE-005 E-Mata-Hari.

ACKNOWLEDGMENT

The authors would like to thank Nicolas Mauran from LAAS-CNRS for his support on the near-field scan test bench.