Ariane Pinche

teaching and research institutions in France or abroad, or from public or private research centers.

CREMMA MEDIEVAL DATASET

Medieval manuscripts corpus for training HTR models from the 12 th to the 15 th century, included in HTR-United catalog: https://htr-united.github.io Dataset has been built with eScriptorium, an interface for HTR ground truth production, and Kraken, an HTR and layout segmentation engine. It is composed of forteen Old French manuscripts written between the 13 th and 15 th centuries.

Transcriptions have been standardized to strengthen HTR models. We chose a graphemic transcription method.

Transcription guidelines avalaible here: https://hal.archives-ouvertes.fr/hal-03697382 MANUSCRIPT DATE TRANSCRIBED LINES BnF, ms fr. 412 13th 6324 BnF, Arsenal 3516 13th 1991 Cologny, bodmer, 168 13th 1976 BnF, ms fr. 24428 13th 1328 BnF, ms fr. 25516 13th 717 BnF, ms fr. 844 13th 224 BnF, ms fr. 17229 13th 164 BnF, ms fr. 13496 13th 161 BnF, Arsenal 3516 13th 105 BnF, ms fr. 22549 14th 2682 Vaticane, Reg. Lat., 1616 14th 1772 Univ. of Pennsylvania, 660 14th 368 BnF, ms fr. 411 14th 179 Univ. of Pennsylvania, 909 15th 2513 All 21656 To ensure data quality, continuous integration workflow has been put in place: • HTRVX: XML schema validator checking Alto files and segmentation vocabulary (segmOnto); • Choco-Mufin: checking the homogeneity of the characters used in the dataset.

Abstract: https://hal.archives-ouvertes.fr/hal-

03719504 HTR MODELS FOR MEDIEVAL MANUSCRIPTS TRAINED WITH KRAKEN https://kraken.re Test scores in-domain 1.0.1 Bicerin, accuracy 89.19%, model based on CREMMA Medieval first dataset, specialized on 13 th and 14t h c. manuscripts. 1.1.0 Bicerin, accuracy 95.30%, model based on CREMMA Medieval extended to 15 th c. manuscripts. 2.0.0 Cortado, accuracy 95.54%, model mixing CREMMA Medieval dataset with early prints (15 th c.) from Gallic(orpor)a Project. Models avalaible here: https://github.com/HTR-United/cremma-medieval/ releases Test scores out-of-domain

		BnF, ms,	BnF, ms,	BnF, NAF
		fr. 17229,	fr. 185,	6213,	ALL
		13th c.	14th c.	15th c.
	Cortado	92.71%	92.07%	87.48%	90.95%
	1.1.0 Bicerin	91.64%	91.34%	83.40%	89.23%

1.0.1 Bicerin 90.66% 88.45% 79.67% 86.50% Example, Cortado model on BnF, NAF 6213 Conclusion: a specialized model per script isn't always necessary, but the variety of the training set increases its robustness, even in our case with early prints. See complete experience here: https://cremmalab.hypotheses.org/modeles-htr CREMMALAB Project Handwritten Text Recognition (HTR) for medieval manuscripts Within the infrastructure of the CREMMA project (Consortium for Handwriting Recognition of Ancient Materials) supported by the DIM (research funded by the Île-de-France Region) MAP (Ancient and Heritage Materials), the CREMMALab project (2021-2022) combines research questions, creation and release of data from French medieval literary manuscripts for HTR. Project Manager: Ariane PINCHE (postdoctoral researcher) EVENTS Conference Ancient documents and automatic recognition of handwriting, 23-24 June 2022 WEBSITES https://cremmalab.hypotheses.org https://github.com/HTR-United/ cremma-medieval


Ibrahim Almuslimani

Philippe Chartier

Mohammed Lemou

Florian Méhats

UNIFORMLY ACCURATE SCHEMES FOR DRIFT-OSCILLATORY STOCHASTIC DIFFERENTIAL EQUATIONS

Keywords: highly-oscillatory, stochastic differential equations, micro-macro decomposition, uniform accuracy. AMS subject classification (2010): 65L20, 74Q10, 35K15

come

Uniformly accurate schemes for drift-oscillatory stochastic differential equations

Introduction

In this paper, we aim at constructing uniformly accurate numerical schemes for solving Itô stochastic differential equations (SDEs) with a (possibly highly) oscillatory drift term of the form dX(t) =f t/ (X(t))dt + σ(X(t))dW (t), X(0) = X 0 ∈ R d , (

where X(t) is a stochastic process with values in R d . A standard assumption in the literature of oscillatory problems and averaging theory is that the drift function (θ, x) ∈ T×R d → f θ (x) is assumed to be periodic 1 with respect to θ (we shall denote accordingly the torus by T = R/Z ≡ [0, 1]), see e.g. [START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF][START_REF] Chartier | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF][START_REF] Castella | Stroboscopic averaging for the nonlinear schrödinger equation[END_REF]. For simplicity and ease of the presentation, and without loss of generality, we assume in the rest of the paper that (θ, x)

∈ T × R d → f θ (x) is 1-periodic.
The diffusion function σ is defined as a smooth function x ∈ R d → σ(x) ∈ R d×m . Finally, W (t) = (W 1 (t), ..., W m (t)) T is an array of m independent one-dimensional Weiner processes. Precise regularity assumptions on f and σ are made in Assumption 1.2. Let us emphasize that is here a parameter whose value can freely vary in the interval ]0, 1] and that equation (1.1) is not restricted to its asymptotic regime where tends to zero , i.e., the oscillations may a Univ Rennes, INRIA, IRMAR -UMR 6625, F-35000 Rennes, France. Ibrahim.Almuslimani@irisa.fr b Univ Rennes, INRIA, CNRS, IRMAR -UMR 6625, F-35000 Rennes, France. Philippe.Chartier@inria.fr c Univ Rennes, CNRS, IRMAR -UMR 6625, F-35000 Rennes, France. Mohammed.Lemou@univ-rennes1.fr d Univ Rennes, CNRS, IRMAR -UMR 6625, F-35000 Rennes, France. Florian.Mehats@univ-rennes1.fr be slow or fast without any restriction on their frequency. This model may appear in many interesting applications including the perturbation of deterministic oscillatory problems by adding a random noise term for modeling purposes. For example, in Section 4, we consider a stochastic Hénon-Heiles model, one might also think of stochastic Fermi-Pasta-Ulam problem and the stochastic nonlinear Schrödinger equation [START_REF] Vilmart | Weak second order multirevolution composition methods for highly oscillatory stochastic differential equations with additive or multiplicative noise[END_REF].

Analogously to the case of deterministic differential equations, in the highly-oscillatory (stiff) regime [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF], standard numerical methods for SDEs, such as Euler-Maruyama method, face a severe time step restriction (∆t = O( )) when applied directly to (1.1). This issue is well documented in the literature for deterministic ODEs [START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF], and several classes of methods were introduced in order to deal with stiffness [START_REF] Bao | Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF][START_REF] Castella | Stroboscopic averaging for the nonlinear schrödinger equation[END_REF][START_REF] Chartier | Multi-revolution composition methods for highlyoscillatory differential equations[END_REF][START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory kinetic equations[END_REF][START_REF] Faou | Asymptotic preserving schemes for the klein-gordon equation in the nonrelativistic limit regime[END_REF]. However, none of the methods introduced therein qualify as uniformly accurate methods as they do not produce numerical approximations with an accuracy and at a cost both independent of the value of ∈]0, 1]. This motivated the introduction of a new methodology based on averaging techniques and exposed in [START_REF] Chartier | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF]: there, the authors elaborate a new technique enabling standard numerical methods to retain their non-stiff order with uniform accuracy for all ∈]0, 1]. In the case of stiff SDEs, as for ODEs, we can differentiate between two kinds of stiffness: stiff dissipative SDEs and highly-oscillatory SDEs. While for stiff dissipative SDEs, many interesting integrators were introduced in the last two decades with nice stability and convergence properties [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF][START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF][START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF], the numerical solution of highly-oscillatory (drift) SDEs has not received so much attention except, for example, [START_REF] Vilmart | Weak second order multirevolution composition methods for highly oscillatory stochastic differential equations with additive or multiplicative noise[END_REF] in which the author derives weak second order multirevolution composition methods which are accurate only for very small values of and fail for 1. This contribution is, up to our knowledge, the first attempt to adapt the technique of micro-macro decomposition to the SDE context in order to construct a uniformly accurate integrator that works well for every ∈]0, 1].

