Hamza Ennaji
email: hamza.ennaji@unilim.fr

Noureddine Igbida
email: noureddine.igbida@unilim.fr.

NOUREDDINE IGBIDA † Van Thanh Nguyen
email: nguyenvanthanh@qnu.edu.vn.

Thanh Van

AUGMENTED LAGRANGIAN METHODS FOR DEGENERATE Hamilton-Jacobi Nguyen

2021 Equations

Hamilton-Jacobi Equations

AUGMENTED LAGRANGIAN METHODS FOR DEGENERATE HAMILTON-JACOBI EQUATIONS

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Here, H is continuous in both arguments, coercive and convex in the second, but not enjoying any property of existence of a smooth strict subsolution. Our study here takes place in the framework of viscosity solution theory, so that the term (sub, super) solution must be understood in this sense.

Recall that the nonexistence of smooth strict subsolution is connected to the fact that the socalled Aubry set, noted by A, is not empty. It is known that in this case the problem has many solutions. Nevertheless, one needs to select the pertinent one from the physical point of view as well as the optimal control one, such as, the maximal subsolution. Such kind of equations are relevant, for instance, in the study of Shape from Shading [START_REF] Durou | Numerical methods for shape-from-shading: A new survey with benchmarks[END_REF][START_REF] Lions | Shape-from-shading, viscosity solutions and edges[END_REF][START_REF] Prados | A viscosity solution method for shape-from-shading without image boundary data[END_REF][START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF][START_REF] Zhang | Shape-from-shading: a survey[END_REF], granular material [START_REF] Igbida | Metric character for the sub-Hamilton-Jacobi obstacle equation[END_REF] and many other problems ranging from mechanics, geometry, seismic analysis, optics to image processing, etc.

Typical situation appears in the study of the so-called Eikonal equation subject to Dirichlet boundary condition, that is

|∇u| = f in Ω u = 0 on ∂Ω, (1.1)
where f ≥ 0 is a continuous function. In this case the Aubry set corresponds to the zero set [f = 0]. In the case where this set is empty, i.e. f > 0 in Ω, the equation (1.1) falls into the scope of Eikonal equation of the type k |∇u| = 1, where k is a continuous positive function on Ω. Among famous numerical methods (especially for Eikonal equations), remain the Fast Sweeping Method (FSM) and the Fast Marching Method (FMM). The FSM (see e.g. [START_REF] Luo | Fast sweeping methods for factored anisotropic Eikonal equations: multiplicative and additive factors[END_REF][START_REF] Danielsson | Euclidean distance mapping[END_REF][START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF]) is based on an upwind difference discretization solved via Gauss-Seidel iterations with alternating sweeping ordering. As to the FMM (see e.g. [START_REF] Albert | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF]), the Eikonal equation is again written via an upwind finite difference approximation. The grid points are divided into three categories: accepted nodes which are points where the values of the solution are already known, narrow band nodes where the computations take place and finally far nodes which will be computed in the next iterations. For other methods and commented references, the interested reader can check [START_REF] Falcone | Numerical methods for Hamilton-Jacobi type equations[END_REF][START_REF] Albert | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF] for a general overview of numerical methods for Hamilton-Jacobi equations and their applications.

The case where A is nonempty (degenerate HJ equation) is different. Recall that the main role of this set concerns the uniqueness of the solution. Roughly speaking a solution is characterized in a unique way by its values on the set of degeneracy, i.e. on A. So, some additional condition is imposed to single out numerically the reasonable solution among infinitely many ones: for example, it is assumed that the solution is given on the Aubry set, which becomes then a part of the boundary of the domain where the problem is discretized (see e.g. [START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF]). Nevertheless, in practice we do not, in general, have this knowledge. To overcome this obstruction, some methods suggest to regularize the problem by cutting from below f at some level > 0, and then let → 0. This approach was exploited in [START_REF] Camilli | Numerical approximation of the maximal solutions for a class of degenerate Hamilton-Jacobi equations[END_REF] using deeply control theoretic interpretation for the approximation of the maximal subsolution. Their scheme was based on a two-step discretization of the control problem associated with the regularized problem: first in the time variable, and then in the space variable. Under some appropriate condition, the approximated solution is shown to be convergent to the maximal solution using stability results. For comprehensive exposition of Hamilton-Jacobi equations and optimal control one can see, for example, [START_REF] Cannarsa | Semiconcave functions, Hamilton-Jacobi equations, and optimal control[END_REF][START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Lions | Generalized solutions of Hamilton-Jacobi equations[END_REF].

1.2. Main contributions. In this paper, we suggest a new approach to tackle general Hamilton-Jacobi problem without requiring any additional assumptions on the Aubry set. The way is to use variational formulation for the characterization of the maximal subsolution. Then using duality method, we give a direct PDE-constrained optimization algorithm to approximate the solution. In particular, this approach empowers the use of fast numerical optimization methods, employed in the recent development of mathematical programming, as alternatives to the commonly used methods we presented previously. We can not cover all the possibilities with this approach in this paper, but we propose to try out its prospects to tackle some examples of degenerate Hamilton-Jacobi equations. We show how one can use methods like augmented Lagrangian one to solve Hamilton-Jacobi equations and present some numerical results. Coming back to the Eikonal equation, let us outline an overview on our approach and main contributions. As a starting point, we observe that problem (1.1) aims to characterize the maximal subsolution; i.e. the maximal continuous functions subject to a Dirichlet condition and the constraint on the gradient (1.2)

|∇u| ≤ f,
For the non-degenerate case, i.e., f (x) > 0 on Ω, the Fenchel-Rockafellar strong duality result (cf. [20, Chap. III]) implies that the optimal value of (1.2) coincides with

min φ∈M b (Ω) N Ω f d|φ| : -div(φ) = ρ in D (Ω) ,
where M b (Ω) N denotes the set of vector valued measures concentrated in Ω. This duality plays an important role in our approach since wherefore the solution of (1.1) is given by the saddle-points of a problem of the type min u,q max φ L(u, q; φ)

that will be precise later. In the degenerate case; i.e. f ≥ 0, the previous duality falls out of the scope of standard Fenchel-Rockafellar duality result since the so called qualification conditions are not fulfilled. Using perturbation techniques, we prove that the duality and the saddle-points characterization still hold true in the degenerate case. As far as the degeneracy of the Hamiltonian is involved the boundary condition takes part definitely in the proofs of duality. Indeed, working with non homogeneous boundary condition like

u = g on ∂Ω,
requires to handle rigorously the trace of vector fields Φ involved both into the dual problem and the saddle point formulation. To overcome these difficulties, we make use of a trace-like operator for the so called divergence-measure field (cf. [START_REF] Chen | Some recent methods for partial differential equations of divergence form[END_REF][START_REF] Chen | Divergence-measure fields and hyperbolic conservation laws[END_REF][START_REF] Chen | Extended divergence-measure fields and the euler equations for gas dynamics[END_REF]). For numerical computation, we operate thus the saddle-point structure in the augmented Lagrangian (cf. [START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF])

L r (u, q; φ) = L(u, q; φ) + r 2 |Λu -q| 2 , for any (u, q, φ)

for a given arbitrary r > 0. Notice here that the addition of the quadratic term r 2 |Λu -q| 2 has the advantage of improving the convergence speed. Then, for numerical experiments we use the so called ALG2 algorithm 1 (also known as ADMM: alternating direction method of multipliers [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]) to construct a sequence {u i } i , {q i } i , {φ i } i by optimizing alternatively in u and q, for i ≥ 1 to lay out numerical approximations of the solutions of (1.1) in miscellaneous context that are meaningful from modeling viewpoint (degenerate, non-degenerate and/or anisotropic Eikonal equation, Riemannian case as well as polyhedral one). For further details about the ALG2 algorithm, we refer the reader to [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF][START_REF] Fortin | Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems[END_REF][START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF].

To unroll the foregoing machinery we describe above from the Eikonal equation to general Hamilton-Jacobi equations H(x, ∇u) = 0, we exploit the metric character of these equations and address similar problems. More precisely, we establish the equivalence between being a viscosity subsolution of HJ equations, i.e. H(x, ∇u) ≤ 0, and ∇u ∈ B σ * where σ * is the dual of 1 ALG2 is an abbreviation of algorithm 2 in [START_REF] Fortin | Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems[END_REF] to calculate efficiently saddle points of functionals of the form Lr. It is essentially based on relaxation of Uzawa's algorithm.

the support function σ of the sublevel sets of H, which will be given later. So, in order to solve the HJ equations we will focus on the problem of the form max Ω udx : ∇u ∈ B σ * and u = g on ∂Ω (1.3) to end up in some sense with the study of problem of maximizing Ω udρ for a given measure

ρ (ρ ≡ 1 in (1.
3)), among 1-Lipschitz functions with respect to d σ . We show thereby that the aforementioned strategy for Eikonal equation works and offers further issue for very general (even degenerate) quasiconvex Hamiltonians. Up to our knowledge, most of the existing methods to approximate solutions of HJ equations of the form H(x, p) = F (p) -f (x) assuming that f > 0.

At last, let us notice that though the approach proposed in [START_REF] Camilli | Numerical approximation of the maximal solutions for a class of degenerate Hamilton-Jacobi equations[END_REF] is more concerned with the discretization of the intrinsic metric d σ , the FMM and FSM repose on a direct discretization of the equation, our strategy is more in the sprite of "optimize then discretize". Besides the degenerate case (see Tests 5.3-5.4-5.5), our approach presents a general method to solve numerically Hamilton-Jacobi equations by performing simple operations like projections into Finsler balls. Actually, it works even for general quasicovex Hamiltonians, the challenge is to be able to project efficiently onto the 0-sublevels of the Hamiltonians. These can be seen in the nontrivial Examples 5.6-5.7-5.8. In addition, let us point out that the other examples, especially the Riemannian, anispotropic and crystalline cases (see tests Examples 5.6-5.7-5.8) are not widely present in the literature. Some papers are exclusively devoted to each case (e.g. [START_REF] Luo | Fast sweeping methods for factored anisotropic Eikonal equations: multiplicative and additive factors[END_REF][START_REF] Mirebeau | Efficient fast marching with Finsler metrics[END_REF][START_REF] Mirebeau | Anisotropic fast-marching on Cartesian grids using lattice basis reduction[END_REF][START_REF] Mirebeau | Hamiltonian Fast Marching: A Numerical Solver for Anisotropic and Non-Holonomic Eikonal PDEs[END_REF] for the anisotropic and non-holonomic Eikonal equations).

Thanks to the structure of the novel variational formulation developed in this paper for (even degenerate) quasiconvex HJ equations, it is clear that other numerical optimization methods may also have potential benefits than the proposed one in this paper. One can proceed for instance also with the Chambolle-Pock method (see for instance [START_REF] Hamza Ennaji | Continuous Lambertian Shape from Shading: A Primal-Dual Algorithm[END_REF]). Preconditioning techniques should furthermore improve the numerical experiments as well.

Let us mention also that our approach is close to the one in [START_REF] Caboussat | On the numerical solution of some Eikonal equations: an elliptic solver approach[END_REF] as they consider a maximization problem. However, they deal with nonconvex constraint contrary to the convex problem we are considering.

1.3. Plan of the paper. This paper is organized as follows. In Section 2, we recall the metric character of HJ equations. More precisely, we present a distance d σ of Finsler type, constructed from the sublevel sets of the Hamiltonian, which is essential to characterize the subsolutions of HJ equations. We also show that the unique maximal subsolution is characterized by the maximization problem of the form max Ω udx : ∇u ∈ B σ * and u = g on D ⊂ Ω .

(1.4)

Dealing with the problem (1.4), it is natural to derive its dual problem. This will be the main goal of Section 3 where we discuss the duality issue for a general Hamiltonian. This will be done via perturbation techniques. Finally we present the augmented Lagrangian, and we see that thanks to the duality result, solving (1.4) amounts to finding saddle points of a Lagrangian L.

We give several numerical tests to illustrate our results.

Metric character of Hamilton-Jacobi equations

2.1. Preliminaries. In this paper Ω ⊂ R N is a regular connected bounded open domain. We consider the following HJ equation

H(x, ∇u) = 0, x ∈ Ω, (2.5)
where The properties and notions introduced in this subsection are more or less known. Proofs and more details can be found in [START_REF] Michael | User's guide to viscosity solutions of second order partial differential equations[END_REF][START_REF] Lions | Generalized solutions of Hamilton-Jacobi equations[END_REF][START_REF] Siconolfi | Hamilton-Jacobi equations and Weak KAM Theory[END_REF][START_REF] Michael | Viscosity solutions of Hamilton-Jacobi equations[END_REF]. We recall briefly some definitions related to viscosity theory. Given two continuous functions φ and ψ, one says that ψ is a strict supertangent (respectively subtangent) to φ at some point x ∈ Ω if x is a strict local maximizer (respectively minimizer) of φ -ψ. Definition 2.1. A continuous function u : Ω → R is said to be a viscosity subsolution of (2.5) if H(x, ∇φ(x)) ≤ 0 for any x ∈ Ω and any C 1 supertangent function φ to u at x. It is said to be a viscosity supersolution of (2.5) if H(x, ∇ψ(x)) ≥ 0 for any x ∈ Ω and any C 1 subtangent function to u at x. Finally, u is a viscosity solution of (2.5) if it is both a subsolution and a supersolution.

H : Ω × R N → R
We denote by S - H (Ω) the family of subsolutions of (2.5). One pertinent property of this family is its stability with respect to uniform convergence, which is not the case for a.e. solutions as we can see by considering a sawtooth function. Moreover, whenever we consider a family C ⊂ S - H (Ω) of locally equibounded functions, then inf u∈C u(x) and sup u∈C u(x) are still subsolutions to (2.5). We similarly denote by S + H (Ω) (resp. S H (Ω)) the family of supersolutions (resp. solutions) of (2.5). For x ∈ Ω, we define the support function of the 0-sublevel set Z(x) by σ(x, q) := sup q • Z(x) = sup{ q, p : p ∈ Z(x)} for q ∈ R N .

The assumptions (H1) -(H2) ensure that σ is a continuous nonnegative function in Ω × R N , convex and positively homogeneous with respect to q. Due to the assumption (H3), σ(x, q) is possible to equal to 0 for q = 0, which leads to the degeneracy and its dual σ * , as given below, may take the value +∞. Here, the dual σ * (also called polar) is defined by

σ * (x, p) := sup q∈R N { p, q : σ(x, q) ≤ 1}. For example, in the case of Eikonal equation |∇u(x)| = f (x), one has σ(x, q) = f (x)|q| and σ * (x, p) = 1 f (x)
|p|. We see that σ * (x, p) take the value +∞ for p = 0 on the zero set of f .

We denote by Γ(x, y) the set of Lipschitz curves defined on [0, 1] joining x, y in Ω. We then define the intrinsic distance by

d σ (x, y) := inf ζ∈Γ(x,y) 1 0 σ(ζ(t), ζ(t))dt,
which is a quasi-distance, i.e. satisfying d σ (x, x) = 0 and the triangular inequality, but not necessarily symmetric. We summarize the basic properties of d σ in the following.

Proposition 2.2. ([22])

1) d σ is a quasi-metric, in the sense that for any x, y ∈ Ω d σ (x, y) ≥ 0 and d σ (x, x) = 0.

