Adrien Laurent

Gilles Vilmart

Multirevolution integrators for differential equations with fast stochastic oscillations

Keywords: highly-oscillatory stochastic differential equations, nonlinear Schrödinger equation, white noise dispersion, geometric integration, quadratic first integral AMS subject classification (2010): 60H35, 35Q55, 34E13

We introduce a new methodology based on the multirevolution idea for constructing integrators for stochastic differential equations in the situation where the fast oscillations themselves are driven by a Stratonovich noise. Applications include in particular highly-oscillatory Kubo oscillators and spatial discretizations of the nonlinear Schrödinger equation with fast white noise dispersion. We construct a method of weak order two with computational cost and accuracy both independent of the stiffness of the oscillations. A geometric modification that conserves exactly quadratic invariants is also presented.

Introduction

This article aims at developing invariant-preserving integrators of second weak order that are robust with respect to the stiffness ε both in accuracy and cost for the following class of highly-oscillatory d-dimensional SDEs driven by a one-dimensional Stratonovich noise dXptq 1c ε AXptq ¥ dW ptq F pXptqqdt, t ¡ 0, Xp0q X 0 , (

where W is a standard one-dimensional Wiener process, the function F : R d Ñ R d is a smooth non-linear map, the stiff parameter ε ¡ 0 is fixed and assumed small, and A R d¢d is a given matrix satisfying e A Id (equivalently A is diagonalizable and has all its eigenvalues in 2iπZ). In the deterministic setting, this last property yields that the solution xptq exppεAtqx 0 of dx dt ε ¡1 Ax is ε-periodic. For stochastic oscillations, it means that the solution Xptq exppε ¡1{2 AW ptqqX 0 of dX ε ¡1{2 AX ¥ dW satisfies XpT q Xp0q for a random time T inftt ¡ 0, § § ε ¡1{2 W ptq § § 1u of mean ε. The class of SDEs (1.1) includes in particular highly-oscillatory Kubo oscillators (see [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF])

dX 2π c ε ¢ 0 ¡1 1 0 X ¥ dW ¢ 0 ¡a a 0 Xdt, a R, (1.2)
or equivalently, dY 2iπε ¡1{2 Y ¥dW iaY dt in the complex setting where Y X 1 iX 2 .

Applying standard SDE integrators to solve equation (1.1) requires in general a time stepsize h ¤ ε to be accurate, which makes these methods dramatically expensive when ε is small. The goal of this paper is to create robust numerical methods, i.e., numerical integrators whose cost and accuracy do not deteriorate when ε becomes small. Several classes of methods have already been developed for highly-oscillatory SDEs with a deterministic fast oscillation (see for instance [START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF][START_REF] Vilmart | Weak second order multirevolution composition methods for highly oscillatory stochastic differential equations with additive or multiplicative noise[END_REF]), but not in the case where the stiff oscillatory part is applied to the noise itself. To numerically face this challenge, we introduce in this paper a new methodology to develop robust methods of any high weak order to approximate the solution of equation (1.1). In particular, we propose a method of weak order two, and a geometric modification of this algorithm that preserves quadratic invariants. Stochastic oscillations as defined in (1.1) typically arise in fiber optics models (see [START_REF]Nonlinear Fiber Optics[END_REF][START_REF]Applications of Nonlinear Fiber Optics[END_REF][START_REF] Garnier | Stabilization of dispersion-managed solitons in random optical fibers by strong dispersion management[END_REF]) with a spatial discretization of the highly-oscillatory nonlinear Schrödinger equation (NLS) with white noise dispersion duptq i c ε ∆uptq ¥ dW ptq F puptqqdt, upt 0q u 0 .

(1.3)

As described for instance in [START_REF] Garnier | Stabilization of dispersion-managed solitons in random optical fibers by strong dispersion management[END_REF], in the case ε 1, the NLS equation (1.3) with a cubic nonlinearity F puq |u| 2 u is a model in dimension d 1 describing the propagation of a signal in optical fibers where x corresponds to the retarded time, while t corresponds to the distance along the fiber. Taking into account the inevitable chromatic dispersion effects of the signal, modeled by a random centered stationary process m with a coefficient ν ¡ 0, yields the following random PDE, fv fx px, tq νimpxq f 2 v ft 2 px, tq ν 2 F pvpx, tqq, vpx 0, tq u 0 ptq. The perfect fiber would satisfy m 0, but in practice, engineers build fibers with a small varying dispersion coefficient. To limit the pulse broadening induced by random dispersion, specialists use a wide range of dispersion management techniques (see for instance [START_REF] Garnier | Stabilization of dispersion-managed solitons in random optical fibers by strong dispersion management[END_REF] and references therein). In [START_REF] Marty | On a splitting scheme for the nonlinear Schrödinger equation in a random medium[END_REF][START_REF] De Bouard | The nonlinear Schrödinger equation with white noise dispersion[END_REF], the authors show that if we denote u ν px, tq vpx{ν 2 , tq, then as ν tends to 0 and under some ergodicity assumptions on m, u ν converges to the solution u of equation (1.3) with ε 1. The non-stiff counterpart of equation (1.3), i.e., for ε 1, has also been studied theorically in [START_REF] Debussche | 1D quintic nonlinear Schrödinger equation with white noise dispersion[END_REF] for a particular nonlinearity. The highly-oscillatory behavior (ε 3 1) appears naturally when observing the propagation in long time with a small nonlinearity (via the change of variable t Ð εt) or the propagation of a small initial data in an optical fiber with a polynomial nonlinearity (via the change of variable u Ð u{ε).

A goal of this article is to develop efficient and cheap numerical methods that can model the propagation of pulses in this context, in order to observe some specific behaviors and, ultimately, to build enhanced fibers. Models of the form (1.3) also appear in the recent work [START_REF] Faou | Linearized wave turbulence convergence results for three-wave systems[END_REF] in the context of stochastic three-wave semi-linear systems. We emphasize that there is a growing interest in the recent litterature for stochastic models involving a fast Stratonovitch noise in the context of ergodic stochastic dynamics. In [START_REF] Abdulle | Accelerated convergence to equilibrium and reduced asymptotic variance for Langevin dynamics using Stratonovich perturbations[END_REF], it is shown for a class of overdamped Langevin equations that adding an appropriate fast Stratnovitch noise permits to increase the convergence rate to equilibrium, while reducing the asymptotic variance at infinity. This suggests that new efficient samplers for the invariant distribution of Langevin type models in context of large dimensional molecular dynamics models could be developed. We also mention the recent homogenization results on stochastic dynamics with fast Stratonovitch noises in [START_REF] Li | Homogenisation on homogeneous spaces[END_REF] where our periodicity assumption is replaced by an ergodicity assumption on the fast component of the dynamics posed on manifolds.

Numerous possibilities exist for numerically integrating equations (1.1) or (1.3). We highlight in particular the exponential integrators [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF][START_REF] Erdoğan | A new class of exponential integrators for sdes with multiplicative noise[END_REF] for the SDE (1.1), and the exponential integrators [START_REF] Cohen | Exponential integrators for nonlinear Schrödinger equations with white noise dispersion[END_REF], the Fourier split-step method [START_REF] Marty | On a splitting scheme for the nonlinear Schrödinger equation in a random medium[END_REF] or the Crank-Nicholson scheme [START_REF] Belaouar | Numerical analysis of the nonlinear Schrödinger equation with white noise dispersion[END_REF] for the SPDE (1.3). These methods have the advantage that they preserve the L 2 invariant of the equation (that is, }uptq} L 2 }u 0 } L 2 for all t ¥ 0) for a class of polynomial nonlinearities. However, they face a severe timestep restriction h ¤ ε when the stiff parameter ε is small. Even in the case of deterministic oscillations, there are restrictions in general, though some robust algorithms exist (see [START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF] for instance). The methods presented in this paper solve this issue of stepsize restriction. The idea is to approximate the solution of equation (1.1) at random times called revolution times because they correspond to complete revolutions of the oscillatory part dX ε ¡1{2 AX ¥ dW . This is in the spirit of [START_REF] Hofmann | Optimal approximation of stochastic differential equations by adaptive step-size control[END_REF] which also approximates the solution of SDEs at random times.

The article is organized as follows. Section 2 is devoted to the presentation of the new integrators. In section 3, we build an asymptotic expansion of the solution of (1.1) and evaluate it at revolution times to derive the new integrators and a limit model for equation (1.1). Section 4 is devoted to the weak convergence theorems and their proofs. In section 5, we present numerical experiments to confirm our theoretical error estimates, and we apply the new methods to solve numerically the Schrödinger equation (1.3).

Multirevolution integrators for stochastic oscillators

Initially created in [START_REF] Melendo | A new approach to the construction of multirevolution methods and their implementation[END_REF][START_REF] Calvo | Approximate compositions of a near identity map by multi-revolution Runge-Kutta methods[END_REF] in the context of celestial mechanics and later extended using geometric integration (see for instance [START_REF] Murua | Order conditions for numerical integrators obtained by composing simpler integrators[END_REF][START_REF] Calvo | On explicit multi-revolution Runge-Kutta schemes[END_REF][START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF]), multirevolution methods represent a class of numerical methods used for solving highly-oscillatory differential equations while reducing the cost of computation. In particular, they can approximate the solution of highly-oscillatory ODEs of the following form at stroboscopic times εN T , where T 1 is the period of dx dt Ax, and N is an integer,

dx dt 1 ε Ax F pxq, xp0q x 0 .
(2.1) The solution x of this equation at times εN T is a perturbation of identity, that is, x satifies xpεtq x 0 Opεtq, thus the solution loses its highly-oscillatory feature when evaluated at stroboscopic times, as shown in Figure 1 for the first component of the solution of equation (2.1) with F pxq ix (respectively F pyq ¢ 0 ¡1 1 0 y in the real setting). The idea of multirevolution is to approximate xpεN q with N Opε ¡1 q with a cost independent of ε.

For stochastic oscillations, the solution Xptq e ε ¡1{2 AW ptq X 0 of dX ε ¡1{2 AX ¥ dW is not periodic, but satisfies XpεT N q X 0 where the T N are random variables called revolution times and defined by

T 0 0, (2.2) T N 1 inf 3 t ¡ T N , ε ¡1{2 |Wpεtq ¡ W pεT N q| ¥ 1 A , N 0, 1, 2, . . .
If X is the solution of (1.1), we show in section 3.1 that X evaluated at times εT N is a perturbation of identity (in a strong and weak sense). Figure 2 illustrates the definition of revolution times and shows the perturbation of identity property on the first component of a Kubo oscillator (1.2) with a 1. We highlight that the revolutions times T N can be simulated without simulating the exact path W . Also we emphasize that the proposed algorithms do not require to simulate W thanks to the use of appropriate discrete random variables. This will be detailed in section 3.4.

We show in section 3.3 that the solution X of (1.1) evaluated at times εT tε ¡1 (when t{ε N is an integer) converges weakly when ε Ñ 0 to the solution y t of the deterministic ODE dy t dt xg 0 ypy t q, y 0 X 0 , where g 0 θ pyq e ¡Aθ F pe Aθ yq and xg 0 y : ³ 1 0 g 0 θ dθ. This ODE is exactly the same one as the asymptotic model for deterministic oscillators of the form (2.1). This asymptotic model naturally yields a weak order 1 deterministic integrator. We propose the two following new multirevolution methods of second weak order for integrating equation (

1 2 1 3N if p k 0 1 2π 2 k 2 N if p 0, k $ 0 ¡1 2π 2 p 2 N if p $ 0, k 0 1 2π 2 p 2 N if p k 0, p, k $ 0 0 else Erp α N p p α N k s 6 9 8
9 7

1 2 3N if p k 0 1 π 2 p 2 N if p k 0, p, k $ 0 0 else p r β N p,k 6 9 8 9 7 1 2π 2 k 2 N if p 0, k $ 0 ¡1 2π 2 p 2 N if p $ 0, k 0 0 else
The definition and construction of these random variables are further discussed in sections 3.2 and 3.4.

Method A (Explicit integrator of weak order two in H N ε to approximate the solution of equation (1.1) at times εT N m for m 0, 1, 2, . . .)

Y 0 X 0 for m ¥ 0 do Y m 1 Y m H Kt{2¡1 ķ¡Kt{2 c 0 k pY m qp α N k H 2 Kt{2¡1 p,k¡Kt{2 c 1 p pY m qpc 0 k pY m qq p β N p,k (2.4)

end for

Method B (Geometric integrator of weak order two in H N ε to approximate the solution of equation (1.1) at times εT N m for m 0, 1, 2, . . . while preserving quadratic invariants)

Y 0 X 0 for m ¥ 0 do Y m 1 Y m H Kt{2¡1 ķ¡Kt{2 c 0 k ¢ Y m Y m 1 2 p α N k (2.5) H 2 Kt{2¡1 p,k¡Kt{2 c 1 p ¢ Y m Y m 1 2 ¢ c 0 k ¢ Y m Y m 1 2 p r β N p,

Analysis and asymptotic expansion of the exact solution

In this section, we first obtain a local expansion of the solution of (1.1) and then evaluate it at particular random times to deal with the highly-oscillatory patterns of the exact solution.

Finally we derive from this expansion an asymptotic limit for equation (1.1) when ε Ñ 0.

Asymptotic expansion of the exact solution

Instead of studying directly equation (1.1), we apply the change of variable t Ð ε ¡1 t to obtain the following equation, whose solution satisfies Y ptq Xpεtq with X solution of (1.1),

dY ptq AY ptq ¥ d W ptq εF pY ptqqdt, Y p0q X 0 , ( 3.1)
where we denote for simplicity the Brownian motion W ptq ε ¡1{2 W pεtq again by W . We introduce the following assumption which guarantees in particular global existence and uniqueness of the solution.

Assumption 3.1. The function F is globally Lipschitz continuous and lies in C 3 P , i.e., there exists constants L, C, K ¡ 0 such that for all y, y 1 ,

y 2 R d |Fpy 1 q ¡ F py 2 q| ¤ L |y 1 ¡ y 2 | § § §F piq pyq § § § ¤ Cp1 |y| K q, i t0, 1, 2, 3u. (3.2)
Also the initial condition X 0 has bounded moments, that is, Er|X 0 | p s V for p ¥ 0.

Therefore we denote ϕ ε,t pX 0 q Y ptq the solution of equation (3.1) and focus in the rest of the paper on the approximation of ϕ ε,t pyq at times t Opε ¡1 q. The variation of constants formula yields ϕ ε,t pyq e AW ptq y ε » t 0 e ApW ptq¡Wpsqq F pϕ ε,s pyqqds.

(3.3)
We deduce the following regularity properties. Lemma 3.2. Under Assumption 3.1, the following estimates hold for all y, y 1 , y 2 R d , where C and K are independent of ε and t,

1. |ϕ ε,t py 1 q ¡ ϕ ε,t py 2 q| ¤ C |y 1 ¡ y 2 | e Cεt , 2. |ϕ ε,t pyq| ¤ Cp1 |y|qe Cεt , 3. ϕ ε,t pyq is C 3 in y and § § §ϕ piq ε,t pyq § § § ¤ Cpεtq i¡1 p1 |y| K qe Cεt for i 1, 2, 3.
The proof is postponed to the appendix. It mainly relies on the Gronwall theorem and the boundedness of the one-periodic function θ Þ Ñ e θA . Using a local expansion of the solution of (3.1) in ε, we define the following first and second order approximations of

ϕ ε,t pyq, ψ 1 ε,
§ § §ϕ ε,t pyq ¡ ψ j ε,t pyq § § § ¤ Cp1 |y| K qe Cεt pεtq j 1 .
The functions ψ j ε,t satisfy the following straightforward inequalities proved with similar arguments as for Lemma 3.2.

Lemma 3.4. With the assumptions and notations of Proposition 3.3, the following esti-

mates hold for all y R d , where C and K are independent of ε and t

, § § ψ 1 ε,t pyq § § ¤ Cp1 |y|qe Cεt , (3.5) § § ψ 2 ε,t pyq § § ¤ Cp1 |y| K qe Cεt , (3.6) § § §ψ 2 ε,t pyq ¡ e AW ptq y § § § ¤ Cp1 |y| K qpεtqe Cεt . (3.7) Proof of Proposition 3.3. Using Assumption 3.1, we get § § ϕ ε,t pyq ¡ ψ 1 ε,t pyq § § ¤ Lε » t 0 § § §ϕ ε,s pyq ¡ e AW psq y § § § ds. Then Lemma 3.2 yields § § §ϕ ε,s pyq ¡ e AW psq y § § § ¤ Cε » s 0 |Fpϕ ε,r pyqq| dr ¤ Cε » s 0 p1 |ϕ ε,r pyq|qdr ¤ Cε » s 0 p1 Cp1 |y|qe Cεr qdr ¤ Cp1 |y|qe Cεs pεsq. We deduce § § ϕ ε,t pyq ¡ ψ 1 ε,t pyq § § ¤ Cp1 |y|qe Cεt pεtq 2 . For j 2, we first denote r ψ 2 ε,t pyq e AW ptq y εe AW ptq » t 0 e ¡AWpsq F pψ 1
ε,s pyqqds. With the same arguments we used for j 1 and inequality (3.5), we have § § §ϕ ε,t pyq ¡ r

ψ 2 ε,t pyq § § § ¤ Cp1 |y| K qe Cεt pεtq 3 . It is sufficient to prove that § § ψ 2 ε,t pyq ¡ r ψ 2 ε,t pyq § § ¤ Cp1 |y| K qe Cεt pεtq 3 . A Taylor expansion of F pψ 1 ε,s pyqq in ε gives F pψ 1 ε,s pyqq F pe AW psq yq εF I pe AW psq yq ¢ e AW psq » s 0 e ¡AWprq F pe AW prq yqdr R ε,s . The remainder R ε,s satisfies |R ε,s | ¤ Cε 2 sup xre AW psq y,ψ 1 ε,s pyqs F P pxq § § § § e AW psq » s 0 e ¡AWprq F pe AW prq yqdr § § § § 2 .
Then, using the polynomial growth of F P and inequality (3.5), we get

|R ε,s | ¤ Cp1 § § §e AW psq y § § § K § § ψ 1 ε,s pyq § § K qpεsq 2 e Cεs ¤ Cp1 |y| K qpεsq 2 e Cεs .
Hence the result. ¥

We shall prove in section 3.3 that the function ψ 2 ε,t in (3.4) evaluated at the revolution times T N (defined in (2.2)) yields a strong order 2 approximation in H εN .

W τ ¢ 1 ¡ t τ i W i ¢ t τ ¡ i W i 1 for iτ ¤ t ¤ pi 1qτ, (3.8)
where W 0 0 and W i 1 W i c τ ξ i with pξ i q i a family of independent standard Gaussian random variables, then it can be shown that we obtain an integrator of strong order two in εt. However, the cost of a standard method computing an approximation of the integrals of equation (3.4) by replacing W with W τ is in Opt 2 {τ 2 q, which makes this method tremendously expensive for t Opε ¡1 q. This is why we develop in section 3.4 weak integrators based on a weak approximation of equation (3.4) with a cost independent of t. We shall replace stochastic integrals by appropriate discrete random variables in order not to simulate any expensive Brownian path W .

Properties of the revolution times

In this section, we study some properties linked to the revolution times T N that will be useful for the analysis. Proposition 3.6. The revolution times T N defined in (2.2) are positive and finite almost surely. Their differences pT N 1 ¡ T N q N ¥0 are independent identically distributed random variables (same law as T 1 ). The Laplace transform Ere zT 1 s of T 1 exists and is analytic for (3.9)

In particular, ErT 1 s 1, ErT 2 1 s 5 3 , and VarpT 1 q 2 3 . Finally, for a fixed ε 0 $ 0, π 2

16

, for all ε s0, ε 0 s and p ¥ 0, we have the estimate

Ere εT N pεT N q p s ¤ Ce CεN pεNq p . (3.10)
The law of the first revolution time T 1 has an analytic density, but there is no closed formula for it. It can be numerically approximated accurately by inverting the Laplace transform. In Figure 3, we observe the convergence in law of T N to a Gaussian variable according to the central limit theorem.

Proof of Proposition 3.6. The first properties can be deduced from [24, Chap.

Ere εT N pεT N q p s ¤ Ere 2εT N s 1{2 ErpεT N q 2p s 1{2 Ere 2εT 1 s N {2 ε p ErT 2p N s 1{2 ¤ Ere 2ε 0 T 1 s εN {2ε 0 pεNq p ErT 2p 1 s 1{2 ¤ Ce CεN pεNq p ,
where we used first the Cauchy-Schwarz inequality and then twice the Jensen inequality. ¥

For developing an algorithm for the weak error, it is useful to know the moments of the random variables appearing in the discretization, that are costly to simulate numerically, in order to replace them with cheap discrete approximations with the same first and second moments. This is the goal of the following proposition.

Proposition 3.7. The random variables

α N k 1 N ³ T N 0 e 2iπkW psq ds, β N p,k 1 N 2 ³ T N 0 e 2iπpW psq ³ s 0 e 2iπkW prq drds, r β N p,k 1 N 2 ³ T N 0 e 2iπpW psq ¡ ³ s 0 e 2iπkW prq ¡ 1 2 ³ T N 0 e 2iπkW prq dr © ds β N p,k ¡ α N p α N k 2 satisfy Erα N k s δ k 4 1 if k 0 0 else
Erβ N p,k s 6 9 9 9 9 9 9 8 9 9 9 9 9 9 7

1 2 1 3N if p k 0 1 2π 2 k 2 N if p 0, k $ 0 ¡1 2π 2 p 2 N if p $ 0, k 0 1 2π 2 p 2 N if p k 0, p, k $ 0 0 else Erα N p α N k s 6 9 8
9 7 ³ T N 0 e 2iπkW psq ³ T N s e 2iπpW prq drds, so that we get if p $ 0,

1 2 3N if p k 0 1 π 2 p 2 N if p k 0, p, k $ 0 0 else Er r β N p,k s 6 9 8 9 7 1 2π 2 k 2 N if p 0, k $ 0 ¡1 2π 2 p 2 N if p $ 0, k 0 0 else Proof. Let k $ 0 (
Erβ N p,k s δ p k ¡ δ k 2π 2 p 2 N .
The case p 0 is obtained by integrating by parts and using the same arguments. Indeed, we find

β N 0,k T N N α N k ¡ β N k,0 i πkN 2 » T N 0 T N e 2iπkW psq dW psq ¡ β N k,0 ,
and Erβ N 0,k s ¡Erβ N k,0 s. Finally, the moments Erα N p α N k s are computed via the equality

β N p,k β N k,p α N p α N k .
Then we obtain Er r

β N p,k s from the formula r β N p,k β N p,k ¡ α N p α N k 2
. ¥ Remark 3.8 (stochastic Fourier series). Let f be a L 2 function on s0, 1r extended on R by 1-periodicity, whose Fourier coefficients are denoted as pc k q kZ . Then we deduce from Proposition 3.7 the following equalities, where the second one is the stochastic version of the Bessel-Parseval theorem,

E » T 1 0 f pWpsqqds & c 0 » 1 0 f pθqdθ, E » T 1 0 |fpWpsqq| 2 ds & ķZ |c k | 2 , E § § § § » T 1 0 f pWpsqqds § § § § 2 ' 5 |c 0 | 2 3 ķ Z ¦ |c k | 2 π 2 k 2 .

Asymptotic expansion at revolution times and limit model

With the results of subsection 3. Notice that c 0 k pyq C d and c 1 p pyq pc 0 k q I pyq C d¢d but ψ ε,t pyq R d . We now evaluate this function ψ ε,t pyq at time t T N to get a second order strong approximation. Proposition 3.9. We define the following quantity

ψ ε,N pyq y H ķZ c 0 k pyqα N k H 2 p ,kZ c 1 p pyqpc 0 k pyqqβ N p,k ,
where pc 0 k pyqq kZ and pc

|ϕ ε,T N pyq ¡ ψ ε,N pyq| 2 % 1{2 ¤ Cp1 |y| K qH 3 , (3.12) |Erφpϕ ε,T N pyqqs ¡ Erφpψ ε,N pyqqs| ¤ Cp1 |y| K qH 3 , (3.13)
that is, ψ ε,N pyq is a numerical approximation of ϕ ε,T N pyq of strong/weak local order two.

Proof. Inequality (3.12) is a straightforward consequence of Proposition 3.6 when evaluating the estimates of Proposition 3.3 at time T N . For the weak local estimate (3.13), using inequality (3.12), the mean value inequality, Lemma 3.2 and equations (3.5) and (3.6), we get

E § § §φpϕ ε,T N pyqq ¡ φpψ j ε,T N pyqq § § § % ¤ E ! sup xrϕ ε,T N pyq,ψ j ε,T N pyqs § § φ I pxq § § § § §ϕ ε,T N pyq ¡ ψ j ε,T N pyq § § § ( ) ¤ Cp1 |y| K qE ! e CεT N pεT N q j 1 sup xrϕ ε,T N pyq,ψ j ε,T N pyqs p1 |x| p q ( ) ¤ Cp1 |y| K qE pεT N q j 1 e CεT N ¡ 1 |ϕ ε,T N pyq| p § § §ψ j ε,T N pyq § § § p ©% ¤ Cp1 |y| K qE pεT N q j 1 e CεT N $
.

Finally we obtain inequality (3.13) by taking H small enough so that we can apply Proposition 3.6. ¥

For a fixed T N ε, when ε Ñ 0 (or equivalently N Ñ V), the solution of (3.1) evaluated at stroboscopic times T N T T ε ¡1 converges weakly to the solution of a deterministic ODE that involves only the first mode c 0 0 xg 0 y ³ 1 0 g 0 θ dθ of g 0 . This asymptotic model is the same one as for the deterministic equation (2.1). The proof is postponed to subsection 4.3. Proposition 3.10 (asymptotic model). Under Assumption 3.1, for T ¡ 0, the solution ϕ ε,T T ε ¡1 pX 0 q (for ε such that T ε ¡1 is an integer) of equation (1.1) converges weakly when ε Ñ 0 to the solution at time T of dy t dt xg 0 ypy t q, y 0 X 0 , (3.14) that is, for all test function φ C 3 P , lim εÑ0 § § Erφpϕ ε,T T ε ¡1 pX 0 qqs ¡ Erφpy T qs § § 0.

Remark 3.11. It can be proven using the results of section 4 that the solution of the asymptotic model (Proposition 3.10) is an order one weak approximation of

XpεT N m q for m ¥ 0 and X solution of equation (1.1). We deduce the following simple one-step explicit deterministic integrator that corresponds to the Euler method applied to equation (3.14), y 0 X 0 , y m 1 y m Hc 0 0 py m q.

(3.15)

Its cost is independent of ε and N , and it has weak order one w.r.t. H, that is, for all m ¥ 0, Erφpϕ ε,T N m pX 0 qqs ¡ Erφpy m qs OpHq.

Construction of the second order integrators

To obtain an integrator of weak order two with a cost independent of ε and N , we truncate the local expansion of Proposition 3.9. We also replace the involved random variables with cheap discrete random variables with the same first and second moments. To simulate the random variable α N k with discrete random variables p α N k with the same first and second moments, we introduce a set pξ k q kN of independent random variables such that Ppξ k ¨1q 1 2 , the covariance matrix pC N α q p,k such that pC N α q 2p¡1:2p,2k¡1:2k ¡k for k 0 (so that the solution stays real while still having the good moments). We also define p Remark 3.12. The methodology presented in section 3 can be generalized to any order. Thus, under more regularity assumptions on F , it is possible to build algorithms similar to Method A of any weak order and that are still robust with respect to the stiffness ε. For order 3, Method A becomes l,p,k . It is also possible to generalize Method B up to any order in the spirit of the middle-point scheme, but the construction of discrete random variables allowing the preservation of quadratic invariants is not obvious for higher orders (although backward error analysis guarantees the preservation of quadratic invariants for the exact random variables based on W ).

£ CovpRepα N p q, Repα N k qq CovpRepα N p q, Impα N k qq CovpImpα N p q, Repα N k qq CovpImpα N p q, Impα N k qq
Y m 1 Y m H ķZ c 0 k pY m qp α N k H 2 p ,kZ c 1 p pY m qpc 0 k pY m qq p β N p,k H 3 ļ,p,kZ c 1 l pc 1 p pc 0 k qqpY m qp γ p1q,N l,p,k

Weak convergence analysis

This section focuses on the proofs of the following two theorems, showing the order two convergence of Methods A and B. satisfies QpY m 1 q QpY m q.

These two theorems focus on approximating the exact solution of equation ( 1 Using the truncation estimates that we previously discussed and the Lipschitz property of g 0 θ , one gets § § § r

ϕ ε,t ¡ r ϕ pKtq ε,t § § § pyq ¤ ε » t 0 § § §g 0 W psq pr ϕ ε,s pyqq ¡ g 0 W psq pr ϕ pKtq ε,s pyqq § § § ds Cεte ¡cKt sup r0,1s § § g 0 pyq § § ¤ Cε » t 0 § § § r ϕ ε,s pyq ¡ r ϕ pKtq ε,s pyq § § § ds Cεtp1 |y|qe ¡cKt .
The Gronwall lemma and Proposition 3.6 yield for mN ε ¤ T ,

E § § § r ϕ ε,T N m pyq ¡ r ϕ pKtq ε,T N m pyq § § § 2 & 1{2 ¤ Cp1 |y|qe ¡cKt .
The structure of the convergence proof is similar to the one for standard SDE integrators, see e.g. [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Chap. 2], but one has to be cautious because our solution is evaluated at stochastic times and the error constants should not depend on ε or N .

Boundedness of the numerical moments

p ψ ε,N pyq y H ķZ c 0 k pyqp α N k H 2 p ,kZ c 1 p pyqpc 0 k pyqq p β N p,k , (4.2)
where p

α N k , p β N
p,k are random variables defined such that for all q ¡ 0,

E ¢ °k |p α N k | 2 k 2 q & and E £ °p,k § § § p β N p,k § § § 2 k 2 q '
are bounded uniformly in N . Then, under Assumption 3.1 and if |y| has bounded moments, for any T ¡ 0, for all m, H such that mεN mH ¤ T , for all q ¡ 0, we have Er § § p ψ m ε,N pyq § § 2q s ¤ C q p1 Er|y| 2q sq, where C q is independent of m, ε and N .

Proof. We first prove § § § p

ψ ε,N pyq ¡ y § § § ¤ CHp1 |y|qM N , (4.3)
where ErpM N q 2q s ¤ C q for all q ¡ 0. We have § § § p

ψ ε,N pyq ¡ y § § § H § § § § § ķZ c 0 k pyqp α N k H p ,kZ c 1 p pyqpc 0 k pyqq p β N p,k § § § § § ¤ CH ¤ ¥ M p0q N d ķZ k 2 § § c 0 k pyq § § 2 M p1q N d pZ § § c 1 p pyq § § 2 d ķZ k 2 § § c 0 k pyq § § 2
where M p0q

N °k |p α N k | 2 k 2 and M plq N °p,k § § § p β N p,k § § § 2 k 2
have moments bounded uniformly in N . Then using the Bessel-Parseval theorem, we get °k

k 2 § § c 0 k pyq § § 2 ³ 1 0 § § f θ g θ pyq 0 § § 2 dθ. Assumption 3.1 yields § § f θ g 0 θ pyq § § ¤ Cp1 |y|q. Then, the Bessel-Parseval theorem applied on g 1 θ gives °p § § c 1 p pyq § § 2 ¤ C, hence the result. We define ∆ψ m p ψ m 1 ε,N pyq ¡ p ψ m ε,N pyq p p ψ ε,N ¡ Idqp p ψ m ε,N pyqq, then p p ψ m 1 ε,N pyqq 2q p p ψ m ε,N pyqq 2q 2q j1 ¢ 2l j p p ψ m ε,N pyqq 2q¡j ∆ψ j m .
Equation (4.3) and the bounded moments of

M N give § § §E p p ψ m ε,N pyqq 2q¡j ∆ψ j m % § § § ¤ E § § § p ψ m ε,N pyq § § § 2q¡j CH j p1 § § § p ψ m ε,N pyq § § § j qM j N & ¤ C q H ¢ 1 E § § § p ψ m ε,N pyq § § § 2q & . We deduce 1 E § § § p ψ m 1 ε,N pyq § § § 2q & ¤ e CqH ¢ 1 E § § § p ψ m ε,N pyq § § § 2q & and by induction E § § § p ψ m ε,N pyq § § § 2q &
¤ e CqmH p1 Er|y| 2q sq ¤ e CqT p1 Er|y| 2q sq.

¥ Proposition 4.5 (bounded moments for the integrator (2.5)). Assume that for y R d , the numerical scheme p

ψ ε,N pyq satisfies p ψ ε,N pyq y H ķZ c 0 k £ y p ψ ε,N pyq 2 p α N k (4.4) H 2 p ,kZ c 1 p £ y p ψ ε,N pyq 2 £ c 0 k £ y p ψ ε,N pyq 2 p r β N p,k , where p α N k , p r β N
p,k are random variables defined such that for all q ¡ 0,

°k § § p α N k § § , °p,k § § § § p r β N p,k § § § § , E ¢ °k |p α N k | 2 k 2 q & and E ¢ °p,k § § § § p r β N p,k § § § § 2 k 2 q &
are bounded uniformly in N . Then, under Assumption 3.1 and if |y| has bounded moments, for H 0 small enough and any T ¡ 0, for all m, H such that mεN mH ¤ T and H ¤ H 0 , for all q ¡ 0, we have Er § § § p

ψ m ε,N pyq § § § 2q s ¤
C q p1 Er|y| 2q sq, where C q is independent of m, ε and N .

Proof. We prove an equivalent of the estimate (4.3) for p ψ ε,N pyq. The growth properties of the Fourier coefficients yield § § § p

ψ ε,N pyq ¡ y § § § ¤ CH £ ķ § § c 0 k § § pyq § § p α N k § § p,k § § c 0 k § § pyq § § § § p r β N p,k § § § § CH £ ķ § § p α N k § § p,k § § § § p r β N p,k § § § § § § § p ψ ε,N pyq ¡ y § § § , As °k § § p α N k § § °p,k § § § § p r β N p,k § § §
§ is bounded, using the same estimates as in the proof of Proposition 4.4, we get for all H ¤ H 0 small enough, § § § p

ψ ε,N pyq ¡ y § § § ¤ CHp1 |y|qM N , (4.5)
where M N has bounded moments. The remaining of the proof is the same as in the proof of Proposition 4.4. ¥

Local weak error

Proposition 4.6 (local error estimate). Assume that for y R d deterministic, the numerical scheme can be written as

p ψ ε,N pyq y H ķZ c 0 k pyqp α N k H 2 p ,kZ c 1 p pyqpc 0 k pyqq p β N p,k R,
where Er|R|s ¤ Cp1 |y| K qH 3 and p

α N k C, p β N p,k R are random variables such that p α N k p α N ¡k and Erp α N k s Erα N k s, Er p β N p,k s Erβ N p,k s, Erp α N k 1 p α N k 2 s Erα N k 1 α N k 2 s.
Under Assumption 3.1, if p ψ ε,N pyq satisfies the assumptions of Proposition 4.4 (or Proposition 4.5), for all test function φ C 3 P , there exists H 0 ¡ 0 such that for all H N ε ¤ H 0 , the following estimate holds, where C and K are independent of ε and

N , § § §Erφpϕ ε,T N pyqqs ¡ Erφp p ψ ε,N pyqqs § § § ¤ Cp1 |y| K qH 3 ,
that is, the numerical scheme has weak local order two.

Proof. Using Proposition 3.9 and its notation

ψ ε,N pyq, it is enough to prove that § § §Erφpψ ε,N pyqqs ¡ Erφp p ψ ε,N pyqqs § § § ¤ Cp1 |y| K qH 3 .
A local expansion gives

φpψ ε,N pyqq φpyq φ I pyqpψ ε,N pyq ¡ yq φ P pyqpψ ε,N pyq ¡ y, ψ ε,N pyq ¡ yq R 1 . As ψ ε,N pyq ψ 2 ε,T N pyq (see equation (3. 4
)), using inequalities (3.6), (3.7) evaluated at T N and Proposition 3.6 yields

Er|R 1 |s ¤ E sup xry,ψ ε,N pyqs § § §φ p3q pxq § § § |ψ ε,N pyq ¡ y| 3 ' ¤ E Cp1 |y| K |ψ ε,N pyq| K qp1 |y| K qpεT N q 3 e CεT N % ¤ E Cp1 |y| K qpεT N q 3 e CεT N % ¤ Cp1 |y| K qH 3 .
We obtain a similar expansion for φp p ψ ε,N pyqq:

φp p ψ ε,N pyqq φpyq φ I pyqp p ψ ε,N pyq ¡ yq φ P pyqp p ψ ε,N pyq ¡ y, p ψ ε,N pyq ¡ yq x R 1 ,
where, using inequality (4.3) (or (4.5)),

Er| x R 1 |s ¤ E sup xry, p ψ ε,N pyqs § § §φ p3q pxq § § § § § § p ψ ε,N pyq ¡ y § § § 3 ' ¤ CE p1 |y| K § § § p ψ ε,N pyq § § § K qp1 |y| K qH 3 M 3 N & ¤ Cp1 |y| K qH 3 E p1 M K N q $ ¤ Cp1 |y| K qH 3 .
Making the difference of both equations gives

φpψ ε,N pyqq ¡ φp p ψ ε,N pyqq φ I pyqpψ ε,N pyq ¡ p ψ ε,N pyqq ¡ φ P pyqpψ ε,N pyq ¡ yq 2 (4.6)
φ P pyqp p ψ ε,N pyq ¡ yq 2 R, where Er|R|s ¤ Cp1 |y| K qH 3 . For the first term of (4.6), we have

Erφ I pyqpψ ε,N pyq ¡ p ψ ε,N pyqqs H ķZ Erφ I pyqpc 0 k pyqpα N k ¡ p α N k qqs H 2 p ,kZ Erφ I pyqpc 1 p pyqpc 0 k pyqqpβ N p,k ¡ p β N p,k qqs.
Then, we get

Erφ I pyqpc 0 k pyqpα N k ¡ p α N k qqs Erα N k ¡ p α N k sφ I pyqpc 0 k pyqq 0.
We can do the same thing with the term in β N p,k and obtain Erφ I pyqpψ ε,N pyq ¡ p ψ ε,N pyqqs 0. Let us now study the second order term Z φ P pyqp p ψ ε,N pyq¡yq 2 ¡φ P pyqpψ ε,N pyq¡yq 2 that appears in (4.6). We develop this expression and keep only the order one and two terms to obtain Z H 2 Y R where Er|R|s ¤ Cp1 |y| K qH 3 (by the same arguments as before)

and

Y ķ 1 ,k 2 φ P pyqpc 0 k 1 pyqp α N k 1 , c 0 k 2 pyqp α N k 2 q ¡ φ P pyqpc 0 k 1 pyqα N k 1 , c 0 k 2 pyqα N k 2 q $ ķ 1 ,k 2 pp α N k 1 p α N k 2 ¡ α N k 1 α N k 2 qφ P pyqpc 0 k 1 pyq, c 0 k 2 pyqq
The condition on the moments of the p α N k yields ErY s 0.

Putting all these arguments together in (4.6), we finally get that § § §Erφpψ ε,N pyqqs ¡ Erφp p ψ ε,N pyqqs § § § ¤ Cp1 |y| K qH 3 .

We deduce the local order two of the proposed numerical scheme. ¥

Remark.

The constant H 0 in Proposition 4.6 depends on F , but also on the polynomial growth power of φ and its first three derivatives. This dependence is expected when trying to evaluate the solution of SDEs at random times. To make H 0 independent of the test functions, one can consider the following sets of test functions

C 3 P,K tφ C 3 , hC ¡ 0, hk ¤ K, dy, § § §φ piq pyq § § § ¤ Cp1 |y| k q, i t0, 1, 2, 3uu.

Global error

¡ 0, E ¢ °k |p α N k | 2 k 2 q & and E £ °p,k § § § p β N p,k § § § 2 k 2 q ' are bounded uniformly in N (respectively °k § § p α N k § § , °p,k § § § § p r β N p,k § § § § , E ¢ °k |p α N k | 2 k 2 q & and E ¢ °p,k § § § § p r β N p,k § § § § 2 k 2 q &
are bounded uniformly in N ), for all T ¡ 0, for all test function φ C 3 P , there exists H 0 ¡ 0 such that for all H ¤ H 0 , for all M ¥ 0 such that M N ε M H ¤ T , there exist two positive constants K and C both independent of ε and N such that § § §Erφpϕ ε,T N M pX 0 qqs ¡ Erφp p

ψ M ε,N pX 0 qqs § § § ¤ CH 2 p1 Er|X 0 | K sq.
Proof. We denote

¤ i ¤ 3, § § § r φ piq m pyq § § § ¤ Ce CεT N m p1 |y| K q.
Thus, knowing the hitting times involved, r

| ¤ M m1 CE e CεT N m $ H 3 ¢ 1 E § § § p ψ M ¡m ε,N pX 0 q § § § K & .
Finally, the moments of p ψ m ε,N pX 0 q are all bounded uniformly in ε, N and m according to Proposition 4.4 (respectively 4.5). Thus The weak order two is obtained using Theorem 4.7. Indeed the use of discrete random variables and Proposition 3.7 gives the convergence of the involved series.

|e M | ¤ M m1 CH 3 p1 Er|X 0 | K sq ¤ CH 2 p1 Er|X 0 | K sq.
For showing that Method B preserves quadratic invariants, it is sufficient to prove that

Q I pyqp °k c 0 k pyqp α N k q 0 and Q I pyqp °p,k c 1 p pyqpc 0 k pyqq p r β N
p,k q 0 (see [START_REF] Hairer | Geometric numerical integration[END_REF]Chap. IV]). The preservation of Q by equation (1.1) yields Q I pyqpAyq 0 and Q I pyqpFpyqq 0. We deduce the following two equations, valid for all y R d ,

y T Sg 0 θ pyq 0, (4.7) y T Sg 1 θ pyqpg 0 ν pyqq ¡pg 0 ν pyqq T Sg 0 θ pyq, (4.8)
where equation (4.8) is obtained by differentiating equation (4.7) in the direction g 0 ν . Using equation (4.7), we have

Q I pyq £ ķ c 0 k pyqp α N k » 1 0 Q I pyqpg 0 θ pyqq ķ e ¡2iπkθ p α N k dθ 0.
For the second order term, equation (4.8) and the values of Proposition 3.7 yield

Q I pyq £ p,k c 1 p pyqpc 0 k pyqq p r β N p,k » 1 0 » 1 0 y T Sg 1 θ pyqpg 0 ν pyqq p,k e ¡2iπpθ e ¡2iπkν p r β N p,k dνdθ ¡ » 1 0 » 1 0 pg 0 ν pyqq T Sg 0 θ pyq p,k e ¡2iπpθ e ¡2iπkν p r β N p,k dνdθ ¡ 1 2 » 1 0 » 1 0 pg 0 ν pyqq T Sg 0 θ pyq p,k
re ¡2iπpθ e ¡2iπkν e ¡2iπpν e ¡2iπkθ s p r

β N p,k dνdθ 0.
Hence Method B is well-posed, has weak order 2 and preserves the invariant Q. ¥

Numerical experiments

In this section, we first illustrate numerically the weak order two of Methods A and B with convergence curves. Then, we apply the new algorithms to solve the nonlinear Schrödinger equation with highly-oscillatory white noise dispersion (1.3).

Weak order of convergence

To confirm the results of Theorem 4.1 and Theorem 4.2, we check numerically if Methods A and B have weak order two of accuracy w.r.t. H uniformly in ε and N . As the Euler-Maruyama method and the algorithms presented in [START_REF] Cohen | On the numerical discretisation of stochastic oscillators[END_REF][START_REF] Belaouar | Numerical analysis of the nonlinear Schrödinger equation with white noise dispersion[END_REF][START_REF] Cohen | Exponential integrators for nonlinear Schrödinger equations with white noise dispersion[END_REF] are completely inaccurate if they do not satisfy the severe timestep restriction h 3 ε, we compare the performance of Methods A and B to the performance of the Euler method (3.15). We first apply the algorithms on equation (1.1) with the linearity F pyq iy, A 2iπ, X 0 1 and ε 10 ¡3 .

Equivalently we can write it in the real setting as

dX 2π c ε ¢ 0 ¡1 1 0 X ¥ dW ¢ 0 ¡1 1 0 Xdt, X 0 ¢ 1 0 .
We plot on a logarithmic scale an estimate of the weak error for approximating X at time T 10 ¡3 T 2 8 where ErT s 0.256. The exact solution XpT q is approximated by the output of Method B for H ε. The parameters N and m are varying under the condition that N m 2 8 . The test function is φpyq 2y 1 4y 2 and the average is taken over 10 7

trajectories. We choose the tolerance 10 ¡13 for the fixed point. On the right picture of We take 8 modes for the Fourier decomposition and the same other parameters as before.

The average is taken over 10 6 trajectories. In both cases, we observe the weak order two of Methods A and B. The irregularities of the curve for a small H come from Monte-Carlo errors. We repeated the same experiment on many other examples and we always observe the desired order two as long as H is small enough.

Numerical experiments on NLS equation with white noise dispersion

We now apply the algorithms to solve on the torus T r¡π, πs the following SPDE of the form (1.3), with a polynomial linearity and the stiffness parameter ε 10 ¡2 , where the unknown u is a random process depending on x T and t ¥ 0. We consider a spectral discretization in space of this equation with K x 2 7 modes upx, tq °|l|¤Kx Y l ptqe ilx .

We obtain an equation of the desired form (1.1) with a truncated nonlinearity and the block-diagonal matrix

A diagp¡2πl 2 ¢ 0 ¡1 1 0 , |l| ¤ K x q.
Beginning with the initial condition u 0 pxq expp¡3x 4 x 2 q on T that decreases fast enough, we apply Methods A and B in the two cases σ 2 and σ 4 with K t 2 6 modes, N 10 revolutions, m 150 iterations and a tolerance of 10 ¡13 for the fixed point iteration. Figure 5 shows the evolution in time of one trajectory given by Method B (with a 300 points evaluation grid in space).

In Figure 6, we observe the discrete L 2 and H 1 norms behavior of one trajectory given by our two algorithms and the Euler method (3.15) (the simulated pα k q k are the same for Methods A and B). The Euler method quickly blows up in both norms. The L 2 norm of Method A is not conserved. In contrast, Method B preserves the L 2 norm according to Theorem 4.2. When σ 4, numerical simulations hint that a blow-up in the H 1 norm always happens for all considered methods at a certain time that increases as ε goes to zero. We recall that in the optic fiber model (1.3), t represents the distance along the optic fiber and a cubic nonlinearity (σ 2) is typically considered [START_REF] Garnier | Stabilization of dispersion-managed solitons in random optical fibers by strong dispersion management[END_REF]. For σ 2, we do not observe any blow-up in the H 1 norm in Figure 6, suggesting the well-posedness of the model for all optic fiber distance. Also, the larger σ is, the sooner the blow-up happens. These behaviors agree with the blow-up conjecture for ε 1 and σ ¥ 4 presented in [START_REF] Belaouar | Numerical analysis of the nonlinear Schrödinger equation with white noise dispersion[END_REF], and suggest that the conjecture persists in the highly-oscillatory regime ε 3 1. F P pϕ ε,s pyqqpf y ϕ ε,s pyqphq, f y ϕ ε,s pyqpkqqsds.

Then § § f Then the Gronwall lemma allows to conclude. The proof is similar for the third derivative.

¥

Figure 1 :

 1 Figure 1: Exact solution evaluated at revolution times for the deterministic oscillator (2.1) with F pyq iy and ε 10 ¡1 .

Figure 2 :

 2 Figure 2: Revolution times (2.2) of a Brownian path (top) and exact solution evaluated at revolution times for the Kubo oscillator (1.2) with a 1 and ε 10 ¡1 (bottom).

Remark 3 . 5 .

 35 If we replace the Brownian motion W in (3.4) by a piecewise linear function W τ defined by

Figure 3 :

 3 Figure 3: Convergence in law of c N V arpT1q T N N ¡ ErT 1 s ¨to a standard Gaussian random variable.

Repzq

π 2 8 . 8 ,

 88 In addition, for x 0, π 2 Ere xT 1 s 1 cosp c 2xq . The variable T 1 has bounded moments, and they are given by ErT k 1 s p¡2q k k! p2kq! p j1 p¡1q j ņ1 ¤¤¤ n j p n i N ¦ ¢ 2p 2n 1 , . . . , 2n j .

  2.3], where the Laplace transform formula is obtained with an analytic continuation of the equality Ere ¡xT 1 s 1 coshp c 2xq for x ¡ 0. Comparing the Taylor expansions of Ere xT 1 s and 1 cosp c 2xq yields (3.9). The estimate (3.10) is proved as follows

» T N s e 2iπpW prq dr. Then, multiplying by 1 N 2 e 2 » T N 0 e 2iπkW psq » T N s e 2iπpW prq dW prqds ¡ 2π 2 p 2 N 2 »e

 12222 2iπkW psq and integrating from 0 to T N yieldsα N k ¡ α N p k N 2iπp N 2iπkW psq dsdW prq,which has zero average by the same arguments as before. Also by the Fubini theorem for stochastic integrals, β N p,k 1 N 2

,qξ l with δ k 4 1

 4 and Γ N its square root. Then, p α N k is defined for k ¥ 0 as p

  Erβ N p,k s with the values of Proposition 3.7. Doing so yields Method A.For Method B, we adapt Method A in the spirit of the middle point scheme for ODEs (see[START_REF] Hairer | Geometric numerical integration[END_REF] Chap. IV]) so that it preserves any quadratic invariant. We also replace p Er r β N p,k s, using the values of Proposition 3.7.

c 2 l pc 0 p , c 0 k 1 N 2 »e

 0012 qpY m qp γ p2q,N l,p,k with the new random variables γ p1q,N l,p,k 2iπkW pqq dqdrds, and where the discrete random variables share the same moments up to order 3 for the p , and order 1 for the p γ piq,N

k

  pr ϕ pKtq ε,s pyqqe 2iπkW psq ds.

Proposition 4 . 4 (

 44 bounded moments for the integrator (2.4)). Assume that for y R d , the numerical integrator p ψ ε,N pyq is given by

Figure 4 ,

 4 Figure 4, we use a modification of a Kubo oscillator introduced in [8] with the nonlinearity F pyq ip1 Repyq 3 Impyq 5 qy. In the real setting, it yields the following two-dimensional SDE

Figure 4 :

 4 Figure 4: Weak error versus the stepsize H N ε for approximating the solution of equation (1.1) at time εT 2 8 for the linear F pyq iy (left) and the non-linear F pyq ip1 Repyq 3 Impyq 5 qy (right) with A 2iπ, X 0 1, ε 10 ¡3 and the test function φpyq 2 Repyq 4 Impyq.

Figure 5 :

 5 Figure 5: Approximation by Method B of |u| and |f x u| with u solution of a spatial discretization with K x 2 7 modes of the nonlinear Schrödinger equation with white noise dispersion (5.1) on the torus T r¡π, πs with the parameters ε 10 ¡2 , σ 2 (top) and σ 4 (bottom).

Figure 6 :

 6 Figure 6: Evolution in long time of the quantities }U t } L 2 ¡}u 0 } L 2 (left) and }U t } H 1 ¡}u 0 } H 1 (right) with U t the approximation computed with Euler method and Methods A and B for one trajectory of equation (5.1) with ε 10 ¡2 , σ 2 (top) and σ 4 (bottom).

k end for Remark 2.1. One

	could apply a Newton iteration to solve the implicit equation (2.5)
	in Method B. However, a few fixed point iterations are sufficient (see discussion in [16,
	Chap. VIII] for non-stiff implicit methods). Indeed, since the Lipschitz constant of the
	iterated map has size OpHq, the convergence rate of the fixed point iterations is independent
	of the smallness of ε.
	Remark 2.2. We observe that p β N p,k and p r β N p,k are always zero except when p 0, k 0 or
	p k 0. Thus the computational cost of one step of Methods A and B grows linearly in
	the number of modes in (2.3).

  The dominated convergence theorem for stochastic integrals allows us to take the limit t Ñ V and yields Erα N k s 0. , let pp, kq $ p0, 0q (the case p k 0 is obtained straightforwardly using Proposition 3.6), we use the Itô formula on e 2iπpW psq and we integrate from s

	α N k i πkN	» T N
	For the coefficients β N p,k to T N , 1 ¡ e 2iπpW psq 2iπp	» T N

the case k 0 is straightforward using Proposition 3.6), then the Itô formula applied to e 2iπkW psq gives 1 2π 2 k 2 N pe 2iπkW ptq ¡ 1q i πkN » t 0 e 2iπkW psq dW psq ¡ 1 N » t 0 e 2iπkW psq ds, which yields at time t T N , 0 e 2iπkW psq dW psq.

Then t Þ Ñ ³ t 0 e 2iπkW psq dW psq is a martingale, so by the Doob theorem, as t T N is finite, Er ³ tT N 0 e 2iπkW psq dW psqs 0 for all t.

s e 2iπpW prq dW prq ¡ 2π 2 p 2

  2, it is now possible to evaluate the local expansions (3.4) at revolution times. To approximate numerically the integrals appearing in equation(3.4) without evaluating F and F I too many times, we first replace the 1-periodic functions g 0 pyqpzq defined in (2.3) by their associated Fourier series with Fourier coefficients pc 0 k pyqq kZ and pc 1 p pyqq pZ . We define the following approximation of ϕ ε,t pyq, ψ ε,t pyq e AW ptq y ε

	and g 1								θ pyq
			ķZ	e AW ptq c 0 k pyq	» t	e 2iπkW psq ds	(3.11)
	ε 2	p	e AW ptq c 1 p pyq	¢ c 0 k pyq	» t	0 » s	e 2iπpW psq e 2iπkW prq drds	.
		,kZ					0		0

θ

  1 p pyqq pZ are the Fourier coefficients of the 1-periodic functions g 0 θ pyq

	and g 1 θ pyq defined in (2.3), α N k , β N p,k are the random variables defined in Proposition 3.7 and
	y R d is deterministic. Under Assumption 3.1, for all test functions φ C 3

P , there exists H 0 ¡ 0 such that for all H N ε ¤ H 0 , the following estimates hold, where C and K are independent of ε and N , E

  P , there exists H 0 ¡ 0 such that for all H ¤ H 0 , for all m ¥ 0 such that mN ε mH ¤ T , there exist two positive constants K and C both independent of ε, N and K t such that|Erφpϕ ε,T N m pX 0 qqs ¡ ErφpY m qs| ¤ CH 2 p1 Er|X 0 | K sq.

	Theorem 4.1. Assume that the Fourier coefficients c 0 k , c 1 p of g 0 θ and g 1 θ in (2.3) are non-zero only for ¡K

t {2 ¤ k, p K t {2. Then, under Assumption 3.1, Method A has weak order two, that is, for all T ¡ 0, for all test functions φ C 3 (4.1) Theorem 4.2. Assume that the Fourier coefficients c 0 k , c 1 p of g 0 θ and g 1 θ in (2.3) are nonzero only for ¡K t {2 ¤ k, p K t {2. Then, under Assumption 3.1, if c 1 0 pc 0 k q and c 1 p pc 0 0 q are Lipschitz continuous uniformly in k and p, Method B is well defined and has weak order two (i.e., it satisfies an estimate of the form (4.1)). In addition, if for a fixed symmetric matrix S R d¢d , the quantity Qpyq 1 2 y T Sy is preserved by equation (1.1), then Method B also preserves the invariant Qpyq 1 2 y T Sy, that is, the solution Y m 1 of equation (2.5)

  Opp1 |y|qe ¡cKt q. For simplicity of the analysis, we thus assumed in Theorems 4.1 and 4.2 that g 0 θ and g 1 θ have a finite number K t of non-zero Fourier modes in (2.3). If this assumption does not hold, the truncation errors Opp1 |y| K qK ¡s t q or Opp1 |y|qe ¡cKt q should be added in the right-hand side of the error estimate (4.1). Let us prove it in the analytic case. We first apply the change of variable r ϕ ε,t pyq e ¡AWptq ϕ ε,t pyq that has no effect at time t T N m . We now have to compare the two solutions of the following integral

	formulations
	r ϕ ε,t pyq y ε	» t	g 0
		0

.1) at the revolution times εT N m , m 0, 1, . . . , but one could compute an approximation at different times by composing with other methods at the end of the integration. Remark 4.3. Since the error constant C in (4.1) is independent of the number K t of Fourier modes, we emphasize that Theorems 4.1 and 4.2 remain valid for infinitely many modes (K t Ñ V). In addition, assuming that F is of class C s 1 P yields a truncation error of the Fourier series in (2.3) of size Opp1 |y| K qK ¡s t q (see, e.g., [19, Sect. III.1.3]), and if g 0 θ is assumed analytic in θ (for example if F is a polynomial), this error becomes exponentially small as W psq pr ϕ ε,s pyqqds,

e

  M Erφpϕ ε,T N M pX 0 qqs ¡ Erφp p m¡1 φ ¥ ϕ ε,T N pm¡1q . Using Lemma 3.2 and φ C 3 P , we obtain for 0

				ψ M ε,N pX 0 qqs
	and rewrite it with a telescopic sum
	e M	M m1	Erφpϕ ε,T N pm¡1q p p ψ M ¡m 1 ε,N	pX 0 qqqs ¡ Erφpϕ ε,T N m p p ψ M ¡m ε,N pX 0 qqqs
		M m1	Er r φ m¡1 p p ψ ε,N p p ψ M ¡m ε,N pX 0 qqqs ¡ Er r φ m¡1 pϕ ε,T N p p ψ M ¡m ε,N pX 0 qqqs

where r φ

  We deduce the global weak order two. ¥With the help of Theorem 4.7, we prove Proposition 3.10 and the convergence of Methods A and B.Proof of Proposition 3.10. Rewriting Theorem 4.7 for order one yields for all H N ε small enough and all M ¥ 0,|Erφpϕ ε,T N M pX 0 qqs ¡ Erφpy εN M qs| ¤ CpεN q 2 p1 Er|X 0 | K sq.Proof of Theorem 4.2. The regularity assumptions yield the Lipschitzness of the c 0

	α N k ¤ C and °p,k 4.7 applies and concludes the proof.	§ § §Erβ k 2 N p,k s	§ § 2 §	converges by Proposition 3.7, Theorem ¥
	the involved c 1 p pyqpc 0 k pyqq with constants independent of k and p. As	°k § § p α N k	k pyq and § and °p,k § § § § § p r β N p,k § § § §

Evaluating in N 1, M T ε and taking the limit ε Ñ 0 yield the result.

¥

Proof of Theorem 4.1. As p are bounded, the right-hand side of equation (2.5) is a contraction for all H ¤ H 0 small enough and the constant does not depend on Y m , so H 0 depends only on F and F I . Thus, the integrator is well-posed for all H ¤ H 0 .

  2 y ϕ ε,t pyqph, kq § § ¤ Cε » t 0 rp1 |ϕ ε,s pyq| K q |f y ϕ ε,s pyqphq| |f y ϕ ε,s pyqpkq| § § f 2 y ϕ ε,s pyqph, kq § § sds ¤ Cεtp1 |y| K qe Cεt |h| |k| Cε

	» t	§ § f 2 y ϕ ε,s pyqph, kq	§ § .
	0

Université de Genève, Section de mathématiques,

2-4 rue du Lièvre, CP 64, CH-1211 Genève 4, Switzerland. Adrien.Laurent@unige.ch, Gilles.Vilmart@unige.ch.

Acknowledgments

The authors would like to thank Georg Gottwald for helpful and stimulating discussions. This work was partially supported by the Swiss National Science Foundation, grants No. 200020_184614, No. 200021_162404 and No. 200020_178752. The computations were performed at the University of Geneva on the Baobab cluster using the Julia programming language.

Appendix

Proof of Lemma 3.2.

1. First, ϕ ε,t pyq is the solution of ϕ ε,t pyq e AW ptq y εe AW ptq » t 0 e ¡AWpsq F pϕ ε,s pyqqds.

Using the boundedness of the continuous periodic function θ Ñ e θA and Assumption 3.1, we get

The Gronwall lemma yields the desired bound.


Paul Lerner
email: paul.lerner@lisn.upsaclay.fr

Olivier Ferret
email: olivier.ferret@cea.fr

Camille Guinaudeau
email: camille.guinaudeau@lisn.upsaclay.fr

Multimodal Inverse Cloze Task for Knowledge-based Visual Question Answering

Keywords: Visual Question Answering, Pre-training, Multimodal Fusion

We present a new pre-training method, Multimodal Inverse Cloze Task, for Knowledge-based Visual Question Answering about named Entities (KVQAE). KVQAE is a recently introduced task that consists in answering questions about named entities grounded in a visual context using a Knowledge Base. Therefore, the interaction between the modalities is paramount to retrieve information and must be captured with complex fusion models. As these models require a lot of training data, we design this pre-training task, which leverages contextualized images in multimodal documents to generate visual pseudo-questions. Our method is applicable to different neural network architectures and leads to a 9% relative-MRR and 15% relative-F1 gain for retrieval and reading comprehension, respectively, over a no-pre-training baseline.

Introduction

Knowledge-based Visual Question Answering about named Entities (KVQAE) is a challenging task recently introduced in [START_REF] Shah | KVQA: Knowledge-Aware Visual Question Answering[END_REF]. It consists in answering questions about named entities grounded in a visual context using a Knowledge Base (KB). Figure 1 provides two examples of visual questions along with relevant visual passages from a KB. To address the task, one must thus retrieve relevant information from a KB. This contrasts with standard Visual Question Answering (VQA [START_REF] Antol | VQA: Visual Question Answering[END_REF]), where questions target the content of the image (e.g. the color of an object or the number of objects) or Knowledge-based VQA (about coarse-grained object categories) [START_REF] Marino | OK-VQA: A visual question answering benchmark requiring external knowledge[END_REF], where one can rely on off-the-shelf object detection [START_REF] Gardères | ConceptBert: Concept-Aware Representation for Visual Question Answering[END_REF]. In KVQAE, both text and image modalities bring useful information that must be

Visual Question (input)

Relevant visual passage in the Knowledge Base "Which constituency did this man represent when he was Prime Minister?" "In which year did this ocean liner make her maiden voyage?"

Queen Elizabeth 2, often referred to simply as QE2, is a floating hotel and retired ocean liner built for the Cunard Line which was operated by Cunard as both a transatlantic liner and a cruise ship from 1969 to 2008.

Macmillan indeed lost Stockton in the landslide Labour victory of 1945, but returned to Parliament in the November 1945 by-election in Bromley.

Fig. 1: Example of visual questions about named entities from the ViQuAE dataset along with relevant visual passages from its Knowledge Base [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF].

combined. Therefore, the task is more broadly related to Multimodal Information Retrieval (IR) and Multimodal Fusion.

There are two paradigms for multimodal IR and for multimodal learning more generally: early fusion (data-and feature-level) and late fusion (score-and decision-level) [START_REF] Kludas | Information fusion in multimedia information retrieval[END_REF]. On the one hand, late fusion is more straightforward as both Natural Language Processing and Computer Vision techniques can be applied independently. However, it neglects interactions between modalities. For example, in Figure 1, attempting to recognize Harold Macmillan without accounting for him being a Prime Minister is suboptimal. On the other hand, the richness of early fusion often comes at the cost of increasing complexity and model parameters. This adds an extra challenge for KVQAE, where the two existing datasets are either small, because of a costly annotation process (ViQuAE [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF]), or generated automatically (KVQA [START_REF] Shah | KVQA: Knowledge-Aware Visual Question Answering[END_REF]), which leads to several limitations discussed in [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF]. To address this challenge, we propose the multimodal Inverse Cloze Task (ICT) for pre-training the two early fusion models we define in this article for tackling KVQAE. Multimodal ICT consists in considering a sentence paired with a nearby image as a visual pseudo-question and its multimodal context as a relevant visual passage. It is related to the visual cloze task proposed by [START_REF] Zhang | Modeling temporal-modal entity graph for procedural multimodal machine comprehension[END_REF], a downstream task that requires modeling temporal events. Textual ICT was first introduced in [START_REF] Lee | Latent Retrieval for Weakly Supervised Open Domain Question Answering[END_REF] to pre-train a neural retriever for textual Question Answering (QA) and can be seen as a generalization of the skip-gram objective [START_REF] Mikolov | Distributed Representations of Words and Phrases and their Compositionality[END_REF].

Our main contributions are: (i) Multimodal ICT, a new pre-training method that allows tackling KVQAE, even for small datasets as ViQuAE; (ii) a multimodal IR framework for KVQAE; (iii) experiments with different neural network architectures, including recently proposed multimodal BERTs.

Related Work

Dense Retrieval. Dense Retrieval is a rapidly evolving field, surveyed in [START_REF] Lin | Pretrained transformers for text ranking: Bert and beyond[END_REF][START_REF] Fan | Pretraining Methods in Information Retrieval[END_REF], with new pre-training tasks, optimizing methods, and variants of the Transformer architecture emerging [START_REF] Ram | Learning to Retrieve Passages without Supervision[END_REF][START_REF] Hofstätter | Efficiently teaching an effective dense retriever with balanced topic aware sampling[END_REF][START_REF] Gao | Unsupervised corpus aware language model pre-training for dense passage retrieval[END_REF][START_REF] Gao | Condenser: a Pre-training Architecture for Dense Retrieval[END_REF]. [START_REF] Lee | Latent Retrieval for Weakly Supervised Open Domain Question Answering[END_REF] were the first to outperform sparse bagof-words representations such as BM25 with dense representations for QA. Their approach relies on three components: (i) pre-trained language models such as BERT [START_REF] Devlin | BERT: Pre-training of deep bidirectional transformers for language understanding[END_REF], which allow to encode the semantic of a sentence in a dense vector; (ii) a contrastive learning objective that optimizes the similarities between questions' and text passages' embeddings (see Section 3); (iii) an unsupervised training task, ICT (see Section 1). [START_REF] Karpukhin | Dense passage retrieval for open-domain question answering[END_REF] criticize the latter for being computationally intensive 3 and argue that regular sentences are not good surrogates of questions. Instead, they propose DPR, which takes advantage of (i) the heuristic of whether the passage contains the answer to the question to deem it relevant; (ii) unsupervised IR methods such as BM25 to mine hard negatives examples, which proved to be the key of their method's success. We aim to take advantage of both approaches by (i) pre-training our model on text QA datasets like DPR; (ii) incorporating multimodality into this hopefully-well-initialized model by adapting the ICT of [START_REF] Lee | Latent Retrieval for Weakly Supervised Open Domain Question Answering[END_REF] to multimodal documents.

Multimodal Fusion and Pre-Training. The success of BERT in NLP [START_REF] Devlin | BERT: Pre-training of deep bidirectional transformers for language understanding[END_REF], which relies on the easily-parallelizable Transformer architecture [START_REF] Vaswani | Attention is all you need[END_REF], an unsupervised training objective, and a task-agnostic architecture, has concurrently inspired many works in the VQA and cross-modal retrieval fields [START_REF] Tan | LXMERT: Learning Cross-Modality Encoder Representations from Transformers[END_REF][START_REF] Lu | ViLBERT: Pretraining Task-Agnostic Visiolinguistic Representations for Vision-and-Language Tasks[END_REF][START_REF] Li | Visualbert: A simple and performant baseline for vision and language[END_REF][START_REF] Su | Vl-bert: Pre-training of generic visual-linguistic representations[END_REF][START_REF] Li | Unicoder-VL: A Universal Encoder for Vision and Language by Cross-Modal Pre-Training[END_REF][START_REF] Chen | Uniter: Universal image-text representation learning[END_REF]. These models are unified under a single framework in [START_REF] Bugliarello | Multimodal Pretraining Unmasked: A Meta-Analysis and a Unified Framework of Vision-and-Language BERTs[END_REF] and partly reviewed in [START_REF] Khan | Transformers in vision: A survey[END_REF]. All of these models rely on the Transformer architecture, often initialized with a pre-trained BERT, in order to fuse image and text. The training is weakly supervised, based upon image caption datasets such as COCO [START_REF] Lin | Microsoft COCO: Common Objects in Context[END_REF] or Conceptual Captions [START_REF] Sharma | Conceptual captions: A cleaned, hypernymed, image alt-text dataset for automatic image captioning[END_REF], and pre-trained object detectors like Faster R-CNN [START_REF] Ren | Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks[END_REF]. [START_REF] Hessel | Does my multimodal model learn cross-modal interactions? it's harder to tell than you might think![END_REF] show that these models learn nontrivial interactions between the modalities for VQA. Multimodal BERTs can be broadly categorized into single-stream and multi-stream. Single-stream models feed both text tokens' embeddings and image regions' embeddings to the same Transformer model, relying on the self-attention mechanism to fuse them. Instead, in the multi-stream architecture, text and image are first processed by two independent Transformers before using crossattention to fuse the modalities. Both architectures have been shown to perform equally well in [START_REF] Bugliarello | Multimodal Pretraining Unmasked: A Meta-Analysis and a Unified Framework of Vision-and-Language BERTs[END_REF]. In this work, we use a single-stream model to take advantage of pre-training on text-only (on QA datasets). Also note that, while inspired by these work, we do not use the same training objectives or data, which are arguably unsuited for named entities' representations, as explained in the next section.

Multimodal Information Retrieval and KVQAE. Multimodal IR has largely been addressed using late fusion techniques (see [START_REF] Depeursinge | Fusion Techniques for Combining Textual and Visual Information Retrieval[END_REF] for a survey) but we are mostly interested in early fusion techniques in this work. [START_REF] Depeursinge | Fusion Techniques for Combining Textual and Visual Information Retrieval[END_REF] review first attempts at early fusion. It was then systematically done by concatenating the features of both modalities in a single vector, with a focus on the feature weighting scheme. Concatenation is confronted with the curse of dimensionality as the resulting feature space equals the sum of the dimensions of each modality's features.

[44] and [START_REF] Luo | Weakly-supervised visual-retriever-reader for knowledge-based question answering[END_REF] concurrently proposed an approach quite similar to ours for Knowledge-based VQA. They adapt DPR [START_REF] Karpukhin | Dense passage retrieval for open-domain question answering[END_REF] by replacing the question encoder with LXMERT [START_REF] Tan | LXMERT: Learning Cross-Modality Encoder Representations from Transformers[END_REF], which allows to fuse the question and image. However, unlike us, they keep the passage encoder based on text-only and use the same pretraining objectives as [START_REF] Tan | LXMERT: Learning Cross-Modality Encoder Representations from Transformers[END_REF], namely Masked Language Modeling, Masked Region Modeling, and Image-Text Matching. We expect that these objectives are suited to learn representations of coarse-grained object categories but not named entities. In other words, they are suited for standard VQA but not KVQAE. For example, Masked Region Modeling relies on an object detector, which is not applicable to KVQAE. While both [START_REF] Qu | Passage retrieval for outside-knowledge visual question answering[END_REF] and [START_REF] Luo | Weakly-supervised visual-retriever-reader for knowledge-based question answering[END_REF] experiment on OK-VQA [START_REF] Marino | OK-VQA: A visual question answering benchmark requiring external knowledge[END_REF], their results are inconsistent: [START_REF] Qu | Passage retrieval for outside-knowledge visual question answering[END_REF] show that their model is competitive with a BM25 baseline that takes as input the question and the human-written caption of the image while the model of [START_REF] Luo | Weakly-supervised visual-retriever-reader for knowledge-based question answering[END_REF] is outperformed by BM25 with an automatically-generated caption. The discrepancies between these works can be explained because they use neither the same KB nor the same evaluation metrics. [START_REF] Guo | A unified end-to-end retriever-reader framework for knowledge-based vqa[END_REF] also experiment with different multimodal BERTs but dispense passage-level annotation for an end-to-end training of the retriever and answer classifier 4 .

Although they experiment with KVQA [START_REF] Shah | KVQA: Knowledge-Aware Visual Question Answering[END_REF], we do not consider the work of [START_REF] Garcia-Olano | Improving and diagnosing knowledgebased visual question answering via entity enhanced knowledge injection[END_REF][START_REF] Heo | Hypergraph Transformer: Weaklysupervised multi-hop reasoning for knowledge-based visual question answering[END_REF] as their systems take a human-written caption as input, which makes the role of the image content unclear. [START_REF] Shah | KVQA: Knowledge-Aware Visual Question Answering[END_REF] follow a late fusion approach at the decision-level. First, they detect and disambiguate the named entity mentions in the question. Then, they rely on a face recognition step as their dataset, KVQA, is restricted to questions about person named entities. Facts from both textuallyand visually-detected entities are retrieved from Wikidata5 and processed by a memory network [START_REF] Weston | Memory networks[END_REF]. In contrast, our work is in line with [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF], who use unstructured text from Wikipedia as KB. Unlike the late fusion approach of [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF], which considers the question and the image independently, we aim at a unified representation of the text and image, both on the visual question and KB sides.

Methods

In this section, we first formalize our KVQAE framework, then describe the models before diving into the three training stages: (i) DPR for textual Question Answering; (ii) Multimodal Inverse Cloze Task, our main contribution;

(iii) Fine-tuning for KVQAE. Finally, we discuss the inference mechanism and implementation details.

Information Retrieval Framework

In our multimodal setting, both visual questions (from the dataset) and visual passages (from the KB) consist of a text-image pair (t, i), as in Figure 1. Our goal is to find the optimal model E to encode adequate representations q = E(t q , i q ) and p = E(t p , i p ) such that they are close if (t p , i p ) is relevant for (t q , i q ) (denoted with the superscripts ( + ) and ( -)). Search then boils down to retrieving the K closest visual passages to the visual question. When computing the similarity between two vectors, here with the dot product, the objective used throughout all the training stages ( §3.3) is to minimize the following negative log-likelihood loss for all visual questions in the dataset (see Figure 2) [START_REF] Lee | Latent Retrieval for Weakly Supervised Open Domain Question Answering[END_REF][START_REF] Karpukhin | Dense passage retrieval for open-domain question answering[END_REF].

-log exp (q • p + ) exp (q • p + ) + j exp (q • p - j ) (1)
This contrastive objective allows to efficiently utilize passages relevant to other questions in the batch as in-batch negatives, since computing the similarity between two vectors is rather inexpensive compared to the forward pass of the whole model. We present two different models E in the next section according to their fusion mechanism.

Models

All of the models described in this section take advantage of CLIP6 [START_REF] Radford | Learning transferable visual models from natural language supervision[END_REF] to represent images and BERT7 [START_REF] Devlin | BERT: Pre-training of deep bidirectional transformers for language understanding[END_REF] to represent either text or multimodal data.

BERT is trained for masked language modeling and next sentence prediction on Wikipedia and BooksCorpus [START_REF] Zhu | Aligning books and movies: Towards story-like visual explanations by watching movies and reading books[END_REF]. CLIP has been trained with a contrastive objective in a weakly-supervised manner over 400M image and caption pairs. It has demonstrated better generalization capacities than fully-supervised models and is efficient for KVQAE, as empirically demonstrated in [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF]. We experiment with two different fusion techniques: ECA and ILF.

Early Cross-Attention fusion (ECA) is carried out by a single-stream Transformer model like the multimodal BERTs described above (e.g. UNITER [START_REF] Chen | Uniter: Universal image-text representation learning[END_REF]). However, instead of relying on a fixed object detector such as Faster R-CNN, we take advantage of CLIP, as motivated above. To enable early fusion, the visual embedding produced by CLIP is projected in the same space as the text using a linear layer with W c ∈ R c×d parameters trained from scratch: e c = CLIP(i) • W c . The resulting embedding is then concatenated with the word embeddings of the text, acting as an additional "visual token".

Cross-entropy

Barack Obama

… …

U.S. Senate (2005-2008)

A sentence about Obama being a senator. In the context we find related sentences.

Believed to be relevant. Intermediate Linear Fusion (ILF) introduces an additional W t ∈ R d×d parameters trained from scratch used to simply project the representation of the [CLS] token in the same space as the CLIP embedding before summing the two8 :

ILF(t, i) = BERT(t) [CLS] • W t + e c .
Because both ECA and ILF produce multimodal representations q and p, ranking is done directly using their similarity q • p. As a baseline, we follow [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF] and linearly combine text and image similarities after zero-mean and unit-variance normalization (omitted in the following equation):

α × BERT(t q ) [CLS] • BERT(t p ) [CLS] + (1 -α) × cos(CLIP(i q ), CLIP(i p )) (2)
The interpolation hyperparameter α is optimized on the validation set using grid search to maximize Mean Reciprocal Rank. The left term (text similarity) is referred to as DPR in the rest of the paper.

Training stages

The models are trained sequentially in three stages. The first two stages correspond to pre-training stages: the first one is dedicated to IR for Question Answering while the second one focuses on multimodal IR. The last stage corresponds to the training for our target task of IR for KVQAE.

Stage 1: DPR for textual Question Answering. Leaving visual representations aside, a DPR model is trained starting from the BERT initialization [START_REF] Karpukhin | Dense passage retrieval for open-domain question answering[END_REF]. DPR consists of two BERT encoders: one for the question t q and one for the text passage t p . We use the model pre-trained by [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF] on TriviaQA, filtered of all questions used in their dataset, ViQuAE. They use the KILT [START_REF] Petroni | KILT: a benchmark for knowledge intensive language tasks[END_REF] version of TriviaQA and Wikipedia, which serves as KB at this stage. Each article is then split into disjoint passages of 100 words for text retrieval, while preserving sentence boundaries, and the title of the article is appended to the beginning of each passage. This yields 32M passages, that is ≈ 5.4 passages per article. Following [START_REF] Karpukhin | Dense passage retrieval for open-domain question answering[END_REF], irrelevant passages (i.e. hard negatives) are mined using BM25 [START_REF] Robertson | Okapi at TREC-3[END_REF].

Stage 2: Multimodal Inverse Cloze Task. This is the main contribution of the paper. We propose to extend the ICT of [START_REF] Lee | Latent Retrieval for Weakly Supervised Open Domain Question Answering[END_REF] to multimodal documents. ICT consists in considering a sentence as a pseudo-question t q and its context as a relevant passage t + p . Note that the title of the article is appended to the beginning of each passage t p (as in Stage 1). We extend it using the contextual images of Wikipedia paragraphs for the pseudo-question and the infobox image for the passage (see Figure 2). [START_REF] Lee | Latent Retrieval for Weakly Supervised Open Domain Question Answering[END_REF] empirically demonstrated that a key success of their approach was to leave the pseudo-question in the relevant passage in 10% of the training samples so that the model will learn to perform word matching, as lexical overlap is ultimately a very useful feature for retrieval. In our case, however, we argue that it is neither necessary, as the model should be strongly initialized from Stage 1 training on TriviaQA, nor beneficial, as the model could then ignore the image modality. Question and passage encoders pre-trained in Stage 1 are used to initialize the visual question and visual passage encoders, respectively.

The process is eased thanks to the WIT dataset [START_REF] Srinivasan | Wit: Wikipediabased image text dataset for multimodal multilingual machine learning[END_REF]. WIT consists of millions of images with associated text from Wikipedia and Wikimedia Commons in 108 different languages. We are, however, only interested in English for this work. While [START_REF] Srinivasan | Wit: Wikipediabased image text dataset for multimodal multilingual machine learning[END_REF] have multiple strategies to find text related to a given Wikipedia image, such as its Commons' caption, we use only the contextual paragraph as text source in order to mimic the downstream KVQAE setting. The resulting English subset of WIT yields 400K infobox images/articles that correspond to 1.2M paragraphs/images. Those 1.2M paragraphs consist of 13.6M sentences, i.e. potential pseudo-questions, which are 26 words long on average. Therefore, to stick as close as possible to stages 1 and 3, where passages are up to 100 words long, passages consist of four sentences. This slightly differs from [START_REF] Lee | Latent Retrieval for Weakly Supervised Open Domain Question Answering[END_REF] who consider passages of up to 288 wordpieces, prior to the pseudo-question masking.

Because both ViQuAE and WIT images are taken from Wikimedia Commons9 , we can estimate from the image URLs that 14% of ViQuAE images overlap with WIT. This might lead to a bias that we analyze in Section 4.1.

Inspired by [START_REF] Aytar | Cross-modal scene networks[END_REF], to prevent catastrophic forgetting and enforce a modalityinvariant representation of the entities, the last l layers of BERT are frozen during this stage. In this way, we tune only the first, modality-specific layers of ECA, Stage 3: Knowledge-based Visual Question Answering about Named Entities. This stage consists in fine-tuning the model on a downstream KVQAE dataset, which provides visual questions (t q , i q ) and relevant visual passages (t + p , i + p ). Following [START_REF] Aytar | Cross-modal scene networks[END_REF], all layers of the model are tuned during this stage. A subtlety of this stage is the selection of irrelevant visual passages (t - p , i - p ). As mentioned in Section 2, it was shown to be essential to DPR [START_REF] Karpukhin | Dense passage retrieval for open-domain question answering[END_REF], and it is more generally important for contrastive learning [START_REF] Kalantidis | Hard Negative Mixing for Contrastive Learning[END_REF]. In [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF], irrelevant passages are mined with BM25 to train DPR. However, we suppose that this is suboptimal for ECA and ILF as BM25 will only mine textually-plausible passages but not visually-plausible ones. Therefore, we use the system provided by [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF] to mine irrelevant passages. It is a late-fusion of DPR, ArcFace [START_REF] Deng | Arcface: Additive angular margin loss for deep face recognition[END_REF], CLIP, and ImageNet-ResNet [START_REF] He | Deep residual learning for image recognition[END_REF]. This leads to different training setups between DPR (used as a baseline) and our models. However, we have experimented both for DPR and found no significant differences 10 .

We use the same KB as [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF], which is based upon KILT's Wikipedia and Wikidata images of the corresponding entities. It consists of 1.5M articles (thus images/entities) split into 12M passages of at most 100 words as in Stage 1.

Visual questions in ViQuAE are split into train [START_REF] Antol | VQA: Visual Question Answering[END_REF]190), validation (1,250), and test [START_REF] Antol | VQA: Visual Question Answering[END_REF]257) without overlap between images' URLs [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF]. We do not experiment with KVQA [START_REF] Shah | KVQA: Knowledge-Aware Visual Question Answering[END_REF] for the following reasons: (i) it is generated automatically from Wikidata so our text-based KB has a poor coverage of the answers; (ii) it comprises yes/no questions for which passage relevance cannot be assessed automatically.

Inference

For efficient retrieval, every passage in the KB is embedded along with its corresponding image by the visual passage encoder beforehand. Given a question grounded in an image, both are embedded by the visual question encoder. Search is then carried out with maximum inner product search using Faiss [START_REF] Johnson | Billion-scale similarity search with GPUs[END_REF].

Implementation Details

Our code is built upon PyTorch [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF], Hugging Face's transformers [START_REF] Wolf | HuggingFace's Transformers: State-of-the-art Natural Language Processing[END_REF] and datasets [START_REF] Lhoest | Datasets: A Community Library for Natural Language Processing[END_REF] (itself wrapping Faiss). It is freely available along with the data and trained models 11 .

To train ECA, we use the same hyperparameters as [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF] for DPR, themselves based upon [START_REF] Karpukhin | Dense passage retrieval for open-domain question answering[END_REF]. In particular, we use a learning rate of 2 × 10 -5 along with the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]. It is scheduled linearly with 100 and 4 warm-up steps for stages 2 and 3, respectively. However, for ILF, we found, based on the validation set, that it converged faster with a learning rate of 2 × 10 -3 and a constant scheduler during Stage 2. We believe this is because ILF fully freezes BERT in Stage 2, so it does not require careful scheduling or a small learning rate. Dropout [START_REF] Srivastava | Dropout: a simple way to prevent neural networks from overfitting[END_REF] is applied in BERT and after projecting embedding with W c and W t with a probability of 0.1 (as in the standard BERT configuration). Likewise, layer normalization [START_REF] Ba | Layer Normalization[END_REF] is applied in BERT and after summing the two embeddings in ILF. Gradients' norms are clipped at 2. Models in stages 2 and 3 are trained with a batch size of 512 and 298 visual questions, respectively. The success of contrastive learning partly relies on a large number of in-batch negatives and, therefore, a large batch size [START_REF] Qu | RocketQA: An Optimized Training Approach to Dense Passage Retrieval for Open-Domain Question Answering[END_REF]. We found that gradient checkpointing [START_REF] Chen | Training Deep Nets with Sublinear Memory Cost[END_REF] enables the use of much larger batch sizes. Instead of [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF] who use four NVIDIA V100 GPUs with 32GB of RAM each for a total batch size of 128 questions, we are able to fit a batch of 298 questions (as stated above) in a single V100 GPU. Stage 2 takes most of the compute budget, with most models converging after ≈ 8K steps, which takes around three days 12 . Checkpoint selection is made based on the validation in-batch Mean Reciprocal Rank, for all stages. In-batch means that only the other visual passages in the batch are ranked and that each visual question is paired with only one relevant visual passage (as during training).

Results

The retrieval models are evaluated in two different ways: (i) by computing standard IR metrics on visual passage retrieval; (ii) by feeding retrieved visual passages to a reader module that is tasked with extracting the concise answer to the question, thus achieving KVQAE. Put differently, either evaluate whether the system is able to retrieve a relevant passage for the question or whether it is able to answer the question. We find both metrics to correlate. Ablation studies are carried out with IR metrics.

ViQuAE is based upon TriviaQA, so it is only distantly supervised: the answer is considered correct if it string-matches the ground truth and, likewise, a passage Table 1: IR evaluation on ViQuAE. l: Number of frozen layers during Multimodal ICT. Superscripts denote significant differences in Fisher's randomization test with p ≤ 0.01. Hits@1 is omitted as it is equivalent to P@1.

# Model

Multimodal ICT MRR@100 P@1 P@20 Hits@20 is deemed relevant if it contains the ground truth 13 . Moreover, Wikipedia aliases of the ground truth are considered to be valid answers.

Information Retrieval

Because of the setting of ViQuAE, it is impossible to get complete coverage of relevant passages. Therefore we do not use any metric based on recall (e.g. R-Precision, mAP, etc.). Instead, we evaluate the models with Precision@K (P@K), Mean Reciprocal Rank (MRR), and Hits@K. Hits@K is the proportion of questions for which IR retrieves at least one relevant passage in top-K. Statistical significance tests are conducted using Fisher's randomization test [START_REF] Fisher | The design of experiments[END_REF][START_REF] Smucker | A comparison of statistical significance tests for information retrieval evaluation[END_REF]. Metrics and statistical tests are computed with ranx [START_REF] Bassani | ranx: A Blazing-Fast Python Library for Ranking Evaluation and Comparison[END_REF] and are reported in Table 1.

The best models pre-trained with Multimodal ICT (d and f) outperform the text-only (a) and late-fusion (b) baselines on all metrics. Some qualitative examples are shown in Figure 3. In the first row, we can see evidence of cross-input cross-modal interactions between the image depicting Winston Churchill and the passage that mentions him (while being illustrated by a totally different image). In contrast, the late fusion baseline exhibits textual bias by returning a passage that mentions several English palaces (highlighted in red). The same observation can be made for the second row, where St Paul's Cathedral is only mentioned in the relevant passage but not depicted in the contextual image. Cross-modal interactions prove useful in this case because of the heterogeneity of visual depictions: Winston Churchill is depicted by a statue in the visual question but by a photograph in the KB.

We can see that Multimodal ICT is essential to ECA (c vs. d). Without it, it performs on par with late fusion. We believe this is because of overfitting on the small training set of ViQuAE. However, we find that fine-tuning on ViQuAE is also essential to ECA, which exhibits catastrophic forgetting because of the sequential learning setup: indeed, after Stage 2, it falls behind DPR (not shown in the table). We see that the freezing technique of [START_REF] Aytar | Cross-modal scene networks[END_REF] helps to prevent catastrophic forgetting to some extent (d vs. e). It is also visible in the upstream WIT pre-training where ECA achieves 91.6 and 92.9 in-batch MRR on WIT's test set with l = 6 and l = 0, respectively: fitting WIT better leads to further forgetting (cf. Table 2). Unlike what is suggested by related work ( §2), we find that the linear fusion model performs on par with the more early, cross-attention based, fusion model (f vs. d in Table 1). This suggests that the improvement over the late fusion baseline indeed comes from the Multimodal ICT pre-training, which is not very sensitive to the model's architecture. Interestingly, cross-input cross-modal interactions (as shown in Figure 3) are possible with both ECA and ILF. So they may be the primary reason for performance improvement. Moreover, the architecture of ILF allows to fully freeze BERT during Stage 2, which circumvents catastrophic forgetting 14 . We leave other training strategies (e.g. multi-tasking, using adapters [START_REF] Houlsby | Parameter-Efficient Transfer Learning for NLP[END_REF]) for future work. Nothing suggests that ECA is better on the 14% of ViQuAE images that overlap with WIT. ECA is better on the out-of-WIT subset (38.0 vs. 36.5 MRR), but it is the other way around for DPR and late fusion.

Reading Comprehension

To extract the answers from the retrieved passages, we keep the same model as [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF]. It uses the Multi-passage BERT architecture [START_REF] Wang | Multimodal Data Enhanced Representation Learning for Knowledge Graphs[END_REF] and is thus based on text only because once the relevant passage has been retrieved, the question may be answered without looking at the image. To limit the variations due to training and the number of experiments, we use the model trained by [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF] off-the-shelf and simply change its input passages. It takes the top-24 passages as input. The model was first trained on TriviaQA (filtered of all questions used in ViQuAE), then fine-tuned on ViQuAE, much like stages 1 and 3. The authors provide five different versions of the model that correspond to different random seeds.

We use Exact Match and F1-score (at the bag-of-words level) to evaluate the extracted answers. In Table 3 we can verify that more relevant passages indeed lead to better downstream answers. The only difference with the IR evaluation is the role of the freezing technique of [START_REF] Aytar | Cross-modal scene networks[END_REF] (d vs. e), which is less clear here.

Generic vs. Specialized Image Representations

Numbers reported in the previous section are actually on par with the best results of [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF]. This is because the latter is based on ArcFace and ImageNet-ResNet, in addition to DPR and CLIP. In particular, [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF] have a heuristic for taking advantage of the face representations provided by ArcFace: they use ArcFace if faces are detected and a combination of CLIP and ImageNet-ResNet otherwise. They show that this method improves retrieval precision for questions about persons (for which face representations are relevant). However, this approach is not scalable for two reasons: (i) there are near 1,000 different entity types in ViQuAE (according to Wikidata's ontology), and not all can benefit from specialized representations; (ii) combining several representations (e.g. CLIP and ImageNet-ResNet) for the same entity type is computationally expensive and quickly saturates. To provide a comparable system to the late fusion of [START_REF] Lerner | ViQuAE, a dataset for knowledge-based visual question answering about named entities[END_REF], we have tried integrating ArcFace and ImageNet-ResNet in ECA. However, we have failed to outperform the CLIP-only version of ECA. Intuitively, we think that ECA dilutes the specialized representations of ArcFace and is unable to preserve them throughout all twelve layers of BERT. Therefore, in this setting, ECA is overall on par with late fusion (37.7 vs. 37.9 MRR, not significant) but better on questions about non-persons (39.3 vs. 35.7 MRR), which again suggests that it is unable to exploit ArcFace's representations.

Conclusion and Perspectives

We have presented a new pre-training method, Multimodal Inverse Cloze Task, for Knowledge-based Visual Question Answering about Named Entities. Multimodal ICT leverages contextual images in multimodal documents to generate visual pseudo-questions. It enables the use of more complex multimodal fusion models than previously proposed late fusion methods. Consequently, our method improves retrieval accuracy over the latter by 10% relative-MRR, leading to a 9% relative-F1 improvement in downstream reading comprehension (i.e. answer extraction), on the recently introduced ViQuAE dataset. We believe it is thanks to cross-modal interactions, which are prohibited by late fusion. More precisely, we qualitatively observed that cross-input interactions occurred between the image of the visual question and the text of the KB, which counteracts the heterogeneity of visual depictions.

We have experimented our pre-training method with two different neural networks architectures: (i) ECA, which follows recently proposed Multimodal BERTs by fusing modalities Early via Cross-Attention; (ii) ILF, a more standard model that fuses modalities through a linear projection. We found that both perform equally well, unlike in standard VQA and cross-modal retrieval. We argue that it might be because cross-input cross-modal interactions are the most important and may be captured by both models. However, further investigations are required, it might also be because of their different training settings, which leads ECA to catastrophic forgetting.

While aiming for generic multimodal representations of named entities, we found that integrating specialized representations in our models, such as ArcFace for faces, was not beneficial. We hypothesize that the studied architectures may be inappropriate but we leave this issue for future studies.

For future work, we think that generalizing Multimodal ICT for re-ranking (processing (t q , i q ) and (t p , i p ) simultaneously) and reading comprehension (generating or extracting the answer from (t p , i p )) is an exciting research lead. Indeed, there is evidence that sharing the same model for IR and reading comprehension, or IR and re-ranking, is beneficial for textual QA [START_REF] Fun | Efficient Retrieval Optimized Multi-task Learning[END_REF] and cross-modal retrieval [START_REF] Geigle | Retrieve Fast, Rerank Smart: Cooperative and Joint Approaches for Improved Cross-Modal Retrieval[END_REF], respectively: two tasks that closely relate to KVQAE.

  Those embeddings, noted [t; e c ], are then fed to the Transformer model, where the attention mechanism should

Fig. 2 :

 2 Fig. 2: Overview of Multimodal Inverse Cloze Task via Wikipedia/WIT.

1 "

 1 In which English palace was this man born?" "Who designed this cathedral?" He was appointed [...] Surveyor of the Fabric of St Paul's Cathedral, and Windsor Castle. St James's Palace was retained for official use.

Fig. 3 :

 3 Fig. 3: Qualitative examples where ECA (l = 6) finds a relevant visual passage in top-1 but late fusion falls behind (column on the right). We can see evidence of cross-input cross-modal interactions with ECA.

Table 3 :

 3 Reading Comprehension evaluation on ViQuAE, averaged over 5 runs of the reader. l: Number of frozen layers during Multimodal ICT.

	# IR Model	Multimodal ICT Exact Match	F1
	a DPR	NA	16.9 ± 0.4 20.1 ± 0.5
	b DPR + CLIP	NA	19.0 ± 0.4 22.3 ± 0.4
	c ECA (l = NA)		17.7 ± 0.6 21.2 ± 0.8
	d ECA (l = 6)		20.6 ± 0.3 24.4 ± 0.2
	e ECA (l = 0)		20.8 ± 0.8 24.3 ± 0.9
	f ILF (l = 12)		21.3 ± 0.6 25.4 ± 0.3

[START_REF] Lee | Latent Retrieval for Weakly Supervised Open Domain Question Answering[END_REF] use a batch size of over

4K questions.

Standard (Knowledge-based) VQA is often treated as a classification task.

https://www.wikidata.org/

With a ResNet-50 backbone[START_REF] He | Deep residual learning for image recognition[END_REF].

Uncased "base" 12-layers version available at https://huggingface.co.

Note that this is equivalent to concatenating both before projecting like [BERT(t) [CLS] ; CLIP(i)] • [Wt; Wc].

https://commons.wikimedia.org/

https://github.com/PaulLerner/ViQuAE

Jean Zay GPUs consume 0.482kW (or 0.259kW after heat recovery) in France, which has an average grid emission factor of 0.0569 kgCO2e/kWh according to https://bilans-ges.ademe.fr/en.

After standard preprocessing (lowercasing, stripping articles, and punctuation).

ILF only achieves 87.1 in-batch MRR on WIT's test set because of the freezing.

This work was supported by the ANR-19-CE23-0028 MEERQAT project. This work was granted access to the HPC resources of IDRIS under the allocation 2021-AD011012846 made by GENCI. We thank the anonymous reviewers for their helpful comments.