Eduardo Hernández-Montero

Andrés Fraguela-Collar

Jacques Henry

AN OPTIMAL QUASI SOLUTION FOR THE CAUCHY PROBLEM FOR LAPLACE EQUATION IN THE FRAMEWORK OF INVERSE ECG

Keywords: Mathematics Subject Classification. 31B20/31A25/35R30/65J22/65N21. March 9, 2018 Cauchy problem, optimal regularization, quasi solution, factorization method, invariant embedding, ECG inverse problem

The inverse ECG problem is set as a boundary data completion for the Laplace equation: at each time the potential is measured on the torso and its normal derivative is null. One aims at reconstructing the potential on the heart. A new regularization scheme is applied to obtain an optimal regularization strategy for the boundary data completion problem. We consider the R n+1 domain Ω. The piecewise regular boundary of Ω is defined as the union ∂Ω = Γ1 ∪ Γ0 ∪ Σ, where Γ1 and Γ0 are disjoint, regular, and n-dimensional surfaces. Cauchy boundary data is given in Γ0, and null Dirichlet data in Σ, while no data is given in Γ1. This scheme is based on two concepts: admissible output data for an ill-posed inverse problem, and the conditionally well-posed approach of an inverse problem. An admissible data is the Cauchy data in Γ0 corresponding to an harmonic function in C 2 (Ω) ∩ H 1 (Ω). The methodology roughly consists of first characterizing the admissible Cauchy data, then finding the minimum distance projection in the L 2 -norm from the measured Cauchy data to the subset of admissible data characterized by given a priori information, and finally solving the Cauchy problem with the aforementioned projection instead of the original measurement.

Introduction

For a century, the main tool for cardiologists to assess the electric wave triggering the contraction of a patient's heart has been the ECG. It provides the effect of the heart potential on the potential at various points of the body surface with a good time resolution. Cardiologists are trained to interpret it, but the resulting information is qualitative. It is thought that building a potential map on the heart from the measurement of potential on the torso will allow a major improvement for the diagnosis. Such devices already exist but the stage of solving numerically a Cauchy problem for the Laplace equation remains the key point to improve.

Given a domain Ω and two separate parts Γ 0 and Γ 1 of its boundary ∂Ω, this problem, understood as the boundary value completion on Γ 1 , from given Cauchy data on Γ 0 , of an harmonic potential defined over Ω, has been widely studied because of its various medical and engineering applications [START_REF] Abda | Boundary data completion: The method of boundary value problem factorization[END_REF][START_REF] Azaez | On cauchy's problem: Ii. completion, regularization and approximation[END_REF][START_REF] Faker | On cauchy's problem: I. a variational steklov-poincaré theory[END_REF][START_REF] Cluitmans | Noninvasive reconstruction of cardiac electrical activity: Update on current methods, applications and challenges[END_REF]. Many examples in the literature show that this Cauchy problem can be approached as a linear inverse problem, severely ill-posed in the Hadamard sense [START_REF] Faker | Why is the cauchy problem severely ill-posed?[END_REF][START_REF] Hadamard | Lecture The Cauchy problem in linear partial differential equations[END_REF]. Improving accuracy of the regularized solution still is a relevant issue [START_REF] Bouyssier | Inverse problem in electrocardography via the factorization method of boundary value problems[END_REF], even though there are regularization strategies of optimal order and asymptotically optimal for this kind of problem in the framework of general theory [START_REF] Sergei | On the a posteriori choice of a regularization parameter in the solution of severely ill-posed problems[END_REF][START_REF] Solodky | Morozov's discrepancy principle for the tikhonov regularization of exponentially ill -posed problems[END_REF][START_REF] Solodky | On optimization of projection methods for solving some classes of severely ill-posed problems[END_REF][START_REF] Hohage | Regularization of exponentially ill-posed problems[END_REF]. Throughout these pages, a new methodology (or scheme) is applied to get an optimal approximation to the solution of the previous boundary data completion problem, referred to as the Admissible Data Methodology (AD). The Cauchy problem for Laplace equation that will be considered is as follows:

∆u ≡ 0 in Ω, ( 1
)
u| Σ = 0, (2)
u| Γ0 = φ, (3)
∂u ∂ν Γ0 = ψ, (4)
where Ω is a bounded domain in R n (n ≥ 2), ν is the unit outward normal to ∂Ω, ∂u ∂ν Γ0 the normal derivative of u over Γ 0 , and ∂Ω is composed of the smooth, and open surfaces Γ 1 , Γ 0 and Σ = ∂Ω \ (Γ 0 ∪ Γ 1 ) with Γ 0 ∩ Γ 1 = ∅, Σ a surface of strictly positive measure in ∂Ω. For us, problem (1)-( 4) will be the theoretical model describing some real phenomenon. The model clearly requires that the pair (φ, ψ) belong to a certain class M of functions where problem (1)-( 4) has a solution, and it makes sense to extend that solution over Γ 1 in the same sense in which boundary conditions over Γ 0 are given.

The methodology described in following pages provides a regularization strategy based on a model testing idea: if for model ( 1)-( 4) a given measured Cauchy data ( φ, ψ) is in a function class containing M , with noise level δ > 0, it must exist an element of M which is close enough to ( φ, ψ) ( of order δ); otherwise, there would be evidence pointing that the model does not describe the phenomenon properly. Based on this thought, a regularization scheme consists in regularizing the ill-posed problem of minimizing the distance from the measurement to the admissible data set M . For that purpose we find a subclass of M where projection is conditionally well posed, and we use that projection to solve problem (1)- [START_REF] Denisov | Elements of the Theory of Inverse Problems[END_REF].

In practice, the solution of the boundary data completion problem using the projection as the new Cauchy data must approximate well enough the solution of the real boundary data completion problem for real phenomenon, whether no evidence to reject the model is observed. This is a theoretical work, and so only the regularization strategy of the boundary data completion problem concerning the Cauchy problem (1)-( 4) is considered, no statistical approach toward a hypothesis test is considered.

In the sense of this modeling philosophy, an exact data is understood as the best theoretical modeling data; this way, (φ † , ψ † ) will denote the exact Cauchy data, and Φ † the exact solution of the boundary data completion related to (1)-( 4) and (φ † , ψ † ). Solution of (1)-(4) will be considered in a weak sense and recovering Dirichlet boundary data means finding a function Φ in the Sobolev space H 1 2 (Γ 1 ), and a harmonic function u in H 1 (Ω) that satisfy (2)-( 4) and

u| Γ1 = Φ, (5)
where Dirichlet data φ belongs to H 1 2 (Γ 0 ) and Neumann boundary condition ψ is taken in L 2 (Γ 0 ) as a priori information, although it is usually considered in the Sobolev space H -1 2 (Γ 0 ). Numerical examples are given in section 4 when Ω is a cylindrical domain.

A common way to set boundary completion problem is via the operational equation AΦ = ρ, where the right hand side is defined by ρ ≡ φ -Bψ, A being the Dirichlet to Dirichlet operator from Γ 1 to Γ 0 defined by auxiliary boundary problem given by (1), (2), (5) and a null Neumann condition over Γ 0 , and B is the Neumann to Dirichlet operator on Γ 0 defined by the auxiliary boundary problem given by ( 1), ( 2), (4) and a null Dirichlet condition on Γ 1 .

When both components of Cauchy data are given with noise of order δ at most, then, the continuous dependency of ρ allows to set the inverse problem as a first kind ill-posed one:

AΦ † = ρ † , ρ † = φ † -Bψ † φ † -φ ≤ δ ψ † -ψ ≤ δ,
Φ † being the exact solution for the unknown exact Cauchy data (φ † , ψ † ), and ( φ, ψ) the known perturbation of (φ † , ψ † ). In this approach, all methods use a priori information (some source condition or noise level) to regularize the Moore-Penrose pseudo-inverse of A i.e. to look for the Φ reg in the A domain minimizing the discrepancy AΦ -ρ . Hence, regularization strategies such as Tikhonov, Landweber, Truncated Singular Value Decomposition (TSVD), conjugate gradient or discretization methods are applied. All these strategies provide order optimal or asymptotically optimal solutions under the a priori information that the exact solution Φ † is bounded by a given constant in a stronger norm than the one of L 2 (Γ 1 ); however, they require a high smoothness condition over the exact solution in order to obtain optimal approximations in the case of severely ill-posed problems [START_REF] Kirsch | An introduction to the Mathematical Theory of Inverse Problems[END_REF].

What we do here is to focus on the admissible data, which will be understood as the Cauchy data such that the weak solution of problem ( 1)-( 4) exists. Set of admissible data is characterized via the operational equation defining the inverse problem, and given a priori information is used to continuously project in some sense the perturbed data over a subset of admissible data and in a such way that the inverse problem is conditionally well-posed over the aforementioned set. Our regularization scheme consists of two simple steps approximately described as follows:

(1) Finding (φ δ , ψ δ ), the minimum L 2 -distance projection from the perturbed Cauchy data ( φ, ψ) over the subset of admissible data defined by given a priori information such that the worst case error defined in section 1.3 of [START_REF] Kirsch | An introduction to the Mathematical Theory of Inverse Problems[END_REF] is bounded (the conditionally well-posed problem). ( 2) Solving the inverse problem with admissible data (φ δ , ψ δ ) instead of the first given Cauchy data ( φ, ψ).

In this work the presented methodology will be called the admissible data methodology (AD), and its optimality will be proved as well. In section 3 it is proved that the AD solution is a quasi-solution according with definition in [START_REF] Denisov | Elements of the Theory of Inverse Problems[END_REF], and that it does not require a so strong a priori information to be optimal, as occurs with Tikhonov regularization strategy, the most popular one, in the framework of the severely ill-posed problems [START_REF] Kirsch | An introduction to the Mathematical Theory of Inverse Problems[END_REF].

Through this scheme, considering a stronger a priori information than before, and applying the classic optimization theory, we show that in such cases there is α ≥ 0, depending on the stronger a priori information, such that Φ α defined by the following optimization problem is an optimal approximation of the exact solution for the inverse problem,

(Φ α , ψ α ) = arg min AΦ + Bψ -φ 2 + ψ -ψ 2 + α (Φ, ψ) 2 0 (Φ,ψ)
where A is the linear compact operator defining the inverse problem, B the Neumann to Dirichlet operator on Γ 0 for the auxiliary boundary data problem defined by ( 1)-( 2) and ( 4)-( 5) with Φ ≡ 0, and (•, •) 0 is a stronger norm than the H

1 2 (Γ 1 ) × L 2 (Γ 0 ) one.
A similar procedure based on the method of fundamental solutions (MFS) has been presented by T. Wei and Y.G. Chen in [START_REF] Wei | A regularization method for a cauchy problem of laplace's equation in an annular domain[END_REF]; the main difference is that their approach assumes that the exact solution is close enough to a finite dimensional subspace of H = {u ∈ H p (Ω) : ∆u ≡ 0}(p > 3/2), which actually imposes a high smoothness a priori condition over the exact solution.

The paper is organized as follows. Section 2 is devoted to formulating the inverse problem in the framework of this methodology. In section 3, the regularization strategy is formally presented, and optimal accuracy is proven. The link with the factorization method in [START_REF] Henry | Factorization of Boundary Value Problems Using the Invariant Embedding Method[END_REF] is shown in section 3 as well, suggesting discretization methods for the implementation of the AD solution in complex geometries of Ω. Finally, in section 4, the methodology is applied when Ω is a cylinder, and a numerical example is given, comparing with the Tikhonov regularization strategy, the most common method employed to solve linear ill-posed problems, which behaves well in less ill-posed problems, but requires strong regularity on the solution to be optimal in severely ill-posed cases.