Moreover, for all x, y, z ∈ Ω one has

d σ (x, y) ≤ d σ (x, z) + d σ (z, y). 2) For any x ∈ Ω d σ (x, .) ∈ S - H (Ω) ∪ S + H (Ω \ {x}). 3) Compatibility condition: v ∈ S - H (Ω) if and only if v(x)-v(y) ≤ d σ (y, x)
for any x, y ∈ Ω. The so called Aubry set is defined as the set where the quasi-metric d σ degenerates. Prescribing a boundary value on ∂Ω does not guarantee the uniqueness of viscosity solutions to (2.5) unless A = ∅, which is not the case in our situation due to the assumption (H3). The Aubry set A appears then to be a uniqueness set for (2.5).

Definition 2.3. We define the Aubry set A as the set of points x ∈ Ω such that there exists

(ζ n ) n ∈ Γ(x, x) with l(ζ n ) ≥ δ > 0 and inf n 1 0 σ(ζ n (t), ζn (t))dt = 0, where l(ζ n ) is the Euclidean length of the curve ζ n .
Proposition 2.4. ( [START_REF] Fathi | PDE aspects of Aubry-Mather theory for quasiconvex Hamiltonians[END_REF], [START_REF] Camilli | Systems of convex Hamilton-Jacobi equations with implicit obstacles and the obstacle problem[END_REF])

1) The Aubry set A is a closed subset of Ω.

2) If x ∈ A then d σ (x, .) ∈ S H (Ω). Moreover, x / ∈ A if and only if (2.5) admits a strict subsolution around x.

3) If

g : A ∪ ∂Ω → R is a continuous function satisfying the compatibility condition g(x) - g(y) ≤ d σ (y, x) on A ∪ ∂Ω, then u(x) = min y∈A∪∂Ω {d σ (y, x) + g(y)}
is the unique viscosity solution of the equation (2.5) such that u = g on A ∪ ∂Ω.

2.2.

Equivalence between HJ and maximization problem. Given a closed subset D ⊂ Ω (typically D = ∂Ω or D = {x} for some x ∈ Ω), we consider the following HJ equation

H(x, ∇u) = 0 Ω \ D u = g D (2.6)
where g : D → R is continuous function satisfying the compatibility condition. Then the unique maximal viscosity subsolution of the equation (2.6) such that u(x) = g(x) for any x ∈ D, is given by u

(x) = min y∈D d σ (y, x) + g(y) (2.7)
with d σ (., .) being the intrinsic metric associated to the Hamiltonian H.

The considerations given in the introduction lead us to look for the maximal subsolutions of the HJ equation, i.e., H(x, ∇u) ≤ 0 or equivalently u(x) -u(y) ≤ d σ (y, x). We will show that such functions are precisely the ones with gradient in the unit ball of σ * . More precisely, we prove the following Proposition 2.5. S -

H (Ω) = {u ∈ W 1,∞ (Ω) and σ * (x, ∇u(x)) ≤ 1 f or a.e x ∈ Ω} := B σ * .
Then, we transform the problem into a question of maximization of the volume Ω udx among the subsolutions u. This leads to the following theorem which is an important step to treat the equation (2.6) via augmented Lagrangian methods.

Theorem 2.6. The maximal viscosity subsolution of (2.6), given by (2.7), is the unique solution of the problem

max Ω z(x)dx, σ * (x, ∇z(x)) ≤ 1 and z = g on D .
(2.8)

For the proof of Proposition 2.5, we recall that the result is more or less known in the case where H(x, 0) < 0 which corresponds to A = ∅ (see [START_REF] Igbida | Sub-gradient diffusion operator[END_REF][START_REF] Igbida | Augmented Lagrangian method for optimal partial transportation[END_REF] for example). Here, under the general condition (H3), we need the following Cauchy-Schwartz-like lemma.

Lemma 2.7. For any q ∈ R N , and for any x ∈ Ω we have

σ * (x, q) ≤ 1 ⇔ p, q ≤ σ(x, p) for any p ∈ R N .
Proof. Assume first that σ * (x, q) ≤ 1. If p, q = 0, it is obvious. For the case p, q > 0, if σ(x, p) = 0 then by homogeneity of σ, σ(x, λp) = 0 for every λ ≥ 0. Consequently, σ * (x, q) ≥ λ p, q → ∞ as λ → ∞, this contradicts σ * (x, q) ≤ 1. This implies that σ(x, p) > 0 and p σ(x, p) , q ≤ σ * (x, q) ≤ 1, as desired. Conversely, if p, q ≤ σ(x, p), by definition of σ * , we take the sup over all p such that σ(x, p) ≤ 1 in the previous inequality to obtain that σ * (x, q) ≤ 1.

