Philippe Chartier
email: philippe.chartier@inria.fr

Nicolas Crouseilles
email: nicolas.crouseilles@inria.fr

Mohammed Lemou
email: mohammed.lemou@univ-rennes1.fr

Florian Méhats
email: florian.mehats@univ-rennes1.fr

Florian M Éhats

Xiaofei Zhao

UNIFORMLY ACCURATE METHODS FOR THREE DIMENSIONAL VLASOV EQUATIONS UNDER STRONG MAGNETIC FIELD WITH VARYING DIRECTION

Keywords: Vlasov-Poisson equation, Three dimensions, Strong magnetic field, Varying direction, Uniformly accurate method, Particle-In-Cell AMS Subject Classification: 65L05, 65L20, 65L70. [35] E. Sonnendrücker, Numerical Methods for Vlasov Equations, Lecture notes, 2016

come     Uniformly accurate methods for three dimensional Vlasov equations under strong magnetic field with varying direction

Introduction

Vlasov models have been widely considered for modelling the dynamics of plasmas as encountered in magnetic fusion devices known as a tokamaks, where a strong external magnetic field is applied so as to confine the charged particles. In this paper, we consider the three dimensional Vlasov-Poisson equation with a strong non-homogeneous magnetic field whose direction may vary [START_REF] Degond | On the asymptotic limit of the three dimensional Vlasov-Poisson system for large magnetic field: formal derivation[END_REF][START_REF] Golse | The Vlasov-Poisson system with strong magnetic field[END_REF]35]

∂ t f ε (t, x, v) + v • ∇ x f ε (t, x, v) + E(t, x) + 1 ε v × B(x) • ∇ v f ε (t, x, v) = 0, (1.1a)
∇ x • E(t, x) = R 3 f ε (t, x, v)dv -n i , (1.1b)
f ε (0, x, v) = f 0 (x, v), (1.1c)
where, for a given T > 0,

f ε : (t, x, v) ∈ [0, T ] × R 3 × R 3 → f ε (t, x, v) ∈ R is the unknown, f 0 : (x, v) ∈ R 3 × R 3 → f 0 (x, v) ∈ R a given initial distribution, where B : x ∈ R 3 → B(x) ∈ R 3
denotes the external magnetic field,

E : (t, x) ∈ R + × R 3 → E(t, x) ∈ R 3
the self-consistent electric-field function, 0 < ε ≤ 1 a dimensionless parameter inversely proportional to the strength of the magnetic field and n i ≥ 0 the ion density of the background. The system (1.1) has a lot of invariants and we will be interested in particular in the Hamiltonian defined by

H(t) := R 3 R 3 1 2 |v| 2 f ε (t, x, v)dxdv + 1 2 R 3 |E(t, x)| 2 dx.
(1.2)

1
The above Vlasov-Poisson model (1.1) is derived from the three dimensional Vlasov-Maxwell equations by considering the electrostatic approximation. Unlike some asymptotically reduced models such as the gyrokinetic equations [START_REF] Latu | Field-Aligned Interpolation for Semi-Lagrangian Gyrokinetic Simulations[END_REF][START_REF] Saint-Raymond | The gyro-kinetic approximation for the Vlasov-Poisson system[END_REF] or the drift-kinetic limit equations [START_REF] Bostan | The Vlasov-Maxwell system with strong initial magnetic field. Guiding-center approximation[END_REF][START_REF] Frénod | Long time behavior of the two-dimensional Vlasov equation with a strong external magnetic field[END_REF][START_REF] Filbet | On the Vlasov-Maxwell system with a strong magnetic field[END_REF], model (1.1) is set at the kinetic scale and is of paramount importance for studying the plasma dynamics in the tokamak device.

In the strong magnetic field limit regime, the charged particles exhibit very fast rotations with cyclotron period proportional to ε, while remaining confined along the magnetic line. In such a case, the small parameter 0 < ε 1 renders the solution f ε (t, x, v) of (1.1) highly-oscillatory in time. Classical numerical integrators such as splitting or finite-difference schemes thus require time steps smaller than the cyclotron period in order to accurately capture the dynamics, thus implying severe computational burden. Recent efforts have aimed at designing numerical schemes which allow step-sizes much larger than the cyclotron period. Upon assuming that the magnetic field has a fixed direction in space, i.e. that B(x) = (0, 0, b(x)) T , b(x) > 0 (a popular choice both for formal and rigorous analyses [START_REF] Bostan | The Vlasov-Maxwell system with strong initial magnetic field. Guiding-center approximation[END_REF][START_REF] Degond | On the asymptotic limit of the three dimensional Vlasov-Poisson system for large magnetic field: formal derivation[END_REF][START_REF] Golse | The Vlasov-Poisson system with strong magnetic field[END_REF]), several multiscale numerical methods have been proposed [START_REF] Ph | Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time two-dimensional Vlasov-Poisson equation with uniform strong magnetic field[END_REF][START_REF] Filbet | Numerical simulations to the Vlasov-Poisson system with a strong magnetic field[END_REF][START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF][START_REF] Filbet | Asymptotically preserving particle-in-cell methods for inhomogeneous strongly magnetized plasmas[END_REF][START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF]. Among them, Filbet et al. constructed Particle-In-Cell schemes in the spirit of asymptotic preserving techniques [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF], which, as ε → 0, are consistent with the drift-limit model [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF] or the gyrokinetic model [START_REF] Filbet | Numerical simulations to the Vlasov-Poisson system with a strong magnetic field[END_REF][START_REF] Filbet | Asymptotically preserving particle-in-cell methods for inhomogeneous strongly magnetized plasmas[END_REF]. These schemes are simple and highly accurate, but the gyro-motion is lost in the limit regime. In contrast, the schemes proposed in [START_REF] Ph | Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time two-dimensional Vlasov-Poisson equation with uniform strong magnetic field[END_REF] capture all the information of the kinetic models with an accuracy uniform with respect to 0 < ε ≤ 1. These uniformly accurate (UA) schemes have computational cost as well as accuracy totally independent of ε (we refer to [START_REF] Crouseilles | Multiscale Particle-In-Cell methods and comparisons for long time two-dimensional Vlasov-Poisson equation with strong magnetic field[END_REF] for a comparison of UA scheme with other multiscale methods). In order to design UA schemes for kinetic models, different numerical approaches may be used: (i) The two-scale formulation technique relies upon an explicit separation of the fast and slow times and allows to smooth out the oscillations [START_REF] Ph | Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field[END_REF][START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time two-dimensional Vlasov-Poisson equation with uniform strong magnetic field[END_REF]. (ii) The multirevolution composition methods, in the spirit of heterogeneous multiscale method [START_REF] Abdulle | The heterogeneous multiscale method[END_REF], are also UA, as confirmed in the recent paper [START_REF] Ph | Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime[END_REF]. Both approaches exploit the periodicity of the solution of the stiff part of the equation. For instance, our recent work [START_REF] Ph | Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field[END_REF] isolates the dominant oscillation frequency owing to a confining property in two dimensions. However, the general case of a strong magnetic field with varying direction, i.e.

B(x) = (b 1 (x), b 2 (x), b 3 (x)) T ,
has been barely considered so far for the Vlasov-Poisson equation (1.1) due to its complicated highly-oscillatory behaviour in three dimensions. Let us also mention recent developments around symplectic Particle-In-Cell method, which allows for good preservation of invariants for very long time (see [START_REF] Hairer | Long-term analysis of a variational integrator for charged-particle dynamics in a strong magnetic field[END_REF][START_REF] Qin | Canonical symplectic particle-in-cell method for long-term large-scale simulations of the Vlasov-Maxwell system[END_REF][START_REF] Kraus | GEMPIC: geometric electromagnetic Particle In Cell methods[END_REF]).

In this work, we propose efficient numerical schemes for solving the three dimensional Vlasov-Poisson equation (1.1) in the strong magnetic field regime by combining multiscale strategies with the Particle-In-Cell (PIC) discretisation. First, we consider the case of a magnetic field with constant intensity |B(x)| = const, for which, as already pointed out in [START_REF] Ph | An averaging technique for transport equations[END_REF][START_REF] Saint-Raymond | The gyro-kinetic approximation for the Vlasov-Poisson system[END_REF], the motion induced by the stiff Lorentz term

1 ε v × B(x) in (1.1) is periodic in time.
Taking advantage of this observation, we derive the limit model of (1.1) by using averaging methods [START_REF] Ph | An averaging technique for transport equations[END_REF], and then introduce three UA schemes, namely (i) the multi-revolution composition (MRC) method, (ii) the two-scale formulation (TSF) method and (iii) the micro-macro (MM) method. All three are of uniform second order in time for all ε ∈]0, 1], though have specific pros and cons: for instance, MRC methods are phase-space volume preserving, while MM easily allows for the full recovery of the gyro-motion. To the best of our knowledge, this key-feature is new and paves the way for an extension to the case of a magnetic field with varying intensity. In this situation, we indeed introduce, under the PIC discretisation, a reparametrization of time to re-normalise the magnetic field. Within this framework, each particle carries its own fictitious time. Hence, and in order to avoid the occurrence of multiple frequencies, we drop in this situation the Poisson part of (1.1) and consider instead the case of an external electric field E(t, x) (this somehow simplifying assumption is relevant as it marks an important first step towards the solution of the full problem). In order to re-synchronise all particles (a necessary step in order to provide an approximation of f ε (t, x, v)), we then use the interpolation strategy of MM which ensures uniform second order except for the angular variable. Eventually, numerical experiments are presented in order to validate uniform accuracy and to compare the various methods. In particular, we simulate the dynamics of equation (1.1) in a three dimensional screw-pinch setup [START_REF] Latu | Field-Aligned Interpolation for Semi-Lagrangian Gyrokinetic Simulations[END_REF].

The remaining of the paper is now organized as follows. Section 2 considers the limit model of (1.1) with a constant intensity B(x) and Section 3 introduces the three aforementioned UA schemes in this situation: Subsection 3.1 is concerned with MRC method, Subsection 3.2 with TSF method and Subsection 3.3 with MM method. Extensions to the case of a varying intensity are presented in Section 4. Finally, numerical results with concluding remarks are exposed in Section 5.

Averaging

A general assumption throughout this paper is that the magnetic field is bounded from below, i.e. that |B(x)| ≥ c 0 for all x ∈ R 3 for some c 0 > 0 independent of ε. In this section, we further assume that the external magnetic field has constant norm

|B(x)| ≡ const > 0, x ∈ R 3 ,
so that the stiff part of equation (1.1) generates periodic motion. This setup has also been considered in [START_REF] Golse | The Vlasov-Poisson system with strong magnetic field[END_REF]. A possible instance of such a B is given by B

(x) = (B 1 (x 1 , x 2 ), B 2 (x 1 , x 2 ), B 3 (x 1 , x 2 ))
where

B 3 = B 2 1 + B 2 2 L ∞ (R 2 ) -B 2 1 -B 2 2 , with (B 1 (x 1 , x 2 ), B 2 (x 1 , x 2 )) ∈ L ∞ (R 2 ) and ∂ x1 B 1 + ∂ x2 B 2 = 0. It can be verified that |B(x)| 2 ≡ B 2 1 +B 2 2 L ∞ (R 2 ) and ∇ x •B ≡ 0.
In such a case, we are able to apply a recently developed averaging method to quickly obtain the limit model of (1.1) as ε → 0.

∂ t f ε (t, x, v) + 1 ε v × B(x) • ∇ v f ε (t, x, v) = 0, (2.1)
is 2π/b-periodic with respect to the fast time-variable t/ε.

Proof. The characteristics of (2.1)

ẋ(t) = 0, v(t) = 1 ε v(t) × B(x(t)), t > 0,
have a periodic solution in t/ε which can be obtained, for instance, by Rodrigues' formula

x(t) = x(0), v(t) = cos(bt/ε)v(0) + (1 -cos(bt/ε))(B(x(0)) • v(0))B(x(0)) -sin(bt/ε)v(0) × B(x(0)). (2.
2)

The statement of the lemma is now an immediate consequence.

Using the observation above, we may apply the following theorem from [START_REF] Ph | An averaging technique for transport equations[END_REF]:

Theorem 2.2. Consider a transport equation of form

∂ t f ε (t, y) + G(y) ε + K(y) • ∇ y f ε (t, y) = 0, f ε (0, y) = f 0 (y),
where the flow map Φ t of ẏ(t) = G(y(t)) is assumed to be 2π-periodic. There exist two formal vector fields G ε (y) and K ε (y) satisfying