More precisely, we focus , in the present paper, on constructing uniformly accurate methods for highly-oscillatory SDEs in the spirit of the methodology explained in [START_REF] Chartier | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF]. We derive a micro-macro system by introducing a change of variable that leads us to treat the average decay and the fast oscillations separately. We show that applying Euler-Maruyama method to the micro-macro system gives, under appropriate assumptions, an approximation of uniform weak order 1 and strong order 1/2 for any value of ∈ (0, 1] and with no restriction on the time step. In more mathematical terms, we prove that, under Assumption 1.2 below, the Euler-Maruyama scheme for solving the micro-macro system derived from (1.1), provides approximations X n on a uniform grid {t n } N n=0 such that

∀n = 0, . . . , N, |E(φ(X(t n ))) -E(φ(X n ))| ≤ C h and E(|X(t n ) -X n | 2 ) 1 2 ≤ Ch 1 2 ,
where h = t n+1 -t n (for simplicity we consider constant time step size) and the constant C is independent of n, h and . We prove as well that the same result can be obtained using the integral scheme introduced in Section 2. After Noticing the simplicity of implementation and the uniform accuracy of the integral scheme, a natural question arises: why do we consider the micro-macro decomposition? The answer is, interestingly, deriving the micro-macro method gives insight about possible generalizations to higher-order schemes inspired from deterministic averaging theory developed in [START_REF] Chartier | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF].

The rest of the paper is organized as follows. In the next subsection, we fix the main notations used throughout the paper and we make assumptions. In Section 2 we introduce the integral scheme and we prove its uniform convergence (weak and strong) with respect to the parameter . In Section 3 we derive the micro-macro scheme, and we show its uniform convergence properties. Finally, in Section 4 we present some numerical experiments that illustrate the efficiency of the schemes.

Notations and Assumptions. Noteworthy, our results are obtained under quite standard assumptions that we recall below. In particular, we shall constantly suppose that Assumption 1.2 below holds true in the sequel. Definition 1.1. Let q ∈ N, we define C q poly as the set of q-times differentiable functions from R d → R whose all derivatives up to order q have at most polynomial growth. We define as well C q lin in the same way but with linear growth instead of polynomial. Assumption 1.2. The functions x → σ i,j (x), i = 1, . . . , d, j = 1, . . . , m are Lipschitz functions of the set C 2 poly . The function (θ, x) → f θ (x) ∈ R d is defined on T × R d and it is uniformly Lipschitz continuous with respect to x (i.e with Lipschitz constant independent of θ). In addition, its components (f θ (•)) i , i = 1, . . . , d, are functions of C 4 poly ∩ C 2 lin . Furthermore, f is continuously differentiable with respect to θ and x → ∂ θ f θ (x) is uniformly Lipschitz continuous with respect to the variable x.

Remark 1.3. The Lipschitz continuity of f θ and σ ensures the existence and uniqueness of the solution of (1.1) in C([0, T ]; L 2 (Ω)). It ensures in addition a linear growth, that is to say

∀(θ, x) ∈ T × R d , |f θ (x)| + |σ(x)| ≤ K(1 + |x|),
where K is a constant independent of θ.

Notations. The derivative with respect to the space variable x will be denoted by a prime. For example, f θ (x) = ∂ x f θ (x) (and gradient or Jacobien when needed), and f θ (x) = ∂ 2

x f θ (x) (Hessian matrix or second derivative tensor when needed). The letter L will be used throughout the paper as the maximum Lipschitz constant of the concerned functions. The constants represented by the capital letter C (indexed or not) are generic constants. We denote by | • | the Euclidean norm of R d .

We will conduct all our analysis for m = 1 in (1.1). This comes for the sake of simplicity and ease of presentation, and to avoid any confusion due to the heavy computations done in the proofs. However, we emphasize that our results apply straightforwardly to SDEs of the form (1.1) with any number of Weiner processes. Hence, σ will denote a function from R d → R d and W (t) will denote a standard one-dimensional Weiner process.

Integral scheme

Let T > 0 be a given final time, N ∈ N\{0}, and h = T /N , and let t n = nh, n = 0, . . . , N . We consider the integral scheme

X n+1 = X n + h f (X n ) + t n+1 / tn/ (f θ (X n ) -f (X n ))dθ + σ(X n )∆W n , (2.1)
where f (x) is the value at x of the average of f defined as

f (x) = T f θ (x)dθ,
and

∆W n = W (t n+1 ) -W (t n ) ∼ N (0, h)
which can be replaced by √ hξ n where ξ n ∼ N (0, 1). Note that the integral term in (2.1) can be rewritten as

(F t n+1 / (X n ) -F tn/ (X n )), where F θ (x) = θ 0 (f τ (x) -f (x)) dτ.
Theorem 2.1. Consider the integral scheme (2.1), where t n = nh for n = 0, . . . , N , and N ∈ N\{0}. Let X(t) be the solution of (1.1). Then, for all φ ∈ C 4 poly , there exists C > 0 independent of h, n and ∈ (0, 1] such that,

∀n ∈ {0, . . . , N }, |E(φ(X(t n ))) -E(φ(X n ))| ≤ Ch, (2.2) ∀n ∈ {0, . . . , N }, E(|X(t n ) -X n | 2 ) 1 2 ≤ Ch 1 2 . (2.3) Lemma 2.2.
Under the assumptions of Theorem 2.1, one has for all n = 0, . . . , N -1

E t n+1 tn [f t/ (X(t)) -f t/ (X(t n ))]dt ≤ Ch 2 . (2.4)
Proof. We first rewrite the integrand

f t/ (X(t)) -f t/ (X(t n )) of (2.4) as t tn d(f t/ (X(s))) = t tn [f t/ (X(s))f t/ (X(s)) + 1 2 f t/ (X(s))(σ(X(s)), σ(X(s)))]ds + t tn f t/ (X(s))σ(X(s))dW (s),
so that

E t n+1 tn [f t/ (X(t)) -f t/ (X(t n ))]dt ≤ E t n+1 tn t tn f t/ (X(s))f t/ (X(s)) + 1 2 f t/ (X(s))(σ(X(s)), σ(X(s))) dsdt + E t n+1 tn t tn f t/ (X(s))σ(X(s))dW (s)dt ,
and using Fubini's theorem, we get the following upper-bound of the left-hand side of (2.4)

t n+1 tn E t tn f t/ (X(s))σ(X(s))dW (s) dt + O(h 2 ),
where O(h 2 ) is bounded thanks to the polynomial growth assumptions and the boundedness of the moments. The expectation in previous integral is an expectation of a stochastic integral, hence it is null since X(s) is independent of the increments of W for times above s.

Lemma 2.3. Under the assumptions of Theorem 2.1, we have for all n = 0, . . . , N -1

E(|X(t) -X(t n )| 2 ; X(t n ) = x) ≤ C(x)(t -t n ).
Proof. From the integral form of equation (1.1)