Proof of Proposition 2.5. We divide the proof into two parts. Firstly, S - H (Ω) ⊂ B σ * . Take a point x ∈ Ω such that u is differentiable at x. For every vector v ∈ S N -1 on the unit sphere, we take ζ h (t) = tx + (1 -t)(x -hv) for every h > 0 and t ∈ [0, 1] and we notice that ζ h ∈ Γ(x -hv, x) and for h small enough, ζ h is close to x. We then have

∇u(x), v ≤ lim h→0 + h -1 (u(x) -u(x -hv)) ≤ lim inf h→0 + h -1 d σ (x -hv, x) ≤ lim inf h→0 + h -1 1 0 σ(ζ h , ζh )dt = lim inf h→0 + 1 0 σ(ζ h , v)dt ≤ σ(x, v),
where we have used the continuity of x → σ(., v) and Lebesgue theorem. Using the definition of σ * , we deduce that σ * (x, ∇u(x)) ≤ 1 as desired. Secondly, B σ * ⊂ S - H (Ω). Assume now that σ * (x, ∇u(x)) ≤ 1 at any differentiable point x of u, i.e., ∇u(x), p ≤ σ(x, p) for all p ∈ R N (Lemma 2.7). For smooth function u, the argument is simply given by

u(y) -u(x) = 1 0 ∇u(ξ(s)) • ξ(s)ds ≤ 1 0 σ(ξ(s), ξ(s))ds,
for Lipschitz curves ξ joining x to y in Ω. This implies that u(y) -u(x) ≤ d σ (x, y). For general Lipschitz function u, one can make use of smooth approximation (see [START_REF] Van | Monge-Kantorovich equation for degenerate Finsler metrics[END_REF]Lemma 3.1]).

Proof of Theorem 2.6. Let us see that for any subsolution v (2.6), we have v ≤ u. Pick any subsolution v of (2.6) satisfying the condition v = g on D. Let x ∈ Ω, we have for any

y ∈ D v(x) = v(x) -v(y) + g(y) ≤ d σ (y, x) + g(y).
This gives that v ≤ u in Ω and then

Ω v(x)dx ≤ Ω u(x)dx. Clearly u is 1-Lipschitz with respect to d σ ,

HJ and duality results

3.1. Generalities. As we said in the introduction among our main interests in this paper is to use augmented Lagrangian methods to give a direct algorithm to approximate the solution of the Hamilton-Jacobi equations.

Actually, the problem (2.8) falls into the scope of the following class of optimization problem

(P) inf u∈X F(u) + G(Λu)
where X and Y are two Banach spaces with the topological dual spaces X * and Y * , F : X → (-∞, +∞], G : Y → (-∞, +∞] are proper, l.s.c., convex functions and Λ : X → Y is a linear operator. The main future for this kind of formulation is to look at its dual problem, which is given by

(D) sup v∈Y * -F * (-Λ * v) -G * (v) .
Here Λ * is the adjoint operator of Λ, while F * , G * are the Legendre-Fenchel transformations of F and G given by

F * (f ) = sup u∈X ( f, u -F(u)) for any f ∈ X * G * (g) = sup q∈Y ( g, q -G(q)) for any g ∈ Y * .
Using the definitions of F * , G * and Λ * , it is not difficult to get the following inequality

sup v∈Y * -F * (-Λ * v) -G * (v) ≤ inf u∈X F(u) + G(Λu),
usually called weak duality. The reverse inequality, the so-called strong duality, plays an important role since the point is that the solutions of (P) and (D) are given by the saddle-points of the problem min u,q max φ L(u, q; φ)

where L(u, q; φ) = F(u) + G(q) + φ, Λu -q . Yet, the strong duality is not true in general. Some qualification conditions are needed. For instance, one has

max v∈Y * -F * (-Λ * v) -G * (v) = inf u∈X F(u) + G(Λu).
Moreover, u solves (P) if and only if

-Λ * u ∈ ∂F(u) and v ∈ ∂G(Λu).
As we pointed out in the introduction, because of the degeneracy of the Hamiltonian, we can not directly use Theorem 3.8 to show duality between the maximization problem and its dual problem. To see this, let us restrict ourselves to the following particular case: the Hamilton-Jacobi equation of Eikonal type, coupled with a zero Dirichlet boundary condition,

|∇u(x)| = f (x) in Ω u = 0 on ∂Ω,
where f is a continuous, nonnegative function on Ω. In this particular case, the problem (M D ) can be rewritten as where

V = C 1 (Ω) ∩ H 1 0 (Ω), F(u) = - Ω udx, Λu = ∇u and G(q) = 0 if |q| ≤ f +∞ otherwise.
For the case of non-degeneracy, i.e., f (x) > 0 on Ω, we have for u 0 = 0 that G(.) is continuous at Λu 0 = 0 in C(Ω) N equipped with the uniform convergence. It follows directly from Theorem 3.8 that

sup Ω udx : |∇u| ≤ f and u = 0 on ∂Ω = min φ∈M b (Ω) N Ω f d|φ| : -div(φ) = ρ in D (Ω) .
However, if f ≥ 0, it is not true that the functional G is continuous at Λu 0 . Therefore the sufficient conditions of Thereom 3.8 are not satisfied.

The main goal of this section is to show rigorously that the duality still holds true for such a degenerate Hamiltionian. Since the values of u are prescribed on D and the solution of (2.6) is given by the distance to D, we can consider the following problem

(M D ) : max Ω D udρ : σ * (x, ∇u(x)) ≤ 1 and u = g on Γ D ,
where

Ω D := Ω \ D, Γ D := ∂Ω D ∩ ∂D and Γ = ∂Ω D \ Γ D . In particular, if D = ∂Ω then Ω D = Ω and Γ D = ∂Ω, Γ = ∅.
As we will see, our dual formulation challenges some kind of trace-like operator for the so called divergence-measure field. To begin with let us sort out formally and briefly our approach, at least in the case where D = ∂Ω. Taking X to be the Banach space W 1,∞ (Ω) and Y to be the space L ∞ (Ω) N , we consider simply

F(u) =      - Ω u dρ if u ∈ W 1,∞ (Ω) and u = g on ∂Ω +∞ otherwise, , Λu = ∇u, and
G(η) =    0 if σ * (x, η) ≤ 1
+∞ otherwise, so that the problem (M D ) reads as (P). For a formal computation of the dual problem, let us notice that G * : Y * → (-∞, +∞],

G * (Φ) = Ω σ(x, Φ) for any Φ ∈ Y * .
The operator Λ * : Y * → X * is given by

Λ * Φ, ξ = Φ, ∇ξ for any ξ ∈ X .
Also,

F * (-Λ * Φ) = sup u∈X ,u=g on ∂Ω Ω -Φ • ∇u dx -F(u) = sup u∈X ,u=g on ∂Ω Ω -Φ • ∇u dx + Ω u dρ = sup u∈X ,u=g on ∂Ω Ω udiv(Φ) - ∂Ω Φ • ν u + Ω u dρ =    - ∂Ω Φ • ν g if -div(Φ) = ρ +∞ otherwise.
In other words, the dual problem reads

inf Ω σ(x, Φ) - ∂Ω Φ • ν g : Φ ∈ Y * , -div(Φ) = ρ in D (Ω) . (3.9)
Note that the above computation is very formal by the appearance of the trace-like term Φ • ν which is not well defined for all Φ.

To handle rigorously the normal trace of the vector-valued dual variable Φ in the dual problem of the type (3.9), we recall the trace-like operator for the so called divergence-measure field (cf. [START_REF] Chen | Some recent methods for partial differential equations of divergence form[END_REF][START_REF] Chen | Divergence-measure fields and hyperbolic conservation laws[END_REF][START_REF] Chen | Extended divergence-measure fields and the euler equations for gas dynamics[END_REF]). To this aim, we assume in this section that

Ω D = Ω \ D is a regular domain with a deformable Lipschitz boundary ∂Ω D .
This is achieved for instance in the case where Ω is a regular domain and D = ∂Ω or D = ω with a regular domain ω ⊂⊂ Ω.

For any 1 ≤ p ≤ ∞, we define the set

DM p (Ω D ) := F ∈ L p (Ω D ) N : div F =: µ ∈ M b (Ω D ) ,
where µ = div F is taken in D (Ω D ). See here, that for any F ∈ DM p (Ω D ), the total variation of div F is given by

|div F |(Ω D ) := sup ϕ dµ : ϕ ∈ C 0 (Ω D ), |ϕ(x)| ≤ 1 for any x ∈ Ω D = sup ϕ dµ : ϕ ∈ C 1 0 (Ω D ), |ϕ(x)| ≤ 1 for any x ∈ Ω D = sup F • ∇ϕ : ϕ ∈ C 1 0 (Ω D ), |ϕ(x)| ≤ 1 for any x ∈ Ω D .
In particular, the space DM p (Ω D ) endowed with the norm

F DM p (Ω D ) := F L p (Ω D ) + |div F |(Ω D )
is a Banach space.

Thanks to [START_REF] Chen | Some recent methods for partial differential equations of divergence form[END_REF][START_REF] Chen | Divergence-measure fields and hyperbolic conservation laws[END_REF][START_REF] Chen | Extended divergence-measure fields and the euler equations for gas dynamics[END_REF], for any 1 < p ≤ ∞, it is possible to define a trace-like operator on the set DM p (Ω D ). Actually, for any F ∈ DM p (Ω D ), we define F • ν the normal trace of F on ∂Ω D , given by F • ν : Lip(∂Ω D ) → R the continuous linear functional such that

F • ν, ξ /∂Ω D = Ω D ξ divF + Ω D ∇ξ • F, for any ξ ∈ C 1 (Ω D ).
Moreover, since Γ D is a deformable Lipschitz boundary, then the restriction of F • ν to Γ D is well defined, this will be denoted by F • ν D (cf. [START_REF] Chen | Some recent methods for partial differential equations of divergence form[END_REF][START_REF] Chen | Divergence-measure fields and hyperbolic conservation laws[END_REF][START_REF] Chen | Extended divergence-measure fields and the euler equations for gas dynamics[END_REF], see also Remark 3.9). Remark 3.9.