G(y) ε + K(y) = G ε (y) ε + K ε (y) and [G ε , K ε ] = 0, such that the system      ∂ τ g(t, τ, y) + G ε (y) • ∇ y g(t, τ, y) = 0, (2.3a) ∂ t g(t, τ, y) + K ε (y) • ∇ y g(t, τ, y) = 0, (2.3b) g(0, 0, y) = f 0 (y). (2.3c)
has a unique formal solution independently of the order in which the equations are solved. Moreover, for all positive time we have f ε (t, y) = g(t, t/ε, y) and the first two terms of K ε = K [1] +εK [2] +O(ε 2 ) may be computed as follows

K [1] = ΠK τ , K [2] = - 1 2 Π τ 0 [K s , K τ ]ds, with K τ (y) := (D y Φ τ (y)) -1 (K • Φ τ )(y),
with Πh := 1/(2π) 2π 0 h(τ )dτ . Remark 2.3. If K ε is truncated at order k in ε and G ε = G(y)+εK(y)-εK ε , then the order in which the equations in (2.3) are solved does matter. However, the difference between the corresponding two solutions is also of order ε k .

Without loss of generality, we assume in the rest of this section that |B(x)| ≡ 1 to derive the averaged model of equation (1.1) as obtained from Theorem 2.2. Details of the derivation may be found in [START_REF] Ph | An averaging technique for transport equations[END_REF] and we thus content ourselves with a sketch of the computations: from (1.1), we have

G = 0 v × B , K = v E ,
and the flow map generated by G is given by

Φ τ (y) = x cos(τ )v + sin(τ )v × B + (1 -cos(τ ))(B • v)B , y = x v .
Then

D y Φ τ (y) = I 3 0 N 1 N 2 ,
where

N 1 = sin(τ ) (v × ∇ x B) + (1 -cos(τ ))[(B • v)∇ x B + B(v T ∇ x B)], N 2 = cos(τ )I 3 -sin(τ )(B × I 3 ) + (1 -cos(τ ))BB T
where we have denoted

v × ∇ x B = [v × ∂ x1 B, v × ∂ x2 B, v × ∂ x3 B]
and similarly for v × I 3 . A straightforward calculation then leads to

Π (D y Φ τ (y)) -1 (K • Φ τ ) = (B • v)B B v
where

B v =BB T E - 1 2 (v × ∇ x B)(v × B) + M v - 5 2 M (B • v)B - 1 2 M [v -2(B • v)B] - 1 2 B × I 3 [M (v × B) + (v × ∇ x B)(B • v)B] , with M = (B • v)∇ x B + B(v T ∇ x B). Eventually, K ε = (B • v)B B v + O(ε),
and limit model at leading order is

∂ t g(t, τ, x, v) + (B • v)B • ∇ x g(t, τ, x, v) + B v • ∇ v g(t, τ, x, v) = 0, t > 0.
By taking τ = 0 and f (t, x, v) = g(t, 0, x, v), we get the leading order averaged model of (1.1) for stroboscopic times t ∈ 2πεN,

∂ t f (t, x, v) + (B • v)B • ∇ x f (t, x, v) + B v • ∇ v f (t, x, v) = 0, (2.4a) ∇ x • E(t, x) = R 3 f (t, x, v)dv -n i , (2.4b)
f (0, x, v) = f 0 (x, v). (2.4c)
As we shall verify later numerically, we have

f ε (t, x, v) -f (t, x, v) = O(ε), 0 < ε 1, t ∈ 2πεN.

Numerical method

In this section, we introduce numerical schemes for equation (1.1) under the assumption |B(x)| ≡ 1. Taking advantage of the periodicity the solution of the stiff part, we apply state-of-art multiscale approaches in combination with PIC discretisation. In this way, we obtain schemes whose accuracy and computational cost are both independent of ε ∈]0, 1]. Our starting point is the following PICrepresentation of f ε as used in, e.g., [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF][START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF][START_REF] Jacobs | Implicit-Explicit time integration of a high-order particle-in-cell method with hyperbolic divergence cleaning[END_REF]35],

f ε (t, x, v) ≈ Np k=1 ω k δ(x -x k (t))δ(v -v k (t)), t ≥ 0, x, v ∈ R 2 . (3.1)
The characteristic equations of model (1.1) for 1 ≤ k ≤ N p are then of the form

ẋk (t) = v k (t), (3.2a) vk (t) = E(t, x k (t)) + 1 ε v k (t) × B(x k (t)), t > 0, (3.2b) x k (0) = x k,0 , v k (0) = v k,0 . (3.2c)
Noticing that

∇ x • E(t, x) = Np k=1 w k δ(x -x k (t)),
we observe that the electric field E in (3.2) has in fact no explicit dependence on time, i.e. E(t, x) = E [X(t)] (x) where X(t) = (x 1 (t), . . . , x Np (t)). We are in a position to briefly present three different UA methods.

3.1. Multi-revolution composition method. For a general exposition of multi-revolution composition (MRC), we refer to [START_REF] Ph | Multi-revolution composition methods for highly oscillatory differential equations[END_REF]. Here, we focus on a uniformly accurate second-order method. MRC framework. Suppose that we wish to solve equation (1.1) on [0, T f ] for some T f > 0. Rescaling time in (3.2) leads to (we omit the particle index for brevity)

ẋ(t) = εv(t), (3.3a)
v(t) = εE [X(t)] (x(t)) + v(t) × B(x(t)), 0 < t ≤ T f ε , (3.3b)
x(0) = x 0 , v(0) = v 0 , (3.3c)
Since the stiff part of (3.3) generates a 2π-periodic motion, equation (3.3) is amenable to MRC [START_REF] Ph | Multi-revolution composition methods for highly oscillatory differential equations[END_REF][START_REF] Ph | Convergence analysis of multi-revolution composition time-splitting pseudo-spectral methods for highly oscillatory differential equations of Schrödinger equations[END_REF]. To do so, we write

T f ε = 2πM f + T r , M f = T f 2πε ∈ N, 0 ≤ T r < 2π. (3.4)
The 2nd order MRC method begins by choosing an integer 0 < M 0 ≤ M f and defining

α = 1 2 1 + 1 M 0 , β = 1 2 1 - 1 M 0 , M = M f M 0 , H = εM 0 . (3.5) Denoting x n ≈ x(2πnM 0 ), v n ≈ v(2πnM 0 )
, the MRC scheme proceeds as follows

x n+1 v n+1 = E β (-2π)E α (2π) x n v n , 0 ≤ n ≤ M -1, (3.6)
where E α (2π) denotes the value at time 2π of the flow of

ẋ(t) = αHv(t), v(t) = αHE [X(t)] (x(t)) + v(t) × B(x(t)), (3.7)
and E β (-2π) the value at time (-2π) of the flow of ẋ(t) = -βHv(t),

v(t) = -βHE [X(t)] (x(t)) + v(t) × B(x(t)). (3.8)
The solution at final time T f is then obtained by applying to

x M v M the flow E r (T r ) at time T r of ẋ(t) = εv(t), v(t) = εE [X(t)] (x(t)) + v(t) × B(x(t)).
(3.9)

Splitting scheme. The full MRC scheme calls for the numerical evaluation of the sub-flows E α (2π), E β (2π) and E r (T r ). This is done here through a splitting, for instance of E α , in

E x α (t) : ẋ(s) = αHv(s), v(s) = 0, 0 < s ≤ t, and E v α (t) : ẋ(s) = 0, 0 < s ≤ t, v(s) = αHE [X(s)] (x(s)) + v(s) × B(x(s)). (3.10) Note that both E x α (t) and E v α (t) can be exactly integrated. The exact flow of E x α (t) is clearly x(t) = x(0) + tαHv(0), v(t) = v(0), t ≥ 0,
while the exact flow of E v α (t), by using the Rodrigues' rotation formula, can also be written explicitly

x(t) =x(0), v(t) = cos(t)v(0) + sin(t)v(0) × B + αH sin(t)E + αH(t -sin(t))(B • E)B + αH(1 -cos(t))E × B + (1 -cos(t))(B • v(0))B, t ≥ 0, where E = E [X(0)] (x(0)) and B = B(x(0)).
In our experiments, we shall take the value of the (micro) time step h = 2π/M , so that

E α (2π) ≈ (E x α (h/2)E v α (h)E x α (h/2)) M .
Approximations for E β (2π) and E r (T r ) are obtained in a similar way. It may then be proved (see (3.6)) that the error of MRC is of size O(M -2 ) for a computational cost of size M 2 , making the overall scheme of order one 1 . It remains to comment on what happens when the user-controlled M increases to the limit where M 0 reaches the critical value M 0 = 1, for which α = 1, β = 0 (3.5) and E β (-2π) ≡ id. In this case, the full MRC scheme may be regarded as just the discretisation of equation (3.3) by Strang's method with time step h. Therefore, as soon as M > 0 implies M 0 = T f /ε/(2π)/M < 1, we replace MRC method by Strang splitting with time step h = 2π/M . Finally, note that all vector fields involved in MRC are divergence free so that their exact flows are phase-space volume preserving as is the MRC method itself.

3.2. Two-scale formulation method. Two-scale formulation (TSF) methods have been developed in [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF]. Their underlying rationale is to consider the fast time as an additional variable. In order to isolate the fast time, we apply the change of unknowns (x(t), v(t)) → (x(t), y(t)) where

y(t) = cos(t/ε)v(t) + (1 -cos(t/ε))(B(x(t)) • v(t))B(x(t)) -sin(t/ε)v(t) × B(x(t)). (3.11) This leads to      ẋ(t) = F x (t/ε, x(t), y(t)), ẏ(t) = F y (t/ε, x(t), y(t)), t > 0, (3.12) x(0) = x 0 , y(0) = v 0 , where F x (τ, x, y) := cos(τ )y + (1 -cos(τ ))(B(x) • y)B(x) + sin(τ )y × B(x), 1 Under the assumption that E(t, x) ∈ C 2 (R + × R 3 ; R 3 ) and B(x) ∈ C 2 (R 3 ; R 3 ) in (3.2).
and

F y (τ, x, y) = cos(τ )E [X] (x) + (1 -cos(τ ))(B(x) • E [X] (x))B(x) -sin(τ )E [X] (x) × B(x) - 1 2 sin(2τ )q τ (y) - 1 2 (2 sin(τ ) -sin(2τ ))q τ ((B(x) • y)B(x)) - 1 2 (1 -cos(2τ ))q τ (y × B(x)) + 1 2 (2 cos(τ ) -cos(2τ ) -1)p τ (y) + 1 2 (3 -4 cos(τ ) + cos(2τ ))p τ ((B(x) • y)B(x)) + 1 2 (2 sin(τ ) -sin(2τ ))p τ (y × B(x)),
with the vector fields

p τ (z) := ((∇ x B(x)F x (τ, x, y)) • z)B(x) + (B(x) • z)(∇ x B(x)F x (τ, x, y)), q τ (z) := z × (∇ x B(x)F x (τ, x, y)), z ∈ R 3 .
Denoting u(t) = x(t) y(t) and F (τ, u) = Fx(τ,x,y) Fy(τ,x,y) , the two-scale formulation of system (3.12) now reads

∂ t U (t, τ ) + 1 ε ∂ τ U (t, τ ) = F (τ, U (t, τ )), t > 0, τ ∈ T, (3.13)
U (0, 0) = u(0),
where T = [0, 2π], and one recovers the solution of (3.12) by taking the diagonal, i.e.

U (t, t/ε) = u(t), t ≥ 0.
It remains to prescribe an appropriate initial data U (0, τ ) to (3.13) so that the solution U has its derivatives uniformly bounded up to some order.

Initial data. In order to derive U (0, τ ), we follow the Chapman-Enskog procedure. From the decomposition

U (t) = ΠU (t, •), h(t, τ ) = U (t, τ ) -U (t), with ΠU (t, •) = 1 2π 2π 0 U (t, τ )dτ, we split (3.13) into    U (t) = ΠF (•, U (t) + h(t, •)), t > 0, ∂ t h(t, τ ) + 1 ε ∂ τ h(t, τ ) = (I -Π)F (τ, U (t) + h(t, τ )), t > 0, τ ∈ T. Denote L = ∂ τ , A = L -1 (I -Π) and we have h(t, τ ) = εAF (τ, U (t) + h(t, τ )) -εL -1 ∂ t h(t, τ ).
Differentiate the above with respect to t on both sides:

∂ t h(t, τ ) = εA∇F (τ, U + h)( U + ∂ t h) -εL -1 ∂ 2 t h(t, τ ). By assuming that ∂ 2 t h = O(1) for ε ∈]0, 1], one gets ∂ t h = O(ε) and h(t, τ ) has the first order asymptotic expansion: h(t, τ ) = εAF (τ, U (t)) + O(ε 2 ). Using the fact U (0) = u(0) -h(0, 0), one gets at initial time h(0, τ ) = h 1st (τ ) + O(ε 2 ), with h 1st (τ ) := εAF (τ, u(0)),
and we denotes the first order initial data as:

U 1st (τ ) := u(0) + h 1st (τ ) -h 1st (0). (3.14)
In fact, one can show rigorously that the equation (3.13) with the well-prepared initial data

U (0, τ ) = U 1st (τ ) offers ∂ t U (t, τ ), ∂ 2 t U (t, τ ) = O(1), ε ∈]0, 1].
(3.15) We refer the readers to [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF] for the mathematical justification. The boundedness of the time derivatives (3.15) is the key to design UA schemes.

Exponential integrator. Thanks to the two-scale formulation (3.13) with the well-prepared initial data U (0, τ ) = U 1st (τ ) from (3.14), we can now directly apply the second order exponential integrator scheme proposed in [START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time two-dimensional Vlasov-Poisson equation with uniform strong magnetic field[END_REF] for integrating (3.13): Choose N τ > 0 an even integer to uniformly discretize τ on T and take a ∆t > 0 to define t n = n∆t. Denote U n (τ ) ≈ U (t n , τ ) for n ≥ 0 and let U 0 (τ ) = U (0, τ ). We update the U n for n ≥ 1 as (U )

1 l = e -il∆t ε (U ) 0 l + p l (F ) 0 l + q l 1 ∆t (F ) * ,1 l -(F ) 0 l , (3.16a)
(U ) n+1 l = e -il∆t ε (U ) n l + p l (F ) n l + q l 1 ∆t (F ) n l -(F ) n-1 l , n ≥ 1, (3.16b)
where for n ≥ 0,

U n (τ ) = Nτ /2-1 l=-Nτ /2 (U ) n l e ilτ , F n (τ ) = Nτ /2-1 l=-Nτ /2 (F ) n l e ilτ , F * ,1 (τ ) = Nτ /2-1 l=-Nτ /2 (F ) * ,1 l e ilτ ,
and

F n (τ ) = F (τ, U n (τ )), F * ,1 (τ ) = F (τ, U * ,1 (τ )) with (U ) * ,1 l = e -il∆t ε (U ) 0 l + p l (F ) 0 l , U * ,1 (τ ) = Nτ /2-1 l=-Nτ /2 (U ) * ,1 l e ilτ ,
and

p l =    iε l e -il∆t ε -1 , l = 0, ∆t, l = 0, q l =      ε l 2 ε -εe -il∆t ε -il∆t , l = 0, ∆t 2 2 , l = 0. Suppose we have the numerical solution U n (τ ) = X n (τ )
Y n (τ ) from the above scheme, then the numerical solution x n ≈ x(t n ), v n ≈ v(t n ) of the original characteristics (3.2) reads:

x n = X n (t n /ε), n ≥ 1, v n = cos(t n /ε)Y n (t n /ε) + (1 -cos(t n /ε))(B(x n ) • Y n (t n /ε))B(x n ) + sin(t n /ε)Y n (t n /ε) × B(x n ).
The derivation and convergence analysis of the above scheme can be found in [START_REF] Crouseilles | Uniformly accurate Particle-in-Cell method for the long time two-dimensional Vlasov-Poisson equation with uniform strong magnetic field[END_REF]. Since the filter (3.11) involves the magnetic field B(x), the filtered system (3.12) which is less smooth than the original form (3.2), needs more regularity for optimal convergence of the algorithm. Assuming that

E(t, x) ∈ C 2 (R + × R 3 ; R 3 ) and B(x) ∈ C 3 (R 3 ; R 3 ) in (3.
2), the two-scale formulation (TSF) exponential integrator (3.16) gives uniform second order accuracy in terms of ∆t for all ε ∈]0, 1] and uniform spectral accuracy in terms of N τ (due to periodicity):

O(∆t 2 + N -m0 τ ).
The total cost of the TSF method is O(∆t -1 N τ log N τ ).

3.3.

Micro-macro method. Now, we present the main new method of this work. It is based on the micro-macro (MM) decomposition that has been proposed very recently in [START_REF] Ph | A new class of uniformly accurate methods for highly oscillatory evolution equations[END_REF]. We shall for the first time consider this approach for the Vlasov-Poisson equation and propose a second order UA scheme. The same notations will be adopted from the previous subsection.

MM decomposition. By the averaging theory [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF], it is known that for general oscillatory problem u(t) = F (t/ε, u(t)), t > 0, (3.17) with 2π-periodicity in τ of F (τ, u), the solution can be written as a composition

u(t) = Φ t/ε • Ψ t • Φ -1 0 (u(0)), (3.18)
where Φ τ (v) is a change of variable with 2π-periodicity in τ for some v, and Ψ t (v) is the flow map of the autonomous equation with initial value v:

Ψt (v) = F 0 (Ψ t (v)), Ψ 0 (v) = v,
for some field F 0 . Though (3.18) is known to hold theoretically for some Φ τ and Ψ t , the explicit formulas of Φ τ and Ψ t are not available. In fact, by plugging (3.18) back to the equation (3.17), the change of variable Φ τ , the flow map Ψ t and the averaged field F 0 can be seen to satisfy the relation

1 ε ∂ τ Φ τ (v) + D v Φ τ (v)F 0 (v) = F (τ, Φ τ (v)), (3.19)
and moreover by taking averaging with respect to τ ∈ [0, 2π] on both sides of (3.19), one can define

F 0 with Φ τ , F 0 = (ΠD v Φ τ ) -1 ΠF (•, Φ τ ).
The above two equalities cannot completely determine Φ τ and Ψ t , so the standard averaging method imposes an extra condition ΠΦ τ = id, which uniquely defines the change of variable Φ τ in an implicit way through (3.19). In general, it is not possible to solve (3.19) to find out the exact Φ τ . However, we can define an approximated Φ τ through a kth-order iteration

Φ [k+1] τ = id + εA F (τ, Φ [k] τ ) -D v Φ [k] τ F [k] 0 , k ∈ N, (3.20)
with initially

Φ [0] τ = id, F [0] 0 = ΠF, which asymptotically gives Φ τ = Φ [k] τ + O(ε k+1
). As a compensation to the composition (3.18) by using the approximated function Φ

[k]
τ , a defect w [k] needs to be introduced:

u(t) = Φ [k] t/ε • Ψ [k] t • (Φ [k] 0 ) -1 (u(0)) + w [k] (t), (3.21)
to ensure that there are no asymptotical truncations made to the exact solution. The decomposition (3.21) is referred as the micro-macro decomposition of the solution of (3.17).

As a matter of fact, the first order approximation Φ [START_REF] Abdulle | The heterogeneous multiscale method[END_REF] τ given by the micro-macro decomposition, i.e. k = 0 in (3.20), coincides with the first order Chapman-Enskog expansion that we introduced in the previous subsection:

Φ [1] τ (r) = r + εAF (τ, r) =: Θ(τ, r). Thus, (3.12) or (3.17) has the first order micro-macro decomposition:

u(t) = Θ(t/ε, r(t)) + w(t), t ≥ 0, (3.22)
where the macro part representing the averaged equation reads

ṙ(t) = ΠF (•, Θ(•, r(t))) , t > 0, (3.23) r(0) = u(0) -εAF (τ, u(0))| τ =0 ,
and the micro part representing the equation for the defect reads

ẇ(t) = G(t/ε, r(t), w(t)), t > 0, (3.24) w(0) = εA [F (τ, u(0)) -F (τ, r(0))] | τ =0 , with G(τ, r, w) := F (τ, Θ(τ, r) + w) -(I -Π)F (τ, r) - d dt Θ(τ, r(t)).
In the first order micro-macro decomposition, the macro part (3.23) is smooth containing no highfrequencies. As for the micro part, it can be shown that [START_REF] Ph | A new class of uniformly accurate methods for highly oscillatory evolution equations[END_REF]

w(t) = O(ε 2 ), ∂ t w(t) = O(ε), ∂ 2 t w(t) = O(1)
. Thanks to the reformulation (3.12), we are able to consider this micro-macro approach for the characteristics (3.2).

An integration scheme. Now based on the micro-macro decomposed systems (3.23) and (3.24), we are going to propose a second order integration for solving (3.17) which is a compact formulation of (3.12) with u(t) = x(t) y(t) and F (τ, u) = Fx(τ,x,y) Fy(τ,x,y) . We solve the macro part (3.23) by a leap-frog finite difference scheme:

r n+1 = r n-1 + 2∆tΠF (•, Θ(•, r n )) , n ≥ 1, r 1 = r 0 + ∆tΠF •, Θ(•, r 0 ) .
For the micro part (3.24), we integrate the equation to have

w(t n+1 ) -w(t n ) = tn+1 tn G(t/ε, r(t), w(t))dt = tn+1 tn H(t/ε, r(t), w(t))dt -Θ(t n+1 /ε, r(t n+1 )) + Θ(t n /ε, r(t n )), (3.25)
where

H(τ, r, w) := F (τ, Θ(τ, r) + w) .
Since H(τ, r, w) is periodic in τ ∈ T, so we have a Fourier expansion

H(τ, r, w) = l∈Z H l (r, w)e ilτ ,
and the integration in (3.25) can be approximated as

tn+1 tn H(t/ε, r(t), w(t))dt = l∈Z tn+1 tn H l (r(t), w(t))e ilt/ε dt ≈ l∈Z tn+1 tn H l (r(t n ), w(t n )) + (t -t n ) d dt H l (r(t n ), w(t n )) e ilt/ε dt ≈ l∈Z tn+1 tn H l (r(t n ), w(t n )) + t -t n ∆t H l (r(t n ), w(t n )) -H l (r(t n-1 ), w(t n-1
)) e ilt/ε dt.

Therefore, for n ≥ 1,

w(t n+1 ) ≈ w(t n ) + l∈Z e iltn/ε α l H l (r(t n ), w(t n )) + β l ∆t H l (r(t n ), w(t n )) -H l (r(t n-1 ), w(t n-1 )) -Θ(t n+1 /ε, r(t n+1 )) + Θ(t n /ε, r(t n )),
and as for n = 0,

w(t 1 ) ≈ w(0) + l∈Z α l H l (r(0), w(0)) -Θ(t 1 /ε, r(t 1 )) + Θ(0, r(0)),
where

α l = ∆t 0 e ilt/ε dt =    iε l 1 -e il∆t ε , l = 0, ∆t, l = 0, , β l = ∆t 0 te ilt/ε dt =      ε l 2 (ε -il∆t)e il∆t ε -ε , l = 0, ∆t 2 2 , l = 0.
In total, the detailed exponential integration scheme based on the micro-macro method reads:

u n+1 =Θ(t n+1 /ε, r n+1 ) + w n+1 , n ≥ 0, (3.26a) r n+1 =r n-1 + 2∆tΠF (•, Θ(•, r n )) , n ≥ 1, (3.26b)
w n+1 =w n + Nτ /2-1 l=-Nτ /2 e iltn/ε α l H l (r n , w n ) + β l ∆t H l (r n , w n ) -H l (r n-1 , w n-1 ) (3.26c) -Θ(t n+1 /ε, r n+1 ) + Θ(t n /ε, r n ), n ≥ 1, r 1 =r 0 + ∆tΠF •, Θ(•, r 0 ) , w 1 = w 0 + Nτ /2-1 l=-Nτ /2 α l H l (r 0 , w n ) -Θ(t 1 /ε, r 1 ) + Θ(0, r 0 ), (3.26d)
r 0 =u(0) -εAF (τ, u(0))| τ =0 , w 0 = εA [F (τ, u(0)) -F (τ, r(0))] | τ =0 , (3.26e)
where N τ is an even integer to truncate the Fourier series. Suppose the numerical solution of MM is obtained as u n = x n y n , then the numerical velocity of (3.2) at t n is given as