X(t) = X(t n ) + t tn f s/ (X(s)ds + t tn σ(X(s))dW (s).
we may write

E(|X(t) -X(t n )| 2 ; X(t n ) = x) = E t tn f s/ (X(s))ds + t tn σ(X(s))dW (s) 2 ≤ 2(t -t n ) t tn E(f s/ (X(s)) 2 )ds + 2 t tn E(σ(X(s)) 2 )ds ≤ 2(t -t n ) t tn K(1 + E(|X(t n )| 2 ))ds + 2 t tn K(1 + E(|X(t n )| 2 ))ds ≤ C(x)(t -t n ),
where we have used, between the first and second lines, the Young's inequality, then the Cauchy-Schwarz inequality for the first integral, and Itô isometry for the second integral. Finally, we have used the bound

E(|X(s)| 2 ) ≤ K(1 + |x| 2 )
, borrowed from [START_REF] Gīhman | Stochastic differential equations[END_REF], between the third and the fourth lines.

Proof of Theorem 2.1. We hereby follow the methodology introduced in [START_REF] Talay | Discrétisation d'une équation différentielle stochastique et calcul approché d'espérances de fonctionnelles de la solution[END_REF], which consists in bounding the moments (first step), then proving the weak convergence of the scheme (second step) and finally establishing its strong convergence (third step).

Step 1. In order to bound the moments of X n of arbitrary order, we shall resort to [18, Lemma 2.2, p. 102], which requires the following estimates

|E(X n+1 -X n |X n = x)| = E t n+1 / tn/ f θ (x)dθ ≤ t n+1 / tn/ K(1 + |x|)dθ = K(1 + |x|)h,
and

|X n+1 -X n | = t n+1 / tn/ f θ (X n )dθ + σ(X n )∆W n ≤ K(1 + |X n |)h + K(1 + |X n |) |∆W n | √ h √ h ≤ M n (1 + |X n |) √ h,
where

M n = K √ h + |∆Wn| √ h
which is of bounded moments since ∆Wn √ h ∼ N (0, 1).

Step 2. Let φ ∈ C 4 poly be a test function. For X(t n ) = X n = x, by performing Taylor expansion of φ with integral remainder, it can be shown that

|E(φ(X(t n+1 )) -φ(X n+1 ))| ≤ E φ (x) t n+1 tn (f t/ (X(t)) -f t/ (X(t n )))dt + E φ (x) t n+1 tn (σ(X(t)) -σ(X(t n )))dW (t) t n+1 tn (σ(X(t)) + σ(X(t n )))dW (t) + 1 6 E 1 0 (1 -τ ) 3 φ (4) (x + τ (X(t n+1 ) -x))dτ (X(t n+1 ) -x) 4 + 1 6 E 1 0 (1 -τ ) 3 φ (4) (x + τ (X n+1 -x))dτ (X n+1 -x) 4 + Ch 2 := I+II+III+IV + Ch 2 ,
where the term Ch 2 comes from the remaining expectations of the second and third derivatives of φ applied to the integrals, which are zero for odd number of stochastic integrals, and bounded by Ch 2 otherwise. The constant C is independent of thanks to the boundedness of the moments of the exact and numerical solutions, and to the regularity assumptions made in Assumption 1.2. By Lemma 2.2, Lemma 2.3 and the boundedness of the moments, the terms I, III and IV are O(h 2 ). As for the integrand of the second term, we have using Itô formula

σ(X(t)) -σ(X(t n )) = t tn σ (X(s))f s/ (X(s)) + 1 2 σ (X(s))(σ(X(s)), σ(X(s))) ds + t tn σ (X(s))σ(X(s))dW (s).
Hence, using the Lipschitz-continuity of σ, the polynomial growth of φ , Lemma 2.3, and the following consequence of Itô isometry

∀g, h ∈ L 2 ad ([a, b] × Ω), E b a g(t)dW (t) b a h(t)dW (t) = E b a g(t)h(t)dt
we have

II ≤ C(x) E t n+1 tn (σ(X(t)) -σ(X(t n )))dW (t) t n+1 tn (σ(X(t)) -σ(X(t n )))dW (t) +2 t n+1 tn σ(X(t n ))dW (t) = C(x) E t n+1 tn (σ(X(t)) -σ(X(t n ))) 2 dt + 2C(x) E t n+1 tn σ(X(t n ))(σ(X(t)) -σ(X(t n )))dt ≤ C(x)L t n+1 tn C(t -t n )dt + 2 E σ(X(t n )) t n+1 tn t tn dσ(X(s))dt = 2 E σ(X(t n )) t n+1 tn t tn σ (X s )f s/ (X s ) + σ (X s )(σ(X s ), σ(X s )) dsdt + 2 E σ(X(t n )) t n+1 tn t tn σ (X s )σ(X s )dW (s)dt + C 1 h 2 .
The first expectation after the last equal sign is clearly O(h 2 ). The second one is equal to zero since X(t n ) is independent of the increments W (r) -W (s) for t n ≤ s ≤ r. The boundedness of the moments and the local weak order 2 imply the global weak convergence of order 1 by a theorem from [START_REF] Milstein | Numerical integration of stochastic differential equations[END_REF] (see also [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Chapter 2.2]).

Step 3. We first derive an upper-bound of E(|X(t n+1 ) -X n+1 | 2 ; X(t n ) = x) as follows

E t n+1 tn (f t/ (X(t)) -f t/ (x))dt + t n+1 tn (σ(X(t)) -σ(x))dW (t) 2 ≤ 2E t n+1 tn (f t/ (X(t)) -f t/ (x))dt 2 + 2E t n+1 tn (σ(X(t)) -σ(x))dW (t) 2 ≤ 2hE t n+1 tn f t/ (X(t)) -f t/ (X(t n )) 2 dt + 2E t n+1 tn σ(X(t)) -σ(X(t n )) 2 dt ≤ 2L 2 (h + 1) t n+1 tn E(|X(t) -X(t n )| 2 )dt ≤ 2L 2 (h + 1) t n+1 tn C(t -t n )dt ≤ Ch 2 ,
where we have used Lemma (2.3). Finally, we conclude that

E(|X(t n+1 ) -X n+1 | 2 ) 1 2 ≤ Ch.
A well-known theorem by Milstein [START_REF] Mil Shteȋn | A theorem on the order of convergence of mean-square approximations of solutions of systems of stochastic differential equations[END_REF] then allows to establish inequality (2.3).

Remark that for the integral scheme we can relax the assumption (f θ (.)

) i , i = 1, . . . , d ∈ C 4 poly ∩ C 2 lin to (f θ (.)) i , i = 1, . . . , d ∈ C 2 poly .

Micro-Macro method

The integral scheme, despite its uniform accuracy, does not generalize to higher order methods. For this reason, and inspired by [START_REF] Castella | Stroboscopic averaging for the nonlinear schrödinger equation[END_REF][START_REF] Chartier | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF], we will introduce and analyze a micro-macro separation of scales that leads the standard Euler-Maruyama method to retain its usual weak and strong orders of convergence. We believe that this will be an excellent starting point and a good insight to try to develop higher order micro-macro methods in future works.

As in general averaging theory, the purpose is to to find a periodic, near-identity and smooth change of variable Φ θ , together with a flow Ψ t , the flow map of an autonomous nonstiff differential equation on R d , such that the solution of the original equation (1.1) takes the composed form [START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF][START_REF] Lochak | Multiphase Averaging for Classical Systems[END_REF][START_REF] Castella | Stroboscopic averaging for the nonlinear schrödinger equation[END_REF][START_REF] Chartier | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF]]

X(t) = Φ t/ • Ψ t (X 0 ). (3.1)
In the current work, we restrict ourselves to first order averaging, and we show that the first order change of variable Φ θ defined, for all 0 < ≤ 1, by the formula

Φ θ (x) = x + F θ (x) = x + θ 0 (f τ (x) -f (x)) dτ, (3.2)
and the first order flow Ψ t satisfying the autonomous equation

dΨ t (x) = f (Ψ t (x))dt + σ(Ψ t (x))dW (t), Ψ 0 (x) = x, (3.3)
both borrowed from deterministic averaging, work well to construct a uniform accurate micromacro scheme of weak order 1 and strong order 1/2 for (1.1). For the derivation of (3.2) and (3.3) from (3.1), we refer to [START_REF] Chartier | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF]. Let us denote by X the solution of (3.3), separating slow and fast scales as follows

X = Φ t/ (X) + Y, (3.4)
leads, using Itô formula, to the micro-macro system of the form

dX = f (X)dt + σ(X)dW, X(0) = X 0 , (3.5
)

dY = f t/ Φ t/ (X) + Y -f t/ (X) -F t/ (X) f (X) (3.6) - 2 F t/ (X) (σ(X), σ(X)) dt + σ(Φ t/ (X) + Y ) -σ(X) -F t/ (X)σ(X) dW, Y (0) = 0.
Here Y represents the averaging error that is characterized by the equation (3.6) and added after applying the change of variable in order to recover the solution X(t) of (1.1). We will show later that the expected value of Y (t) is of size (see Lemma 3.3), which is analogous to the averaging error for deterministic oscillatory equations (σ ≡ 0).

Remark 3.1. The existence and the uniqueness of the solution Y (t) of (3.6) is an immediate consequence of the existence and the uniqueness of the solutions X(t) and X(t) of (1.1) and (3.5), and the decomposition (3.4).

We use the same uniform discretization as for the integral scheme. Our aim is now to prove a uniform (in ) convergence result for the following micro-macro scheme, which is nothing but the Euler-Maruyama method applied to (3.5, 3.6)

X n+1 = X n + h f (X n ) + √ hσ(X n )ξ n , X 0 = X 0 , (3.7)
Y n+1 = Y n + h f t/ Φ tn/ (X n ) + Y -f tn/ (X n ) (3.8) -h F tn/ (X n ) f (X n ) + 1 2 F tn/ (X n ) (σ(X n ), σ(X n )) + √ h σ(Φ tn/ (X n ) + Y n ) -σ(X n ) -F tn/ (X n )σ(X n ) ξ n , Y 0 = 0,
where the increment ξ n is a random quantity sampled from a normalized Gaussian centered at zero and with variance 1.

Main result.

Theorem 3.2. Consider the Euler-Maruyama scheme (3.7, 3.8) for solving the micro-macro system (3.5, 3.6) and let

X n = Φ tn/ (X n ) + Y n
for n = 0, . . . , N . Let X(t) be the solution of (1.1). Then, for all φ ∈ C 4 poly , there exists C > 0 independent of h, n, and , such that, one has