(1) Thanks to [START_REF] Chen | Some recent methods for partial differential equations of divergence form[END_REF][START_REF] Chen | Divergence-measure fields and hyperbolic conservation laws[END_REF][START_REF] Chen | Extended divergence-measure fields and the euler equations for gas dynamics[END_REF], for any F ∈ DM p (Ω D ), with 1 < p ≤ ∞, it is possible to define the normal trace F • ν locally by using Lipschitz deformation of the boundary. This formulation is very useful in the case where the boundary is partitioned into disjoint deformable Lipschitz patches.

(2) In the case where p = 1 as well as the case where L p (Ω D ) N is replaced by the space M b (Ω D ) N , the trace may be defined as well, but only as continuous linear form on a subset of Lip(∂Ω D ) (cf. [START_REF] Chen | Some recent methods for partial differential equations of divergence form[END_REF][START_REF] Chen | Divergence-measure fields and hyperbolic conservation laws[END_REF][START_REF] Chen | Extended divergence-measure fields and the euler equations for gas dynamics[END_REF]). Now, combining this consideration with the formal computation for (3.9), we introduce the following optimization problem

(OF D ) : inf φ∈DM p (Ω D ) Ω D σ(x, φ(x))dx -φ • ν D , g : -div (φ) = ρ in D (Ω D \ Γ D ) ,
where the divergence constraint is understood as follows:

Ω D ∇ξ • φ dx -φ • ν D , ξ = Ω D \Γ D ξdρ, for any ξ ∈ Lip(Ω D ),
with φ•ν D being trace-like term on Γ D as defined above. In other words, we impose that φ•ν Γ = 0 on Γ.

Our main result in this section is the following duality result.

F (h) = inf φ∈DM p (Ω D ) Ω D σ(x, φ(x))dx+ Γ D gdh-φ•ν D , g : -div(φ) = ρ+h in D (Ω D \ Γ D ) ,
for any h ∈ M b (Ω D ). Then F is convex and l.s.c. Indeed, take h 1 , h 2 ∈ M b (Ω D ) and set

h := th 1 + (1 -t)h 2 for t ∈ [0, 1]. Let φ 1,n , φ
) = ρ + h i , in D (Ω D \ Γ D ) and F (h i ) = lim n Ω D σ(x, φ i,n (x))dx + Γ D g dh i -φ i,n • ν D , g for i = 1, 2.
Set φ n = tφ 1,n + (1 -t)φ 2,n . We clearly see that φ n are admissible for h and

F (h) ≤ lim n Ω D σ(x, φ n (x))dx + Γ D g dh -φ n • ν D , g = lim n Ω D σ(x, tφ 1,n + (1 -t)φ 2,n )dx + Γ D g dh -(tφ 1,n + (1 -t)φ 2,n ) • ν D , g ≤ lim n t Ω D σ(x, φ 1,n )dx + Γ D g dh -φ 1,n • ν D , g + (1 -t) Ω D σ(x, φ 2,n )dx + Γ D g dh -φ 2,n • ν D , g ≤ tF (h 1 ) + (1 -t)F (h 2 )
and this proves convexity. For the lower semicontinuity, take a sequence h n h in M b (Ω D ). For every n ∈ N, we consider a sequence (φ k n ) k∈N of DM p (Ω D ) such that

F (h n ) = lim k→∞ Ω D σ(x, φ k n (x))dx + Γ D g dh n -φ k n • ν D , g .
We may find some

ψ n ∈ L 1 (Ω D , R N ) such that -div(ψ n ) = h -h n in D (Ω D \ Γ D ), ψ n L 1 → 0 and ψ n • ν D , g → 0. In fact, we have h -h n ∈ M b (Ω D ) → W -1,p (Ω D ) for p > N and p := p p -1
. We consider the following p-Laplace equation

           -∆ p u n = h -h n in Ω D u n = 0 on Γ D |∇u n | p-2 ∇u n • ν Γ = 0 on Γ.
(3.10)

The system (3.10) admits a unique solution u n ∈ W 1,p (Ω D ) such that u n = 0 on Γ D . Hence, if we set ψ n = |∇u n | p-2 ∇u n , we see that ψ n ∈ L p (Ω D ), and then in L 1 (Ω D ). Moreover, we have

-div(ψ n ) = h -h n in D (Ω D ). Since h -h n is bounded in W -1,p (Ω D )
, it is not difficult to prove that u n is bounded in W 1,p (Ω D ). So, by taking a subsequence if necessary, we have u n u in W 1,p (Ω D ), and uniformly in Ω D . On the other hand, we have

Ω D |ψ n | p dx = Ω D |∇u n | p dx = h -h n , u n W -1,p (Ω D ),W 1,p (Ω D ) -→ n→∞ 0.
In particular, this implies that |ψ n | -→ n→∞ 0 in L 1 (Ω D ). Moreover, taking g ∈ Lip(Ω D ) be such that g = g on Γ D , we have

ψ n • ν D , g = Ω D ψ n • ∇g dx - Ω D g d(h -h n ) -→ n→∞ 0.
This being said, we clearly have -div(φ k n + ψ n ) = ρ + h, i.e. φ k n + ψ n are admissible fluxes for h. By semicontinuity of the integral, we have

F (h) ≤ Ω D σ(x, (φ k n + ψ n )(x))dx + Γ D g dh -(φ k n + ψ n ) • ν D , g ≤ Ω D σ(x, φ k n (x))dx + Γ D g dh n -φ k n • ν D , g + Ω D σ(x, ψ n (x))dx + Γ D g d(h -h n ) -ψ n • ν D , g .
Letting k → ∞ we get

F (h) ≤ F (h n ) + Ω D σ(x, ψ n (x))dx + Γ D g d(h -h n ) -ψ n • ν D , g .
Now, letting n → ∞, and using the fact that ψ n → 0 in L 1 (Ω D ) N , and h n h in M b (Ω D ), as n → ∞, we obtain the lower semicontinuity, i.e.

F (h) ≤ lim inf n F (h n ).
Next let us compute F * . For any u ∈ C(Ω D ), we have

F * (u) = sup h∈M b (Ω D ) Ω D udh -F (h) = sup h∈M b (Ω D ) φ∈DM p (Ω D ) Ω D udh - Ω D σ(x, φ(x))dx - Γ D g dh + φ • ν D , g : -div(φ) = ρ + h in D (Ω D \ Γ D ) = I 1 (u) + I 2 (u),
where I 1 (u) := -

Ω D
udρ and

I 2 (u) := sup h∈M b (Ω D ) φ∈DM p (Ω D ) Ω D \Γ D ud(ρ + h) - Ω D σ(x, φ(x))dx + Γ D (u -g) dh + φ • ν D , g : -div(φ) = ρ + h in D (Ω D \ Γ D ) .
Using Lemma 3.11 below, we deduce that, for any u ∈ Lip(Ω D ), we have

F * (u) =          - Ω D udρ if σ * (x, ∇u) ≤ 1 and u = g on Γ D ∞ otherwise.
Finally, using the fact

inf(OF D ) = F (0) = F * * (0) = sup u∈Lip(Ω D ) -F * (u) = max(M D ),
we deduce the result.

Ω D \Γ D ud(ρ + h) - Ω D σ(x, φ(x))dx + Γ D (u -g) dh + φ • ν D , g : -div(φ) = ρ + h in D (Ω D \ Γ D ) =          0 if σ * (x, ∇u) ≤ 1 and u = g on Γ D ∞ otherwise.
Proof. Take u as a test function in the divergence constraint -div(φ

) = ρ + h in D (Ω D \ Γ D ),
we get

I(h, φ) := Ω D \Γ D ud(ρ + h) - Ω D σ(x, φ(x))dx + Γ D (u -g) dh + φ • ν D , g = Ω D ∇u • φdx - Ω D σ(x, φ(x))dx + Γ D (u -g) dh + φ • ν D , g -u .
If σ * (x, ∇u) ≤ 1 and u = g on Γ D , then following Lemma 2.7, we obtain sup I(h, φ) ≤ 0.

Actually, sup I(h, φ) = 0 in this case by taking h ≡ -ρ and φ ≡ 0. If u(x 0 ) = g(x 0 ) for some x 0 ∈ Γ D , then we consider Dirac mass at x 0 , h = n sign(u(x 0 ) -g(x 0 ))δ x 0 for n ∈ N, and fix

Φ 0 ∈ DM p (Ω D ) such that -div Φ 0 = ρ in D (Ω D \ Γ D ), we have sup I(h, φ) ≥ Ω D Φ 0 • ∇u - Ω D σ(x, Φ 0 (x))dx + n|u(x 0 ) -g(x 0 )| + Φ 0 • ν D , g -u .
Letting n → ∞, we get the result. For the remaining case, i.e. u = g on Γ D and σ * (x, ∇u) > 1 on a non negligible set, we see first that, for any u ∈ Lip(Ω D ), there exists a measurable function that we denote by q u : Ω D → R N , such that q u (x) • ∇u(x) = σ * (x, ∇u(x)), a.e. in Ω D .

Indeed, recall that σ * (x, ∇u(x)) = max p { p, ∇u(x) : σ(x, p) = 1} and the function x → p, ∇u(x) + II [σ(x,.)=1] is measurable. Then, q u is given by the measurable representation in the set arg max

p { p, ∇u : σ(x, p) = 1} . Now, if u = g on Γ D and σ * (x, ∇u) > 1 in a subset A ⊂ Ω D such that |A| = 0, we consider Φ n = n q u |q u | χ A * η , where η is a sequence of mollifiers. It is clear that there exists h ∈ M b (Ω D ),
such that -div Φ n = ρ + h. For any n, we have

sup I(h, φ) ≥ Ω D Φ n • ∇u - Ω D σ(x, Φ n (x))dx.
Letting → 0, we get

sup I(h, φ) ≥ n A 1 |q u | (q u • ∇u -σ(x, q u )) dx ≥ n A 1 |q u | (σ * (x, ∇u(x)) -1) dx.
Then, letting n → ∞, we get the result.