v n = cos(t n /ε)y n + (1 -cos(t n /ε))(B(x n ) • y n )B(x n ) + sin(t n /ε)y n × B(x n ).
The micro-macro (MM) scheme (3.26) is uniformly second order accurate. In practical programming, one only needs a subroutine to evaluate F (τ, u). When the electric and magnetic field E and B in particle system (3.2) are given external functions such as polynomials, the dependence of the fast time scale t/ε (or τ ) in F and Θ can be found out explicitly and the averaging with respect to τ (through the operator Π) in the MM scheme can be pre-computed exactly. Then the MM method will have a discretisation error in time of O(∆t 2 ) with optimal computational cost O(∆t -1 ). In case that the exact evaluation of t/ε is impossible or too costly, one can always perform those computations of the fast time scale with the additional variable τ by FFT with uniform spectral accuracy thanks to the periodicity. In such case, the error bound of MM is

O(∆t 2 + N -m0 τ ),
and the total cost is O(∆t -1 N τ log N τ ), which are the same as TSF.

Full recovery of oscillation. Since MM method finds out the dependence of the fast scale in a rather explicit way, it can easily recover the complete gyro-motion of the particles, i.e. the full oscillatory trajectory of the solution of (3.2), by interpolating respectively the macro part and micro part.

Let r n and w n be the numerical solutions obtained from MM under a step size ∆t > 0. For an arbitrary t > 0, if t n < t < t n+1 , then we can use the linear interpolation to get

r n I (t) = t n+1 -t ∆t r n + t -t n ∆t r n+1 , w n I (t) = t n+1 -t ∆t w n + t -t n ∆t w n+1 .
Noting that r(t) is the averaged part and w(t) satisfies ∂ 2 t w(t) = O(1), together with the accuracy order of r n and w n from the MM scheme, it is clear that the above linear interpolation gives uniform second accuracy for approximating r(t) and w(t). Then with micro-macro decomposition (3.22), we get the interpolated numerical solution of (3.17) as

u n I (t) = Θ(t/ε, r n I (t)) + w n I (t), t n ≤ t ≤ t n+1 , (3.27)
which fully recovers the oscillation information with ease. It is direct to see

|u(t) -u n I (t)| = O(∆t 2 + N -m0 τ ).
Restart strategy. In practical long time computing, we observe that the MM scheme (3.26) could have numerical instability issue. The instability is developed from the micro part (3.24) in MM decomposition (3.22) as time evolves, since w(t) = O(ε) does not hold for arbitrary long time in general. Here, we propose a restart strategy to improve its long time performance.

Choose T 0 > 0 as the period to restart the MM decomposition. For some m ∈ N, we consider the oscillatory problem (3.17) for u m (t) = u(mT 0 + t) as um (t) = F (mT 0 /ε + t/ε, u m (t)), 0 < t ≤ T 0 .

Then we apply the proposed MM strategy on the above, which leads to MM decomposition as

u m (t) = Θ(mT 0 /ε + t/ε, r(t)) + w(t), 0 ≤ t ≤ T 0 , (3.28) with ṙ(t) = ΠF (mT 0 /ε + •, Θ(mT 0 /ε + •, r(t))) , 0 < t ≤ T 0 , r(0) = u m (0) -εAF (mT 0 /ε + τ, u m (0))| τ =0 , and
ẇ(t) = G(mT 0 /ε + t/ε, r(t), w(t)), 0 < t ≤ T 0 , w(0) = εA [F (mT 0 /ε + τ, u m (0)) -F (mT 0 /ε + τ, r(0))] | τ =0 .
The integration scheme (3.26) is then applied to solve the above two systems.

As can be seen in the numerical results later, this restart strategy for solving (3.2) is stable in long time computing. Its accuracy and computational cost are essentially the same as the direct scheme without restart.

Remark 3.1. In the case that B(x) = B 0 (x)+O(ε) with |B 0 (x)| ≡ const, all the proposed algorithms in this section can be extended to such case without any essential difficulties. Remark 3.2. Although physically the magnetic field should be divergence free, i.e. ∇ x • B = 0, all the algorithms we proposed in this section do not rely on the divergence free property of B(x) to offer the uniform accuracy.

Extension to varying intensity magnetic field

In this section, we extend previous methods to the case of Vlasov equation with a general magnetic field (whose intensity may vary), i.e.

|B(x)|

= b(x) = const, x ∈ R 3 .
We start by commenting on the difficulties to be encountered in this situation. As soon as |B(x)| = b(x) varies with x while remaining bounded from below by some c 0 independent of ε, the characteristic equation for each particle

ẋk (t) = v k (t), vk (t) = E(t, x k (t)) + 1 ε v k (t) × B(x k (t)), t > 0,
generate high oscillations. However, the dynamics of the linear part of the equation is non-periodic and thus does not allow for the application of averaging techniques. A possible remedy consists in time-reparametrisation. However, another difficulty then arises from the Poisson equation itself

∇ x • E(t, x) = Np k=1 w k δ(x -x k (t)),
which couples a huge number (N p 1) of particles with different frequencies.

Rescaling the time for each particle. Each particle has a periodic oscillation with respect to its own time s k = s k (t), given by

ṡk (t) = b(x k (t)), s k (0) = 0. (4.1)
Note that s k (t) is strictly increasing and that

s k (t) → ∞, as t → ∞, since b(x) ≥ c 0 > 0. Denoting xk (s k ) := x k (t), ṽk (s k ) := v k (t), we indeed have              d ds k xk (s k ) = ṽk (s k ) b(x k (s k )) , d ds k ṽk (s k ) = E(t(s k ), xk (s k )) b(x k (s k )) + 1 ε ṽk (s k ) × B(x k (s k )) b(x k (s k )) , s k > 0, (4.2) xk (0) = x k,0 , ṽk (0) = v k,0 ,
where the intensity of the magnetic field is scaled to one. Assuming that the electric field E(t, x) is a given external field with no ε-dependent oscillation in t, then the particle system (4.2) is decoupled for each k. Therefore, the numerical methods introduced in the previous section can all be applied to (4.2) for each particle in its own time s k with uniform accuracy. In order to build up an approximation of the function f ε (t, x, v) through (3.1), it is then necessary to re-synchronise for all particles. However, reverting s k to the physical time t is not straightforward, as it requires to numerically solve the nonlinear equation (4.1) or its equivalent for the inverse map

ṫ(s k ) = 1/b(x k (s k )), t(0) = 0. (4.3)
Given a physical time t = T > 0 (or conversely S k > 0), the best we can hope for is to determine s k (T ) (or conversely t(S k )) up to an error of size O(∆t p ) if a pth-order numerical method is applied. This source of possible error needs to be properly controlled. Here we illustrate how it can be done for the micro-macro (MM) method.

Interpolating to synchronise. For the sake of brevity, we omit k and denote

B(x) = B(x) b(x) , Ẽ(t, x) = E(t, x) b(x) .
We filter (4.2) and ( 4.3) as before by introducing ỹ(s) := cos(s/ε)ṽ(s) + (1 -cos(s/ε))( B(x(s)) • ṽ(s)) B(x(s)) -sin(s/ε)ṽ(s) × B(x(s)), (4.4) and obtain

                   d ds x(s) = Fx (s/ε, x(s), ỹ(s)), d ds ỹ(s) = Fy (s/ε, x(s), ỹ(s)), s > 0, (4.5) d ds t(s) = 1 b(x(s)) , x(0) = x 0 , ỹ(0) = v 0 , t(0) = 0,
which has now the appropriate format (3.17). Here Fy is defined similarly as in (3.12) (see Section 3.2) with the scaled vector fields Ẽ, B and Fx = F x /b(x). We then solve system (4.5) with the MM scheme (3.26) with time step ∆s > 0 and denote t n ≈ t(s n ) the numerical solution of t(s n ) at s n = n∆s. Then, using the notations r := (r x , r y ) and w = (w x , w y ) for the macro and micro parts (see Section 3.3), the numerical solution of x and ỹ at s n = n∆s is

xn := r n x + εA Fx (s n /ε, r n x , r n y ) + w n x ≈ x(s n ), ỹn := r n y + εA Fy (s n /ε, r n x , r n y ) + w n y ≈ ỹ(s n ). Assume that b(•) ∈ C 1 (R 3 ) and 0 < c 0 ≤ b(x) ≤ C b for all x ∈ R 3 for some C b > 0. Then, from t n ≥ t(s n ) -|t n -t(s n )| ≥ s n C b -C∆s 2
we see that whenever ∆s > 0 is small enough, the value of t n will eventually become greater than any arbitrary positive value. For a given final time T > 0, we thus stop the algorithm when t n ≤ T ≤ t n+1 . Note that the function t(s) satisfies

d 2 ds 2 t(s) = O(1), 0 < ε ≤ 1,
so that we can interpolate the value of t(s) from t n and t n+1 with second order uniform accuracy:

θ := T -t n+1 t n -t n+1 , T = θt n + (1 -θ)t n+1 , s * = θs n + (1 -θ)s n+1 .
Interpolation is further used to obtain

r * x = θr n x + (1 -θ)r n+1 x , r * y = θr n y + (1 -θ)r n+1 y , w * x = θw n x + (1 -θ)w n+1 x , w * y = θw n y + (1 -θ)w n+1 y .
As stated in Section 3.3, all functions used above in the interpolation have uniformly bounded second order derivative. As a consequence, the so-obtained approximations are uniformly second order. Eventually, the numerical solutions of (4.5) at s = s(T ) are given by

x(s(T )) ≈ x * , ỹ(s(T )) ≈ ỹ * , with x * := r * x + εA Fx (s * /ε, r * x , r * y ) + w * x , ỹ * := r * y + εA Fy (s * /ε, r * x , r * y ) + w * y . (4.6)
Note that the dependence in the fast-time s/ε within the MM method only appears in the O(ε)terms. As a consequence, the approximation errors of x(s(T )) and ỹ(s(T )) by (4.6) are still of uniform second order, although an error is introduced on s * /ε owing to |s(T ) -s * |/ε = O(∆s 2 /ε).

To reconstruct an approximation of the distribution function f ε (T, x, v), we need x(T ) and v(T ). For the position variable, we directly have x(T ) = x(s(T )) due to the definition. As for the velocity variable v(T ), we need to invert the change of variable (4.4), where the fast scale s/ε occurs in some O(1)-terms. However, the parallel component v := (B • v)B/ B 2 of the velocity as well as |v| do not suffer from the same problem, thanks to the following observations (let us recall that B( Therefore, the strategy proposed in this section is of overall uniform second order for the computation of

x(t), v (t), |v(t)|, t ≥ 0.
This, in turn, allows for a uniformly accurate approximation of macroscopic quantities such as the density or the kinetic energy

ρ ε (t, x) := R 3 f ε (t, x, v)dv, ρ ε v (t, x) := R 3 |v| 2 f ε (t, x, v)dv,
as well as the magnetic moment [START_REF] Hairer | Long-term analysis of a variational integrator for charged-particle dynamics in a strong magnetic field[END_REF][START_REF] Northrop | The adiabatic motion of charged particles[END_REF] µ ε (t) :=

R 3 R 3 f ε (t, x, v) |v ⊥ | 2 |B(x)| dxdv, (4.7)
with v ⊥ := vv .

Numerical results

This section is devoted to present the numerical results from the proposed numerical schemes. We shall firstly test and compare the accuracy, efficiency and long time performance of the schemes considering a single test particle for some three dimensional simulations in the two following cases: constant intensity and varying intensity magnetic field. Then, we shall focus on the nonlinear Vlasov-Poisson case under the influence of a constant intensity magnetic field.

5.1. Accuracy study. We investigate the performance of the proposed numerical methods by considering a single particle system in three dimensions:

ẋ(t) = v(t), v(t) = E(x(t)) + 1 ε v(t) × B(x(t)), t > 0, (5.1)
x(0) = x 0 , v(0) = v 0 .
We take two 3D vector fields B(x) : x ∈ R 3 → R 3 and E(x) = -∇ x φ(x) with some φ(x) : R 3 → R, then (5.1) is an Hamiltonian system with the energy conserved as

H s (t) := 1 2 |v(t)| 2 + φ(x(t)) = H s (0), t ≥ 0. (5.2)
Note here we do not require B being divergence free for these accuracy tests, since all the presented properties of the proposed schemes hold in general as long as B is a smooth enough vector field.

Hence, we first focus on a constant intensity magnetic field before considering the general case to test the proposed methods MRC, TSF and MM.

Example 5.1. (Constant intensity). We take the two external fields in the system (5.1)

E(x) =    cos(x 1 /2) sin(x 2 ) sin(x 3 )/2 sin(x 1 /2) cos(x 2 ) sin(x 3 ) sin(x 1 /2) sin(x 2 ) cos(x 3 )    , B(x) =    sin(x 1 + x 2 ) cos(x 1 + x 2 ) sin(x 3 ) cos(x 1 + x 2 ) cos(x 3 )    , x = (x 1 , x 2 , x 3 ),
where |B(x)| = 1 and E(x) derives from the potential φ(x):

E(x) = -∇ x φ(x), φ(x) = -sin(x 1 /2) sin(x 2 ) sin(x 3 ).
We choose the initial data for (5.1) as

x 0 = (1/3, -1/2, √ π/2), v 0 = (1/2, e/4, -1/3).
A reference solution is obtained by using the fourth order Runge-Kutta method with small step size ∆t = 10 -5 . We firstly study the convergence of the three proposed methods (MRC, TSF and MM) aiming to illustrate their uniform accuracy for all ε ∈]0, 1]. To do so, we solve the system under different ε till T = π/2 and compute the error