∀n ∈ {0, . . . , N }, |E(φ(X(t n ))) -E(φ(X n ))| ≤ C h, (3.9) ∀n ∈ {0, . . . , N }, E(|X(t n ) -X n | 2 ) 1 2 ≤ Ch 1 2 . (3.10)
Again, the proof of the theorem follows the usual steps from [START_REF] Talay | Discrétisation d'une équation différentielle stochastique et calcul approché d'espérances de fonctionnelles de la solution[END_REF]. Note that the main novelty of our result lies in the fact that estimate (3.9) is uniform w.r.t.

∈ (0, 1]. The following lemma will be needed in the proof of Theorem 3.2.

Lemma 3.3.

There exists c > 0 such that for all ∈ (0, 1] and all t ∈ [0, T ], one has

E(|Y (t)|) ≤ c .
Proof. First, note that we have for all t ∈ [0, T ] and all positive integer m (see [START_REF] Gīhman | Stochastic differential equations[END_REF])

E(|X(t)| 2m ) ≤ K(1 + |X 0 | 2m ).
Moreover, we have

E(|Y (t)| 2 ) ≤ 2tE t 0 f s/ (Φ s/ (X(s)) + Y (s)) -f s/ (X(s)) -F s/ (X(s)) f (X(s)) + 1 2 F s/ (X(s))(σ(X(s)), σ(X(s))) 2 ds + 2E t 0 σ(Φ s/ (X(s)) + Y (s)) -σ(X(s)) -F s/ (X(s))σ(X(s)) 2 ds ≤ 4t t 0 E f s/ (Φ s/ (X(s)) + Y (s)) -f s/ (X(s)) 2 ds + 4t 2 t 0 E F s/ (X(s)) f (X(s)) + 1 2 F s/ (X(s))(σ(X(s)), σ(X(s))) 2 ds + 4 t 0 E σ(Φ s/ (X(s)) + Y (s)) -σ(X(s)) 2 ds + 4 2 t 0 E F s/ (X(s))σ(X(s)) 2 ds := 4tA + 4t 2 B + 4C + 4 2 D.
In the first inequality, we have used the Cauchy-Schwarz inequality for the first term and Itô isometry for the second one. We have also used the triangular inequality as well as Young's inequality. Given that E(|X(s)| 2m ) < +∞ and that the components of F θ are in C 2 poly by Assumption 1.2, we can conclude that the terms B and D are uniformly bounded. Now, by definition of the change of variables (3.2) and the uniform Lipschitz continuity of f θ and σ, we get

4tA + 4C ≤ 4tL 2 t 0 E(| F s/ (X(s)) + Y (s)| 2 )ds + 4L 2 t 0 E(| F s/ (X(s)) + Y (s)| 2 )ds ≤ 8L 2 (T + 1) 2 t 0 E(|F s/ (X(s))| 2 )ds + 8L 2 (T + 1) t 0 E(|Y (s)| 2 )ds ≤ C 1 2 + C 2 t 0 E(|Y (s)| 2 )ds.
Therefore, E(|Y (t)| 2 ) satisfies the following inequality

E(|Y (t)| 2 ) ≤ C 3 2 + C 4 t 0 E(|Y (s)| 2 )ds. Gronwall's lemma imply that E(|Y (t)| 2 ) ≤ c 2 2 . Finally, E(|Y (t)|) 2 ≤ E(|Y (t)| 2 ) = c 2 2 , thus E(|Y (t)|) ≤ c
where c is a positive constant independent of and t.

Proof of Theorem 3.2. We begin by the proof of the weak convergence: Collecting the variables X and Y into Z = (X, Y ) and similarly z = (x, y), we may rewrite the micro-macro system (3.5, 3.6) as dZ = g t/ (Z)dt + Σ t/ (Z)dW where

g θ (z) = f (x) f θ (Φ θ (x) + y) -f θ (x) -F θ (x) f (x) -2 F θ (x) (σ(x), σ(x)) and Σ θ (z) = σ(x) σ(Φ θ (x) + y) -σ(x) -F θ (x)σ(x)
and accordingly equations (3.7), (3.8) as

Z n+1 = Z n + hg tn/ (Z n ) + √ hΣ tn/ (Z n )ξ n
where the ξ n ∼ N (0, I m ) are independent Gaussian random variables.

Step 1. Since the vector fields g t/ and Σ t/ have uniform linear growth in z (this follows from the linear growth of f , σ, Φ, F , F , and F ), the proof of the boundedness of the moments of the numerical solution can be obtained following the same arguments as in Milstein's lemma for the non-oscillating case ( = 1). As a matter of fact, we have that

(i) E(|Z 0 | 2m ) < +∞ since Z 0 = (X 0 , 0) ∈ R 2d (deterministic); (ii) |E(Z n+1 -Z n | Z n = z)| = h|g tn/ (z)| ≤ K(1 + |z|)h; (iii) |Z n+1 -Z n | ≤ h|g tn/ (Z n )| + |Σ tn/ (Z n )||∆W n | ≤ M n (1 + |Z n |) √ h, with ∆W n = W (t n+1 ) -W (t n ) ∼ N (0, h) and M n = K| ∆Wn √ h | has clearly bounded moments uni- formly in n.
We have used the fact that g and Σ grow linearly in Z thanks to Assumption 1.2. Under the above conditions, Lemma 2.2 from [18, p. 102] implies the boundedness of the moments of arbitrary order of Z n .

Step 2. We define, for a function ψ that depends explicitly on θ and z,

(L θ ψ)(z) = ∂ θ ψ(θ, z) + ∂ z ψ(θ, z)g θ (z) + 1 2 ∂ 2 z ψ(θ, z) (Σ θ (z), Σ θ (z)) . (3.11)
We have for any test function φ ∈ C 4 poly (independent of theta)

E(φ(Z n+1 )|Z n = z) ≤ φ(z) + h(L tn/ φ)(z) + Ch 2 + 1 6 E 1 0 (1 -τ ) 3 φ (4) (z + τ (Z n+1 -z))dτ (Z n+1 -z) 4 Z n = z .
where, as in the proof of Theorem 2.1, the constant C is independent of and where the term Ch 2 comes from the remaining expectations of the second and third derivatives of φ applied repeatedly to hg + √ hΣξ n , which are zero for odd moments of ξ n , and bounded by Ch 2 otherwise. For the part of the remainder coming from the above first order Taylor expansion of the expectation of φ(Z n+1 ), it is clearly bounded by Ch 2 , where the constant C is independent of . This stems from the polynomial growth of the test function φ and its derivatives up to order 4, as well as from the bound

|Z n+1 -Z n | ≤ M n (1 + |Z n |) √ h.
Performing the Taylor expansion of the expectation of the test function φ applied to the exact solution leads to

E(φ(Z(t n+1 ))|Z n = z) = φ(z) + h(L tn/ φ)(z) + R,
where, according to [17, P. 26

] 2 R = E t n+1 tn s tn (L 2 τ / φ)(Z(τ ))dτ ds|Z(t n ) = z .
2 There are three other terms with vanishing expectations.

UNIFORMLY ACCURATE SCHEMES FOR DRIFT-OSCILLATORY STOCHASTIC DIFFERENTIAL EQUATIONS 11

Let us notice that the derivative with respect to θ will not appear in L tn/ φ because our test function φ does not depend explicitly on time. However, it will appear in L 2 tn/ φ owing to the explicit dependence of g t/ and Σ t/ on t, and only because of this derivation with respect to θ, we will have terms of order O( 1 ) that need to be bounded.

Let g 1 θ and g 2 θ be the two components of g θ (both belong to R d ), then

∂ θ g θ (z) = (∂ θ g 1 θ (z), ∂ θ g 2 θ (z)) T = (0, ∂ θ g 2 θ (z)) T , with ∂ θ (g 2 θ (z)) = ∂ θ (g 2 θ (x, y)) = 1 ∂ θ f t/ (Φ t/ (x) + y) + 1 f t/ (Φ t/ (x) + y)∂ θ Φ t/ (x) - 1 ∂ θ f t/ (x) -∂ θ F t/ (x) f (x) - 1 2 ∂ θ F t/ (x)(σ(x), σ(x)).
Using the same notations for Σ θ (z), we have

∂ θ (Σ 2 t/ (z)) = 1 σ (Φ t/ (x) + y)∂ θ Φ t/ (x) -∂ θ F t/ (x)σ(x).
It remains to bound the difference

1 ∂ θ f τ / (Φ τ / (X(τ )) + Y (τ )) - 1 ∂ θ f τ / (X(τ )), (3.12)
since all other terms are bounded independently of thanks to the smoothnes of f θ , σ and their derivatives with respect to the space variable, and using the fact that ∂ θ Φ τ / (X) = O( ). Now, taking into account the uniform Lipschitz continuity of ∂ θ f θ , the polynomial growth of the test function φ, and the above lemma, we have