Remark 3.12. Going over the duality inferred by Theorem 3.10, we have inf

Ω D σ(x, φ(x))dx -φ • ν D , g : φ ∈ DM p (Ω D ), -div (φ) = ρ in D (Ω D \ Γ D ) = max Ω D u(x)dρ, σ * (x, ∇u(x)) ≤ 1 a.e. x in Ω D and u = g on Γ D .
It is not clear if the inf is a min . This is closely connected to the regularity of the trace of divergence-measure field. However, one sees that if this is true, i.e. the inf is a min, then the respective extremums u and φ satisfy the following PDE:

                           -div(φ) = ρ in D (Ω D \ Γ D ) φ(x) • ∇u(x) = σ(x, φ(x)) in Ω D u = g on Γ D σ * (x, ∇u(x)) ≤ 1 in Ω D φ • ν Γ = 0 on Γ.

The augmented Lagrangian technique

4.1. Formulation of the problem. We set again

X = W 1,∞ (Ω D ) and Y = L ∞ (Ω D ) N .
For any u ∈ X and η ∈ Y , we define

F(u) =        - Ω D u dρ if u = g on Γ D +∞ otherwise , G(η) =    0 if σ * (x, η) ≤ 1 +∞ otherwise
, and Λu = ∇u.

Thus, the problem (M D ) can be rewritten in the form

-min u∈X F(u) + G(Λu).
Introducing a new primal variable q ∈ Y we can write (M D ) in the following alternative form min

(u,q)∈X ×Y Λu=q F(u) + G(q)
so that (M D ) and (OF D ) can be recast in a saddle point form

min (u,q)∈X ×Y sup φ∈DM p (Ω D ) L(u, q; φ)
where

L(u, q; φ) = F(u) + G(q) + Ω D φ • (Λu -q) dx, for any (u, q, φ) ∈ X × Y × DM p (Ω D ).
More precisely, we have Proposition 4.13. u is a solution of (M D ) if and only if the couple (u, q) ∈ X × Y with q := Λu is a solution of sup

φ∈DM p (Ω D ) min (u,q)∈X ×Y L(u, q; φ)
which is equal to min

(u,q)∈X ×Y sup φ∈DM p (Ω D )
L(u, q; φ).

Proof. The duality of the previous section takes part of this result. For completeness let us give here the main arguments. It is not difficult to see that

min F(u) + G(Λu) : u ∈ X = min (u,q)∈X ×Y :Λu=q F(u) + G(q) = min (u,q)∈X ×Y sup φ∈DM p (Ω D ) F(u) + G(q) + Ω D φ (Λu -q) = min (u,q)∈X ×Y sup φ∈DM p (Ω D )
L(u, q; φ).

In a similar way, using the definition of F * , G * and Λ * , one sees that

sup -F * (-Λ * φ) -G * (φ) : φ ∈ DM p (Ω D ) = sup φ∈DM p (Ω D ) min (u,q)∈X ×Y L(u, q; φ).
On the other hand, thanks to the strong duality given in Theorem 3.10 and direct computation for the last line below, we get

-min F(u) + G(Λu) : u ∈ X = max(M D ) = inf Ω D σ(x, φ(x))dx -φ • ν D , g : φ ∈ DM p (Ω D ), -div (φ) = ρ in D (Ω D \ Γ D ) = -sup -F * (-Λ * φ) -G * (φ) : φ ∈ DM p (Ω D ) .
Thus, sup

φ∈DM p (Ω D ) min (u,q)∈X ×Y L(u, q; φ) = min (u,q)∈X ×Y sup φ∈DM p (Ω D )
L(u, q; φ).

For a given r > 0, we recall that the augmented Lagrangian (cf. [START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF]) is given by

L r (u, q; φ) = F(u) + G(q) + φ, Λu -q + r 2 |Λu -q| 2 , for any (u, q, φ) ∈ X × Y × DM p (Ω D ).
In the same way, it can easily be proved that Proposition 4.14. Let r > 0. Then, u is a solution of (M D ) if and only if the couple (u, q) ∈ X × Y with q := Λu is a solution of

sup φ∈DM p (Ω D ) min (u,q)∈X ×Y L r (u, q; φ)
which is equal to (S r ) : min

(u,q)∈X ×Y sup φ∈DM p (Ω D )
L r (u, q; φ). Now for the numerical computation concerning the problem (M D ), we will focus on the saddle point problem (S r ). Recall that the addition of the quadratic term r 2 |Λu -q| 2 has the advantage of improving the convergence of the dual approach (one can see [START_REF] Ekeland | Convex analysis and variational problems[END_REF][START_REF] Glowinski | Numerical methods for nonlinear variational problems[END_REF]).

4.2. Application of ALG2. We approximate the domain Ω D via a triangulation T h . For k ≥ 1, we denote by P k the space of polynomials with real coefficients and of degree at most k. We define

X h ⊂ W 1,∞ (Ω D )
as the space of continuous functions on Ω D belonging to P k on each triangle. Similarly, Y h is the space of vector valued functions belonging to (P k-1 ) d on each triangle. Then the problem (M D ) is discretized by the following finite-dimensional optimization problem:

inf u∈X h F(u) + G(Λu).
Then, for a given q 0 ∈ Y h , φ 0 ∈ Y * h , using ALG2 algorithm (also called ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]) we construct a sequence {u i } i , {q i } i , {φ i } i by optimizing alternatively in u and q, for i ≥ 1. Convergence of the algorithm in the finite-dimensional setting is well known (cf. [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF]). The convergence can be also established in the case of infinite-dimensional Hilbert spaces. For further details in this direction and about the ALG2, we refer the reader to [START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF][START_REF] Fortin | Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems[END_REF][START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF].

Algorithm 1 ALG2 iterations 1st step: u i+1 ∈ arg min u∈X h F(u) + φ i , Λ(u) + r 2 |Λ(u) -q i | 2 .
2nd step: q i+1 ∈ arg min

q∈Y h G(q) -φ i , q + r 2 |Λ(u i+1 ) -q| 2 .
3rd step: We update the multiplier φ by

φ i+1 = φ i + r(∇u i+1 -q i+1 ).
Remark 4.15.

• The first step amounts to solve a Laplace equation with mixed boundary conditions. Indeed, we have for every z ∈ X h with z = 0 on Γ D r ∇u i+1 , ∇z = ρ, z + (rq i -φ i ), ∇z which is equivalent to solve • The second step is a pointwise projection. Indeed, we use P 1 finite element for q and φ, we have at each vertex

-r∆u = ρ + div(φ -rq) in Ω D ,
x k 0 ∈ ∂II B * (x k ,.) (q i+1 (x k )) -φ i (x k ) -r(∇u i+1 (x k ) -q i+1 (x k )),
which is equivalent to perform pointwise projections:

q i+1 (x k ) = Proj B * (x k ,.) φ i (x k ) r + ∇u i+1 (x k )
where B * (x, .) = {p ∈ R N : σ * (x, p) ≤ 1}.

Error criterion.

Basing on the primal-dual optimality conditions, we use the following stopping criterion

1. MaxLip := sup x σ * (x, ∇u(x)), 2. Div := -div(φ) -ρ L 2 , 3. Dual := σ(x, φ(x)) -φ(x) • ∇u L 2 , 4. NBD φ := Γ(φ • ν) 2 1/2
.

We expect MaxLip ≤ 1 and Div, Dual, NBD φ to be small. In addition, we compute u -u exact for different norms where u is the computed solution and u exact is the exact solution, whenever the latter is easily found. Let us mention that NBD φ will concern only the Test 2 where we prescribe data on a closed set D, so that Γ = ∅.

To implement the algorithm we use FreeFem++ [START_REF] Hecht | New development in Freefem++[END_REF], which is particularly adapted to solve the Laplace equation in the first step of ALG2. We use P 2 finite element for u and P 1 finite element for φ and q (see e.g. [START_REF] Benamou | Augmented Lagrangian methods for transport optimization, mean field games and degenerate elliptic equations[END_REF]). All the tests are executed on a macOs Mojave 10.14.4. The exact solution is u exact (x, y) = d((x, y), ∂Ω) = min(x, 1 -x, y, 1 -y). We measure u -u exact in different norms, with different mesh sizes N h and for 1000 iterations. This test shows us that ALG2 iteration converges to the accurate solution u as the mesh is refined.

Numerical experiments

N h Time execution u -u exact L 2 u -u exact L 1 u -u exact L ∞ 30

Test 2:

We consider the same equation |∇u| = 1 on Ω = (-1, 1) 2 \ D where D is the euclidean ball centered at (0.05, 0.09) and of radius 0.25, and we set u to be equal to 0 on D.

(a) (a) In this case B * is an ellipse and the projection can be computed as in the Riemannian case [START_REF] Yu | Algorithms of projection of a point onto an ellipsoid[END_REF]:

N h Time execution u -u exact L 2 u -u exact L 1 u -u exact L ∞ 30
Proj B * (q) = q if q ∈ B * (ζM + I n ) -1 q otherwise,
where ζ is the unique positive root of the function

F (ζ) = q † ζ M qζ -1 with qζ = (ζM + I n ) -1
q, which can be found with a dichotomy method. We perform a test as in [START_REF] Luo | Fast sweeping methods for factored anisotropic Eikonal equations: multiplicative and additive factors[END_REF]Example 2] by taking Given a compact convex set K of R N such that 0 ∈ int(K), then its gauge function reads g K (p) = inf{λ ≥ 0; λp ∈ K}.

M = l 1 (x, y) -l 3 (x, y) -l 3 (x, y) l 2 (x, y) with l 1 (x, y) = 1 e -2
We consider the following HJ equation H(x, ∇u) = g K * (∇u) -1, in Ω and u = 0 on ∂Ω, where K * is the polar set of K defined through K * = {p ∈ R N : p, q ≤ 1, ∀q ∈ K}.

It is well known (see e.g. [START_REF] Rockafellar | Convex analysis[END_REF][START_REF] Hiriart-Urrut | Convex analysis and minimization algorithms[END_REF]) that g K * = σ K . In the last two examples we take some explicit forms of the support function σ.