error = |x(T ) -x num | |x(T )| + |v(T ) -v num | |v(T )| , (5.3)
where x num and v num are the numerical values obtained by the different schemes. For the TSF and MM methods, we define the time step ∆t = T /M with M ∈ N and we fix the grid points for τ -direction as N τ = 32. For MRC, we define the numerical parameters from a given M ∈ N as follows:

H = εM f /M and h = 2π/M if M f /M ≥ 1 and ∆t = 2π/M if M f /M < 1 (M f being defined by (3.4)).
The error (defined by (5.3)) produced by the three methods at T = π/2 with respect to the number of (macro) grid points M or with respect to ε is given in Figure 1. As expected, the three methods enjoy the uniform second order accuracy property since the rate of convergence is essentially insensitive to the value ε ∈]0, 1]. The typical behavior of uniformly accurate methods can be observed on the error as a function of ε: the curves obtained for different M are almost parallel. Note that the results obtained by TSF and MM are very close whereas the error produced by MRC becomes smaller when ε decreases.

In Figure 2 we look at the error of TSF and MM with respect to the number of grid points N τ in the auxiliary variable τ (the time step is fixed to ∆t = 10 -5 ). This error is important to study since these two methods involve an additional variable τ which may make them less competitive. We can see in Figure 2 that the error decreases spectrally as the number of grid points N τ increases. Moreover, for small values of ε, a very small number of N τ is needed to reach high accuracy: ε ≤ 2 -7 , N τ = 16 is enough for machine precision. Finally, let us remark that the results obtained for MM is much less sensitive than TSF: when ε = 1/2, N τ = 32 enables to reach machine precision for MM whereas TSF requires N τ = 128.

We now intend to compare the efficiency of TSF, MM and MRC in different regimes (ε = 1/2 and 1/2 14 ). Let us first fix the numerical parameters. According to the previous comments, in the regime ε = 1/2 we take N τ = 128 for TSF and N τ = 32 for MM whereas in the regime ε = 1/2 14 , we take N τ = 8 for both TSF and MM. We test the long time behavior of the three methods by investigating the relative error on the numerical total energy defined by

error(t n ) = |H n s -H s (0)| |H s (0)| , H n s = 1 2 |v n | 2 + φ(x n ), (5.4)
where H n s is the numerical approximation of H s (t n ) given by (5.2). We plot in Figure 3 the error (considering the maximum of (5.4) among all the iterations) against the computational time of the three methods for ε = 1/2 or 1/2 14 (different time steps have been chosen). For a given error, when ε = 1/2 the MRC method is more efficient than TSF or MM, but it is no longer true when ε becomes smaller. This is explained by the fact that N τ can be chosen smaller in the asymptotic regime, making TSF and MM more competitive. MRC for ε = 1/2 reads as the Strang splitting, while for ε = 1/2 14 the convergence of MRC becomes first order in terms of total computational cost. Then, in Figure 4, the time history of (5.4) is plotted for the three methods till T = 32π, for ε = 1/2 14 . The TSF and MM methods run with N τ = 32 and M = 1024 (∆t = 0.098) or M = 2048 (∆t = 0.049). We report that, for this test, MM becomes unstable in large time so that the restart strategy is used every T 0 = 8π. For MRC, we used M = 64 or M = 128. Figure 4 clearly shows that MRC has the best long time behavior among the three methods. Indeed, TSF and MM has a linear drift in the energy error as time evolves whereas for MRC, it remains of the same order (about 10 -5 ) for large time. Let us remark that the energy error converges quadratically for the three methods with respect to number of time grid points M .

Finally, we consider the scheme MM (M = 32, N τ = 32) to illustrate the reconstruction of the whole trajectory for t ∈ [0, π] (so that ∆t = 0.0982). To do so, we still consider the system (5.1) with example 5.1, with ε = 1/2 5 . In Figure 5, we plot a reference trajectory (obtained with a very small time step) and the numerical solution obtained by MM using the strategy proposed in Section 3.3 (i.e. with a coarse time grid and using the linear interpolation strategy in (3.27)). Using a few grid points, we can see that the MM method is able to fully restore the complex trajectory (highly oscillatory confined behavior around magnetic field line) of the particle under trivial computational cost.

Example 5.2. (Varying intensity) Secondly, we investigate the numerical performance of the strategy proposed in Section 4 for a magnetic field with varying direction and varying intensity on the particle system (5.1).

We shall consider the particle system (5.1) with the same electric field E(x) as before, but here the magnetic field is

B(x) =    1 -sin(x 2 )/2 1 + cos(x 3 )/2 1 + cos(x 1 )/2    , x = (x 1 , x 2 , x 3 ), which satisfies ∇ x • B = 0 but has a varying intensity in x ∈ R 3 since |B(x)| 2 = 3 + cos(x 1 ) + cos(x 3 ) -sin(x 2 ) + cos(x 1 ) 2 /4 + sin(x 2 ) 2 /4 + cos(x 3 ) 2 /4.
We choose the same initial data as before for example 5.1 and solve the problem via the new time formulation (4.2) with the MM method (see Section 4). The reference solution is again obtained by directly solving (5.1) with the fourth order Runge-Kutta method under small step size (∆t = 10 -5 ).

First, we are interested in the error (defined by (5.3)) against the number of grid points M for the quantities x(t), v (t) := v(t) • B(x(t))B(x(t)/ B(x(t) 2 and |v(t)| at T = 1. In Figure 6, we can observe that the proposed MM scheme converges as number grid points M increases (∆s decreases) with uniform second order accurate rate for all ε]0, 1].

Then, in Figure 7, the time history of the energy error (defined by (5.4)) of the method with ∆s = 1/8 and ∆s = 1/16 till a physical time T = T (s) = 100 and under three different ε is shown. Let us remark that the restart strategy is used at every time step. In Figure 7, we observe that the scheme computes the energy (5.2) with uniform second order accuracy for ε ∈]0, 1]. Under a rather large step size (∆s ε), the scheme is stable in long time computing, and even if a slight linear drift in the energy error is observed, the energy error (about 10 -3 ) is rather good for all ε considered. In Figure 7, the relation between the new time s and the physical time t(s) is also plotted to illustrate that the physical time t(s) is a monotone increasing function.

Finally, we study the dynamics of the magnetic moment defined by

I(t) = 1 2 |v ⊥ (t)| 2 |B(x(t))| , (5.5)
which is an analogy of (4.7) at the particle level (5.1). We use MM with ∆s = 1/16 (so that it is accurate enough) to solve (5.1) till t = 100 with three different values of ε (ε = 2 -9 , 2 -10 , 2 -11 ). in the simulation, which is consistent with the results obtained in [START_REF] Hairer | Long-term analysis of a variational integrator for charged-particle dynamics in a strong magnetic field[END_REF]. Our scheme captures this adiabatic quantity even when ∆s ε whereas the scheme used in [START_REF] Hairer | Long-term analysis of a variational integrator for charged-particle dynamics in a strong magnetic field[END_REF] needs ∆s < ε.

5.2.

Simulation of the Vlasov-Poisson system. In this last part, we focus on the numerical simulation of the full 3D Vlasov-Poisson equation (1.1) using the MRC method. The chosen initial data is a Maxwellian in velocity and a ring-shape distribution in space with a perturbation in angle [START_REF] Filbet | Numerical simulations to the Vlasov-Poisson system with a strong magnetic field[END_REF]:

f 0 (x, v) = n 0 2π
(1 + η cos(kθ)) e -5(r-5) 2 e - x 3 -direction. As a diagnostic, we consider the following quantity:

ρ ε (t, x) = R 3 f ε (t, x, v)dv, x ∈ Ω.
For the PIC method, we choose N p = 100 × N x1 N x2 N x3 particles and the projection of the particles on the spatial grid is done by cubic splines. In Figures 9, the density ρ ε is displayed at different times for ε = 1/2 5 with M = 256 whereas α = 0 in the magnetic field, so that B is homogeneous and aligned with the x 3 direction. There is two different dynamics which can be seen in the results: an instability develops in the direction orthogonal to the magnetic field (one can see four vortices at time t = 64π) and a slight parallel dynamics in the plane parallel to the magnetic field. In Figures 10 and11, a non-homogeneous magnetic field is considered (α = 0.003). We can observe that the dynamics is different dynamics from the homogeneous case. Indeed, the instability leading to the formation of four vortices is different and one can see stronger non homogeneous phenomena in the x 3 direction due to the expression of the magnetic field.

Finally, in Figure 12, we plot the energy error for both configurations (α = 0 and α = 0.003). Very good conservations are obtained for long time. Moreover, we consider the relative error between the Vlasov-Poisson system (1.1) and the asymptotic model (2.4) as a function of ε, for α = 0.003. To do so, we compute the L ∞ norm (in space) of |ρ ε (t = π, x) -ρ(t = π, x)|/|ρ ε (t = π, x)| at the final time t = π. We can see that when ε decreases, the error is O(ε), as predicted by the theory.

Lemma 2 . 1 .

 21 If |B(x)| ≡ b for some constant b > 0, then the solution of

  ) = ( B(x(s)) • ỹ(s)) B(x(s)), |ṽ(s)| = |ỹ(s)|, which allow to get v (T ) ≈ ( B(x(s(T ))) • ỹ(s(T ))) B(x(s(T ))), |v(T )| = |ỹ(s(T ))|.

Figure 1 .

 1 Figure 1. Errors of MRC, TSF and MM with respect to time steps M under different ε (left) or with respect to ε under different M (right) for example 5.1.

Figure 2 .

 2 Figure 2. Error of TSF and MM with respect to N τ under different ε for example 5.1.

Figure 3 .Figure 4 .Figure 5 .

 345 Figure 3. Efficiency comparison of TSF, MM and MRC in classical (left) or asymptotic regime (right) of ε for example 5.1: error versus computational time.

Figure 6 .

 6 Figure 6. Error of MM under different M = ∆s -1 in x, |v| and v at T = 1 in example 5.2 of varying intensity.

1 2 |v| 2 , ( 5 . 6 ) 1 1 + α 2 x 2 1 + α 2 x 2 14 Figure 7 . 11 Figure 8 .

 256112147118 Figure 7. Energy error of MM (restart each step) for ε = 1/2, 1/2 5 , 1/2 14 till t = 100 and the evolution of t(s) in example 5.2.

Figure 9 .

 9 Figure 9. Vlasov-Poisson case: pseudo-color snapshots of ρ ε under ε = 1/2 5 at t = 0, 16π, 32π, 64π with initial condition 5.6 with α = 0.

Figure 10 .

 10 Figure 10. Vlasov-Poisson case: pseudo-color snapshots of ρ ε under ε = 1/2 5 at t = 0, 16π, 32π, 64π with initial condition 5.6 with α = 0.003.

Figure 11 .Figure 12 .

 1112 Figure 11. Vlasov-Poisson case: pseudo-color snapshots of ρ ε under ε = 1/2 5 at t = 88π, 128π with initial condition 5.6 with α = 0.003.

Acknowledgements This work is supported by the French ANR project MOONRISE ANR-14-CE23-0007-01. This work has been carried out within the framework of the French Federation for Magnetic Fusion Studies (FR-FCM) and of the Eurofusion consortium, and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.


Issam Raïs

Daniel Balouek-Thomert

Anne-Cécile Orgerie

Laurent Lefèvre

Manish Parashar

Leveraging energy-efficient non-lossy compression for data-intensive applications

Keywords: Big Data, Compression, Energy efficiency, HPC, Datacenter