E 1 ∂ θ f τ / (Φ τ / (X(τ )) + Y (τ )) - 1 ∂ θ f τ / (X(τ )) ≤ L E Φ τ / (X(τ )) + Y (τ ) -X(τ ) ≤ L E F τ / (X(τ )) + L E |Y (τ )| = LE F τ / (X(τ )) + Lc ≤ L(K(1 + E |X(τ )|) + c) < +∞
independently of the value of . Hence, we have

E t n+1 tn s tn (L 2 τ / φ)(Z(τ ))dτ ds|Z(t n ) = z ≤ C h 2 ,
and the local order 2 is proved.

Step 3. The boundedness of the moments and the local weak order 2 imply the global weak convergence of order 1 by a theorem from [START_REF] Milstein | Numerical integration of stochastic differential equations[END_REF] (see also [18, chap2.2]).

We have proved, for any test function φ ∈ C 4 poly (R 2d , R), for all n = 0, 1, . . . , N , that

|E(φ(Z(t n )) -φ(Z n ))| ≤ Ch,
where Z = (X, Y ) and C is independent of n and . We need to prove that for any test function φ ∈ C 4 poly (R d , R), for all n = 0, 1, . . . , N ,

|E(φ(X(t n )) -φ(X n ))| = E(φ(Φ tn/ (X(t n )) + Y (t n )) -φ(Φ tn/ (X n ) + Y n )) ≤ C 1 h.
where C 1 is independent of n and . Let φ ∈ C 4 poly (R d , R), for each fixed parameter θ ∈ T, we consider the test function

ψ θ (Z) = φ(Φ θ (X) + Y ). Note that ψ θ ∈ C 4 poly (R d , R), since Φ θ ∈ C 4 poly (R d , R d ) (
by assumption on f θ ). Hence, for each n = 0, 1, . . . , N , we have

E(ψ tn/ (Z(t n )) -ψ tn/ (Z n )) = E(φ(Φ tn/ (X(t n )) + Y (t n )) -φ(Φ tn/ (X n ) + Y n )) = |E(φ(X(t n )) -φ(X n ))| ≤ C 1 h, where C 1 = C.
Strong convergence. For the strong convergence with order 1/2 (for Z), in addition to the bounded moments and the local weak order 2, using [START_REF] Mil Shteȋn | A theorem on the order of convergence of mean-square approximations of solutions of systems of stochastic differential equations[END_REF], it is sufficient to show first order strong convergence after one step, i.e,

E(|Z(t n+1 ) -Z n+1 | 2 ; Z(t n ) = z) 1 2 ≤ Ch, (3.13)
where the generic constant C is independent of n, h, and . By Wagner-Platen expansion [17, chap1.2.2], we have

Z(t n+1 ) -Z n+1 = t n+1 tn s tn Λ τ / Σ τ / (Z(τ ))dW (τ )dW (s) + t n+1 tn s tn L τ / Σ τ / (Z(τ ))dτ dW (s) + t n+1 tn s tn Λ τ / g τ / (Z(τ ))dW (τ )d(s) + t n+1 tn s tn L τ / g τ / (Z(τ ))dτ d(s)
where Λϕ θ (z) = ϕ θ (z)Σ θ (z), L θ is defined in (3.11) and

(Λ θ ψ)(z) = ψ (z)Σ θ (z).
The only terms we need to bound are the derivatives with respect to θ of g θ (Z) and Σ θ (Z) at θ = τ / arising from the application of the differential operator L τ / to g τ / (Z) and Σ τ / (Z) respectively. It can be checked that the other terms in the integrals are again uniformly bounded by the regularity assumptions, Lemma 3.3, and the boundedness of the moments of Z(t). Once this is done, we can conclude that (3.13) is satisfied. The uniform boundedness of ∂ θ g τ / (Z) was already proved in the weak convergence case. Now, we have

∂ θ Σ τ / (Z) = (0, ∂ θ Σ 2 θ (Z)), and
∂ θ Σ 2 τ / (Z) = 1 σ (Φ τ / (X) + Y )∂ θ Φ τ / (X) -∂ θ F τ / (X)σ(X)
The right hand side is uniformly bounded since ∂ θ Φ θ (x) = O( ) (thanks to the periodicity of f θ with respect to θ and the definition of f ), σ (x) and ∂ θ F θ (x) have at most polynomial growth, and σ has linear growth. It follows that, after several applications of the Itô isometry, the Cauchy-Schwarz and the Young inequalities,

E(|Z(t n+1 ) -Z n+1 | 2 ) ≤ Ch 2 ,
and thus E(|Z(t n+1 ) -Z n+1 | 2 ) 1 2 ≤ Ch. Now we conclude the local strong order 1 for X. First, note that

E(|Z(t n+1 ) -Z n+1 | 2 ) = E(|X(t n+1 ) -X n+1 | 2 ) + E(|Y (t n+1 ) -Y n+1 | 2 ) ≤ Ch 2 , Now, we consider the SDE dX = f t/ (X)dt + σ(X)dW (t). (4.1)
In all our experiments = 2 -2i , i = 2, 3, 4, 5, final time T = 1; M denotes the number of computed samples and ∆t denotes the time step size.

4.1. Weak convergence. In this section we use X 0 = (0.7, 0.7, 0.7, 0.7), ∆t = 2 -i , i = 1, . . . , 5, M = 10 4 .

4.1.1. Multiplicative noise. We consider the above SDE (4.1) with multiplicative noise where σ(X) = 0.2(0, 0, X 1 , X 2 ) T . We use the test function φ(X) = X 1 to measure the weak convergence. In Figure 1A we plot the weak error with respect to the time step for different values of (left figure) using the micro-macro method (3.7)- (3.8). We can see that the convergence behavior looks almost the same, with weak order one, for all the different values of . The right picture of Figure 1A shows that for a fixed time step, the weak error remains almost constant when varying . The above description applies also to the integral scheme (see Figure 1B).

4.1.2. Additive noise. The weak error of the micro-macro scheme (3.7)-(3.8) applied to the SDE (4.1) with additive noise is shown in Figure 2. We set σ(X) = (0, 0, 0.2, 0.2) T , and we perform the test with two different test functions. We see again the uniform weak order one.

Strong convergence.

In this section we use X 0 = (0.12, 0.12, 0.12, 0.12), ∆t = 2 -i , i = 4, . . . , 8, M = 10 2 .

4.2.1. Multiplicative noise. We consider the above SDE (4.1) with multiplicative noise where σ(X) = 0.5X. Figure 3 shows the uniform strong order 1 2 for both methods.

4.2.2. Additive noise. We consider the above SDE (4.1) with additive noise where σ(X) = (0, 0, 0.5, 0.5) T . In addition to the uniform convergence, Figure 4 shows strong order one for the micro-macro method (3.7)-(3.8) since when the noise is additive, Euler-Maruyama method coincides with Milstein method of strong order one. This applies to uniformly accurate methods too.

4.3.

Inefficiency of Euler-Maruyama method for particular time steps. Although Euler-Maruyama method seems to work quite well, it still fails for some particular choices of time steps, while the Micro-Macro method (3.7)-(3.8) does not. We recall that the importance of uniformly accurate methods appears more when using higher order schemes. See Figure 5. We consider the logistic SDE dX = X(1 -X) + sin t dt + 0.2XdW (t), X(0) = 2.

We set = 0.1 and we plot in Figure 5A the reference solution (in blue) calculated with very small time step using the integral scheme and the solution obtained using EM with time step ∆t = 0.99(2π ). In Figure 5B, we plot the reference solution (in blue) calculated with very small time step using the integral scheme and the solution obtained using the uniformly accurate method (3.7)-(3.8) with time step h = 0.99(2π ).

Acknowledgement. The authors would like to thank Gilles Vilmart for useful discussions and comments.

The work of I. Almuslimani is supported by the Swiss National Science Foundation, project No: P2GEP2_195212.

which implies that each component of Z n+1 = (X n+1 , Y n+1 ) converges strongly with local order 1 to the corresponding component of Z(t n+1 ) = (X(t n+1 ), Y (t n+1 )). Next, for a given X(t n ) = X n , we have

Therefore,

and the local strong order 1 is proved.

Remark 3.4. When applying the Euler-Maruyama method directly to (1.1), a term that contains 1 ∂ θ f t/ (X(t n )) will appear in the remainder. This term causes the local truncation error to be of order O( 1 ) which causes instability in the high frequency regime ( 1) even when Assumption 1.2 is satisfied. In contrast, applying Euler Maruyama method to the micro-macro system (3.5,3.6) leads instead of 1 ∂ θ f t/ (X(t n )) to the term (3.12) which can be bounded independently of using Assumption 1.2.