Operational Formulation of boundary data completion, and admissible data definition

Present section is devoted to formally set the operational equation that defines the inverse problem. Through this section and the next one, Ω ⊂ R n+1 (n ≥ 1) will be a domain with a piecewise smooth boundary ∂Ω, Γ 0 and Γ 1 will denote two disjoint C ∞ surfaces of dimension n in ∂Ω, and Σ = ∂Ω \ (Γ 0 ∪ Γ 1 ). It will also be assumed that the euclidean distance between closures of Γ 0 and Γ 1 is not null, and that the interior of Σ in the ∂Ω is regular too. Let us define D α f as the Sobolev partial derivative of order α = (α 1 , ..., α n ) of f ∈ L 1 loc (Ω), H s (Ω) the Sobolev space W s,2 , and H 1 0 (Ω) the closure in the H 1 -norm of the test function D(Ω) (functions infinitely differentiable with compact support in Ω). The L 2 -norm, and the standard operator norm will be denoted by the symbol • . When X is a normed space, the open ball in X centered at x with radius r will be denoted by B X (x, r) ; notation cl X (M ) is reserved to the closure of the set M in the topology of X. As particular cases, notation conventions B X (0, r) = B X (r), and B X (x, r) = cl X (B X (x, r)) are established.

Following subspaces of H 1 (Ω) and H 1 2 (S) are set forth:

E 0 = v ∈ H 1 (Ω) : v| Σ = 0 , E 00 = v ∈ H 1 (Ω) : v| Σ∪Γ1 = 0 , and
E 1 2 (Γ i ) = φ ∈ H 1 2 (Γ i ) : ∃v ∈ E 0 : v| Γi = φ , i = 0 or 1. Notation E -1
2 (Γ i ) will be reserved for the dual of E 1 2 (Γ i ). By Poincaré inequality [START_REF] Luc | An Introduction to Sobolev Spaces and Interpolation Spaces[END_REF], spaces E 0 (Ω), and E 00 (Ω) will be provided with the equivalent norm induced by the inner product:

(∇u, ∇v) L 2 = Ω ∇u∇vdx = Ω n+1 k=1 du dx k dv dx k dx.
Now, consider boundary value problem ( 6)- [START_REF] Bouarroudj | On a direct study of an operator riccati equation appearing in boundary value problems factorization[END_REF], and definition 2.1, slightly differing from the usual definition of weak solution :

∆u ≡ 0 in Ω, (6)
u| Σ = 0, ( 7
) u| Γ1 = Φ, (8)
∂u ∂ν Γ0 = ψ. (9)
Definition 2.1. For a given couple of boundary conditions (Φ, ψ) in E 1 2 (Γ 1 ) × L 2 (Γ 0 ), function u belonging to E 0 is a weak solution of ( 6)-( 9) if equation [START_REF] Bouyssier | Inverse problem in electrocardography via the factorization method of boundary value problems[END_REF] and boundary condition u| Γ1 = Φ are both fulfilled.

Ω ∇u∇vdx = Γ0 ψvdS, ∀v ∈ E 00 . ( 10
)
The set of the standard test function D(Ω) is contained in E 00 , so that the weak solution is also a solution in the sense of [START_REF] Mijailov | Ecuaciones Diferenciales en Derivadas Parciales[END_REF]. The existence, uniqueness, and continuous dependence on the boundary data of the weak solution to problem (6)-( 9) are proven in an analogous way as Mijailov did in [START_REF] Mijailov | Ecuaciones Diferenciales en Derivadas Parciales[END_REF] for elliptic boundary value problems.

Definition 2.2 (Admissible Cauchy data). Data (φ, ψ) in H 1 2 (Γ 0 ) × L 2 (Γ 0 ) will be called an admissible data for the Cauchy problem (1)-(4), if there exist Φ in E 1 2 (Γ 1 ), and u in H 1 (Ω) such that: u solves eq.( 6)- [START_REF] Bouarroudj | On a direct study of an operator riccati equation appearing in boundary value problems factorization[END_REF] in the sense of Definition 2.1 for the pair (Φ, ψ), and the Dirichlet condition u| Γ0 = φ is also fulfilled . Function Φ will be called a solution of the Dirichlet boundary data completion problem for (1)-( 4).

In order to establish the direct problem that defines the admissible Cauchy data class, functions u 1 and u 2 will denote solutions of problem ( 6)- [START_REF] Bouarroudj | On a direct study of an operator riccati equation appearing in boundary value problems factorization[END_REF] for boundary data (Φ, ψ ≡ 0) and (Φ ≡ 0, ψ), respectively. Function u 1 depends only on Φ, u 2 on ψ, and the solution of the auxiliary boundary problem ( 6)-( 9) is given by u = u 1 + u 2 ; hence, identifying the Dirichlet condition in (1)-( 4) can be approached as a first kind linear inverse problem via the Dirichlet to Dirichlet operator A, from Γ 1 to Γ 0 , and the Neumann to Dirichlet operator B defined in Γ 0 :

A : E 1 2 (Γ 1 ) -→ L 2 (Γ 0 ) Φ → u 1 | Γ0 , B : L 2 (Γ 0 ) -→ L 2 (Γ 0 ) ψ → u 2 | Γ0 .
Each function u 1 depending on Φ, and u 2 depending on ψ solves a well-posed boundary data problem, and the trace operator is compact from H 1 (Ω) into L 2 (Γ 0 ) [START_REF] Lions | Non-Homogeneous Bondary Value Problems and Applications[END_REF][START_REF] Luc | An Introduction to Sobolev Spaces and Interpolation Spaces[END_REF][START_REF] Mijailov | Ecuaciones Diferenciales en Derivadas Parciales[END_REF]; hence the linear operators A and B are compact. We extend the range of both operators A and B to L 2 (Γ 0 ) instead of E 1 2 (Γ 0 ) because, taking into account the noise in the measurements, we consider L 2 -norms in the regularization stage.

It is a well known fact that Cauchy problems have a unique solution [START_REF] Aronszajn | A unique continuation theorem for solutions of elliptic partial differential equations or inequalities of second order[END_REF][START_REF] Isakov | Inverse Problems for Partial Differential Equations[END_REF]. Therefore, AΦ = 0 implies that u 1 solves the Cauchy problem for a null Cauchy data given on Γ 0 , and hence A is injective. On the other hand, let us prove that B is injective, self-adjoint and positive defined: choose v in [START_REF] Bouyssier | Inverse problem in electrocardography via the factorization method of boundary value problems[END_REF] as the solution of problem ( 6)-( 9) for a third boundary data (Φ ≡ 0, ψ); B is self-adjoint follows immediately from Definition 2.1; taking ψ = ψ proves that B is positive defined, and under same equality the assumption Bφ = 0 implies that solution of ( 6)-( 9) is null for (Φ ≡ 0, ψ), which proves that B is injective. Now, by definition of A and B, the set M of all admissible data for the Cauchy problem can be defined equivalently as the linear space:

M = (AΦ + Bψ, ψ) : Φ ∈ E 1 2 (Γ 1 ), ψ ∈ L 2 (Γ 0 ) . (11)
The forward problem from which boundary data completion is established is to determine the trace over Γ 0 of the solution of boundary problem ( 6)-( 9) from given mixed boundary data (Φ, ψ), which in the form of an abstract equation means to evaluate AΦ + Bψ. Then we set the inverse problem of interest as follows:

AΦ = ρ, ρ = φ -Bψ. (12)
Characterizing E 1 2 (Γ 1 ) represents an essential step toward an admissible data characterization. The norm in

E 1 2 (Γ 1 ) is defined by Φ E 1 2 = ∇u 1 .
If P is the Dirichlet to Neumann operator in Γ 1 for the boundary data problem ( 6)-( 9) with ψ ≡ 0: P : Φ → ∂u1 ∂ν Γ1 , then, in virtue of Green's formula:

Φ 2 E 1 2 = Φ, P Φ E 1 2 (Γ1)×E -1 2 (Γ1)
.

All the properties of P are studied in detail in [START_REF] Henry | Factorization of Boundary Value Problems Using the Invariant Embedding Method[END_REF].

Regularization Strategy

In later sections it is shown that the Dirichlet to Dirichlet mapping A can be extended in the L 2 sense for cylindrical domains and be considered as a bounded mapping from L 2 (Γ 1 ) to L 2 (Γ 0 ). However, it could happen that A coud not be extended in the general case. The possibility to extend A is helpful for the regularization stage, so that, it will be considered as an additional a priori information when possible. Using a shorter notation, D (A) will denote the domain of A, where D (A) will be understood as L 2 (Γ 1 ) or E 1 2 (Γ 1 ), depending on whether A can be extended or not: A : D (A) → L 2 (Γ 0 ), meaning that the symbol A will be employed to denote the original operator defined in previous section, or its extension to L 2 (Γ 0 ) as appropriate. Now, let u 1 , and u 2 be as before. By replacing u by u 1 , and v by u 2 in (10) it immediately follows that Ω ∇u 1 ∇u 2 dx = 0, and replacing both u and v by u 2 gives that ∇u 2 2 = Γ0 ψBψdS. Hence, analogously to

how norm • E 1 2 (Γ1)
was defined, the natural norm over the admissible data set M in eq [START_REF] Burger | Convergence rates of convex variational regularization[END_REF] is given by the H 1 -norm of the corresponding weak solution of the BVP ( 1)-( 4):

∇(u 1 + u 2 ) 2 = Φ 2 E 1 2 + Γ0 ψBψdS, with Φ = A -1 (φ -Bψ).
However, B is bounded; then, the really natural norm over the admissible data class in this case is a stronger one defined by:

(φ, ψ) ad = Φ 2 E 1 2 (Γ1) + ψ 2 , Φ = A -1 (φ -Bψ).
For the inverse problem ( 12), a measured (or perturbed in synthetic examples) Cauchy data ( φ, ψ) will be known instead of the exact one (φ † , ψ † ), with a L 2 -noise level of order δ > 0, meaning that the inverse problem to regularize is

AΦ † = φ † -Bψ † ; (φ † , ψ † ) ad < ∞, φ † -φ ≤ δ, ψ † -ψ ≤ δ. (13)
Then, defining

M δ = {(φ, ψ) ∈ M : (φ -φ, ψ -ψ) ≤ √ 2δ}, we have that (13) implies (φ † , ψ † ) ∈ M δ ∩ B ad (K) for a given K > 0. ( 14
)
Due to the compact injection of

E 1 2 (Γ 1 ) in L 2 (Γ 1
), a quasi solution regularization scheme is immediately suggested from [START_REF] Sergei | On the a posteriori choice of a regularization parameter in the solution of severely ill-posed problems[END_REF] to the inverse problem [START_REF] Chamorro-Servent | Improving the spatial solution of electrocardiographic imaging: A new regularization parameter choice technique for the tikhonov method[END_REF] when D (A) = L 2 (Γ 1 ). As it will be shown shortly, the admissible data solution (AD solution) is indeed a quasi solution scheme with a slightly more regular definition of ( 14) as cornerstone. For a given K as a priori information, a natural candidate to be a regularization of Φ † , whether the distance from measurement to cl ad (B(K)) can be reached at some element (φ δ , ψ δ ), is of course Φ δ = A -1 (φ δ -Bψ δ ). What AD methodology is looking for in our case, are sufficient conditions to guarantee the existence of (φ δ , ψ δ ), and simultaneously that Φ δ tends to Φ † when δ tends to 0.