The continuous increase of data volumes poses several challenges to established infrastructures in terms of resource management and expenses. One of the most important challenges is the energy-efficient enactment of data operations in the context of data-intensive applications. Computing, generating and exchanging growing volumes of data are costly operations, both in terms of time and energy. In the late literature, different types of compression mechanisms emerge as a new way to reduce time spent on data-related operations, but the overall energy cost has not been studied. Based on current advances and benefits of compression techniques, we propose a model that leverages nonlossy compression and identifies situations where compression presents an interest from an energy reduction perspective. The proposed model considers sender, receiver, communications costs over various types of files and available bandwidth. This strategy allows us to improve both time and energy required for communications by taking advantage of idle times and power states. Evaluation is performed over HPC, Big Data and datacenter scenarios. Results show significant energy savings for all types of file while avoiding counter performances, resulting in a strong incentive to actively leverage non-lossy compression using our model.

I. INTRODUCTION

Energy consumption is a growing global concern. With the multiplication of connected devices per person around the world, reducing the energy consumption of large scale computing system is a mandatory step to address in order to build a sustainable digital society.

To face this growing concern many solutions have been developed at multiple levels of computing facilities: hardware, middleware, and application. Examples of hardware leverages are Dynamic Voltage and Frequency Scaling (DVFS) [START_REF] Suleiman | Dynamic voltage frequency scaling (dvfs) for microprocessors power and energy reduction[END_REF] or ShutDown Techniques [START_REF] Benoit | Reducing the energy consumption of large-scale computing systems through combined shutdown policies with multiple constraints[END_REF]. At the middleware level, energy efficient policies for task, jobs and resource managers can also be labeled as leverages [START_REF] Georgiou | A scheduler-level incentive mechanism for energy efficiency in hpc[END_REF]. Last, at the application level, the way an application is implemented has important effects on its energy and power consumption [START_REF] Raïs | Exploiting the table of energy and power leverages[END_REF].

Nowadays applications tend to generate very large amount of data. Thus, data centric computing has raised since late years in the literature. In a data driven world, the convergence of problematics like HPC and big data are inevitable. Scientists are now struggling to deal with all these generated data on the fly.

One promising technique that currently rises in literature is the lossy compression, that lightly pre-processes data to know if it is relevant to take it into account in the computing process [START_REF] Calhoun | Exploring the feasibility of lossy compression for pde simulations[END_REF]. This leverage is promising but implies a new way of pre-processing data for every application. Thus, this leverage is not generic and has to be evaluated in terms of energy consumption.

In this paper, we tackle the problem from a different perspective. We aim at underlining the relevance of non lossy compression process in various contexts from an energy efficient perspective. By using such a leverage through our proposed solution, we expect a reduction in the overall energy consumption as well as a reduction in response time.

The HPC -BigData context is interesting for the compression leverage. Indeed, workflows made of applications that generate and more importantly exchange a lot of data, creating contention on the network, are good candidates for such a technique. We propose a model that gives binary answers to the usage of compression and decompression to improve energy efficiency of data exchanges.

We evaluate our approach on different representative simulated scenarios with extracted traces from HPC and BigData applications running on Titan Supercomputer, on datacenter and Edge facilities.

This paper is structured as follow. Section II will present the compression model as an energy leverage. Section III will focus on the experimental platform while Section IV will analyze evaluation results. Section V presents an overview of selected related works. Section VI presents some conclusion and future works.

II. COMPRESSION MODEL AS ENERGY LEVERAGE AND

MODELISATION

The aim of the proposed compression model is to answer the following question: "Is it energy efficient to compress a file on a specific sender node and decompress it on a receiver node for a given size, file format, and available bandwidth?".

Figure 1 explains the comparison of energy costs addressed by the model. Figure 1a presents the cost of the model with compression and decompression on sender and receiver nodes, For our model, on the sender side, sending the file starts only when the compression is done. Thus the reception starts only when the compression is finished, on the receiver side. Decompression starts only if the file is completely received, on the receiver node. The energy consumed by the sender and the receiver during these actions is compared to what would happened during the exchange of the original file without the steps of compression and decompression, as displayed in Figure 1b.

CompressionEnergyCost(Green)

+ EmissionCompressedF ileEnergyCost(Red) + ReceptionCompressedF ileEnergyCost(Purple) + DecompressionEnergyCost(Blue) < EmissionF ileEnergyCost(Orange) + ReceptionF ileEnergyCost(Yellow)
Thus the model is summed up with the previous formula (in parenthesis the corresponding color in Figure 1) where the answer is a binary statement. The left part of the equation corresponds to the model Figure 1a while the right side corresponds to Figure 1b.

A. Modeling energy cost of compression and decompression phases

To understand the implication of compression and decompression process to energy consumption, we used a compression tool on various file formats and for each of these, various file sizes.