Numerical experiments

Throughout this section (except for the last experiment) we will focus on the Hénon-Heiles model (see [START_REF] Burrage | Structure-preserving Runge-Kutta methods for stochastic Hamiltonian equations with additive noise[END_REF][START_REF] Han | High-order stochastic symplectic partitioned Runge-Kutta methods for stochastic Hamiltonian systems with additive noise[END_REF]). We consider the Hamiltonian

It can be checked that the variable X(t) satisfies the following ODE [START_REF] Chartier | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF] dX dt

Conclusion

In this work, we have introduced two uniformly accurate methods for solving numerically stochastic differential equations with oscillatory drift. The first one is the so-called integral scheme (2.1) and can be derived quite straightforwardly, whereas the second one is obtained through a more elaborate transformation, namely a micro-macro decomposition (3.5, 3.6). Both schemes exhibit weak-order 1 and strong-order 1/2, as proved in the corresponding sections and confirmed numerically in Section 4. Given their comparable performance, the first scheme is arguably better for its simplicity. However, it is our belief that the micromacro scheme exposed here could be generalized to higher order methods which would be the stochastic counterpart of existing deterministic uniformly accurate methods [START_REF] Baumstark | Uniformly accurate exponential-type integrators for Klein-Gordon equations with asymptotic convergence to the classical NLS splitting[END_REF][START_REF] Chartier | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF].


domain: spi phys

text:
Maria-Elisavet Chatzitheodoridi
email: maria-elisavet.chatzitheodoridi@onera.fr

Abigael Taylor
email: abigael.taylor@onera.fr

Olivier Rabaste
email: olivier.rabaste@onera.fr

Hélène Oriot
email: helene.oriot@onera.fr

Fast Computation of the Optimal Integrated Sidelobe Level Filter for Large-Scaled Problems

Keywords: Mismatched Filter, Integrated Sidelobe Level, Continuous Phase Frequency Shift Keying, Time Complexity, Fast Algorithm, SARComm System, Waveform Diversity

The growth of spectral demand due to congestion in the electromagnetic spectrum is skyrocketing through the past years. Joint radar-communications systems are thus attracting more attention as they can share their bandwidth to alleviate the frequency spectrum. In such a system, communication codes are used to both transmit information and perform radar missions. The latter notably requires mitigating the sidelobe level energy, which is generally achieved via the use of a mismatched filter. However, such a filter may be costly to optimize, especially for largescale problems. Hence, this paper proposes a fast algorithm for computing the optimal integrated sidelobe level mismatched filter, based on the Lagrangian dual function of the primal problem. Once the Lagrangian function is calculated, the optimum values of its variables are computed, after formulating it as a monodimensional problem, in order to provide the desired filter. The resulting computational time shows a dramatic improvement on the execution time of the proposed method compared to the classical Matlab convex solver CVX, especially for large-scale problems.

I. Introduction

The continuous growth of communication applications gradually congests the electromagnetic spectrum. Concurrently, the need of alleviating the bandwidth occupation, adds up to the spectrum congestion problem [START_REF] Griffiths | Challenge problems in spectrum engineering and waveform diversity[END_REF]. To deal with this problem, it is imperative to combine two different applications, such as radars and communications systems. One of the possible ways of coupling them, without creating interference between them, is to design a joint radar communications framework that uses the same signal and allows simultaneous functioning at the same frequency bandwidth [START_REF] Hayvaci | Spectrum sharing in radar and wireless communication systems: A review[END_REF], [START_REF] Tedesso | Code shift keying based joint radar and communications for emcon applications[END_REF].

The radar application considered in this paper is a synthetic aperture radar (SAR) framework. In such context, the waveform is usually a linear frequency modulated pulse also called chirp [START_REF] Cumming | Digital Processing of Synthetic Aperture Radar Data: Algorithms and Implementation[END_REF]. However, such signals are not able to transmit binary information and thus cannot be used to create a combined SAR-communications system. To do so, communications codes are used, more precisely in our case we chose continuous phase frequency-shift keying (CPFSK) signals [START_REF] Anderson | Digital phase modulation[END_REF], in order to simultaneously generate the radar image and transmit binary information. Moreover, in SAR framework, target detection is essential and it is usually performed by the use of matched filtering, that aims at maximizing the target signal-to-noise ratio (SNR) [START_REF] Levanon | Radar signals[END_REF], [START_REF] Wang | Digital Signal Processing Techniques and Applications in Radar Image Processing[END_REF]. However, potentially strong sidelobes are created while using the matched filter processing in the presence of strong targets, which may prevent the detection of weaker targets [START_REF] Blunt | Overview of radar waveform diversity[END_REF] and potentially deteriorate the image contrast. Thus, a different compression filter can be used in order to mitigate the sidelobe level, which is called mismatched filter (MMF) [START_REF] Blunt | Overview of radar waveform diversity[END_REF]. An extensive literature has already been devoted to the subject of mismatched filter formulation and optimization [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF]- [START_REF] Chatzitheodoridi | A sidelobe level minimization mismatched filter using continuous phase frequency-shift keying codes for the off-grid delay problem[END_REF].

For sampled signals, apart from the particular case of spectral analysis and linear frequency modulation where the pulse compression can be carried out via weighting windows [START_REF] Harris | On the use of windows for harmonic analysis with the discrete fourier transform[END_REF], mismatched filtering is usually expressed as a convex problem under convex constraints. The criterion used in the cost function is usually either the Peak-to-Sidelobe Ratio (PSLR) or the Integrated Sidelobe Level (ISL) [START_REF] Stoica | Transmit codes and receive filters for radar[END_REF]- [START_REF] Baden | Optimal peak sidelobe filters for biphase pulse compression[END_REF], [START_REF] Maio | Design of radar receive filters optimized according to Lp-norm based criteria[END_REF]- [START_REF] Aittomäki | Mismatched filter design for radar waveforms by semidefinite relaxation[END_REF], and the convex constraints are typically considered to control the loss-in-processing gain or the mainlobe shape [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF].

In this paper, we consider the ISL criterion in order to create a joint SAR-communications system for which the sidelobe level should be controlled. Such a system needs to satisfy essential constraints for both applications. More precisely, the signals are over-sampled for the SAR image generation, and the transmitting codes are different at each acquisition in order to maximize the data rate. In such configuration, where the transmitting signals differ from pulse to pulse, the mismatched filter needs to be computed for every transmitted code. Optimizing a large number of over-sampled waveforms leads to solving a large-scale problem of high computational cost.

In many papers that deal with sidelobe suppression [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF], [START_REF] Maio | Design of radar receive filters optimized according to Lp-norm based criteria[END_REF], [START_REF] Jiu | Knowledge-based spatial-temporal hierarchical mimo radar waveform design method for target detection in heterogeneous clutter zone[END_REF], [START_REF] Jia | A novel method to suppress short-range clutter in airborne radar[END_REF], the Matlab toolbox CVX : Matlab Software for Disciplined Convex Programming [START_REF] Grant | CVX: Matlab Software for Disciplined Convex Programming, version 2.1[END_REF] is used to solve proteiform convex optimization problems. However, as mentioned in [START_REF] Grant | CVX: Matlab Software for Disciplined Convex Programming, version 2.1[END_REF], this solver is not recommended for largesized vector variables due to its computational cost. For instance, the problem considered in this paper introduces several thousands of pulses, each one of them containing some thousands of sampled points, and thus the overall computational cost for the generation of all the optimal filters was estimated around 11 years! It is thus essential to provide a different way to solve the convex optimization problem in order to reduce the computational time.

We show in this paper that the original optimization problem with convex cost function and convex constraints can be efficiently solved via the associated Lagrangian dual problem, in the the case of minimizing the ISL criterion. It can be mostly solved analytically, by converting the dual problem, containing several thousands of variables to a simple mono-dimensional optimization problem, which can be easily solved via a classic Newton-Raphson algorithm.

It is shown in the results section that the computational time of the proposed algorithm is dramatically improved, even just for only one long-size filter computation. Thus for some thousand filters that need to be computed for the SAR acquisitions, the gain in computational time is extraordinary. It was then possible to compute a resynthesized SAR image obtained from several thousand of different CPFSK signals and the associated optimal mismatched filters.

This paper is organized in 5 sections. The introduction of the paper is presented in section I and illustrates the state-of-the art, and the paper description and motivation. In section II, the definition of the convex mismatched filter optimization problem is provided, in order to allow in section III to introduce a different way to implement this problem by solving it mathematically using the Lagrangian dual function. Section IV presents the computational performance of the proposed algorithm and the classical Matlab Toolbox CVX, usually used for solving convex problems, and also a generated SAR image using different signal and by extension different mismatched filters. Finally, section V provides a conclusion of the paper.