The AD solution of problem [START_REF] Chamorro-Servent | Improving the spatial solution of electrocardiographic imaging: A new regularization parameter choice technique for the tikhonov method[END_REF] requires the L 2 minimum distance projection from perturbed Cauchy data over the admissible data set, because of the noise level in [START_REF] Cluitmans | Noninvasive reconstruction of cardiac electrical activity: Update on current methods, applications and challenges[END_REF] is given in the L 2 -norm as well. Unfortunately, the problem to find the aforementioned projection is ill-posed, because generally the admissible data set is not a closed space in L 2 (Γ 0 ) × L 2 (Γ 0 ). To get a good enough behaviour of the projection some a priori information is required. In our case, convexity of the compact subset to which the solution belongs is that required information.

Setting the strategy

As a priori information it will be required that the exact couple (φ † , ψ † ) belongs to a convex subset M of admissible data of the form

M = {(AΦ + Bψ, ψ) ∈ M : Φ ∈ Z 1 , ψ ∈ Z 0 } ; ( 15
)
where Z 0 is a convex subset of L 2 (Γ 0 ), and

Z 1 is compact and convex in D (A) (L 2 (Γ 1 ) or E 1 2 (Γ 1 )
), the required a priori for the construction of a quasi solution of the inverse problem AΦ = ρ.

The minimum L 2 -distance projection from ( φ, ψ) to cl (L 2 (Γ 0 )) 2 (M ) exists, and is unique since cl (L 2 (Γ 0 )) 2 (M ) is closed and convex [20, sec. 3.12]. Let (φ δ , ψ δ ) be the aforementioned projection, also defined by

(φ δ , ψ δ ) = arg min φ -φ 2 + ψ -ψ 2 . (φ, ψ) ∈ cl (L 2 (Γ 0 )) 2 (M ) (16)
We shall see that (φ δ , ψ δ ) is in fact an admissible data, and that the corresponding Dirichlet data Φ δ = A -1 (φ δ -Bψ δ ) converges to the exact solution Φ † when δ tends to 0.

Lemma 3.1. The projection (φ δ , ψ δ ) verifies (φ † -φ δ , ψ † -ψ δ ) ≤ δ 1 , where 0 ≤ δ 1 = 2δ 2 -φ δ -φ 2 -ψ δ -ψ 2 ≤ √ 2δ.
Proof. First,

φ † -φ δ 2 + ψ † -ψ δ 2 = φ † -φ 2 + ψ † -ψ 2 -φ δ -φ 2 -ψ δ -ψ 2 +2 φ † -φ δ , φ -φ δ + ψ † -ψ δ , ψ -ψ δ , (17)
but, cl

(L 2 (Γ 0 )) 2 (M ) is convex and (φ δ , ψ δ ) must ensure (see [20, sec. 3.12]) φ -φ δ , φ -φ δ + ψ -ψ δ , ψ -ψ δ ≤ 0, ∀ (φ, ψ) ∈ cl (L 2 (Γ 0 )) 2 (M ) . ( 18
)
The result follows from ( 17), [START_REF] Lions | Non-Homogeneous Bondary Value Problems and Applications[END_REF], and because (φ

† , ψ † ) ∈ cl (L 2 (Γ 0 )) 2 (M ).
In other words, projection (φ δ , ψ δ ) can be thought as a perturbation of (φ † , ψ † ) with smaller or equal noise error than ( φ, ψ). The immediate step is to solve (19) instead of ( 13), with δ 1 as in Lemma 3.1:

AΦ † = φ † -Bψ † ; (φ † -φ δ , ψ † -ψ δ ) ≤ δ 1 . ( 19
)
Lemma 3.2. The projection (φ δ , ψ δ ) in Lemma 3.1 is an admissible data, and

Φ δ = A -1 (φ δ -Bψ δ ) converges to Φ † when δ tends to 0. Proof. Let {(φ n , ψ n )} n∈N be a sequence in M that converges to (φ δ , ψ δ ). Denoting ρ n = φ n -Bψ n (ρ δ = φ δ -Bψ δ ), and Φ n = A -1 ρ n ∈ Z 1 , it immediately follows that we can choose Φ δ = A -1 ρ δ in Z 1 as any accumulation point of {Φ n } n∈N , which exists since Z 1 is compact in D (A)
. By continuity of A and B we have

AΦ δ + Bψ δ = φ δ . So (φ δ , ψ δ
) is admissible and uniqueness of Φ δ holds true in virtue of the injectivity of A.

To finish the proof, for δ > 0 let N δ be the smallest natural such that

∀n ≥ N δ : Φ N δ -Φ δ D(A) ≤ δ, and φ N δ -φ + ψ N δ -ψ ≤ δ 1 + B ,
then, by triangular inequality

Φ δ -Φ † D(A) ≤ δ + Φ N δ -Φ † D(A)
. But, Φ N δ belongs to the intersection of Z 1 , and the following set

Z δ ( φ, ψ) = {Φ ∈ D (A) : AΦ -ρ ≤ δ} , ρ = φ -B ψ. Hence, Φ N δ converges to Φ † in virtue of Theorem 1 in [4, ch.6, sec 1].
Lemmas 3.1-3.2 actually show the heart of the methodology. We are taking full advantage of the available a priori information to regularize the minimum distance projection from measured data to the admissible data with a quasi solution regularization scheme.

There is more than one way to select Z 1 . Mainly, we are going to work with the one which is related to the worst case error definition in [1, sec1.3]. Let X and Y denote Banach spaces, X 1 a dense subspace of X endowed with the norm • 1 , stronger than the norm in X, T : X → Y a bounded linear operator and ω T (δ, K, • 1 ) the worst case error for T corresponding to the noise level δ > 0 and a priori information x † 1 ≤ K in the framework of the inverse problem T x † = y † :

ω T (δ, K, • 1 ) = sup { x X : x ∈ X 1 , T x Y ≤ δ, x 1 ≤ K} .
In our case, D (A) will play the role of X in the previous definition of the worst case error, Y = L 2 (Γ 0 ), and X 1 will be replaced by F , a subspace with compact injection over D (A); for congruence in the notation • F replaces • 1 . When A can be extended and D (A) = L 2 (Γ 1 ) the immediate choice of F is E 1 2 (Γ 1 ). Otherwise, whether D (A) = E 1 2 (Γ 1 ), we can choose F by considering extra a priori information as in section 3.2, or, for example when is also assumed that the exact solution Φ † belongs to a well known finite dimensional subspace of E 1 2 (Γ 1 ). In the same way as norm (•, •) ad , we define over the subspace of admissible data

M F = (AΦ + Bψ, ψ) : Φ ∈ F, ψ ∈ L 2 (Γ 0 ) its own norm, stronger than (•, •) ad : (φ, ψ) 0 = Φ 2 F + ψ 2 ; Φ = A -1 (φ -Bψ).
In previous framework, Z 1 = B F (K), and Z 0 such that M = B ad (K) i.e

First assumption of a priori information (1 st assumption):

An upper bound K of (φ † , ψ † ) 0 will be given as a priori information.

Now, by how it has been built, the AD regularized solution Φ δ must be optimal in the sense of Definition 1.18 in [START_REF] Kirsch | An introduction to the Mathematical Theory of Inverse Problems[END_REF] for the just above assumed a priori information over the exact Cauchy data. Theorem 3.3. Under the 1 st assumption and if ( φ, ψ) does not belong to cl (L 2 (Γ0)) 2 (M ) then, function Φ δ defined in Lemma 3.2 is an optimal regularized solution of the linear inverse problem of the first kind [START_REF] Luc | An Introduction to Sobolev Spaces and Interpolation Spaces[END_REF] in the sense of Definition 1.18 in [START_REF] Kirsch | An introduction to the Mathematical Theory of Inverse Problems[END_REF], for the error √ 2δ max{1, B } in the data, and a priori information Φ † F ≤ 2K. Proof. By Lemmas 3.1 and 3.2 it happens that, for any n ≥ N δ1 , Φ † -Φ n ∈ B F (K), and

A(Φ † -Φ n ) ≤ δ 1 max{1, B } 1 + φ δ -φ n + ψ δ -ψ n ; (20)
then, defining 1,n = φ δ -φ n + ψ δ -ψ n → 0, and

2,n = Φ n -Φ δ D(A) → 0: Φ † -Φ δ D(A) ≤ Φ † -Φ n D(A) + 2,n ≤ ω A ((1 + 1,n )δ 1 max{1, B }, K, • F ) + 2,n . (21)
However, (1 + 1,n )δ 1 < √ 2δ for any large enough n since ( φ, ψ) ∈ cl (L 2 (Γ0)) 2 (M ).

Remark 3.4. Best case for the 1 st assumption is K = (φ † , ψ † ) 0 .

Even if Theorem 3.3 in the present form is not helpful in numerical calculations, it shows that the regularization error under this regularization methodology is optimal under the 1 st assumption. It is well known that inverse problem ( 19) is exponentially ill-posed, and as it shall be seen in section 4, last definition of Z 1 is equivalent to say that Φ † satisfies a logarithmic source condition (Φ † = -ln -p (A A)ζ, ζ ≤ K 1 , p > 0, where A denotes the A adjoint). The result that ω A (δ, K, F ) must be of order O(-ln -p (δ)) is also shown in [START_REF] Tautenhahn | Optimality for ill-posed problems under general source conditions[END_REF][START_REF] Hohage | Regularization of exponentially ill-posed problems[END_REF], and other recent studies on convergence rate estimates can be found in [START_REF] Albani | Optimal convergence rates results for linear inverse problems in hilbert spaces[END_REF][START_REF] Resmerita | Regularization of ill-posed problems in banach spaces: Convergence rates[END_REF].

Additional a priori information, and link with the Tikhonov regularization

First practical issue of the AD solution is that there is not an equivalent Lagrangian formulation of the problem, because M is not closed in the L 2 -norm. In fact, the imposed constraint by the 1 st assumption may be not defined in (φ δ , ψ δ ) i.e the projection (φ δ , ψ δ ) could be a non-regular point for the constraint defining M in the sense of [START_REF] David | Optimization by Vector Spaces Methods[END_REF]. The present section is devoted to provide conditions for which the equivalent Lagrangian formulation of the main optimization problem exists under the 1 st assumption. The main motivation to do this is to get an approach that allows to turn the value of K from a priori information into a regularization parameter when K is unknown, preserving the optimality of the solution for some regularization parameter, as is discussed in [START_REF] Burger | Convergence rates of convex variational regularization[END_REF][START_REF] Resmerita | Regularization of ill-posed problems in banach spaces: Convergence rates[END_REF].

Stronger a priori information over solution Φ † will be required through this section in order to get an equivalent optimization problem to [START_REF] Henry | Factorization of Boundary Value Problems Using the Invariant Embedding Method[END_REF], where the projection (φ δ , ψ δ ) is a regular point for the corresponding constraint.

Consider the Neumann to Dirichlet operator defined in Γ 1 by

Q : Ψ ∈ L 2 (Γ 1 ) → u 3 | Γ1
, where u 3 is the solution of problem ( 22)-( 25) in an analogous sense to Definition 2.1:

∆u 3 ≡ 0 in Ω, ( 22
)
u 3 | Σ = 0, ( 23
)
∂u 3 ∂ν Γ1 = Ψ, ( 24
)
∂u 3 ∂ν Γ0 = 0. ( 25
)
Just like in the case of operator B, Q is a compact mapping from L 2 (Γ 1 ) to itself, and bounded from L 2 (Γ 1 ) to E 1 2 (Γ 1 ). As an operator in L 2 (Γ 1 ) to itself it is also self adjoint and positive defined. Moreover, it is clear that P is the left inverse of Q (I = P Q), and QP Φ = Φ if P Φ is a regular distribution (P Φ ∈ L 2 (Γ 1 )).

Remark 3.5. Uniqueness of solution u 3 is true in this case by the boundary condition [START_REF] Solodky | On optimization of projection methods for solving some classes of severely ill-posed problems[END_REF], without which the solution is unique up to an additive constant. This fact makes necessary to be careful if the methodology needs to be extended to cases where ∂Ω is composed only by the two regular and smooth surfaces Γ 0 , and Γ 1 . The ECG inverse problem, for instance, is so.

Define E 1 2 s (Γ 1 ) = E 1 2 (Γ 1 ) ∩ Q s (L 2 (Γ 1 )) (s > 0) as follows E 1 2 s (Γ 1 ) = Φ ∈ E 1 2 (Γ 1 ) : Φ = Q s ζ, ζ ∈ L 2 (Γ 1 ) , s > 0; Φ E 1 2 s (Γ1) = ζ . ( 26
) Injection E 1 2 s2 (Γ 1 ) → E 1 2
s1 (Γ 1 ) is compact as soon as s 2 > s 1 > 1 2 . After the definition of E 1 2 s (Γ 1 ) (s > 0) we are ready to work with a stronger condition than the 1 st assumption. Now we will consider the following definition of the space F :

F = E 1 2 s (Γ 1 ) where s ≥ 1 2 if D (A) = L 2 (Γ 1 ), and s > 1 2 if D (A) = E 1 2 (Γ 1 ).
The 1 st assumption is rewritten as

Second assumption of a priori information (2 nd assumption)

There exists ζ † ∈ L 2 (Γ 1 ) such that Φ † = Q s ζ † (s depending on D (A))
, and an a priori given

K > 0 such that (φ, ψ) 0 = ζ † 2 + ψ † 2 ≤ K.
Through the 2 nd assumption we have that solving problem replace ( 16) is equivalent to solve

(ζ δ , ψ δ ) = arg min AQ s ζ + Bψ -φ 2 + ψ -ψ 2 ζ∈L 2 (Γ 1 ) , ψ∈L 2 (Γ 0 ) such that ζ 2 + ψ 2 -K 2 ≤ 0 (27)
and define the AD solution by

Φ δ = Q s ζ δ .
Now, solution of ( 27) is a regular point for its constraint because (ζ δ , ψ δ ) satisfies this constraint. The objective functional and the constraint are both convex and have Gateaux derivatives linear in their increments; then, by the generalized Kuhn-Tucker theorem in the classic optimization theory [START_REF] David | Optimization by Vector Spaces Methods[END_REF], there exists α K ≥ 0 (a Lagrange multiplier) such that solving problem ( 27) is equivalent to solving its Lagrangian formulation for α = α k :

(ζ δ,α , ψ δ,α ) = arg min J α (ζ, ψ) (ζ,ψ)∈L 2 (Γ1)×L 2 (Γ0) (28)
where

J α (ζ, ψ) = AQ s ζ + Bψ -φ 2 + ψ -ψ 2 + α ζ 2 + ψ 2 .
Define the operator W by W : ψ → -∂u2 ∂ν Γ1 (u 2 , and u 1 the same solutions of ( 6)-( 9) that define B, and A respectively). Solution u 2 continuously depends on ψ in L 2 (Γ 0 ), and the outward normal derivative operator over the boundary of Γ is continuous from H

1 (Ω) to E -1 2 (Γ 1 ), meaning that W is bounded from L 2 (Γ 0 ) to E -1 2 (Γ 1
). Moreover, with almost no change on the existence, uniqueness, and continuity demonstrations, solution u 3 defining Q can be extended to Ψ in E -1 2 (Γ 1 ). This means that operator Q can be extended to the compact operator Q :

E -1 2 (Γ 1 ) → L 2 (Γ 1 ); then, operator Qs W (s > 0) is compact from L 2 (Γ 0 ) to L 2 (Γ 1 ). Now, when Φ = Q s ζ with ζ in L 2 (Γ 1 )
, again by Green's formula and because u 1 is orthogonal to u 2 with the inner product (∇•, ∇•) L 2 , it happens that:

Γ0 (AQ s ζ) ψdS = Q s ζ, W ψ E 1 2 (Γ1)×E -1 2 (Γ1) = ζ, Qs W ψ E 1 2 (Γ1)×E -1 2 (Γ1) = Γ1 ζ Qs W ψdS, ∀ψ ∈ L 2 (Γ 0 );
(29) So clearly, Qs W is the L 2 adjoint of A s = AQ s (A s = Qs W ). This way, the J α Gateaux derivative in (ζ, ψ), and applied to x = (x 1 , x 2 ) is given by

DJα(ζ, ψ)(x) = Γ 0 A s Asζ + Bψ -φ + αζ x1dS + Γ 0 B Asζ + Bψ -φ + (1 + α)ψ -ψ x2dS,
where I is the identity in L 2 (Γ 0 ). Hence, (ζ δ,α , ψ δ,α ) solves:

(A s A s + αI) ζ + A s Bψ = A s φ, ( 30
)
BA s ζ + B 2 + (1 + α)I ψ = B φ + ψ. ( 31
)
Parameter α is a Lagrange multiplier, meaning that the α value making equivalent the optimization problems is null (α

K = 0) if there is ( ζ, ψ) in B L 2 (Γ1)×L 2 (Γ0) (K) such that φ = A s ζ + B ψ. Otherwise α K must fulfill ζ δ,α K 2 + ψ δ,α K 2 = K 2 , Φ δ,α K E 1 2 s (Γ1) = ζ δ,α K , Φ δ = Φ δ,α K . ( 32
)
Normal equations ( 30)-(31) provide the first corollary to Theorem 3.3.

Corollary 3.6. If the 2 nd assumption is fulfilled, then, for any α > 0 the pair (Φ δ,α , ψ δ,α ) is defined by (33)-(34)

ψ δ,α = B I -A s L -1 s,α A s B + (1 + α)I -1 B I -A s L -1 s,α A s φ + ψ , (33)
Φ δ,α = Q s L -1 s,α A s φ -Bψ δ,α (34)
where L s,α = A s A s + αI. Moreover, under the 2 nd assumption and if ( φ, ψ) does not belong to cl (L 2 (Γ0)) 2 (M ), there exists α K > 0 such that

Φ † -Φ δ,α K ≤ ω A √ 2δ max{1, B }, 2K, E 1 2 s (Γ 1 ) .
Proof. To get (33)-( 34) it is only needed to algebraically solve system (30)-(31). A well known fact is that L s,α is invertible. By a singular value decomposition it is also easy to see that I -A s L -1 s,α A s is positive defined, implying that any positive real belongs to the resolvent of B I -A s L -1 s,α A s B, which finally proves that Φ δ,α and ψ δ,α are well defined by ( 33)-(34).

As usual in Lagrangian formulations for inverse problems, eq (32) says that K = (Φ † , ψ † ) 0 is what really matters to know a priori in this kind of approach. One has to look for the smallest value of K that makes sense, and minimize the worst case error

ω A √ 2δ max{1, B }, 2K, E 1 2
for a given δ > 0. The dependency of α on K actually makes that parameter α can be considered as a priori information when K is known; however, via the Lagrangian formulation, when K is unknown and (Φ † , ψ † ) 0 < ∞ is the only known fact, then α becomes a regularization parameter. In the last mentioned scheme, previous discussion tells that there is a regularization parameter α AD for which the AD solution in its Lagrangian formulation is optimal, where the question of how to choose the parameter α is introduced as a new issue. In fact, system (30)-( 31) can be rewritten as

(T s T s + αI)(ζ, ψ) = T s ( φ, ψ), T s (ζ, ψ) = A s B 0 I ζ ψ .
In other words, under the 2 nd assumption, the AD solution is the evaluation by Q s of the first coordinate of the Tikhonov solution with parameter α = α K for the inverse problem T s x = y. When the a priori information α K becomes a regularization parameter α, it is important to recall that T s is not a compact operator, so that one has to be carefully selecting the strategy to choose it. In some cases the strategy to choose the regularization parameter for Tikhonov solution, or the properties of its convergence rate to the exact solution depends on the property of compactness of the operator defining the inverse problem. Remark 3.7. It is clear that ψ δ converges to ψ † when δ tends to 0, and because of that the pair (Φ δ , ψ δ ) is a quasi solution of the inverse problem T 0 x = y; however, the ill-posedness of the Cauchy problem is dominated by operator A, and this is why the approach AΦ = ρ has been chosen.

Dirichlet boundary data completion problem in cylindrical domains

Let Γ be a bounded domain in R n (n ≥ 1), Ω = (0, a) × Γ, ∂Γ the boundary of Γ, Σ = [0, a] × ∂Γ, and Γ z = {z} × Γ, where z lies in [0, a], Γ a will take the place of Γ 1 in previous sections. The space L 2 (Γ z ) is isometric to L 2 (Γ) for all non negative real z, so we will consider all these spaces to be identical.

In this case, by the separation of variables method, solution of the problem ( 6)- [START_REF] Bouarroudj | On a direct study of an operator riccati equation appearing in boundary value problems factorization[END_REF] in Ω is given by

u(x; z) = ∞ k=1 w k (z)v k (x), w k (z) = Φ k cosh(zλ k ) + ψ k λ k sinh((a -z)λ k ) cosh(aλ k ) , ( 35
)
where {v k } k∈N is an orthonormal and complete system of L 2 (Γ) composed by eigenfunctions of minus Laplace operator defined on H 2 (Γ) ∩ H 1 0 (Γ), {λ k } k∈N is the set of corresponding eigenvalues. Eigenvalues are repeated according with their multiplicity, and their eigenspaces are finite dimensional. Let Φ k and ψ k be the k-th Fourier coefficient of Φ and ψ in L 2 (Γ) with respect to the system {v k }. The analytic form of A, B, P and Q (as the inverse of P ) are obtained from ( 35)

AΦ = ∞ k=1 Φ k cosh(aλ k ) v k (A = W ), Bψ = ∞ k=1 ψ k sinh(aλ k ) λ k cosh(aλ k ) v k , P Φ = ∞ k=1 λ k Φ k sinh(aλ k ) cosh(aλ k ) v k , QΨ = ∞ k=1 Ψ k cosh(aλ k ) λ k sinh(aλ k ) v k . ( 36
)
The characterization of the admissible data set immediately follows from ( 12) and the analytic form of mappings A and B.

Theorem 4.1. The pair (φ, ψ) ∈ L 2 (Γ) × L 2 (Γ) is an admissible data if and only if ∞ k=1 λ k (φ k - ψ k λ k ) 2 e 2aλ k < ∞; (37)
In that case, the function Φ on the left side of ( 12) is defined by:

Φ = ∞ k=1 φ k cosh(aλ k ) - ψ k λ k sinh(aλ k ) v k . ( 38
)
It is clear that recovering Φ from admissible data (φ, ψ) is severely ill-posed because singular values of A tend to zero exponentially. It is important to recall that λ 1 , λ 2 , ... is a non decreasing sequence that goes to infinity. Notice that