Table I presents results of energy consumption of compression and decompression of pbzip tool with various formats and sizes of file on Nova server nodes (from the Grid'5000 experimental testbed described on Table II).

The compression rate is computed by dividing the original size of a file by the size obtained after compression. We noticed that for every experience, on every file structure and size of file generated, the compression rates are very close (the worst case is 1.3185% for random binary files and the best case is 250 % for XML file), showing a high correlation between the file structure and the compression rate.

The proposed energy consumption prediction model is thus based on the size and the structure of the file. This prediction model is a linear regression based on results obtained for various explored sizes. Every file structure has its own regression model for compression and decompression. Thus, for a given size and a file structure, the regression tells us how much energy is consumed to compress or decompress the file. On a production scenario, this regression could be built on the fly and refined every time a file is compressed and decompressed.

Figure 2 shows results of energy consumption for compression in function of size of various file structures. The obtained linear regression and actual measurements are represented for every case as a green and blue lines, respectively. Every energy consumption obtained is very close to the obtained linear reduction.

B. Modeling energy cost of emission and reception phases

The energy cost during a file exchange between two nodes results from the usage of three major components: the sender, For simplicity, we consider that the cost of network is, at full or idle load, almost constant [START_REF] Ye | Analysis of power consumption on switch fabrics in network routers[END_REF]. Thus, network will have a constant energy cost here. We place this leverage in the case of contention in the network. Therefore, available bandwidth is also a parameter for our answer. The idle power consumed by the sender and receiver are also parameters of our model. We only take into account point-to-point communications.

For the sender and the receiver, the best scenario is when the nodes send and receive the expected data while consuming the same energy as being in their idle state. Indeed, it is, after its shutdown state, the lowest reachable power state. Thus, from an energy point of view, to determine if it is beneficial to use compression and decompression or send as generated, we compare the best case (idle for the initial size on both sides) with the energy cost exposed by our regression on compression and decompression for the given file structure and size. Thus, needed information to estimate this best case scenario is the idle power for both sender and receiver, the available bandwidth during the transit of the file, and the size of the raw file to send.

With this information, the model answers yes to the question if the cost of idle during compression sending and decompressing is lower than sending raw generated file over the contended network, from an energy point of view.

III. EXPERIMENTAL SETUP

In order to evaluate our model, we explore parameters that consider different bandwidths, sizes of data exchanges and energy consumption of nodes. In this section, we describe the infrastructure, parameters and traces that were used.

A. Calibration architecture

Our evaluation is carried out on the Grid'5000 infrastructure. Grid'5000 is a large-scale and versatile testbed for experiment-driven research in all areas of computer science, with a focus on parallel and distributed computing including Cloud, HPC and Big Data [START_REF] Balouek | Adding virtualization capabilities to the Grid'5000 testbed[END_REF]. It provides a software stack, including management tools and monitoring frameworks. This includes the Kwapi utility that allows a user to monitor and profile the power consumption of platform nodes (servers and network equipment) [START_REF] Clouet | A Unified Monitoring Framework for Energy Consumption and Network Traffic[END_REF]. Kwapi gives us access to external wattmeters on the Grid'5000 Lyon site, that monitor entire nodes with a 0.125 Watts accuracy. Characteristics of nodes used to calibrate the energy consumption model, named Nova, are given in Table II.

B. Simulation parameters

Bandwidth We define three evaluation scenarios for the bandwidth parameter: Datacenter, SuperComputer and Edge scenario. The Datacenter scenario relies on the Gigabit Ethernet as the network technology as actually implemented on Grid'5000. We simulate a randomly equally distributed contention, varying between 1 to 100 in order to factor the maximal capability of this technology. The SuperComputer scenario is based on the specifications of Cray Titan XP7 1 , using a bandwidth of 100 Gb/s between each computing node. We make the hypothesis that the quality of links between the computing nodes and the data nodes are the same as the Datacenter scenario. For the Edge scenario, a 10 Mbps link is used for edge links, which according to [START_REF] Akamai | State of the internet report[END_REF] corresponds to a higher connectivity than the average connection speed in USA. Due to the high instability of the connection over Internet, we raise the maximal contention factor to 1000. Size of exchanges We define three scenarios for the size of exchanges: Regular, S3D and BigData size. Regular size is based on the settings used to generate the regression models:

1 https://www.top500.org/system/177975 size varies from 250MB to 2GB files, chosen randomly for a fixed number of exchanges. The S3D size replays the data exchanges of the S3D application on Titan Cray XK7. S3D is a direct numerical simulation of combustion flow that runs in production. Following the obtained traces, every exchange between compute and data nodes is around 0.2 MBytes (exactly 238328 Bytes). The BigData size relies on [START_REF] Bao | A data colocation grid framework for big data medical image processing-backend design[END_REF]. In this paper, authors describe pictures exchanges of size between 6MB and 20MB, per chunks of 30 to 160. Thus, we used values comprised between 960 MB and 3.2 GB as parameters.

Power and Energy characteristics The energy consumption of the nodes that send and receive the data is split in two parameters. The Nova parameter corresponds to the idle value of energy consumption of the chosen node, as described in Table II. The average monitored value is equal to 100 W. The Little parameter corresponds to a low-consumption processor present on smartphones or edge devices. We used a value of 4 W, as referred in [START_REF] Raïs | Quantifying the Impact of Shutdown Techniques for Energy-Efficient Data Centers[END_REF] where authors evaluated the specifications of a Big-Medium-Little architecture composed of 3 different processors to adapt computing needs to demand. A summary of the simulation parameters is given in Table III.

C. Simulation Scenarios

The simulation scenarios combine the parameters described in the previous section. Such scenarios are representative of datacenter management, Big data applications and HPC applications. A summary of the simulation scenarios is given in Table IV.

The scenario A describes what could happen during the lifetime of a datacenter, when exchanging files at application level or between nodes in a point-to-point fashion. It represents the exchange between two nodes with homogeneous architectures. The scenarios B and C describe the file exchanges of an application performing large scale medical image processing (as presented in [START_REF] Bao | A data colocation grid framework for big data medical image processing-backend design[END_REF]) where huge amounts of data are exchanged, during every iteration of the application. B represents the exchange between two regular nodes belonging to the same supercomputer, while C considers a regular node and a low-consumption node at the edge of the network. The scenario D describes the file exchanges of the S3D numerical simulation [START_REF] Chen | Terascale direct numerical simulations of turbulent combustion using s3d[END_REF] based on a collected trace where small amounts of data are exchanged, during every iteration of the application. For every file sent, the model takes into account the linear regression to predict the cost of compression and decompression. From this predicted cost, the actual bandwidth and the node power characteristics, the model answers if it is energy efficient to use compression on the sender node and decompression on the receiver for the given format and size of file.

IV. EVALUATION

In this section, we evaluate the different scenarios previously described. For each scenario, we measure three different values over 10,000 data exchanges: (i) No Compression : the baseline energy consumption value without any use of the model, (ii) Active model : using the model as a leverage when it benefits a communication, (iii) Always ON : where the compression is always used without the model regulating the choices.

For all file formats and scenarios, activating the model of compression is always beneficial from an energy perspective. For the scenario A, and for the most complicated file format to compress, the Random binary files, the percentage of gain is around 8, 6%. Although, for very compressible files, like SQL files, we noticed excellent percentage of energy savings (up to 79%). Intermediately structured files also show significant energy savings, 19% and 31% for BP and XML files, respectively.

Regarding the always active leverage on the A scenario, one can observe that for every file format, the proposed model outperforms the energy gains of this configuration. For example, almost 19% and 24% of difference concerning energy gain are seen for XML and BP files, respectively. Note that for BP files, the Always ON results in a negative percentage of gain compared to the No Compression scenario as some files should not have been compressed. In other words, not using compression would have been better than always using it, but using the proposed model permits energy savings.

For all scenarios from SuperComputer to Edge configurations, the model permits to save important amounts of energy. In the worst case, scenarios B and for random binary files, 14.48% of energy is gained using our model. In the best cases, scenarios C with XML files and D with BP file, 97.95% and 97.30% of energy is gained compared to the same exchanges without model. This evaluation over different file formats and extensive scenarios demonstrates the relevancy of nonlossy compression from an energy consumption perspective. As seen, the best case is not always the same. We underline the fact that, under all scenarios and studied file format, the proposed model enables significant energy savings while ensuring no counter performances.

V. RELATED WORK

The impact of lossless compression techniques on energy consumption has been studied in the literature. Barr and al. demonstrate that, using several typical compression algorithms, there is an energy consumption increase when compression is applied before wireless transmission [START_REF] Barr | Energy-aware lossless data compression[END_REF]. On the contrary, in [START_REF] Xu | Impact of data compression on energy consumption of wireless-networked handheld devices[END_REF], authors investigate the need of compression to reduce the energy consumed to save battery on hand-held devices. Their results show that on key cases for devices downloading files, the gain of energy consumption when using lossless compression could be interesting. Welton and al. allocate idle CPU resources to compress network traffic, thus reducing the amount of data transferred over the network and increasing the effectiveness of network bandwidth [START_REF] Welton | Improving i/o forwarding throughput with data compression[END_REF]. In [START_REF] Lee | Improving energy efficiency of gpus through data compression and compressed execution[END_REF], the authors explore register file data compression for GPUs to improve power efficiency. Compression reduces the width of the register file during read and write operations, which in turn reduces dynamic power.

Satish and al. show that current cluster implementations suffer from high latency data communication with large volumes of transfers across nodes, leading to inefficiency in performance and energy consumption [START_REF] Satish | Large-scale energyefficient graph traversal: a path to efficient data-intensive supercomputing[END_REF]. Authors conclude that these constraints can be overcome by using a combination of efficient low-overhead data compression techniques to reduce transfer volumes along with latency-hiding techniques.

Very recent papers [START_REF] Di | Optimization of error-bounded lossy compression for hard-to-compress hpc data[END_REF]- [START_REF] Foster | Computing just what you need: online data analysis and reduction at extreme scales[END_REF] have studied lossy compression as a key leverage to reduce the data to process. It consists in preprocessing data using detectors to detect if data is interesting to be sent to the computing kernel. This differs from the focus of this paper, as we propose a generic model having no input nor strong domain-specific knowledge from the application, resulting in no insights on what is the relevance of the data exchanged to allow losses.

Despite several investigations on using lossy or non lossy compression as leverage to consume less energy, none of the previous papers presented a simple model that could answer, at a given time, if the compression could be beneficial from a current state of the system (size of file currently exchanged, actual available bandwidth, energy consumption of nodes) and applied it to the complexity of datacenter, supercomputer and edge systems. Our approach underlines a model that helps an automated system to determine, at runtime, if compressing and decompressing the needed file could lead to a reduction of the energy consumption. It provides a straightforward binary answer to this problem and proves to be easily integrable into third-party systems.

VI. CONCLUSION AND FUTURE WORK

Evaluating the energy efficiency and its impact on application performance becomes increasingly challenging when computation units are a moving target and data exchange are growing in size and diversity. It is of paramount importance to address data centric applications as the convergence of HPC and Big Data domain to design and leverage energy efficient techniques. In this work, we propose a model that determines the relevancy of using the non-lossy compression leverage on a target setting, for energy efficiency purpose.

Evaluation highlights compression performance of various heterogeneous workloads while considering supercomputer nodes, edge devices and datacenter specifications. We show that besides enabling energy savings without counter performance, our model is also useful to explore the optimization opportunity in infrastructures and system design. Future work includes the extension of this model to stream processing technologies. Compression of machine-generated data from Edge and IoT devices would extend the battery life of sensors. Also, by including the latency criteria into the compression leverage, the model would contribute to the cost-efficient enactment of stream processing topologies under changing data volume.

Fig. 1 :

 1 Fig. 1: Model of energy cost with or without compression and decompression

Fig. 2 :

 2 Fig. 2: Compression energy consumption for various size and linear regression

TABLE I :

 I Energy consumption of pbzip under various sizes' files

		Energy cost (Joules)
	Size (Bytes)	Compression Decompression Compression Rate
		Rand Binary Files
	250MB	286.0	192.0	1.3185
	500MB	826.0	322.0	1.3185
	750MB	988.0	631.0	1.3185
	1.0GB	1468.0	650.0	1.3185
	1.25GB	1810.0	807.0	1.3185
	1.50GB	2160.0	945.0	1.3185
	1.75GB	2402.0	1291.0	1.3185
	2.00GB	2869.0	1324.0	1.3185
			BP files
	250MB	2180.0	160.0	6.5482
	500MB	4290.0	226.0	6.5486
	750MB	6334.0	368.0	6.5496
	1.0GB	8806.0	494.0	6.5497
	1.25GB	10919.0	613.0	6.5517
	1.50GB	12827.0	628.0	6.5516
	1.75GB	15140.0	793.0	6.5511
	2.00GB	17136.0	848.0	6.5509
			SQL files
	250MB	332.0	178.0	11.7680
	500MB	660.0	178.0	11.7621
	750MB	972.0	305.0	11.7592
	1.0GB	1357.0	450.0	11.7592
	1.25GB	1535.0	468.0	11.7623
	1.50GB	1844.0	651.0	11.7619
	1.75GB	2242.0	756.0	11.7608
	2.00GB	2623.0	895.0	11.7600
			XML files
	250MB	2149.0	545.0	250.9874
	500MB	4205.0	138.0	250.9998
	750MB	6298.0	155.0	250.9972
	1.0GB	8394.0	293.0	251.0136
	1.25GB	10404.0	301.0	251.1332
	1.50GB	12698.0	349.0	251.1085
	1.75GB	14886.0	326.0	251.0894
	2.00GB	16846.0	347.0	251.0799

TABLE II :

 II Calibration nodes characteristics

	Features	Nova node
	Server model	Dell PowerEdge R430
	CPU model	Intel Xeon E5-2620
	Number of CPU	2
	Cores per CPU	8
	Memory (GB)	64
	Storage (GB)	2 x 300 (HDD)

TABLE III :

 III Summary of simulation parameters

		SuperComputer	100 Gb/s
	Bandwidth	Datacenter	1 Gb/s
		Edge	10 Mb/s
		Regular	250MB to 2GB
	Size of exchanges	BigData	960MB to 3.2GB
		Edge	0.2MB
		Nova	100W
	Power Consumption (Idle state)	Little	4 W

TABLE IV :

 IV Simulation scenarios

	Scenario Data exchange	Bandwidth	Sender Receiver
	A	Regular	Datacenter	Nova	Nova
	B	Big Data	SuperComputer	Nova	Nova
	C	Big Data	Edge	Nova	Little
	D	S3D	Edge	Nova	Nova

TABLE V :

 V Simulated energy gains of non-lossy compression model over 4 different scenarios. Each value represents an exchange of 10,000 files.

		A		B		C		D
	Scenario	Energy(J)	%gain	Energy(J)	%gain	Energy(J)	%gain	Energy(J)	%gain
				Rand Binary Files
	No Compress.	1.12e+08	0.0	2.45e+08	0.0	6.50e+09	0.0	2.47e+07	0.0
	Active model	1.02e+08	8.760	2.10e+08	14.48	4.95e+09	23.79	1.91e+07	22.88
	Always ON	1.06e+08	4.836	2.10e+08	14.48	4.95e+09	23.79	1.91e+07	22.88
				BP files
	No Compress.	1.12e+08	0.0	2.45e+08	0.0	6.49e+09	0.0	3.00e+07	0.0
	Active model	9.07e+07	18.976	1.47e+08	39.87	1.10e+09	83.019	8.08e+05	97.30
	Always ON	1.17e+08	-4.860	1.47e+08	39.87	1.10e+09	83.018	8.08e+05	97.30
				SQL files
	No Compress.	1.12e+08	0.0	2.47e+08	0.0	6.50e+09	0.0	2.47e+07	0.0
	Active model	2.31e+07	79.437	3.65e+07	85.22	5.69e+08	91.25	2.28e+06	90.76
	Always ON	2.36e+07	79.002	3.65e+07	85.22	5.69e+08	91.25	2.28e+06	90.76
				XML files
	No Compress.	1.13e+08	0.0	2.46e+08	0.0	6.49e+09	0.0	2.47e+07	0.0
	Active model	7.73e+07	32.068	1.07e+08	56.41	1.33e+08	97.95	2.87e+06	88.40
	Always ON	9.88e+07	13.152	1.07e+08	56.41	1.33e+08	97.95	2.87e+06	88.40

ACKNOWLEDGMENTS

This research is partially supported by the French FSN ELCI project and the Inria-Rutgers SUSTAM (Sustainable Ultra Scale compuTing, dAta and energy Management) associate team, by NSF via grants numbers ACI 1339036, ACI 1441376 and by the Research Council of Norway (RCN) IKTPluss program, project number 270672. Experiments were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RE-NATER and several Universities as well as other organisations (https://www.grid5000.fr).


Emmanuel Paradis
email: emmanuel.paradis@ird.fr

Multidimensional scaling with very large data sets

Keywords: dimension reduction, distance data, projection method, random matrices

Multidimensional scaling has a wide range of applications when observations are not continuous but it is possible to define a distance (or dissimilarity) among them. However, standard implementations are limited when analyzing very large data sets because they rely on eigendecomposition of the full distance matrix and require very long computing times and large quantities of memory. Here, a new approach is developed based on projection of the observations in a space defined by a subset of the full data set. The method is easily implemented. A simulation study showed that its performance are satisfactory in different situations and can be run in a short time when the standard method takes a very long time or cannot be run because of memory requirements.

Introduction

Multidimensional scaling (MDS) aims to find a set of coordinates in one or several dimensions, so that the distances derived from these coordinates are the closest possible to the observed distances. By contrast to principal component analysis (PCA), MDS (also called principal coordinates analysis, PCoA) does not require a set of original coordinates. This is a very attractive feature since a distance can naturally be defined for many types of data whereas these data do not easily define a system of coordinates (e.g., DNA sequences, graphs, psychological profiles). A search of the string "multidimensional scaling" in the Web of Science returned 7313 hits (accessed 2017-07-31) with the highest proportions of publications in the fields of environmental sciences & ecology (15.9%), psychology (14.9%), and computer science (14.7%).

MDS is basically done through a decomposition of the symmetric distance matrix among observations. This matrix as thus has many rows and columns than there are observations, and its decomposition with standard eigenvalue analysis algorithm can be a limiting factor when this number is large. In a context where many scientific fields collect large amounts of data, this is clearly a severe limitation. Some examples include genomic analyses from high throughput sequencing technologies which typically handle millions of DNA sequences, or environmental data from high-resolution remote sensing which represent variables for millions of localities on the Earth's surface.

In this paper, I present an approach to avoid the limitations of the standard procedure of MDS. In the next section, I present the different computational procedures considered in this paper. Section 3 reports the results from a simulation study. The last section discusses these results and presents some perspectives.

Computational procedures

Section 2.1 below describes the standard MDS procedure; Section 2.2 explains how this procedure can be extended with random matrix algorithms. Sections 2.3 and 2.4 introduce new procedures.

Standard multidimensional scaling

Let us denote n the number of observations, and ∆ the n×n symmetric matrix of distances with δ ij (= δ ji ) being the distance between observations i and j. MDS proceeds by doing an eigendecomposition of the doubly-centred distance matrix:

- 1 2 J∆ 2 J, (1)
with J = I -1 n 11 T . The matrix of coordinates Z are calculated with:

Z = V Λ 1/2 ,
where V is an n × k matrix with the first k eigenvectors extracted from (1), and Λ is a diagonal k × k matrix with the first k eigenvalues (λ 1 , . . . , λ k ). Z has therefore n rows and k columns. The value of k is the number of dimensions of the projected space and is chosen by the investigator, usually depending on the values of λ.

This procedure involves manipulating matrices with n 2 elements which is very expensive in terms of memory requirements and computing times when n is large.

Eigendecomposition with random matrices

1-D projection

The principle of this method is quite simple. In a first step, m observations are selected (m < n) and a standard MDS is performed on them setting k = 1. We denote as z i the coordinate of the ith observation (i = 1, . . . , m). In a second step, for each observation j not among the m selected in step 1, the coordinate z j is found by minimizing:

f = m i=1 (δ ij -d ij ) 2 , (2) with d ij = |z i -z j |.
We may write the latter as:

d ij = γ i (z i -z j ) γ i =    +1 z j < z i -1 z j ≥ z i ,
A simple algebraic development leads to:

f = m i=1 δ 2 ij -2 m i=1 (δ ij γ i z i -δ ij γ i z j ) + m i=1 z 2 i -2z j m i=1 z i + mz 2 j .
We can now write the partial derivative of f with respect to z j :

∂f ∂z j = 2 m i=1 δ ij γ i -2 m i=1 z i + 2mz j .
Writing z = 1 m m i=1 z i and solving the last expression leads us to find the value z j that minimizes f :

1 m m i=1 δ ij γ i -z + z j = 0.
Given that there are only m distances δ ij to calculate, and that z is calculated only once, this equation can be solved numerically relatively easily. During the first step, m(m -1)/2 distances are computed for the classical MDS. Thus, the memory requirements of this method are very small. Overall, instead of calculating n(n -1)/2 distances for a full MDS, this method requires to calculate nm -m(m + 1)/2 distances.

The choice of the m observations is important for this projection procedure to work as best as possible. They should represent as much as possible the distribution of the n observations. This is not straightforward since the projection of these n points in the MDS space is not known a priori. In order to solve this issue, the following algorithm selects m observations so their distances are representative of the whole set of distances: 1. Select one observation i at random among the n ones; store i.

2.

Compute the distances δ ij with j among the values not yet stored.

3. Find j so that δ ij is max of the distances calculated at step 2; store j . 4. Compute the distances δ j j with j among the values not yet stored. 5. Select i so that δ j i is median of the distances calculated at step 4; store i. 6. Repeat steps 2-5 until m values are stored. This algorithm may be used whatever the number of dimensions k. There could be some variants of this algorithm: for instance, instead of the median in step 5, one could select a random distance δ j i so that the observations will be represented in proportion of their relative frequencies (somehow similar to an MCMC procedure).

Higher dimension projection

The method developed above can be generalized to more than one dimension in a straightforward way. The distances d ij in (2) would then be calculated as Euclidean distances in the projection space:

d ij = k l=1 (z il -z jl ) 2
In this case, the partial derivatives can be derived but are too complicated to be useful, so it is more efficient to rely on a standard minimization procedure to minimize (2). Two procedures were tested here: the classical BFGS method [START_REF] Broyden | The convergence of a class of double-rank minimization algorithms 1. General considerations[END_REF][START_REF] Fletcher | A new approach to variable metric algorithms[END_REF][START_REF] Goldfarb | A family of variable metric updates derived by variational means[END_REF][START_REF] Shanno | Conditioning of quasi-Newton methods for function minimization[END_REF]) and the PORT routine [START_REF] Gay | Usage summary for selected optimization routines[END_REF]).