5.7. Riemannian case. We take σ((x, y), v) = β 1 v 2 1 + β 2 v 2 2 , with β 1 , β 2 > 0. It is not difficult to see that σ * (q) = q 2 1 β 1 + q 2 2 β 2 . As in [START_REF] Benamou | A numerical solution to Monge's problem with a Finsler distance as cost[END_REF], the projection onto the unit ball of σ * , B * = {q, σ * (q) ≤ 1} is given by

Proj B * (q) =    q if q ∈ B * β 1 q 1 β 1 + ζ , β 2 q 2 β 2 + ζ otherwise,
where ζ is the zero of the function

F (ζ) = 1 - β 1 q 1 (β 1 + ζ) 2 + β 2 q 1 (β 2 + ζ) 2 ,
which can be computed with a dichotomy algorithm.

For this test, we take .

(a) The projection onto B * can be performed easily (see [START_REF] Benamou | A numerical solution to Monge's problem with a Finsler distance as cost[END_REF][START_REF] Igbida | Augmented Lagrangian method for optimal partial transportation[END_REF]). We start by determining the vertices s 1 , • • • , s k of B * and the corresponding outward normal vectors ν i to the edges of B * . Afterwards, if v ∈ B * , we distinguish to cases: either v ∈ [s i , s i+1 ] + R + ν i , and in this case we project v onto the segment [s i , s i+1 ], or it belongs to a sector s i + R + ν i + R + ν i+1 and in this case its projection is s i . We perform a test with p 1 = (1, -1), p 2 = (1, -0.8), p 3 = (-0.8, 1), p 4 = (-1, 1), p 5 = (-1, -1). We take N h = 64 and 600 iterations.

1 . Introduction 1 . 1 .

 111 Motivation and related works. Our main interest lies in the study from theoretical and numerical point of view the Hamilton-Jacobi (HJ) equation H(x, ∇u) = 0 in an open bounded domain Ω ⊆ R N coupled, for instance, with a Dirichlet boundary condition on ∂Ω.

  which can be rewritten as ∇u ∈ B |.| (0, f ) where |.| stands for the Euclidean norm. This enables to rewrite the problem into the following sub-gradient constrained optimization max Ω udx : ∇u ∈ B |.| (0, f ) and u = 0 on ∂Ω .

  hence by Proposition 2.5 we have σ * (x, ∇u) ≤ 1. Consequently u solves (2.8). If w is another solution, then Ω u(x)dx = Ω w(x)dx and w ≤ u in Ω by the first step. Consequently u = w, as desired.

Theorem 3 . 8 .

 38 (Strong duality [20, Chap. III]) Assume moreover that there exists u 0 ∈ X such that F(u 0 ) < +∞, G(Λu 0 ) < +∞ and G is continuous at Λu 0 . Then the dual problem admits at least a solution v ∈ Y * and the strong duality holds, i.e.

  sup Ω udx : |∇u| ≤ f and u = 0 ∂Ω , or equivalently (P) : inf u∈V {F(u) + G(Λu)}

3. 2 .

 2 Duality for HJ equation. To the duality, we consider a more general problem by considering, for a nonnegative Radon measure ρ ∈ M + b (Ω) and a closed subset D ⊂ Ω, the problem max Ω udρ : σ * (x, ∇u(x)) ≤ 1 and u = g on D .

Theorem 3 . 10 .

 310 Let ρ ∈ M + b (Ω). Then we have the following strong duality (M D ) := max Ω D udρ : σ * (x, ∇u(x)) ≤ 1 and u = g on Γ D = inf φ∈DM p (Ω D ) Ω D σ(x, φ(x))dx -φ • ν D , g : -div (φ) = ρ in D (Ω D \ Γ D ) := (OF D ). Proof. Consider on M b (Ω D ) the following functional F : M b (Ω D ) →] -∞, ∞] defined by

Lemma 3 . 11 .

 311 Let u ∈ Lip(Ω D ), we have sup h∈M b (Ω D ) φ∈DM p (Ω D )



  together with the following mixed boundary conditions:   Dirichlet boundary condition u = g on Γ D Homogeneous Neumann boundary condition (r∇u + φ -rq) • ν Γ = 0 on Γ.

5. 1 .

 1 Test 1: We first examine the case |∇u| = f (x, y) ≡ 1 in Ω = (0, 1) 2 and u = 0 on ∂Ω.

Figure 1 .

 1 Figure 1. Left to right: 3D plot of the solution u, contour plot of u, the flux φ

Figure 2 .

 2 Figure 2. Error criterion for 1000 iterations with N h = 120

Figure 3 .

 3 Figure 3. Left to right: 3D plot of the solution u, contour plot of u, the flux φ

Figure 4 . 5 . 3 .

 453 Figure 4. Error criterion for 500 iterations with N h = 50

Figure 5 .Figure 6 . 5 . 4 . 4 :

 56544 Figure 5. Left to right: 3D plot of the solution u, contour plot of u, the flux φ

Figure 7 .Figure 8 . 5 . 5 .

 7855 Figure 7. Left to right: 3D plot of the solution u, contour plot of u, the flux φ

Figure 9 .

 9 Figure 9. Left to right: 3D plot of the solution u, contour plot of u, the flux φ

Figure 10 .

 10 Figure 10. Error criterion for 600 iterations with N h = 120

√ 2 . 2 l 2

 222 (x-0.5) 2 +2.(x-0.5).(y-0.5)+(y-0.5) (x, y) = 2l 1 (x, y) l 3 (x, y) = l 1 (x, y).

Figure 11 .

 11 Figure 11. Left to right: 3D plot of the solution u, contour plot of u, the flux φ

Figure 12 .

 12 Figure 12. Error criterion for 300 iterations with N h = 80.

2 e - 2 √ 2 (

 222 -0.5) 2 +2(x-0.5)(y-0.5)+(y-0.5) 2 ) and β 2 = x-0.5) 2 +2(x-0.5)(y-0.5)+(y-0.5) 2 )

Figure 13 .

 13 Figure 13. Left to right: 3D plot of the solution u, contour plot of u, the flux φ

Figure 14 .

 14 Figure 14. Error criterion for 400 iterations with N h = 64

Figure 15 .

 15 Figure 15. Left to right: 3D plot of the solution u, contour plot of u, the flux φ

Figure 16 .

 16 Figure 16. Error criterion for 600 iterations with N h = 64.

Acknowledgements. The authors are grateful to the anonymous referees for carefully reading this paper and for the interesting remarks and suggestions. Research of the last author is funded by Vietnam National Foundation for Science and Technology Development (NAFOS-TED) under grant number 101.01-2018.309. His research has been partially carried out during a visit at VIASM, he would like to thank the institution for supporting the visit and hospitality.


Rami Klaimi

Charbel Stefan Weithoffer

Abdel Nour

Catherine Douillard

Simplified recursion units for Max-Log-MAP: New trade-offs through variants of Local-SOVA

Keywords: Max-Log-MAP decoding, local-SOVA algorithm, low-complexity decoding, high-throughput turbo decoders

The Log-domain BCJR algorithm is broadly used in iterative decoding processes. However, the serial nature of the recursive state metric calculations is a limiting factor for throughput increase. A possible solution resorts to high-radix decoding, which involves decoding several successive symbols at once. Despite several studies aiming at reducing its complexity, high-radix processing remains the most computationally intensive part of the decoder when targeting very high throughput. In this work, we propose a reformulation specifically targeting the complexity reduction of the recursive calculation units by either limiting the required number of operations or by selectively removing unnecessary ones. We report a complexity reduction of the add-compare-select units in the order of 50% compared to the recently proposed local-SOVA algorithm. In addition, our results show that several performance/complexity trade-offs can be achieved thanks to the proposed simplified variants. This represents a promising step forward in order to implement efficient very high throughput convolutional decoders.

I. INTRODUCTION

The Bahl, Cocke, Jelinek, and Raviv (BCJR) algorithm [START_REF] Bahl | Optimal decoding of linear codes for minimizing symbol error rate (corresp.)[END_REF] or its variants are widely used to compute the maximum a posteriori (MAP) estimate of transmitted symbols in iterative decoders or detectors [START_REF] Berrou | Near Shannon limit error-correcting coding and decoding: Turbo-codes[END_REF]. Recently, a number of studies have focused on fully parallel or fully pipelined hardware implementations of the Max-Log-MAP (MLM) algorithm [START_REF] Robertson | A comparison of optimal and sub-optimal MAP decoding algorithms operating in the log domain[END_REF] to achieve a throughput of tens of Gb/s [START_REF] Maunder | A fully-parallel turbo decoding algorithm[END_REF]. The price to pay to reach such high-throughput levels is a significant increase in area due to the computational complexity of MLM decoding [START_REF] Weithoffer | Fully pipelined iteration unrolled decoders-The road to Tb/s turbo decoding[END_REF]. To address this detrimental effect, the local-SOVA algorithm (LSOVA) was introduced in [START_REF] Le | Revisiting the Max-Log-Map algorithm with SOVA update rules: new simplifications for high-radix SISO decoders[END_REF]. It relies on a new low-complexity soft-output calculation unit (SOU) that applies a path-based decoding variant of the MLM algorithm. This algorithm was adapted for a hardware implementation in [START_REF] Weithoffer | Low-complexity computational units for the local-SOVA decoding algorithm[END_REF] and was shown to significantly reduce the hardware complexity compared to a standard MLM decoder. Therefore, the resulting increased area efficiency enables higher throughput implementations of turbo decoders.