-ln -p (W A)ζ = ∞ k=1 --2aλ k -ln 4(1 + e -2aλ k ) -2 -p ζ k , (39)
then, for p = 1 2 , and for K a depending on K and a, (φ † , ψ † ) ∈ M implies the logarithmic source condition:

Φ † ∈ M p,Ka = Φ : Φ = -ln -p (W A)ζ, ζ ≤ K a . (40)
The worst case error w A (δ, K, • 1 ) can be bounded by estimates in [START_REF] Tautenhahn | Optimality for ill-posed problems under general source conditions[END_REF][START_REF] Hohage | Regularization of exponentially ill-posed problems[END_REF]. On the other hand, the solution u must be in H 1 (Ω), and {v k } is also an orthogonal set in H 1 0 (Γ) with inner product (∇•, ∇•) L 2 satisfying ∇v k = λ 2 k . Hence, by a simple calculation from (35)

Φ, P Φ E 1/2 ×E -1/2 = ∞ k=1 λ k Φ 2 k sinh(aλ k ) cosh(aλ k ) v k , (41)
In the above equation, it is clear that ∇u is upper and lower bounded by a scalar multiple of (-∆) 1/4 Φ . In other words, a solution of the boundary data problem (1)-( 4) exists if and only if E 1/2 (Γ) is the Hilbert space H 1/2 00 defined in [START_REF] Lions | Non-Homogeneous Bondary Value Problems and Applications[END_REF] with inner product as in [START_REF] Abda | Boundary data completion: The method of boundary value problem factorization[END_REF], all details can be reviewed at [16, chap 2]:

f, g H 1/2 00 dS = Γ (-∆) 1/4 f (-∆) 1/4 gdS (42)
Then, H 1/2 00 (Γ) is also defined as the range of the compact operator (-∆) -1 4 , it follows that for every Φ in H 1/2 00 (Γ), there exists τ in L 2 (Γ) such that (-∆) -1 4 τ = Φ i.e. operator Q in second assumption of a priori information in section 3.2 can be replaced by (-∆) -1 2 .

Remark 4.2. The sequence