II. Mismatched Filter Definition

Let s be a N -length sampled signal expressed as

s = s 1 s 2 . . . s N T , (1)
where . T is the transpose operator. At reception, the signal is processed by a mismatched filter in order to constrain the sidelobe energy. The length K of the mismatched filter q ∈ C K can be equal to or greater than the length of the signal (K ≥ N ). The output signal can be expressed as [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF]

y = Λ * K (s) q, (2)
where . * is the complex conjugate,

q = q 1 q 2 . . . q K T ,
and where Λ K is the correlation matrix of s of size (K + N -1) × K.

If the filter is chosen as q = s, then the classical matched filter is obtained, that maximizes the SNR at the cost of potentially high sidelobes. But it can be chosen so as to optimize various criteria, dependent on the considered application.

As the SAR application is considered in this paper, the main criterion of interest is the Integrated Sidelobe Level (ISL), due to numerous scatterers present in the scene, that can be added up. When referring to the optimization problem formulation using the ISL criterion, extensive literature is available [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF], [START_REF] Stoica | Transmit codes and receive filters for radar[END_REF]- [START_REF] Baden | Optimal peak sidelobe filters for biphase pulse compression[END_REF], [START_REF] Maio | Design of radar receive filters optimized according to Lp-norm based criteria[END_REF]- [START_REF] Aittomäki | Mismatched filter design for radar waveforms by semidefinite relaxation[END_REF], and the problem is expressed as a L 2 -norm minimization problem [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF]:

min q q H Mq, s.t. s H q = s H s, q H q ≤ 10 β 10 s H s, (3)
where M = Λ T K (s)FΛ * K (s) is a symmetric semi-definite positive matrix, F is a diagonal matrix of ones except for some zero values which correspond to the mainlobe position indices and β is a positive constant in dB that expresses the loss-in-processing gain. This value represents the maximum tolerated loss for the resulting filter with respect to the matched filter [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF].

As mentioned in [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF], the first constraint, which allows to discard the trivial null solution, is linear and thus convex, and the second one, that enables to control the loss-inprocessing gain, is also convex as it is a positive semidefinite quadratic constraint. Globally, the optimization problem is thus convex, as the L 2 -norm is convex, with convex constraints. As a consequence, any local optimum of the optimization problem is necessarily a global optimum.

It is important to note that if the second constraint on the loss-in-processing gain is removed from the problem (3), then the optimization problem can be analytically solved [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF]. Indeed, the solution can be obtained by applying the Lagrange multipliers method and it is expressed as:

q * analytic (s, M) = (s H s) M -1 s s H M -1 s .
where M = Λ T K (s)FΛ * K (s). This analytic solution provides the optimal sidelobe level but a worse loss-inprocessing gain, which is potentially unacceptable.

Considering that the problem ( 3) is convex, it can be solved using a convex solver, such as the CVX: Matlab Software for Disciplined Convex Programming [START_REF] Grant | CVX: Matlab Software for Disciplined Convex Programming, version 2.1[END_REF], as it is performed in [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF], [START_REF] Maio | Design of radar receive filters optimized according to Lp-norm based criteria[END_REF], [START_REF] Jiu | Knowledge-based spatial-temporal hierarchical mimo radar waveform design method for target detection in heterogeneous clutter zone[END_REF], [START_REF] Jia | A novel method to suppress short-range clutter in airborne radar[END_REF]. However, as mentioned in [START_REF] Grant | CVX: Matlab Software for Disciplined Convex Programming, version 2.1[END_REF], CVX is not suitable for large-scale problems due to high computational costs. Especially in the SAR-communications framework, it is infeasible to use this solver, as the need of using different and oversampled pulses creates a very large problem with very high computational time. Thus, it is necessary to bypass the use of CVX and manage to solve the convex optimization problem (3) via a much faster algorithm.

III. Mismatched Filter Optimization using the Dual Cost Function without the CVX Solver

In this section, a solution of the problem (3) is provided without the use of CVX solver. As mentioned, CVX solves generic convex problems and it is not recommended to deal with very large-scale problems, due to computational complexity. More precisely, the variable number, that need to be determined using CVX, is at least of the order of the filter size. As, the size of the filter is a large number in our framework, the computational complexity of the problem skyrockets. Thus, the Lagrangian function is calculated and the solution of the problem is provided by the dual problem.

The associated Lagrangian function combines the cost function and the weighted sum of the constraints of the primal problem [START_REF] Boyd | Convex optimization[END_REF]. It is expressed as

L(q, λ, ν) = q H Mq + λ(q H q -αs H s) + ν(s H q -s H s) = -λαs H s -νs H s independent of q + q H (M + λI)q + νs H q dependent on q . ( 4
)
where α = 10 β 10 , λ ≥ 0 and ν ∈ R. After computing the gradient with respect to q of the Lagrangian function gives,

∇ q L(q, λ, ν) = 2(M + λI)q + νs, (5)
the optimum filter is found by taking the null of the Lagrangian gradient,

∇ q L(q * , λ, ν) = 0 ⇐⇒ q * = - 1 2 ν(M + λI) -1 s. ( 6
)
The dual function of the problem is defined as,

g(λ, ν) = inf q L(q, λ, ν) = L(q * , λ, ν).
Finally, after injecting the solution (6) in the Lagrangian function ( 4), the dual function can be written as

g(λ, ν) = -λαs H s -νs H s - 1 4 ν 2 s H (M + λI) -1 s. (7)
The optimal pair (λ, ν) needs to be found in order to deduce the optimal filter for the dual problem from the equation [START_REF] Levanon | Radar signals[END_REF]. Since the dual function g is concave [START_REF] Boyd | Convex optimization[END_REF], its optimum can be easily found. More precisely, we have,

∂g ∂ν (λ, ν) = -s H s - 1 2 νs H (M + λI) -1 s, ∂g ∂ν (λ, ν) = 0 ⇐⇒ ν * = - 2s H s s H (M + λI) -1 s . ( 8
)
After determining the optimal value of ν, the equation ( 8) is injected in the equation [START_REF] Wang | Digital Signal Processing Techniques and Applications in Radar Image Processing[END_REF], in order to provide a function g that depends only on the λ variable. We led thus our problem to have to solve numerically a monodimensional optimization problem of a concave function.

More precisely, let us denote the mono-dimensional function,

g(λ) = max ν g(λ, ν) = -λαs H s + (s H s) 2 s H (M + λI) -1 s (9)
which is a concave function.

Optimizing g(λ) cannot be performed analytically due to its structure, but the optimal value λ * can be found by using the Newton-Raphson method.

The Newton-Raphson method iterative formula is [START_REF] Ortega | Iterative solution of nonlinear equations in several variables[END_REF]

λ k+1 = λ k - g (λ k ) g (λ k ) . ( 10
)
For the sake of simplicity, let us denote,

h(λ) = s H (M + λI) -1 s (11)
The first and second derivatives of the function g, expressed in [START_REF] Rabaste | Mismatched filter optimization for radar applications using quadratically constrained quadratic programs[END_REF], are computed with respect to the λ variable in order to specify and perform the iterative algorithm using the formula [START_REF] Blunt | Doppler-compensated adaptive pulse compression[END_REF],

g (λ) = -αs H s -(s H s) 2 h (λ) h 2 (λ)
, and

g (λ) = -(s H s) 2 • h (λ)h(λ) -2[h(λ)] 2 h 3 (λ) .
The derivatives of g are dependent on function h, expressed in [START_REF] Stoica | Transmit codes and receive filters for radar[END_REF], and its first and second derivatives. This means that for every iteration, the function h and its associated derivatives need to be computed, that leads to inverse matrix computations, at each iteration. In order to avoid inverting matrices at each step of the gradient algorithm, which are more costly from a computational point of view, an eigenvalue decomposition of M is performed. More precisely, M can be written as M = P H DP, with P a square matrix, D a diagonal one and the conjugate transpose is used because M is symmetric. Then, the function h can be simplified as,

h(λ) = s H P H (D + λI) -1 Ps = s H P H diag 1 d + λ Ps ( 12
)
where d is the vector that determines the diagonal elements of the matrix D, and now the inverse matrix term can be easily calculated [START_REF] Mirsky | An introduction to linear algebra[END_REF].

Similarly, the first and the second derivatives of the h function from the equation ( 12) can be expressed after performing the eigenvalue decomposition as

h (λ) = s H P H diag - 1 (d + λ) 2 Ps, h (λ) = s H P H diag 2 (d + λ) 3 Ps.
We can remark that the advantage of these expressions is that the square matrix P and the diagonal matrix D can be computed only once and for all for the optimization problem and the function h expressed in [START_REF] Cilliers | Pulse compression sidelobe reduction by minimization of Lp-norms[END_REF] and its derivatives are finally computed only by performing multiplications and no matrix inversions.