Simulation study

A simulation study was conducted to answer two questions: What are the computing times of the different procedures described above? and, What are their respective accuracy? The data were simulated from a standard normal distribution with one or two dimensions and with different samples sizes n (1000, 5000, 10,000, and 50,000). The projection procedures were performed with m = 100 points chosen randomly. To assess the accuracy of the results, two quantities were calculated. First, the classical stress was calculated with Kruskal's formula [START_REF] Kruskal | Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis[END_REF]:

S = n i,j (δ ij -d ij ) 2 n i,j δ 2 ij .
This quantity varies between 0 (complete mismatch of the distances) and 1 (perfect match).

Second, a measure of the accuracy of the inferred distances was calculated with:

S = n i,j (δ ij -d ij ) 2 δ ij .
This second quantity is similar to Sammon's criterion [START_REF] Sammon | A nonlinear mapping for data structure analysis[END_REF]) and gives more emphasis on the precision of each distance whereas Kruskal's stress puts emphasis on the overall precision of the distances. In addition, in order to assess whether the observed results were better than simply randomly positioning the points in the MDS space, both quantities were calculated after randomizing the simulated data.

The simulations were run with k = 1 and k = 2. In the second case, the two variables were either independent or with a correlation of 0.7. Additionally, skewed distributions were generated with k = 1 in order to simulate aggregation of points. Two cases were considered: an exponential distribution with rate equal to one, and a mixture of two normal distributions with means -5 and 5 and proportions 0.9 and 0.1, respectively. In these two cases, the analyses were performed with m = 100 observations selected randomly like previously, and using the algorithm described above. All simulations were replicated 100 times and run on a computer equipped with a duo-core, 2.1 GHz processor and 16 GB of RAM, and running Ubuntu 16.04. All computations were implemented in R version 3.4.1 (R Core Team 2017); the code is available as supplementary material with this article.

1-D MDS

With a moderate sample size n = 1000, the three methods considered here completed in less than one second (Table 1). However, with n = 10,000, the standard MDS took almost dour minutes while the same procedure with a random decomposition took slightly more than eleven seconds, and the projection method took less than one second. Most importantly, the computing times of the latter appeared to be linearly related with n whereas the two others seemed to run proportionally to n 3 for the standard MDS, or to n 2 for the random Standard MDS 1.6 × 10 -15 2.7 × 10 -15 3.6 × 10 -15 -Random eigendecomposition 1.7 × 10 -2 7.6 × 10 -3 6.9 × 10 -3 -1-D Projection (m = 100) 9.1 × 10 -6 9.5 × 10 -6 9.3 × 10 -6 9.7 × 10 -6

"random" 9.5 × 10 -1 8.5 × 10 -1 8.5 × 10 -1 8.5 × 10 -1 decomposition. The computing time of the two decomposition-based methods was not 129 assessed with n = 50,000.

Unsurprisingly, the standard MDS performed the best considering either the stress S (Table 2) or the accuracy of the inferred distances S (Table 3). The projection method 132 performed the second best and better than the random decomposition method. All methods performed better than randomly projecting the points. with PORT 4.5 × 10 -10 4.9 × 10 -10 4.5 × 10 -10 4.7 × 10 -10 "random" 6.6 × 10 -1 6.5 × 10 -1 6.5 × 10 -1 6.5 × 10 -1

2-D MDS

The two methods based on matrix decomposition showed similar computing times than in one dimension (Table 4). On the other hand, the projection method was slower but its computing times still scaled proportionally to n. The PORT-based projection method was almost twice faster than the BFGS-based one.

With two independent variables, the accuracy of the projection methods were less than the standard MDS but these methods appeared still accurate (Tables 5 and6). The PORTbased variant was more accurate than the BFGS-based one. By constrast, the random decomposition method performed poorly. When the two variables were correlated, the performance of the standard MDS and the projection methods were similar to the previous situation; however, the performance of the random decomposition method deteriorated considerably and were only slightly better than randomly positioning the data (Tables 7 and8). 2.2 × 10 3 3.5 × 10 3 -2-D Projection (m = 100) with BFGS 5.1 × 10 -9 1.2 × 10 -7 5.0 × 10 -7 1.2 × 10 -5

with PORT 4.0 × 10 -13 1.1 × 10 -11 4.0 × 10 -11 1.1 × 10 -9

"random" 8.8 × 10 5 2.2 × 10 7 8.9 × 10 7 2.2 × 10 9

Table 7: Stress (S) from different 2-D MDS methods with correlated variables.

Method n 1000 5000 10,000 50,000

Standard MDS 1.4 × 10 -15 2.2 × 10 -15 2.7 × 10 -15 -Random eigendecomposition 5.5 × 10 -1 5.5 × 10 -1 5.5 × 10 -1 -2-D Projection (m = 100) with BFGS 5.4 × 10 -8 5.9 × 10 -8 5.5 × 10 -8 6.0 × 10 -8

with PORT 1.9 × 10 -10 2.1 × 10 -10 2.0 × 10 -10 2.1 × 10 -10 "random" 7.4 × 10 -1 7.4 × 10 -1 7.5 × 10 -1 7.4 × 10 -1

Table 8: Accuracy of inferred distances (S ) from different 2-D MDS methods with correlated variables.

Method n 1000 5000 10,000 50,000

Standard MDS 3.3 × 10 -24 2.7 × 10 -22 1.6 × 10 -21 -Random eigendecomposition 2.9 × 10 5 7.0 × 10 6 2.9 × 10 7 -2-D Projection (m = 100) with BFGS 8.9 × 10 -9 2.5 × 10 -7 7.7 × 10 -7 2.3 × 10 -5

with PORT 1.1 × 10 -14 2.8 × 10 -12 1.0 × 10 -11 2.8 × 10 -10

"random" 1.5 × 10 6 3.7 × 10 7 1.5 × 10 8 3.7 × 10 9

Skewed Distributions

For the two methods based on matrix decomposition, the results with the skewed distributions were very similar to the previous ones, and are thus not reported here. For the 1-D projection method, the analyses were performed with a uniform random sample and using the algorithm described above. Both algorithms resulted in similar computing times; however, the uniform sampling resulted in decreased accuracy while the above algorithm yielded performance comparable to the previous ones with non-aggregated data. This difference in accuracy was small for the data simulated from an exponential distribution (Table 9) whereas it was important for the mixture of normal variables where the skewness of the data was much more pronounced (Table 10).

Discussion

This article presents a method to perform MDS on very large data sets in short times and with small memory requirements. The objective of the present study was to develop a method easily implemented and generally applicable. One initial motivation was to avoid the need to perform a matrix decomposition of the full distance matrix as used in standard MDS.

Two approaches were considered in this study: the first one used matrix decomposition algorithms based on random matrices, and the second one used a projection algorithm based on a subset of points. The first approach did not appear as a viable solution to handle large data sets: it was too slow for sample sizes larger than 10,000 and was very inacurrate in some situations. This method performed poorly with correlated variables, which was an unexpected result. It is unclear whether this a pathological specific case or a more general problem with random decomposition of distance matrices. Further tests will be needed to clarify this point.

One issue not treated in depth in the present work is how to select the number of dimensions (k). In standard MDS, this value is selected depending on the eigenvalues extracted from the decomposition of the distance matrix. Typically, in practical applications of MDS two dimensions are selected in order to provide an interpretable graphical display. With the projection method proposed in this paper, since the number of dimensions determines the algorithm used, this number may be selected with respect to the eigenvalues of the standard MDS done on the subset of size m.

Another issue not explored here is the choice of the size of the subset (m). It was found that a value m = 100 is appropriate in the situations considered here: it makes possible to perform the projection easily since a larger value would make this procedure slower and more complicated. An interesting result was the good performance of the selection algorithm presented in this paper, particularly if the data were aggregated. This also deserves further study.

The present approach can have a wide range of practical applications. Many applied researchers need to analyze increasingly larger data sets, such as in ecological habitat modeling based on remote sensing [START_REF] Hansen | Humid tropical forest clearing from 2000 to 2005 quantified by using multitemporal and multiresolution remotely sensed data[END_REF] or in genomic analysis [START_REF] Erlich | A vision for ubiquitous sequencing[END_REF]. It will thus be interesting to see how the present method behaves and performs in practical applications.

Table 1 :

 1 Computing times of different 1-D MDS methods (in seconds).

	Method	n
		1000 5000 10,000 50,000
	Standard MDS	0.33 26.47 220.94	-
	Random eigendecomposition 0.09 1.09 11.66	-
	1-D Projection (m = 100)	0.11 0.41	0.79 3.91
	Table 2: Stress (S) from different 1-D MDS methods.
	Method	n
	1000	5000	10,000	50,000

Table 3 :

 3 Accuracy of inferred distances (S ) from different 1-D MDS methods.

	Method

Table 4 :

 4 Computing times of different 2-D MDS methods (in seconds).

	Method

2-D Projection (m = 100) with BFGS 5.4 × 10 -8 5.7 × 10 -8 5.7 × 10 -8 5.7 × 10 -8

Table 6 :

 6 Accuracy of inferred distances (S ) from different 2-D MDS methods.

	Method

Table 9 :

 9 Results of simulations with an exponential variable. CT: computing times (in seconds). Sampling m = 100 observations with the algorithm presented in this paper (A) or uniformly (U).

	n

Table 10 :

 10 Results of simulations with a mixture of two random variables. CT: computing times (in seconds). Sampling m = 100 observations with the algorithm presented in this paper (A) or uniformly (U).

	n

Acknowledgments

I am grateful to an anonymous reviewer for helpful comments on a previous version of this paper. This is publication ISEM 2017-201.

SUPPLEMENTARY MATERIAL

functions MDS.R: R functions used to perform the methods described in this paper.

sim.R: R functions used to perform the simulations reported in this paper.