λ k cosh(aλ) k is bounded (let us say by a constant C); then, AΦ E 1 2 (Γ) = AΦ, P AΦ E 1/2 (Γ)×E -1/2 (Γ) = ∞ k=1 λ k Φ 2 k sinh(aλ k ) cosh(aλ k ) 3 v k ≤ C Φ 2 , (43)
which means that A is continuous from L 2 (Γ) to E 1/2 (Γ) (compact from L 2 (Γ) to itself) i.e. we are in the case where A can be continuously extended: D(A) = L 2 (Γ). Hence, both choices for 2 nd assumption, ( (-∆) -1 4 or Q 1 2 , are equivalents to the 1 st assumption, and produce an optimal regularized solution in this case, because the worst case error is considered in L 2 (Γ 1 ). In further sections we will choose (-∆) -1 4 instead of Q 1 2 for the regularization framework because the numerical example is set in a cylinder where (-∆) -1/4 is easier to compute.

4.1.

Regularized solution in a semi-discretized scheme Let V (m) be the finite dimensional L 2 subspace generated by {v 1 , v 2 , . . . , v m }, V (m)⊥ the orthogonal complement of V (m) , Pr the minimum distance projection over V (m) , and Pr ⊥ the corresponding minimum distance projection over V (m)⊥ . Since {v 1 , v 2 , . . .} is a complete system of L 2 (Γ), and still in the framework of the 2 nd assumption with (-∆) -1 4 instead of Q s with s = 1 2 as mentioned at the end of the previous section:

Φ † = (-∆) -1 4 τ † with τ † in L 2 (Γ), it holds true that Pr ⊥ Φ † = ∞ k=m+1 τ 2 k λ k ≤ δ whether m = m δ,K is such that max K λ 1/2 m δ,K , Pr ⊥ φ , Pr ⊥ ψ ≤ δ. (44)
System {v 1 , v 2 , ...} is complete and orthonormal in L 2 (Γ), the analityc form of A and B in (4) shows that every v k is an eigen function of both operators, as well. Hence, Pr commutes with A and B, implying that it also holds true that A Pr Φ † = Pr φ † -B Pr ψ † , Pr φ † -Pr φ ≤ δ, and Pr ψ † -Pr ψ ≤ δ. This way, considering (Pr φ, Pr ψ) as the measured Cauchy data instead of ( φ, ψ), and denoting the corresponding AD solution by Φ δ,m δ,K a new corollary of Theorem 3.3 is achieved.

Corollary 4.3. The function Φ δ,m δ,K is an asymptotically optimal regularized solution to the inverse problem (13) whether (Pr φ, Pr ψ) does not belong to cl (L 2 (Γ0)) 2 (M ).

Proof. Result follows immediately from previous discussion, and from Theorem 3.3:

Φ † -Φ δ,m δ,K = Pr Φ † -Φ δ,m δ,K + Pr ⊥ Φ † ≤ w A ((1 + B )δ 1 , 2K) + δ. (45)
The regularized solution (Φ δ,m δ,K , ψ δ,m δ,K ) = (Φ m,δ,α , ψ m,δ,α ) is given by solving the normal equations ( 46)-(47):

A 2 + α(-∆) 1/4 Φ + ABψ = A Pr φ, ( 46
)
BAΦ + B 2 + (1 + α)I ψ = B Pr φ + Pr ψ, (47)
for α = α k such that Φ m,δ,α , ψ m,δ,α 0 = K.

Numerical example

The scheme in subsection 4.1 is applied for m = 550 in the cylinder of revolution Ω = (0, a) × Γ, with Γ = {(x, y)|x 2 +y 2 = 1}. As (38) shows, the regularization error exponentially depends on how high the cylinder is (value of a) i.e the ill-posedness order of the inverse problem grows exponentially in function of cylinder height, making the problem virtually unavailable for larger values of a, so that, the example is developed for a = 1. In the present section the AD solution is compared with the Tikhonov regularization strategy, because the last one is still the most common method. The comparison is made in different regularity cases for exact data Φ † :

• Φ † ∈ R((-∆) -1 4
), where R means range. This is the poorest a priori information that can be provided to the operational formulation of the problem. This is a case of low regularity, where optimality for Tikhonov solution is not guaranteed (Φ † ∈ R((W A) p (-∆) 1 4 ), p < 1 with A as a L 2 continuous operator defined in R((-∆) -1 4 )). Optimality for Tikhonov solution is considered with X 1 = R((W A) p ), and x X1 = (W A) -p x in the worst case error definition in section 3; for a given p ≥ 0. In this first case of low regularity: p = 0.

• Φ † ∈ R((W A) p (-∆) -1 4 ), 1 ≤ p ≤ 2.
High regularity such that Tikhonov solution is optimal in the same sense as before. ( see Theorem 2.12 [1, p. 38]).

• Φ † ∈ R((W A) p (-∆) -1 4 ), p > 2.
High regularity where no regularization parameter exists making optimal the Tikhonov solution (Φ † ∈ R((W A) p ), p > 2. In our particular case, W A < 1, by analogous proof of Theorem 1.21 in [START_REF] Kirsch | An introduction to the Mathematical Theory of Inverse Problems[END_REF] one has that ω A (δ, K,

• 1 ) < O(δ 2
3 ), and Theorem 2.13, also in [START_REF] Kirsch | An introduction to the Mathematical Theory of Inverse Problems[END_REF], tells us that Tikhonov solution cannot be optimal for p > 2.

For a given noise level δ > 0, and regularization parameter α, the AD and Tikhonov solutions will be respectively denoted by Φ m,δ,α AD , and Φ m,δ,α T . The best possible regularized solutions in a synthetic example are defined by the smallest regularization parameters solving

α AD = argmin α>0 Φ † -Φ m,δ,α AD , α T = argmin α>0 Φ † -Φ m,δ,α T .
The relative error of a regularized solution Φ is defined by

RE( Φ) = Φ † -Φ Φ †
, as usual. Corresponding histograms are shown in Figure 1. Mean values and standard deviations (std) of a comparative numerical test between the best possible solutions for AD and Tikhonov solutions are presented in Table 1; the comparative is made for any p ∈ {0, 1, 3} determining regularity of the exact solution (Φ † ∈ R (W A)

p (-∆) -1/4 ), and noise level δ ∈ {1e -3, 1e -6, 1e -9}. In all cases a 300 size sample of synthetic data (Φ † , ψ † , φ, ψ) were built via the analytic form (35):

• ζ † , ψ † , δ (1) , δ (2) ∈ V (550) such that any one of their Fourier coefficients in the {v 1 , v 2 ...} system is a pseudo random number with a normal standard distribution. Besides the error, we build a large sample of pseudo-random exact pairs (ζ † , ψ † ) in L 2 (Γ) × L 2 (Γ) to control the regularity on Φ † through its Fourier series.

• Φ † = (W A) p (-∆) -1 4
ζ † , according with the cases described above.

• φ = AΦ † + Bψ † + δ δ (1)
δ [START_REF] Kirsch | An introduction to the Mathematical Theory of Inverse Problems[END_REF] , and ψ = ψ † + δ δ (2) δ (2) . 1. Mean values and standard deviations for the best possible solutions for AD and Tikhonov solutions. Samples size:300. Height of the cylinder: a = 1. For p = 0 the identity I = (W A) 0 is considered in L 2 (Γ a ). Figure 2 shows that low regularity on the exact solution implies a higher order on the worst case error. Two cases are seen, in the first one the exact solution is not regular enough to belong to a class of functions for which Tikhonov strategy behaves well, while the second example was made applying W A to the exact data in the first example, in order to include an example where best possible Tikhonov solution must be optimal. In both cases the shared Neumann component is the following non differentiable function given in polar coordinates: ψ † (r, θ) = |r -θ|(r sin(3r) cos(θ)) 3 , as long as Dirichlet boundary conditions are

RE(Φ

Φ † (r, θ) = (1 -r) cos(θ) 3 r , and Φ † (r, θ) = W A (1 -r) cos(θ) 3 r ; φ † = AΦ † + Bψ † .
Cauchy data (φ † , ψ † ) were perturbed in the same pseudo-random way than before with a noise level δ = 10 -4 to obtain ( φ, ψ). Remark 4.4. In this case, an eigenfunction of -∆ has the form v(r, θ) = J i (λr) cos(jθ), v(r, θ) = J i (λr) sin(jθ),

where J i is the special Bessel function of first kind and order i, and λ are such that J i (λ) = 0.

Link with the invariant embedding

The boundary data completion problem in Γ a has been approached as the first kind linear inverse problem in eq. [START_REF] Chamorro-Servent | Improving the spatial solution of electrocardiographic imaging: A new regularization parameter choice technique for the tikhonov method[END_REF]. That can also be done in the framework of the invariant embedding in [START_REF] Abda | Boundary data completion: The method of boundary value problem factorization[END_REF], with the slightly difference that the embedding is taken from a to 0 instead of from 0 to a, and state u defined by ( 6)-( 9) is the only one used. The boundary data problem ( 6)-( 9) is embedded in a family of boundary data problem defined in the sub-domains Ω s = (s, a) × Γ. Consider the following families of BVP: . This way, equation (49) holds, and ( 12) is obtained from it, at s = 0:

∆u (s) 1 = 0 in Ω s , u (s) 1 Σ = 0, u (s) 1 Γa = Φ, ∂u (s) 1 ∂ν Γs = 0, ∆u (s) 2 = 0 in Ω s , u (s) 2 Σ = 0, u (s) 2 Γa = 0, ∂u (s) 2 ∂ν Γs = ζ. Let A(s) : H 1/2 00 (Γ a ) → H
u(s; •) = B(s) ∂u 2 ∂ν Γz + A a (s)Φ. ( 49
)
Through the same formal derivation scheme as in [START_REF] Abda | Boundary data completion: The method of boundary value problem factorization[END_REF][START_REF] Bouarroudj | On a direct study of an operator riccati equation appearing in boundary value problems factorization[END_REF][START_REF] Henry | Factorization of Boundary Value Problems Using the Invariant Embedding Method[END_REF], with h = a -z, all the following equations are fulfilled for h ∈ (0, a) , where:

d dh A(h) + B(h)∆ Γ A(h) = 0, A(0) = I, d dh B(h) + B(h)∆ Γ B(h) = I, B(0) = 0.
∆U (s) 1 = 0 in Ω s , U (s) 1 Σ = 0, U (s) 1 Γs = τ, ∂U (s) 1 ∂ν Γ0 = 0, ∆U (s) 2 = 0 in Ω s , U (s) 2 Σ = 0, U (s) 2 Γs = 0, ∂U (s) 2 ∂ν Γ0 = ψ.
According with [START_REF] Abda | Boundary data completion: The method of boundary value problem factorization[END_REF][START_REF] Bouarroudj | On a direct study of an operator riccati equation appearing in boundary value problems factorization[END_REF][START_REF] Henry | Factorization of Boundary Value Problems Using the Invariant Embedding Method[END_REF]:

d ds P + P 2 = -∆ Γ , P (0) = 0, d ds W + P W = 0, W (0) = I.
Hence, the invariant embedding in [START_REF] Henry | Factorization of Boundary Value Problems Using the Invariant Embedding Method[END_REF] allows to set the normal equations of the AD solution depending on operators satisfying the above formal differential equations, also suggesting that it can be implemented in the same discretized and semi-discretized schemes than the factorization method in more complex geometries of Ω, or when the goal is to solve as fast as possible the same inverse problem for a set of measured Cauchy data recorded in a period of time, computing once for all the involved operator and avoiding time consuming numerical methods. The key of those statements is that the same regularization approach developed in sections 2 to 3.1 apply for the following cases:

• Ω a cylinder as before, where null Dirichlet boundary condition [START_REF] Faker | Why is the cauchy problem severely ill-posed?[END_REF] is replaced by the Neumann one ∂u ∂ν Σ = 0, and the spaces E 0 and E 00 are defined as follows

E 0 = v ∈ H 1 (Ω) : ∂u ∂ν Σ = 0 , E 00 = v ∈ E 0 : v| Γ1 = 0 ,
• Ω = Ω 0 \ cl R n (Ω 1 ), where Ω 0 and Ω 1 are bounded domains in R n with regular boundaries, and such that cl R n (Ω 1 ) ⊂ Ω 0 . In this case Γ 0 = ∂Ω 0 , Γ 1 = ∂Ω 1 , boundary condition ( 7) is omitted, and

E 0 = E 00 = v ∈ H 1 (Ω) : v| Γ1 = 0 .
In both of previous modifications of the inverse problem, L 2 (Γ 1 ) ⊥ = Φ ∈ L 2 (Γ 1 ) : Γ1 ΦdS = 0 must be the domain of Q. Both examples are developed in [START_REF] Henry | Factorization of Boundary Value Problems Using the Invariant Embedding Method[END_REF] in the framework o the factorization method when there exists a C 1 diffeomorphism between Γ 0 and Γ 1 .

Discussion

The proposed solution was built following a new regularization scheme, which do not provide a general strategy in the regularization theory; however, it allows to take full advantage on all available a priori information about the solution for a given problem, following an intuitive way to build optimal approximations to the solution. When the measured Cauchy data is projected over the set of admissible data, we are looking for the nearest Cauchy data in a subclass of functions defined by strong enough a priori information. This subclass is such that the weak solution of the Cauchy problem is well defined on it, in other words, the approximation is optimal because we forced it.

The AD solution in this case leads to a constrained optimization problem, were constraints are defined by the available a priori information. If this information guarantees that the solution is a regular point for the given constraints, they can become regularization parameters via the equivalent Lagrangian multipliers. In that case and by the following modification T s (Φ, ψ) = (A s Φ + Bψ, ψ) described in section 3, the AD solution can be thought as the image of a compact operator applied over the first component of the Tikhonov solution for the linear and non compact operator. Another way of looking at it is through eq. ( 34), where one can easily see that (A * s A S + αI)A * s ( φ -Bψ δ,α ) is the Tikhonov solution for the original inverse problem Aφ = ρ were the perturbed output data have been continuously modified by eq. ( 33). Anyhow, it seems that careful modifications of well known strategies can be applied under the 2 nd assumption of a priori information, strategies like discrepancy principle, L-curve [START_REF] Hansen | Regularization tools: A matlab package for analysis and solution of discrete ill-posed problems[END_REF], U -curve, balancing principle [START_REF] Solodky | On optimization of projection methods for solving some classes of severely ill-posed problems[END_REF], or ADP technique [START_REF] Chamorro-Servent | Improving the spatial solution of electrocardiographic imaging: A new regularization parameter choice technique for the tikhonov method[END_REF].

For more complex geometries of Ω, when the goal is to solve as fast as possible the same inverse problem for different output data recorded in a period of time, section 4.2 and [START_REF] Bouyssier | Inverse problem in electrocardography via the factorization method of boundary value problems[END_REF][START_REF] Henry | Factorization of Boundary Value Problems Using the Invariant Embedding Method[END_REF] suggest that the invariant embedding can be an option for numerical implementations of the admissible data methodology.

The last remarkable aspect is that, applying this procedure, a kind of deterministic hypothesis test is given when δ and (Φ † , ψ † ) 0 are known. In practice, if the projection (φ δ , ψ δ ) satisfies the inequality (φ δ , ψ δ ) -( φ, ψ) ≥ √ 2δ, then it is possible to say that there is something wrong with the device which recorded the Cauchy data, or, in the worst of the cases, that the model of the phenomena is not good enough for the desired purposes.

Finally, the obtained regularized solution was compared with Tikhonov regularization strategy, since it is still the most common applied regularization strategy, even though it is well known that it cannot be optimal for high enough smoothness of the input data in linear inverse problems [1, Theorem 2.13 p.39].

Conclusion

In cylindrical domains is enough to known the value of the natural norm of the exact Cauchy data (φ † , ψ † ) in the class of admissible data; in the general assuming slightly stronger a priori information is enough. With the employed methodology, an optimal regularized approximation to the boundary data completion problem is attained. However, in the Lagrangian formulation, the regularized solution depends on the a priori information K = (Φ † , ψ † ) 0 , which is rarely available in practice. The remaining research issue is to obtain a method to choose the parameter α in (28) that only depends on the measurement.

Numerical examples show that the best possible AD solution behaves well in a semi-discretized scheme; however, it is important to explore other discretization schemes for different, and more complex geometries of Ω, where computation of eigenfunctions v k , and eigenvalues λ k is a problem itself.

Figure 1 .

 1 Figure 1. Histograms of RE AD-T ( φ, ψ), in 300 size samples, with noise level error δ, and regularity of the exact data determined by p (Φ † ∈ R((W A) p (-∆) -1/4 )). Height of the cylinder: a = 1. For p = 0 the identity I = (W A) 0 is considered in L 2 (Γ a ).

  s ) be the Dirichlet to Dirichlet operator defined by: A(s)Φ = u s ) the Neumann to Dirichlet operator defined by B(s)ζ = u

Figure 2 .

 2 Figure 2. Comparison of AD and Tikhonov solutions. First line: Φ exact = Φ † is given by ((1r) cos(θ) 3 )/r. Projection of Φ exact over V (550) , P rΦ exact (left), best possible AD solution (center), best possible Tikhonov solution (right). Second line: Φ exact = W A ((1 -r) cos(θ) 3 )/r . Projection of Φ exact over V (550) , P rΦ exact (left), best possible AD solution (center), best possible Tikhonov solution (right). Semi-discretization scheme in section 4.1. Height of the cylinder: a = 1.


Imad Alawe ⇤

Adlen Ksentini

Yassine Hadjadj-Aoul ⇤

Philippe Bertin

Improving traffic forecasting for 5G core network scalability: A Machine Learning approach

5G is expected to provide network connectivity to not only classical devices (i.e. tablets, smartphones, etc) but also to the Internet Of Things (IOT), which will drastically increase the traffic load carried over the network. 5G will mainly rely on Network Function Virtualization (NFV) and Software Defined Network (SDN) to build flexible and on-demand instances of functional networking entities, via Virtual Network Functions (VNF). Indeed, 3GPP is devising a new architecture for the core network, which replaces point to point interfaces used in 3G and 4G, by a producer/consumer-based communication among 5G core network functions, facilitating deployment over a virtual infrastructure. One big advantage of using VNF, is the possibility of dynamically scaling, depending on traffic load (i.e. instantiate new resources to VNF when the traffic load increases, and reduce the number of resources when the traffic load decreases). In this paper, we propose a novel mechanism to scale 5G core network resources by anticipating traffic load changes through forecasting via Machine Learning (ML) techniques. The traffic load forecast is achieved by using and training a Neural Network on a real dataset of traffic arrival in a mobile network. Two techniques were used and compared: (i) Recurrent Neural Network (RNN), more specifically Long Short Term Memory Cell (LSTM); and (ii) Deep Neural Network (DNN). Simulation results showed that the forecast-based scalability mechanism outperforms the threshold-based solutions, in terms of latency to react to traffic change, and delay to have new resources ready to be used by the VNF to react to traffic increase.

I. INTRODUCTION

Network virtualization has paved the way for 5G to meet high flexibility and agility expectations. Network virtualization allows running complex network functions (namely, Virtual Network Function -VNF) on top of a virtualized infrastructure, hosted at central or edge telco cloud. Chained together, the VNFs allow building flexible and on demand virtual networks. Such flexibility is highly required in 5G, where services need to be created on demand and for a given duration. Relying on Network virtualization, standardization groups, like ETSI and 3GPP, are multiplying efforts to build a new environment for 5G. ETSI, with the Network Function Virtualization (NFV) reference architecture 1 , has created a reference model to orchestrate and manage VNF running on top of Virtual Infrastructure. Meanwhile, 3GPP is developing specifications to deploy 3GPP elements, including the core network (CN) functions and some parts of the Radio Access Network (RAN) functions [START_REF] Gpp | System architecture for the 5g system[END_REF] [START_REF]Telecommunication management;study on management and orchestration of network slicing for next generation network[END_REF], on top of a virtualized infrastructure. Particularly, the envisioned 5G core network has been devised with the aim to run the new functions on a virtualized environment. One option of the 5G core architecture is to replace classical 3GPP interfaces (i.e. point to point) by a bus-like communication to interconnect the core network elements, following the concept of web services.

VNF dimensioning is one of the main challenges of running core network elements in a virtualized environment, in terms of resources, which highly impacts performances [START_REF] Alawe | On evaluating different trends for virtualized and sdn-ready mobile network[END_REF]. Indeed, the virtualized core network elements should perform similarly to their hardware version; supporting the same Quality of Service (QoS): supported number of users, provided data rate, etc. To reach the same performance as the hardware version, over-provisioning the allocated resources for a VNF could be envisioned but it might be costly. Alternatively, adaptive provisioning of virtual resources could represent an interesting solution; where the VNF is over sized (by increasing the number of CPUs or by multiplying the number of VMs instances of the VNF) when traffic load is growing, and down sized when traffic is decreasing. Such solution is known as a scale-up/down (if the number of CPU is increased/decreased) or scale-out/in processes (if the number of VMs serving the same VNF is augmented/reduced), and is highly used in cloud computing. To scale VNF resources, a solution to take scaling decisions is needed. While in cloud computing the decision is mainly based on system level information, like consumed CPU or memory, obtained via the virtulization platform, mobile networks require service level information, like the number of users connected to the system or the traffic load, which indicate if the VNF is overloaded or not.

In this paper, we are interested in the scalability of a 5G core network key element, namely the Access and Mobility Management Function (AMF). The latter represents the bottleneck of the mobile network control plane, as it handles UEs attach requests. Accordingly, VMs hosting the AMF should be well dimensioned in order to avoid increasing the attach duration of UEs, and the percentage of rejected requests. The AMF should be over sized (scale-out), if the number of attach requests is high, and down sized (scale-in) if the number of such requests is low. In [START_REF] Alawe | On the scalability of 5g core network: the amf case[END_REF], we proposed a mechanism to scale the AMF using a threshold-based algorithm, which was derived using control theory. The proposed mechanism uses a threshold (the number of attach requests received by the AMF) to take the decision to scale-out or -in. However, due to the time needed to scale-out (i.e. the time needed for a new VM to be operational), the proposed solution suffers from delay to react to an increase in the number of attach requests; which leads to an increase of the attach procedure duration, and hence the rejection of some requests. To overcome this issue, we propose in this paper to use Machine Learning (ML) to forecast requests' arrivals, and hence anticipate the scale-out process; avoiding the rejection of requests and keeping the attach duration low.

The remainder of this paper is organized as follows. Section II introduces the necessary background to understand the paper, focusing on the 5G core network architecture and related works. Section III first introduces our motivations to ensure AMF scalability; then, the machine learning model and its usage to predict the user attach requests, and the proposed ML-based mechanism to scale-out and -in the AMF. Section IV presents the result of using the ML algorithm to predict the user attach based on the Telecom Italia dataset2 . Finally, we conclude the paper in section V.

II. BACKGROUND

A. 5G core network architecture

The 5G core network architecture has been devised with the aims: (i) to keep compliance with the concept of 3GPP communication reference interfaces (point-to-point); (ii) to evolve in near future to a service-oriented architecture, where the 5G core control network functions communicate via a common bus, through the concept of producers and consumers. Indeed, 3GPP SA2 group has defined two architectures, one based on the reference interfaces, and the other one following a service oriented architecture. In this paper, we will focus on the second one as the new components are expected to be run on top of a virtualized environment. For the first architecture readers may refer to [START_REF] Gpp | System architecture for the 5g system[END_REF]. Regarding the 5G core functions, some are similar to 4G EPC ones and some are completely newly designed. Notably, the access control and session management functions are combined in the Mobility Management Element (MME) of EPC, but separated in 5G core to better support fixed access and ensure scalability and flexibility. The most relevant network functions, as defined in 5G Core, are:

• Core Access and Mobility Management Function (AMF), which handles access control and mobility. If fixed access is involved, the mobility management will not be required in the AMF. It is in charge of handling attach requests of UEs from several Radio Access Networks (RAN); as being the entry point to the core network control plane it creates bottleneck issues when undersized; • Session Management Function (SMF), which sets up and handles sessions according to the network's policy;

• Policy Control Function (PCF), which corresponds to the Policy and Charging Rules Function (PCRF) in LTE.

It is worth noting that this function will integrate also a policy framework for network slicing; • User Plane Function (UPF), which can be deployed, based on the service type, in several configurations and locations to be able to handle and forward users' data traffic; • Unified Data Management (UDM), which is similar to the Home Subscribe Server (HSS) in EPS. However, it is envisioned that it integrates subscriber information for both fixed and mobile access in the New Generation (NG) Core; In the following we shortly describe the AMF and SMF functions as this work focuses on the AMF scalability. For more details on the other functions reader may refer to [START_REF] Gpp | System architecture for the 5g system[END_REF]. The AMF ensures UE-based authentication, authorization, mobility management, etc. UEs, even using multiple access technologies, are connected to a single AMF, since the latter is agnostic to the access technology. The SMF is responsible for session management and IP address allocation to UEs. In addition, it selects and controls the UPF for data transfer. When a UE has multiple sessions, several SMFs may be allocated to each session aiming at managing the latter individually and possibly provide different functionalities per session. Based on operator policy or AF provided policy on packet flow, the PCF defines policies about mobility and session management, which will be enforced by the AMF and SMF.

Figure 1 illustrates the service oriented 5G core network architecture. Only the UPF is deployed on the legacy transport network infrastructure; the other core network functions handling the control plane are expected to run in more centralized telco cloud platform. While the communications with the UE, RAN and UPF are using the reference interfaces, the core control plane functions are connected together via a common communication bus. As stated earlier, the CN functions use the concept of producer/consumer, where a CN function (CNF) may register, using the NRF, for specific events provided by another CNF via an API. Thus, either 1:1 or 1:N communication is possible. For instance, the AMF service exposes information to the other CNFs regarding events and statistics related to mobility; or a PCF provides all the operations related to policy rules to other CNFs. The service-oriented 5G core architecture could be easily deployed in a virtualized environment (e.g. VM or container), wherein libraries of functions may be requested from a VNF catalog and composed into virtual Core Network on demand.

B. Related Work

As illustrated in Figure 1, the AMF is the entry point of the 5G core network; which makes it the system bottleneck. Thus, a special focus should be given to the dimensioning of the VNF(s) hosting the AMF to ensure system scalability. Scalability in this case is more a software deployment issue rather than a hardware upgrade as most of the CNF are running as VNFs, on top of a virtualized infrastructure.

One of the most common scaling solutions consists in using static thresholds. In [START_REF] Dutta | Qoe-aware elasticity support in cloud-native 5g systems[END_REF], [START_REF] Carella | An extensible autoscaling engine (ae) for software-based network functions[END_REF] and [START_REF] Alvi | Load based automatic scaling in virtual ip based multimedia subsystem[END_REF] authors have proposed scalability based on thresholds. In the same spirit, we worked also on threshold based scalability in a previous attempt to manage AMF instances in [START_REF] Alawe | On the scalability of 5g core network: the amf case[END_REF]. The main idea is to deploy a new instance (scale-out) when the CPU usage (or other resources) exceeds a certain fixed "scale-out" threshold. Reciprocally, an instance is shut off (scale-in), when a considered threshold exceeds a fixed value. Static thresholds-based solutions are already deployed and widely used in commercial and open-source cloud solutions, like [START_REF] Amazon | Web services auto scaling[END_REF] and [START_REF] Cloudstack | Open source cloud computing[END_REF]. However, static thresholds-based approaches are not suitable for mobile networks, and have several limitations due to the fact that they react to events (i.e. reactive solution). Indeed, deploying a new instance of a VNF, when reacting to events, may take from one minute to a dozen of minutes depending on the data center architecture and the complexity of the VNF instance to scale; leading to delay the scaling-out process, hence increasing the attach duration of UE and rejection of some requests.

An alternative to static-threshold-based solutions is the usage of adaptive techniques, as proposed in [START_REF] Arteaga | An adaptive scaling mechanism for managing performance variations in network functions virtualization: A case study in an nfv-based epc[END_REF]. In this work, the authors have proposed a mechanism that combines Q-Learning with Gaussian Processes-based system models, which allows to adapt to dynamic environments and improve the scaling policy before taking any action. The authors assume that their proposal reacts better than static threshold rules. However, it remains a reactive solution, and hence inherits the same weaknesses. Moreover, using Reinforcement learning, it may take a considerable time before starting to have the correct decisions, which is not acceptable in systems like 5G, due to the time constraint of certain type of applications. Note also the side effects associated to this type of approaches, consisting in the influence of the external environment on performance and therefore the risk of quality degradation.

Proactive solutions for scalability have been also proposed in the literature. The authors, in [START_REF] Bilal | Dynamic cloud resource scheduling in virtualized 5g mobile systems[END_REF], suggested to pro-actively schedule resources for VNFs based on CPU usage forecasting from a historical dataset. Although this technique is interesting for its pro-activity, it does not offer adaptation in case of traffic pattern evolution and is not compliant with mobile networks requirements. Indeed, the proposed solution considers only system level information, like consumed CPU or memory, ignoring service level information. VNF scaling in mobile networks should take into account strict rules as defined in standards and therefore deploying only a new instance by a cloud orchestrator is not sufficient for scaling an end-to-end service. Consequently, mobile networks require service level information, like the number of connected users or traffic load when scaling-out or -in. Moreover, the parameters of the "offline" prediction equation, proposed by the authors, are determined statically and are completely independent from the considered dataset. This implies stability in the weights that are given for each day prediction, and thus a risk of divergence in the prediction with the slightest change of pattern. Concerning the "online" approach, the authors propose to use an average (exponential average) estimator. However following [START_REF] Fu | Using lstm and gru neural network methods for traffic flow prediction[END_REF], when a pattern can be determined, this is clearly not the best solution.

III. ML-BASED SCALABLE MECHANISM FOR 5G CORE NETWORK

A. Scalability issue of the AMF

To illustrate the need of an efficient scalability mechanism for the AMF, we conducted a set of experimentations using OpenAirInterface (OAI) solution 3 . OAI is a set of open source software tools that allows rapid prototyping of 5G mobile network; including the Radio Access Network and the Core Network functions. We used the emulated version of OAI, known as OAI simulator (OAISIM), which enables the emulation of a set of connected UEs and one eNodeB. The AMF was run on top of a VM, with two configurations: 1 CPU and 1 Gigabyte (Gb) of Memory (noted as first configuration); 2 CPU and 2 Gb of memory (noted as second configuration). It is important to note that we refer to AMF, even if we used the MME function of the OAI EPC. Indeed, our focus is mainly to compute the number of accepted UE and the time needed for a UE to attach to the network when varying the number of attach requests per second. Both metrics refer to the performance of authentication procedure, which is one of the common functions between the MME and AMF.

Figure 2(a) represents the accepted UE according to the number of UE attach requests per second. Clearly we observe that from a certain threshold, both configurations start rejecting UE attach. Obviously, stronger the In this case, we observe that from the same threshold 20 UEs/sec, for the first configuration, and 40 UEs/sec second for the second one, the attach duration is exponentially increasing, which is violating the requirements specified by the 3GPP standard 4 . Obviously, there is a need to have a mechanism to solve this issue, which may consist on either: (i) adding more resources (CPU and memory) to the VM, i.e. scaling-up; or (ii) adding a new VM to the AMF function (with an internal load balancer), i.e. scaling-out.

B. Traffic forecasting: Deep learning-based approach

As explained above, there is a need of a mechanism to allow dynamic dimensioning of the AMF resources to avoid degrading the system QoS. The envisioned mechanism should be proactive to avoid the weakness of reactive solutions; hence requiring a traffic prediction solution. For this aim (i.e. traffic prediction), two Neural-networksbased techniques were envisioned and tested aiming at predicting the upcoming traffic load. The first one is based on Deep learning Neural Network (DNN), while the second one relies on Recurrent Neural Networks (RNN), more specifically on Long Short Term Memory cell (LSTM). In a deep neural network (DNN), a particular neuron can be activated by the contribution of several inputs at once, and therefore there is no real recognition of a traffic pattern. On the contrary, a recurrent neural network, and particularly an LSTM network, receives the inputs successively and through a memory and omission mechanism, it manages to memorize traffic patterns, as it has been demonstrated in several works [START_REF] Hochreiter | Long short-term memory[END_REF].

The main procedures for traffic prediction using learning techniques are: getting the dataset, formatting the data, designing the neural network, training the neural network and then starting predicting (or testing your network). Each step is important in order to get the best accuracy of prediction compared to real values.

1) DataSet:

The main key of prediction using deep learning, is the dataset. Nowadays, operators are able to collect those data easily as they are the owner of their infrastructure and have a clear idea of the evolution of the load over their network. In this work, we have used the dataset collected during the Big Data Challenge organized by Telecom Italia 2 . The data set is open source [START_REF] Italia | A multi-source dataset of urban life in the city of milan and the province of trentino dataverse[END_REF] and rich in information. It can be used to predict the evolution of the load in a core network.

The considered data set contains several kinds of data, like sms in/out, call in/out and the Call Detail Records (CDRs) with an interval of ten minutes for two months; including information on both control plane and data plane traffic. We consider in this work that at least 10% of the traffic is control plane traffic, and the remaining 90% is data plane traffic. Therefore, for the remainder of this paper, we will only use the control plane traffic (10% of the load in each period of ten minutes) for prediction.

2) Data Formatting: Once the dataset is collected and 10% of each period is extracted, the data should be formated in order to be injected in the neural network for the training phase. This part requires a high attention as the data format can affect the accuracy of prediction of the neural network. Formatting intelligently the data will allow the neural network to learn more in the training phase, and thus predicting with a high accuracy. However, adding lots of Key Performance Indicators (KPI) may lead to decreasing the prediction accuracy further, particularly for small datasets. Therefore a balance should be found between the dataset size and the KPI injected in the data formatting. In our case, we decided to format the dataset by classifying the load into 10 different classes, which can be interpreted as: class 1 is a load interval where (x) AMF instances are needed; Class 2 (x + y) AMF are needed, and so on. Note that adding a much larger number of output classes would require a much larger dataset. Indeed, the addition of classes implies the increase of the number of parameters of the neural network which has an effect of requiring more data to have a significant precision.