In order to ensure that the optimal pair (λ * , ν * ) of the dual problem coincides with the optimal solution of the primal convex problem defined in (3), strong duality should hold. It consists on proving that the duality gap is equal to zero, which means that the difference between the solution of the primal problem (3) and the solution of the dual problem ( 4) is zero.

Strong duality can be ensured when the primal problem is convex and Slater's constraint qualification is verified [START_REF] Boyd | Convex optimization[END_REF]. The first condition is already justified in Section II and the second qualification is proved in the following.

Lemma III.1 (Slater's condition). [23]

There exists an x 0 ∈ relint (D), where relint (D) the relative interior of D, and D the set of feasible solutions of the primal problem, that verifies the following equality and strict inequality of the problem,

s H x 0 = s H s, (i) x H 0 x 0 < αs H s. (ii)
This lemma III.1 is verified for x 0 = s. Indeed, when x 0 = s, the equality (i) is true. Besides, the strict inequality (ii) clearly holds whenever α > 1, which is always the case since α = 1 corresponds to no loss in processing gain, which means the matched filter solution.

Algorithm 1 illustrates the procedure to obtain one optimal mismatched filter, by providing the optimal pair (λ * , ν * ). It is repeated as many times as the number of transmitted signals.

IV. Simulation and Results

This section is dedicated to the framework explanation and the comparison of the computational cost between the proposed algorithm and the CVX solver. Then, an illustration of a SAR image is provided, which makes use of multiple mismatched filters computed with the Algorithm 1.

A. Computational cost between the proposed method and the CVX solver

The algorithm proposed in this paper can be applied to any sampled signal as long as we consider the integrated sidelobe level criterion. In order to allow simultaneously the generation of a SAR image and the transmission of binary information, communication codes are used. More precisely, continuous phase frequency-shift keying (CPFSK) codes [START_REF] Anderson | Coded modulation systems[END_REF], a family of continuous phase modulated codes, are chosen for their attractive properties, such as constant envelop signals, continuous phase and well-contained spectral energy [START_REF] Tierney | Adaptive waveform design for sar in a crowded spectrum[END_REF].

Algorithm 1: Mismatched Filter Optimization using Gradient Descent Input: Loss-in-progressing gain α, initialization ρ, max iterations max i , precision , Output: Optimal pair (λ * , ν * ), optimal mismatched filter q * ; s ← code generation;

M ← Λ T K (s)FΛ * K (s); [P, D] ← SVD(M); λ 1 ← ρ; for i ∈ [1, max i ] do h i ← s H P H diag 1 diag(D) + λ i Ps; h i ← s H P H diag - 1 (diag(D) + λ i ) 2 Ps; h i ← s H P H diag 2 (diag(D) + λ i ) 3 Ps; g i ← -αs H s -(s H s) 2 h i h 2 i ; g i ← -(s H s) 2 • h i h i -2h 2 i h 3 i ; λ i+1 ← λ i - g i g i ; if |g i | < then break; end end ν * ← - 2s H s s H (M + λ * I) -1 s ; q * ← - 1 2 ν * (M + λ * I) -1 s.
Table I lists some signal characteristics and algorithm parameters used for the SAR image generation while using CPFSK codes. Algorithm 1, displayed in section III, provides the procedure to compute one optimal mismatched filter. The number of different signals, transmitted by the airborne platform, indicates the repetition number of the algorithm execution. This number is determined by the pulse repetition frequency and the total integration time of the platform, which in our framework corresponds to as much as 6522 transmitted signals, as indicated in Table I. For each filter computation, initialization paramters are given to the algorithm, also indicated in Table I. More precisely, for the initialization of the λ variable, we fix λ 1 = 0.5, a maximum iteration parameter is also needed and set to 100 iterations for which the algorithm stops and provides the computed λ, and a precision parameter is set to 1e-10, to ensure the numerical convergence of the algorithm. Moreover, each one of these signals is composed of 3000 samples and the length of the associated mismatched filter is set to three times longer, which leads to 9000 samples. The choice of longer filters provides significantly lower sidelobe level energy than considering Table II: Computational cost of the proposed algorithm and the CVX solver for one optimal mismatched filter.

same signal length filters, as shown in [START_REF] Chatzitheodoridi | A mismatched filter for integrated sidelobe level minimization over a continuous Doppler shift interval[END_REF], [START_REF] Chatzitheodoridi | A sidelobe level minimization mismatched filter using continuous phase frequency-shift keying codes for the off-grid delay problem[END_REF]. Table II shows the computational time of the Algorithm 1 and of the Matlab CVX solver for one and for all the generated mismatched filters. Their computational time difference is enormous. The computational cost of one filter using our method is estimated to be equal to 3.48 minutes, after averaging over 30 different randomly drawn sequences, and it is compared with the 14.5 hours spent on average by the CVX algorithm on the same sequences. Extrapolated to the 6522 different sequences required for the formation of the SAR image, this computational time would be equal to 11 years, as shown in Table II. On the contrary, the proposed solution obtains all 6522 mismatched filters in less than 16 days.

B. SAR Image Generation from Multiple Transmitted CPFSK Signals

This section is dedicated to the creation and representation of SAR images and their associated impulse response, generated with different transmitted CPFSK signals and with either matched or mismatched filters. The proposed Algorithm 1 is implemented in Matlab and it is used in order to compute all 6522 mismatched filters, as indicated in Table I.

In order to compare the performance of the matched and mismatched filters with CPFSK codes on realistic images, we worked here with a real data image acquired with chirp signals. From this real SAR image, the scene reflectivity is obtained by deconvolving the image with the impulse response of the chirp signal. Then it is possible to convolve this scene reflectivity with each of the CPFSK signals, and perform the range compression either with the matched or the mismatched filter in order to obtain synthesized outputs.

In Figure 1, we first represent the impulse responses of chirp and CPFSK codes with the matched or the mismatched filter on a single point scatterer, and then the corresponding re-synthesized SAR images using the SAR image reconstruction algorithm.

Figure 1a represents the impulse response and the SAR image of the chirp and it is used as a comparison reference for the rest of the resulting images. It provides a benchmark in terms of image quality, but as it was obtained using chirp signals, no information transmission is possible. Figures 1b &1c illustrate the impulse responses and the SAR images of the CPFSK codes either using matched or mismatched filtering, respectively. Ideally, we want to provide a performance that converges to the chirp results, but at the same time provide the ability of transmitting information. When comparing the impulse responses from sub-figures 1b & 1c, the one generated with the mismatched filters provides an overall lower and more diffused sidelobe energy compared to the one generated using matched filtering. Due to the large number of transmitted signals, the sidelobe energy is spread in the "2D plan" and not only on the range and azimuth direction of the corner reflector, as seen at the chirp impulse response in sub-figure 1a. This effect may reduce high sidelobes over a direction but can potentially create more blurred zones on the image. The SAR image generated with CPFSK codes and mismatched filters, in sub-figure 1c, is close in contrast to the SAR image generated with the chirp signal, illustrated in sub-figure 1a. Still the contrast in the chirp image is better and the characteristics in azimuth direction are sharper. When comparing the SAR images generated with CPFSK codes and matched filters, subfigure 1b and mismatched filters, sub-figure 1c, we can remark that the energy level of sub-figure 1b is higher.This effect can potentially cover areas with not a lot of signal, or weak targets if the level of the sidelobes in the image increases.

V. Conclusion

Our paper proposes a fast algorithm for the computation of the mismatched filter which minimizes the integrated sidelobe level. It is based on optimization of the dual problem, which is shown to resort mainly to the optimization of a simple concave mono-dimensional problem that can be easily solved via a classical Newton-Raphson algorithm. This method provides the optimal mismatched filter, with a constraint on the loss in processing gain, and it is sufficiently fast to enable the computation of largesized vectors, which is infeasible while using other solvers, such as the convex CVX solver, due to computational time limits.

Table I :

 I Signal Characteristics and Algorithm Parameters.

	Parameters	Value
	Bandwidth B	305 MHz
	Pulse duration Tp	4.8 µs
	Modulation index h	0.4
	Chip duration T	19.2 ns
	Sampling frequency Ts	625 MHz
	Constellation	16
	Signal Length	3000
	Allowed loss in processing gain β	2 dB
	Initialization ρ	0.5
	Maximum iterations max i	100
	Precision		1e-10
	Pulse repetition frequency	5000 Hz
	Number of transmitted signals	6522
	Number of filters	Proposed Algorithm	CVX	Ratio
	1 6522	3.48 minutes 15.8 days	14.5 hours 11 years	249.7