While [START_REF] Le | Revisiting the Max-Log-Map algorithm with SOVA update rules: new simplifications for high-radix SISO decoders[END_REF] and [START_REF] Weithoffer | Low-complexity computational units for the local-SOVA decoding algorithm[END_REF] focused on the SOU, this work targets the serial nature and complexity of the add-compare-select (ACS) units that calculate the forward and the backward recursions. In this sense, the results presented in the remainder of this paper are complementary to [START_REF] Le | Revisiting the Max-Log-Map algorithm with SOVA update rules: new simplifications for high-radix SISO decoders[END_REF], [START_REF] Weithoffer | Low-complexity computational units for the local-SOVA decoding algorithm[END_REF]. Increasing the radix order for the ACS units, i.e. processing several trellis sections at the same time, lowers the latency but is associated with an increase in area complexity and a long critical path [START_REF] Kim | A unified parallel radix-4 turbo decoder for mobile WiMAX and 3GPP-LTE[END_REF], [START_REF] Sánchez | High speed low complexity radix-16 Max-Log-MAP SISO decoder[END_REF].

Hence, reducing the number of ACS units without increasing the complexity of the performed computations remains an open question. This problem was nevertheless addressed in [START_REF] Klaimi | Low-complexity decoders for non-binary turbo codes[END_REF] in the context of non-binary turbo decoders, where a modified bubble sorter was shown to reduce significantly the required number of ACS units. Corresponding results confirmed that the decoder complexity cannot be lowered without reducing the complexity of the recursions.

In this paper, we develop two different LSOVA variants that reduce the complexity of the recursive calculations. While these simplifications do apply on simple trellis-based decoding or detection algorithms, in this work we take turbo decoding as an example to assess the effect of the proposals on the decoding performance. We report significant complexity reductions in comparison with the LSOVA algorithm in [START_REF] Le | Revisiting the Max-Log-Map algorithm with SOVA update rules: new simplifications for high-radix SISO decoders[END_REF], [START_REF] Weithoffer | Low-complexity computational units for the local-SOVA decoding algorithm[END_REF]. Proposed variants enable further reductions in area complexity and in latency of turbo decoder hardware implementations.

This paper is structured as follows: Section II gives a background overview of the reference decoding algorithms. In Section III, we propose two complementary variants of these algorithms: Section III-A explains how to reduce the complexity of the recursive calculations while Section III-B proposes a technique to compute all the soft decisions in a single step, thus reducing the decoding latency. Performance comparisons in terms of error correction capability and computational complexity are provided in Section IV. Finally, Section V concludes the paper.

II. REFERENCE DECODING ALGORITHMS

A. The Max-Log-MAP decoding algorithm

In this section, we briefly recall the equations of the MLM decoding algorithm [START_REF] Robertson | A comparison of optimal and sub-optimal MAP decoding algorithms operating in the log domain[END_REF]. The forward and the backward metrics and for state at trellis step are computed using the branch metrics as:

( ) = max ∈S ( -1 ( ) + -1 ( , )) (1)
( ) = max ∈S ( +1 ( ) + ( , )) (2)
S and S are the sets of trellis states at steps -1 and + 1, respectively, that are connected to state at trellis step .

In addition, the soft output is calculated as a log-likelihood ratio (LLR) at step or time index :

L = max ( , ) ∈{0•••2 } 2 | ( , )=0 ( ( ) + +1 ( ) + ( , )) - max ( , ) ∈{0•••2 } 2 | ( , ) =1 ( ( ) + +1 ( ) + ( , )) (3)
where is the memory of the considered convolutional code and ( , ) is the systematic bit that connects state at trellis step to state at trellis step + 1, if the connection exists. When aiming for high throughput, several trellis sections can be decoded together. Precisely, a radix-2 trellis allows the decoding of bits at each decoding step. Equations ( 1), ( 2) and ( 3) have then to be adapted accordingly. Additionally, for high-throughput decoding, the data frames are in practice decomposed into windows of size , processed in parallel.

B. The local-SOVA decoding algorithm

The LSOVA is a path-based variant of the MLM algorithm introduced in [START_REF] Le | Revisiting the Max-Log-Map algorithm with SOVA update rules: new simplifications for high-radix SISO decoders[END_REF]. It uses the same recursive metrics calculations (1) and ( 2) but the LLRs are computed using path-based local update rules inspired by the SOVA. Conventionally, a path in a trellis diagram is defined as a sequence of states and is associated with an input bit sequence and a path metric. In [START_REF] Le | Revisiting the Max-Log-Map algorithm with SOVA update rules: new simplifications for high-radix SISO decoders[END_REF], an alternative definition of a path, with a more local sense that focuses on a particular trellis section is adopted. When considering a radix-2 trellis, for each trellis section, a path is defined for each state as a 3tuple consisting of the path metric , the hard decisions = { ,0 , . . . , , -1 } labeling this path, and their associated reliability values = { ,0 , . . . , , -1 }. Two paths and can be merged into path according to the following:

= 0 ( , ) = max( , ) (4)
( ) = 1 ( ( ), ( ))∀ ∈ [0, -1] (5)
where 1 selects the hard decision of the winning path resulting from (4). The reliabilities are updated following:

( ) = 2 ( ( ), ( ))∀ ∈ [0, -1] (6)
with 2 being the Hagenauer rule (HR) [START_REF] Hagenauer | A Viterbi algorithm with soft-decision outputs and its applications[END_REF] or the Battail rule (BR) [START_REF] Battail | Pondération des symboles décodés par l'algorithme de Viterbi[END_REF] depending on ( ) = ( ) or ( ) ≠ ( ). Let = arg max( , ), = arg min( , ) and Δ = ( ) -( ):

( ) = min( ( ), Δ + ( )) if ( ) = ( ) (BR) min( ( ), Δ) if ( ) ≠ ( ) (HR) (7)
A merge and update tree is traversed to get the final hard and soft outputs.

III. LOWERING THE COMPLEXITY OF THE FORWARD AND BACKWARD RECURSIONS

Forward and backward recursions are performed serially for the trellis stages within a window. Corresponding computations require a proportional amount of ACS units, depending on the radix order , the number of states 2 and the targeted frame size . For large and , this number becomes prohibitive for fully pipelined/parallel hardware architectures, motivating the study of alternative algorithms with a reduced number and/or simplified ACS units. To address this problem, we propose two low-complexity ACS computation variants.

Fig. 1 shows the conventional computation of the forward and backward recursions launched simultaneously from the edges of a size-window in a radix-4 trellis. The straight lines correspond to the winning transitions in ( 1) and ( 2), while the dotted lines correspond to the losing transitions. Once both recursions meet in the middle (MIM) of the window, i. e. when + /2 and + /2 are computed, two possible strategies can be applied, corresponding to the two proposed simplified variants. The first strategy applies lowcomplexity recursion units beyond + /2 to return to the edges of (down to for and up to + for ). This strategy is denoted by low complexity return (LCR). The second strategy halts the execution of the recursion units and applies LSOVA over a unique merge tree to update the soft output (reliability) values for the /2 trellis steps of the forward and backward recursions jointly. Indeed, the decoding decision of MLM at + /2 is the same regardless from the path direction (forwards from to + /2 or backwards from + to + /2). This second strategy is called modified dualsided (MDS) LSOVA. In the remainder of this paper, since the processing is limited to a window, index is omitted ( = 1) for simpler notations.

A. Meet in the middle and low-complexity return

Following Fig. 2, the metric difference between transition and the winning transition converging to state at trellis section is denoted by ( ) for the forward recursion, whereas the metric difference between this same transition and the one that diverges from state is denoted by Δ ( ) for the backward recursion. For each trellis stage while performing to the MLM algorithm where values of ( ) and ( ) are computed accurately and separately for the trellis stages, an estimate of ( ) = ( ) + ( ) is computed recursively in our LCR proposal. By considering solely winning transitions from the forward and backward recursions before MIM, this proposal largely reduces computational complexity.

Low-complexity backward computation beyond MIM: Continuing with the same example, let us consider the calculation of ( ), ∈ [0, 3] from the four +2 ( ), ∈ [0, 3] values beyond MIM. The corresponding four transitions per state were already considered, with one labeled as winning in the forward metric calculation of +2 , before MIM. This is illustrated in the example of Fig. 3. Now for the same transitions, considering the opposite direction from the viewpoint of the states at trellis step , three labelling configurations can be observed: denoting by ( ) the number of winning transitions emerging from state at position , we have (0 When ( ) > 0, ( ) can be directly deduced exclusively from +2 . When ( ) = 0, a readjustment calling upon memorized ( ) and Δ ( ) values, representing the difference between the metrics of currently considered and the winning transitions, is needed to estimate ( ):

) = 0, (1) = (3) = 1 and (2) = 2.
( ) =        +2 ( 1 ) if ( ) = 1 max( +2 ( 1 )) if ( ) > 1 max( +2 ( 2 )) -( 2 ) if ( ) = 0 (8)
where 1 is the set of winning transitions emerging from state during the calculation of +2 and 2 is the set of all the transitions emerging from state .

Low-complexity forward computation beyond MIM: Similarly, at the right side of the MIM calls upon transitions labeled as winning in the backward metric calculation of -2 :

( ) =        -2 ( 1 ) if ( ) = 1 max( -2 ( 1 )) if ( ) > 1 max( -2 ( 2 )) -Δ ( 2 ) if ( ) = 0 (9)
Using ( 8) and [START_REF] Sánchez | High speed low complexity radix-16 Max-Log-MAP SISO decoder[END_REF], is estimated at each position in the sub-trellis and is used in the LSOVA SOU to compute the soft output values. With this proposal, we slightly increase the memory requirements of the decoder, while reducing the performed operations after MIM: the number of comparisons needed for the metric computations is greatly decreased. Moreover, no branch metric memorization/re-computation is required. These reductions are expected to lead to significant improvements of both the complexity and efficiency of resulting hardware architectures. To sum up, the proposed reformulation of the MLM recursive units noticeably reduces the computational complexity while still suffering from the serial nature of these computations.

B. Meet in the middle and modified dual-sided LSOVA

Contrary to the LCR approach, the proposed MDS-LSOVA algorithm avoids the computation of the recursive metrics after MIM, without significant impact on the error correcting performance. To do so, we compute ( ) and ( ) using paths = { , , } and = { , , } resulting from the backward and the forward recursion computations for each state at time index before MIM. A first approach to carry out these computations in one shot is by considering a radix-2 /2 LSOVA for the computation of + /2 and + /2 . Hence, no path merge occurs before MIM. Then, only one LSOVA SOU merge tree is required to update both sides of MIM while avoiding recursive computations after MIM.