3) Neural networks-based prediction:

The first technique used for predicting the class of load of the next period is DNN. DNN is a neural network composed of three main layers: the Input layer, the middle layer(s) (hidden layer(s)) and finally the output layer. This type of neural network is also known as the Feed Forward Neural Network (FFNN). The data goes in the neural network from the input layer through the middle layer(s), and finally go out by the output layer.

The second technique tested for predicting the average load of the upcoming period of time is RNN, which is another type of neural network. Unlike DNN, RNN have loops. Those loops allow to inject previous events while predicting future events. It is a sort of memory block allowing to take decision depending on the previous event(s). Also, in RNN multiple kinds of cells are used today. In this paper we chose to rely on the LSTM cell. The advantage of the LSTM cell is to allow data patterns to be stored without degradation over time.

In order to compare the accuracy of the DNN and the RNN (LSTM), one DNN and one RNN (LSTM) were designed using the Tensor Flow library 5 . We divided and formatted the data into classes of average load. Then, we split the original dataset into two separate datasets. The first one contains 60% of the data, for training the two neural networks; while the remaining 40% are left to test the robustness of the two networks. Then both networks are trained with the 60% dataset and then asked to predict the remaining 40%. Figure 3 illustrates the prediction of both the DNN network and the RNN network. From Figure 3 we notice that RNN (LSTM) network performs better than DNN network. The RNN (LSTM) network achieves an accuracy of ⇡ 90%, while the DDN network obtains ⇡ 80% of accuracy; hence the RNN is performing more than 10% better than DNN when predicting the class of average load of the upcoming period of time. The accuracy of the prediction is mainly due to the ability of LSTM networks to recognize a traffic pattern, unlike DNN networks. The accuracy of the DNN network can, however, be improved but this would require more neurons and more data, something that is not always available.

Based on these results, we will consider using RNN (LSTM) network for the AMF scalability mechanism.

C. Prediction-based scalability

Having described the prediction model, now we present the prediction-based scalability algorithm. Based on the ETSI NFV reference architecture, the scalability of a VNF is a decision taken by the NFV Orchestrator (NFVO), using the information provided by both the Virtual Infrastructure Manager (VIM) and the VNF Manager (VNFM). The former provides information on the VM execution environment (like CPU usage, Memory usage), while the latter provides information obtained from the Element Manager (EM) hosted in the VNF. The information provided by EM, represents the inside VNF service performances, like the number of attach requests or the size of the queue of the AMF, etc. In [START_REF] Alawe | On the scalability of 5g core network: the amf case[END_REF], we proposed that the EM of the AMF reports this information to the VNFM, hence the NFVO. The latter takes decisions to scale-in or -out using predefined thresholds. However, due the time needed to scale out (run a new VM), the proposed solution suffers from delay to react to the increase of UE attach; which leads to an increase in the time of the users' attach with some losses of requests.

To remedy to this problem, in this work, we will use the RNN to predict the needed number of AMFs that efficiently manages the traffic load, aiming at keeping the attach duration under an acceptable value (using the results shown in section III.A). The NFVO uses the RNN to predict the traffic (i.e. class of traffic), hence derives the optimal number of AMF to pro-actively deploy before the change of traffic load; anticipating any network load changes and removing the time needed to start a new VM instance, in case of scale-out process. Besides, in this Fig. 3: Reality vs NN prediction work there is no need to monitor the VNF performance via the VNFM, as the NFVO relies only on the trained RNN, which permits to reduce the communication load introduced by the exchanged messages between the EM and VNFM.

IV. PERFORMANCE EVALUATION

To evaluate the performances of the proposed solution and compare it with the threshold-based mechanism, we developed a discrete event simulator using Simpy 6 . The choice of Simpy was motivated by its efficiency when it comes to developing a simulator from scratch as well as its integration with other tools we have used, such as tensorflow. The simulator consists of:

• UE traffic generator, which allows to initiate UE attach requests according to the dataset • an eNodeB that dispatches the arriving of control requests upon the deployed instances of AMF using a load balancing mechanism, if more than one AMF instance is used. It is based on the algorithm introduced in [START_REF] Alawe | On the scalability of 5g core network: the amf case[END_REF] • each VNF (hosting an AMF component) can processes up to 20 UE attach requests per second.

• an NFVO that manages the scalability of the AMF. We launched several simulations to compare the performances of the proposed RNN-based scalability mechanism with the threshold-based solution. The former predicts the average load for the next period of time (10 minutes) and deploys the exact number of AMF instances needed to handle the upcoming load; while the latter launches a new instance (scale-out) when the overall arrival average load exceeds 80% of the overall system capacity (depending on the number of active AMF instances at the moment), and removes an instance (scale-in) when the overall average load decreases below 10% of the overall core network capacity. Both mechanisms were tested over 48 hours of arrivals data and with a deployment latency (i.e. time needed for the AMF instance to be operational) of 60 seconds. We chose this value according to [START_REF] Frangoudis | Cdn-as-a-service provision over a telecom operators cloud[END_REF], wherein the authors have studied the deployment latency for different VNF deployment scenarios on top of a private cloud infrastructure.

Figure 4 illustrates the attach requests rates vs the number of AMF instances deployed by both mechanisms; i.e. the threshold-based and the RNN-based. From Figure 4 we notice that both mechanisms are able to follow the arrival flow by deploying more AMF instances when the arrival load increases, and by reducing the number of AMF instances when the traffic load decreases. We remark that the RNN-based mechanism follows perfectly the arrival load. Furthermore, the RNN-based solution always anticipates the traffic load, and the number of AMF instances is always following the traffic trends. In contrast, the threshold-based mechanism takes more time to react. Indeed, we observe that in many cases, the two curves (Traffic load and number of AMF in Figure 4(a)) cross each other, meaning that the threshold-based solution is taking time to react; leading to a degraded network QoS. This observation is confirmed in Figure 5, which illustrates the number of blocked attach requests in both solutions. Indeed, we remark that several requests are rejected in the case of the threshold-based solution compared to the RNN-based solution. The accurate prediction of the RNN-based solution allows to react in advance to the traffic changes by pro-actively instantiating the necessary VNF resources; avoiding overloading the AMF while waiting for the additional resource to be instantiated. Moreover, this result confirms that the RNN-based solution accurately estimates the needed number of AMF instances to manage a specific traffic load. For example when the load is maximal the RNN-based mechanism deploys up to 8 instances, while threshold-based solution deploys up to 7. Reciprocally, for the case when the arrival load is minimum, the threshold-based mechanism deploys at least 2 instances, while the RNN-based mechanism instantiates at least 3 AMF instances. Thus, the RNN calculates more accurately the number of instances needed to manage the load, allowing to stabilize the system and significantly decrease the number of rejected requests. Indeed, as the RNN predicts the arrival load of the upcoming period, it deploys the AMF instances needed before the upcoming load in a way that the AMF instances are up and ready with the increase of the Load.

To sum up, the RNN(LSTM) technique performs better than the threshold technique as firstly it acts in a proactive way unlike the threshold technique. Secondly it allows to deploy accurately the exact needed number of AMF instances to absorb the upcoming load. Finally, it allows to remove the deployment latency by scaling out in advance the needed AMF instances. In that way, the mobile core is ready to absorb the upcoming load increase without suffering from overhead and extra latency for the control procedures.

V. CONCLUSION AND FUTURE WORK

In this paper, we proposed a novel solution, based on machine learning, to scale-out and -in the AMF resources in a virtualized environment. Using a neural network, which has been trained on mobile network traffic dataset to forecast the user attach request rate, allows to predict the exact number of AMF instances needed to handle the upcoming user traffic. By being proactive, the proposed solution allows to avoid the deployment latency when scaling out resources. The latter may degrade the network performances, since the AMF is overloaded while waiting for the additional resource to be instantiated. Simulation results confirmed the efficiency of the RNN-based solution compared to a threshold-based solution. Future work aims to use the trained RNN to estimate the number of UPF (User plane 5G core network function) needed to handle the data plane traffic, in a virtualized environment.

Fig. 1 :

 1 Fig. 1: Service oriented 5G Core network

Fig. 2 :

 2 Fig. 2: Benchmarking result of the AMF

Fig. 4 :

 4 Fig. 4: Evolution of AMF instances number with control traffic Requests (a) threshold-based solution (b) RNN-based solution

Figure 6 Fig. 6 :

 66 Fig. 5: Number of Blocked control traffic requests of AMF instances number with control traffic Requests (a) threshold-based solution (b) RNN-based solution

ETSI, Network Functions Virtualisation (NFV); Architectural Framework, https://bit.ly/2GGPZpN, accessed in march

Telecom Italia, Telecom Italia Big Data Challenge: http://www.telecomitalia.com/tit/en/bigdatachallenge/contest.html, accessed in march 2018

http://www.openairinterface.org/, accessed in march 2018

3GPP, Non-Access-Stratum (NAS) protocol for Evolved Packet System (EPS); Stage 3, https://bit.ly/2H0Tred, accessed in march 2018

TensorFlow, An open-source machine learning framework for everyone: https://www.tensorflow.org/, accessed in march 2018

Team-SimPy, Discrete event simulation for python: https://simpy.readthedocs.io/en/latest/, accessed in march 2018