However, the computational complexity of a recursive metric or increases exponentially with the number of trellis sections within a radix. With = /2, the complexity becomes prohibitive for a practical implementation with typical values of (32 or 64). As an alternative, we propose to use a lower radix order (typically 4, i. e. = 2) for the recursions before MIM, while still applying only one SOU at MIM. This requires modifications to the LSOVA algorithm due to the merges that occur after each lower-order radix during the recursions in this case. Indeed, if we consider for example a forward recursion with radix 4, it happens that a state at position does not lead to any winning transition for the merge at + 2 (i. e. ( ) = 0) as seen in Fig. 3. Then, the forward metric accumulated (before ) up to this state will fail to reach the merge tree of the SOU at MIM and will not be considered for the computation of the final reliability values, entailing large performance penalties. To address this issue, we propose to update the reliability values for all previous trellis steps during the computation of the recursion metrics for each new radix-2 trellis segment. In the case of forward recursions, let 1 ( ), 2 ( ), . . . 2 ( ) be the paths at trellis position to be merged to compute ( + ) at trellis position + . Then, + = is computed from , and + bits are to be updated according to [START_REF] Le | Revisiting the Max-Log-Map algorithm with SOVA update rules: new simplifications for high-radix SISO decoders[END_REF]. The merge operation computes

= 0 ( 1 , . . . , 2 ) = max( 1 , . . . , 2 ) (10)
( ) = 1 ( 1 ( ), . . . , 2 ( ))∀ ∈ [1, + ]. (11)
The reliability values are updated through

( ) = 2 ( 1 ( ), . . . , 2 ( ))∀ ∈ [1, + ] (12)
where 1 and 2 are defined as in ( 5) and [START_REF] Le | Revisiting the Max-Log-Map algorithm with SOVA update rules: new simplifications for high-radix SISO decoders[END_REF]. Similar merge operations are to be performed during the backward recursions. At MIM, the 2 length-/2 winning paths of the forward and backward recursions meet. The result of the addition of their metrics and is provided as input to a LSOVA SOU tree involving 2 -1 merge operations, generally organized in layers [START_REF] Le | Revisiting the Max-Log-Map algorithm with SOVA update rules: new simplifications for high-radix SISO decoders[END_REF]. Each merge operation involves updating the reliability values of the hard decisions labeling the merged paths on both sides of MIM, due to symmetry between the forward and backward processes.

The number of performed updates has a great impact on computational complexity. Aiming to reduce it, we studied the effect of limiting the update depth. This was motivated by the well-known feature of convolutional codes that all trellis paths merge to the maximum likelihood path after some trellis steps. When considering no limitation, the number of reliability updates to be performed for each trellis state when processing trellis position in the forward recursion, ( ), is . Therefore, when considering radix-2 processing for the recursions, the total number of updates performed during the forward recursions is

upd = 2 2 =1 upd ( ) = 2 2 2 2 + 1 (13)
Due to symmetry, the total number of updates performed during the backward recursions upd is also equal to upd . In addition, the total number of updates in the SOU is

SOU upd = (2 -1) (14)
In summary, MDS-LSOVA replaces half of the recursive computations by updates and can thus largely limit latency since the necessary update operations can be highly parallelized and pipelined. Moreover, limiting the number of updates can have a significant impact on complexity and error correction performance. Several trade-offs were assessed when limiting the value of to different fixed values in (13) as well as the number of updates in [START_REF] Bickerstaff | A 24Mb/s radix-4 logMAP turbo decoder for 3GPP-HSDPA mobile wireless[END_REF], depending on the layer , ∈ [1, ] . They are illustrated in the following example.

IV. CASE STUDY: SIMULATION RESULTS AND COMPLEXITY ANALYSIS

In this section, we consider a turbo code using tail-biting = 3 recursive systematic convolutional code with generator polynomials (1, 15/13) 8 as component code. A frame size = 128 bits, window size = 32 and an almost regular permutation designed following [START_REF] Garzón-Bohórquez | Protographbased interleavers for punctured turbo codes[END_REF] are used. Table I shows the different configurations for the truncation of the updates. SOU upd ( ), ∈ [START_REF] Bahl | Optimal decoding of linear codes for minimizing symbol error rate (corresp.)[END_REF][START_REF] Robertson | A comparison of optimal and sub-optimal MAP decoding algorithms operating in the log domain[END_REF] is the number of updates performed by the SOU for each merge operation in layer (instead of ). As explained above, forward and backward recursions stop at MIM for MDS-LSOVA. Therefore, the starting and ending positions of the decoding windows have to be shifted between the successive iterations of the decoding process. This can be easily done with a tail-biting code. The starting and ending positions used for each window in our simulation setup are shown in Table II. For initialization, we used the metric values at MIM as initial values for the next iteration, as in the next iteration initialization technique [START_REF] Bickerstaff | A 24Mb/s radix-4 logMAP turbo decoder for 3GPP-HSDPA mobile wireless[END_REF]. Fig. 4 shows the error correction performance in terms of bit error rate (BER) of the proposed algorithms in AWGN channel using BPSK modulation. Two coding rates, = 1/3 and = 8/9, were simulated. For both rates, the LCR algorithm and the full-complexity MDS-LSOVA perform within 0.2 dB of the reference LSOVA, which is strictly equivalent to the MLM algorithm. Observed penalties are due to the simplified recursion metrics in the case of LCR and to the limited propagation of the forward and backward metrics in a window (only up to /2 instead of ) during one iteration for the fullcomplexity MDS-LSOVA. Regarding the configurations with a limited update length, performance is gradually degrading from 7 to 1 , up to 1 dB compared to full-complexity MDS-LSOVA performance. To choose between configurations, complexity must be taken into account.

The complexity of the different algorithms whose BER performance is shown in Fig. 4 is provided in Table III. It is expressed in number of ACS units and updates in the case of radix-4 processing ( = 2). Before MIM, the LSOVA and LCR perform 2 (2 -1) 2 = 192 ACS operations and updates for each recursion, forward or backward. For the MDS-LSOVA, the same amount of ACS operations is required but the updates are carried out on a certain amount of previous trellis steps, depending on the truncation length. After MIM, there are no more update operations in the recursions and the number of ACS computations is the same as before MIM for the reference LSOVA, while it is statistically halved for the LCR algorithm according to (8) and ( 9). In the LSOVA and LCR algorithms, = 16 SOUs are required, each SOU performing 2 -1 = 7 ACS operations, each comparison being followed by = 2 updates (one for each decoded bit). The MDS-LSOVA no longer needs any recursive computation after MIM and requires only one SOU for the whole window, involving 2 -1 = 7 ACS operations and (2 -1) = 224 updates for the full-complexity algorithm. The complexity of the simplified variants of the MDS-LSOVA is computed according to the parameters given in Table I.

All the proposed algorithms reduce the number of ACS units compared to the reference LSOVA. The LCR and MDS-LSOVA algorithms reduce the number of ACS operations by more than 20% and 50%, respectively, at the price of an increased number of updates in some cases. All the simplified MDS-LSOVA configurations require a lower total number of operations than LSOVA. Therefore, the choice of the appropriate configuration is a compromise between error correction performance and complexity/latency reduction. For example, the LCR decoder offers a significant reduction in number of ACS units while maintaining the same latency as the LSOVA decoder. Differently, MDS-7 performs close to the reference algorithm, with an expected halved latency but with only a slight reduction in complexity.

V. CONCLUSION

We proposed several variants of the LSOVA algorithm with several performance/complexity/latency trade-offs. The LCR algorithm reduces the complexity of the ACS units by more than 20% for a penalty of up to 0.2 dB depending on the code rate. The second variant goes further to halve the number of ACS units at the cost of more update operations for the soft output computation in certain cases. The decoder designer is offered then a much wider panel of trade-offs compared to prior art. In fact, recent works on LSOVA showed that pathbased complexity reductions directly translate to reductions in area complexity for hardware implementation [START_REF] Weithoffer | Low-complexity computational units for the local-SOVA decoding algorithm[END_REF]. The results of this paper confirm that path-based variants of the BCJR algorithm have a high potential in achieving low latency decoding, while guaranteeing high throughput.

Fig. 1 :

 1 Fig. 1: Windowed MLM decoding in a radix-4 trellis ( = 2).

Fig. 2 :

 2 Fig. 2: ( ) and Δ ( ) for a radix-4 trellis. Straight and dotted lines correspond to winning and losing transitions. the recursive calculations until MIM, we propose to memorize ( ) and Δ ( ) instead of the intermediate result of the forward and the backward recursions. Beyond MIM, in contrast

Fig. 3 :

 3 Fig. 3: Example of winning (bold) and loosing (dotted) transitions for the backward computation of beyond MIM.

Fig. 4 :

 4 Fig.4: BER comparison of a turbo code using proposed and reference decoding algorithms.= 128 bits, = 32, code rates = 1/3, = 8/9 and 8 iterations. AWGN channel and BPSK modulation.

TABLE I :

 I Proposed MDS-LSOVA update truncation configurations.

	Conf.	Upd.	upd	SOU upd (1)	SOU upd (2)	SOU upd (3)
	1		4	4	4	4
	2		8	8	8	8
	3		4	4	4	16
	4		8	8	8	16
	5		4	8	8	16
	6		4	16	16	16
	7		8	16	16	16

TABLE II :

 II Window boundaries for odd/even iter. ( = 128, = 32).

	Iteration	Win	Win1	Win2	Win3	Win4
	odd		[0, 31]	[32, 63]	[64, 95]	[96, 127]
	even		[16, 47]	[48, 79]	[80, 111]	[112, 15]

TABLE III :

 III Complexity comparison of the considered decoding algorithms in number of ACS and update operations required to decode one window. Radix 4 ( = 2), = 32.

ACKNOWLEDGMENT

This work was partially funded by the French National Research Agency TurboLEAP project (ANR-20-CE25-0007).