L Brice

D Kalise
email: dkaliseb@ic.ac.uk.

Z Kobeissi
email: zkobeissi@math.univ-paris-diderot.fr.

M Lauri Ère

Á Mateos

Gonz Ález

F J Silva

ON THE IMPLEMENTATION OF A PRIMAL-DUAL ALGORITHM FOR SECOND ORDER TIME-DEPENDENT MEAN FIELD GAMES WITH LOCAL COUPLINGS

, whose core consists in iteratively solving linear systems and applying a proximity operator. We apply that method to timedependent MFG and, for large viscosity parameters, we improve the linear system solution by replacing the direct approach used in [14] by suitable preconditioned iterative algorithms.

Introduction

In this work we consider the following MFG system with local couplings (MFG)

       -∂ t u -ν∆u + H(x, ∇u) = f (x, m(x, t)) in T d × [0, T ], ∂ t m -ν∆m -div(∇ p H(x, Du)m) = 0 in T d × [0, T ], m(•, •) = m 0 (•), u(•, T ) = g(•, m(•, T )) in T d .
In the notation above ν > 0, d ∈ N, T d is the d-dimensional torus, H : T d × R d → R is jointly continuous and convex with respect to its second variable, f , g : T d × R → R are continuous functions and m 0 : T d → R satisfies m 0 ≥ 0 and T d m 0 (x)dx = 1. System (MFG) has been introduced by J.-M. Lasry and P.-L. Lions in [START_REF] Lasry | Jeux à champ moyen II. Horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Mean field games[END_REF] in order to describe the asymptotic behaviour of symmetric stochastic differential games as the number of players tends to infinity. Several analytical techniques can be used to prove the existence of solutions to (MFG) under various assumptions on the data. Even if system (MFG) has been introduced recently, there are already too many references about its theoretical study to be exhaustive in this introduction and we content ourselves in referring the interested reader to the monographs [START_REF] Bensoussan | Mean field games and mean field type control theory[END_REF][START_REF] Gomes | Regularity theory for mean-field game systems[END_REF], the surveys [START_REF] Cardaliaguet | Notes on Mean Field Games: from P.-L. Lions' lectures at Collège de France[END_REF][START_REF] Gomes | Mean field games models -a brief survey[END_REF] and the references therein for the state of the art of the subject.

A useful argument that can be used to establish the existence of solutions to (MFG) is the variational one, already presented in [START_REF] Lasry | Mean field games[END_REF]. The main idea behind is that, at least formally, system (MFG) can be seen as the first order optimality condition associated to minimizers of the following optimization problem where, in the definition of b, H * (x, •) denotes the Legendre-Fenchel conjugate of H(x, •). Under the assumption that f (x, •) and g(x, •) are non-decreasing, problem (P ) is shown to be a convex optimization problem and convex duality techniques can be successfully applied in order to provide existence and uniqueness results to (MFG). This argument has been made rigorous in several articles: let us mention [START_REF] Cardaliaguet | Weak solutions for first order mean field games with local coupling[END_REF][START_REF] Cardaliaguet | Mean field games systems of first order[END_REF] in the context of first order MFGs, the paper [START_REF] Cardaliaguet | Second order mean field games with degenerate diffusion and local coupling[END_REF] for degenerate second order MFGs, and finally [START_REF] Mészáros | A variational approach to second order mean field games with density constraints: the stationary case[END_REF][START_REF] Mészáros | On the variational formulation of some stationary second order mean field games systems[END_REF] for ergodic second order MFGs.

The variational approach described above has also been successful in the numerical resolution of system (MFG). In this direction, we mention the article [START_REF] Lachapelle | Computation of mean field equilibria in economics[END_REF] dealing with applications in economics, the paper [START_REF] Achdou | Mean field games: numerical methods for the planning problem[END_REF] concerned with the so-called planning problem in MFGs, the works [START_REF] Benamou | Augmented Lagrangian methods for transport optimization, mean field games and degenerate elliptic equations[END_REF][START_REF] Andreev | Preconditioning the augmented lagrangian method for instationary mean field games with diffusion[END_REF] focused on the resolution of a discretization of (P) by the Alternating Direction Method of Multipliers (ADMM) and [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF] where several first order methods are implemented and compared for the stationary version of (MFG). Let us mention that the variational approach is closely related to the so-called mean field optimal control problem, for which numerical methods have been studied in [START_REF] Burger | Mean field games with nonlinear mobilities in pedestrian dynamics[END_REF][START_REF] Albi | Mean-field control hierarchy[END_REF], among others.

In this paper we consider a finite difference discretization of problem (P). Assuming that f (x, •) and g(x, •) are non-decreasing, the discretization we consider is such that it preserves the convexity properties of problem (P) and the first order optimality conditions for its solutions, which are shown to exist, coincide with the finite difference scheme for MFGs introduced in [START_REF] Achdou | Mean field games: numerical methods[END_REF]. A very nice feature of this approach is that the solutions of the resulting discretized MFGs are shown to converge to the solutions of (MFG). We refer the reader to [START_REF] Achdou | Mean field games: convergence of a finite difference method[END_REF], where the convergence result is obtained under the assumption that (MFG) admits a unique classical solution, and to [START_REF] Achdou | Convergence of a finite difference scheme to weak solutions of the system of partial differential equations arising in mean field games[END_REF] in the framework of weak solutions (see [START_REF] Porretta | Weak solutions to Fokker-Planck equations and mean field games[END_REF] for the definition of this notion). We solve the discrete convex optimization problem by using the primal-dual algorithm introduced in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. As was pointed out in [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF] (see also [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF] in the context of transport problems), the primal-dual algorithm we consider seems to be faster than the ADMM when ν in (MFG) is small (or null). On the other hand, the efficiency of both methods is arguable when ν is large. This is due to the fact that, in both algorithms, at each iteration one has to invert a matrix whose condition number importantly increases as the viscosity parameter increases. Naturally, preconditioning strategies (see e.g. [START_REF] Benzi | Preconditioning techniques for large linear systems: A survey[END_REF]) can then be used in order to improve the efficiency of both algorithms. This strategy has been already successfully implemented in [START_REF] Andreev | Preconditioning the augmented lagrangian method for instationary mean field games with diffusion[END_REF] for the ADMM.

Our main objective in the present work is to take a closer look at the phenomenon described at the end of previous paragraph when considering the primal-dual algorithm. We have implemented standard indirect methods for solving the linear systems appearing in the computation of the iterates of the primaldual algorithm. As our numerical results suggest, it is very important to design suitable preconditioning strategies in order to be able to find solutions of the discretization of problem (P ) efficiently, and in a robust way with respect to the viscosity. For this, we explore different preconditioning strategies, and in particular, multigrid preconditioning (see also [START_REF] Achdou | Iterative strategies for solving linearized discrete mean field games systems[END_REF][START_REF] Andreev | Preconditioning the augmented lagrangian method for instationary mean field games with diffusion[END_REF], where multigrid strategies have been implemented for other solution methods).

The article is organized as follows. In section 2 we introduce some standard notation and we recall the finite difference scheme for (MFG) introduced in [START_REF] Achdou | Mean field games: numerical methods[END_REF]. The variational interpretation of this finite difference scheme is discussed in section 3. Next, in section 4, we recall the primal-dual algorithm introduced in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] and we consider its application to the discretization of (P). In section 5, we summarize the preconditioning strategies we consider and we discuss a numerical example, which is the time-dependent version of one of the examples treated in [START_REF] Achdou | Mean field games: numerical methods[END_REF][START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF].

2. Preliminaries and the finite difference scheme in [START_REF] Achdou | Mean field games: numerical methods[END_REF] In this section we introduce some basic notation and present the finite difference scheme introduced in [START_REF] Achdou | Mean field games: numerical methods[END_REF], whose efficient resolution will be the subject of this article. For the sake of simplicity, we will assume that d = 2 and that given q > 1, with conjugate exponent denoted by q = q/(q -1), the Hamiltonian H :

T 2 × R 2 → R has the form H(x, p) = 1 q |p| q ∀ x ∈ T 2 , p ∈ R 2 .
In this case the function b defined in (1.1) takes the form

b(x, m, w) =    |w| q qm q-1 if m > 0, 0 if (m, w) = (0, 0), +∞ otherwise.
Let N T , N h be positive integers and set ∆t = T /N T , the time step, and h = 1/N h , the space step. We associate to these steps a time grid T ∆t := {t n = n∆t ; n = 0, . . . , N T } and a space grid T 2 h := {x i,j = (ih, jh) ; i, j ∈ Z}. Since T 2 h intends to discretize T 2 , we impose the identification z i,j = z (i mod N h ),(j mod N h ) , which allows to assume that i, j ∈ {0, . . . , N h -1}. A function y := T 2 ×[0, T ] → R is approximated by its values at (x i,j , t n ) ∈ T 2 h × T ∆t , which we denote by y n i,j := y(x i,j , t n ). Given y : T 2 h → R we define the first order finite difference operators (2.1) (D 1 y) i,j = y i+1,j -y i,j h , and (D 2 y) i,j = y i,j+1 -y i,j h , [D h y] i,j = ((D 1 y) i,j , (D 1 y) i-1,j , (D 2 y) i,j , (D 2 y) i,j-1 ).

For later use, given a ∈ R, we set a + := max(a, 0) and a -:= a + -a, and we define

[D h y] i,j = ((D 1 y) - i,j , -(D 1 y) + i-1,j , (D 2 y) - i,j , -(D 2 y) + i,j-1
). The discrete Laplace operator ∆ h y : T 2 h → R is defined by (∆ h y) i,j = -1 h 2 (4y i,j -y i+1,j -y i-1,j -y i,j+1 -y i,j-1 ). For y : T 2 h × T ∆t → R we define the discrete time derivative

D t y k i,j = y k+1 i,j -y k i,j
∆t .

The Godunov-type finite difference discretization of (MFG) introduced in [START_REF] Achdou | Mean field games: numerical methods[END_REF] is as follows: find u, m :

T 2 h × T ∆ → R such that for all 0 ≤ i, j ≤ N h -1 and 0 ≤ k ≤ N T -1 we have (MFG h,∆t )            -D t u k i,j -ν(∆ h u k ) i,j + 1 q | [D h u k ] i,j | q = f (x i,j , m k+1 i,j ), D t m k i,j -ν(∆ h m k+1 ) i,j -T i,j (u k , m k+1 ) = 0, m 0 i,j = mi,j , u N T i,j = g(x i,j , m N T i,j ), where (2.2) mi,j := |x-xi,j |∞≤ h 2 m 0 (x)dx ≥ 0,
and the operator T (u , m ) :

T 2 h → R, with u , m : T 2 h → R, is defined by T i,j (u , m ) := 1 h -m i,j 1 q | [D h u ] i,j | 2-q q-1 (D 1 u ) - i,j + m i-1,j 1 q | [D h u ] i-1,j | 2-q q-1 (D 1 u ) - i-1,j +m i+1,j 1 q | [D h u ] i+1,j | 2-q q-1 (D 1 u ) + i,j -m i,j 1 q | [D h u ] i,j | 2-q q-1 (D 1 u ) + i-1,j -m i,j 1 q | [D h u ] i,j | 2-q q-1 (D 2 u ) - i,j + m i,j-1 1 q | [D h u ] i,j-1 | 2-q q-1 (D 2 u ) - i,j-1 +m i,j+1 1 q | [D h u ] i,j+1 | 2-q q-1 (D 2 u ) + i,j -m i,j 1 q | [D h u ] i,j | 2-q
q-1 (D 2 u ) + i,j-1 , with the convention:

(2.3) | [D h u ] i,j | 2-q q-1 [D h u ] i,j = 0 if q > 0 and [D h u ] i,j = 0.
The existence of a solution (u h,∆t , m h,∆t ) of system (MFG h,∆t ) is proved in [START_REF] Achdou | Mean field games: numerical methods[END_REF]Theorem 6] as a consequence of Brouwer fixed point theorem. Furthermore, if we assume that f and g are increasing with respect to their second argument, and one of them is strictly increasing, this solution is unique when h is small enough (see [START_REF] Achdou | Mean field games: numerical methods[END_REF]Theorem 7]). As we will see in the next section, these results can also be obtained by variational arguments. The convergence, as h and ∆t tend to 0, of suitable extensions of u h,∆t and m h,∆t to T 2 × [0, T ] to a solution (u, m) of (MFG) is proved in [START_REF] Achdou | Mean field games: convergence of a finite difference method[END_REF] under the assumption that (u, m) is unique and sufficiently regular. The later smoothness assumption has been relaxed in [START_REF] Achdou | Convergence of a finite difference scheme to weak solutions of the system of partial differential equations arising in mean field games[END_REF].

The finite dimensional variational problem and the discrete MFG system

Following [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF] in the stationary case and [START_REF] Achdou | Mean field games: numerical methods for the planning problem[END_REF] for the planning problem, we introduce some finitedimensional operators that will allow us to write easily a finite dimensional version of problem (P). Denoting by R + the set of non-negative real numbers and by R -the set of non-positive real numbers, we define 4 , we denote by y l the l-th coordinate of y, 1 ≤ l ≤ 4. Let A : M → U and B : W → U be defined by

K := R + × R -× R + × R -and we denote by P K (v) the orthogonal projection of v ∈ R 4 onto K. Let M := R (N T +1)×N h ×N h , W := (R 4 ) N T ×N h ×N h and U := R N T ×N h ×N h . If y ∈ R
(3.1) (Am) k i,j = D t m k i,j -ν(∆ h m k+1 ) i,j , (Bw) k i,j = (D 1 w k,( 1
) ) i-1,j + (D 1 w k,( 2
) ) i,j + (D 2 w k,( 3
) ) i,j-1 + (D 2 w k,( 4
) ) i,j , for all 0 ≤ i, j ≤ N h -1 and 0 ≤ k ≤ N T -1.
One can easily check (see e.g. [START_REF] Achdou | Mean field games: numerical methods[END_REF]) that the corresponding dual operators are given by

(3.2) (B * u) k i,j = -[D h u k ] i,j for all 0 ≤ k ≤ N T -1, (A * u) k i,j = -D t u k-1 i,j -ν(∆ h u k-1 ) i,j , if 1 ≤ k ≤ N T -1, (A * u) 0 i,j = - 1 ∆t u 0 i,j , (A * u) N T i,j = 1 ∆t u N T -1 i,j -ν(∆ h u N T -1 ) i,j ,
for all u ∈ U. For latter use, notice that

Ker(B * ) = {u ∈ U | ∀ k = 0, . . . , N T -1 there exists c k ∈ R such that u k i,j = c k ∀ i, j}, and so (3.3) Im(B) = Ker(B * ) ⊥ = u ∈ U i,j u k i,j = 0 ∀ k = 0, . . . , N T -1 . Let us define b : R × R 4 → R ∪ {+∞} (3.4) b(m, w) :=        |w| q qm q-1 , if m > 0, w ∈ K, 0, if (m, w) = (0, 0), +∞, otherwise,
and the functions B, F :

M × W → R, G : M × W → M × R N h ×N h as (3.5) B(m, w) = 1≤k≤N T , 0≤i,j≤N h -1 b(m k i,j , w k-1 i,j ), F(m, w) = 1≤k≤N T , 0≤i,j≤N h -1 F (x i,j , m k i,j ) + 1 ∆t 0≤i,j≤N h -1 G(x i,j , m N T i,j ), G(m, w) := (Am + Bw, m 0 ). Note that if (m, w) ∈ M × W is such that G(m, w) = (0, m) then (3.6) h 2 i,j m k i,j = 1 ∀ k = 0, . . . , N T .
Indeed, by periodicity,i,j (∆ h m k+1 ) i,j = 0 and i,j (Bw) k i,j = 0 for all k = 0, . . . , N T -1. This implies that

0 = i,j (Am + Bw) k i,j = i,j m k+1 i,j ∆t - i,j m k i,j ∆t ,
and so h 2 i,j m k i,j = h 2 i,j mi,j = 1 for all k = 0, . . . , N T . The discretization of the variational problem (P) we consider is

(P h,∆t ) inf (m,w)∈M×W B(m, w) + F(m), subject to G(m, w) = (0, m),
where we recall that F and G in (3.5) are defined in (1.1).

We have the following result Theorem 3.1. For any ν > 0 problem (P h,∆t ) admits at least one solution (m h,∆t , w h,∆t ) and associated to it there exists u h,∆t : M × W → R such that (MFG h,∆t ) holds true. Moreover, (m h,∆t ) k i,j > 0 for all k = 1, . . . , N T , i, j = 0, . . . , N h -1.

In order to prove the result above, let us first show a lemma that implies the feasibility of the constraints in (P h,∆t ). Lemma 3.1. There exists ( m, w) ∈ M × W such that

(3.7) G( m, w) = (0, m), wk i,j ∈ int(K) ∀ i, j = 1, . . . , N h -1, k = 1, . . . , N T -1, mk i,j > 0, ∀ i, j = 1, . . . , N h -1, k = 1, . . . , N T .
Proof. Let us define m0 i,j = mi,j and mk i,j := 1 for all k = 1, . . . , N T and i, j. Since h 2 i,j mk i,j = 1 for all k = 0, . . . , N T , by (3.3) and the definition of A we easily get that A m ∈ Im(B). Therefore, there exists ŵ ∈ W satisfying G( m, ŵ) = (0, m). Then, given δ > 0, we set for all k = 0, . . . , N T -1 and i, j

wk i,j := ŵk,(1) + max i,j ŵk, (1)
i,j

+ δ, ŵk,(2) -max i,j ŵk,( 2
) i,j -δ, ŵk,(3) + max i,j ŵk,( 3
) i,j + δ, ŵk,(4) -max i,j ŵk, (4) i,j
-δ , which satisfies wk i,j ∈ int(K) and (B w) k = (B ŵ) k . The result follows. Now, we prove the existence of solutions to (P h,∆t ).

Lemma 3.2. Problem (P h,∆t ) admits at least one solution (m h,∆t , w h,∆t ) and every such solution satisfies (m h,∆t ) k i,j > 0 for all k = 1, . . . , N T , i, j = 0, . . . , N h -1.

Proof. Let (m n , w n ) be a minimizing sequence for (P h,∆t ). Lemma 3.1 implies that B( m, w) + F( m) < +∞. Therefore, there exists a constant C 1 > 0 such that

(3.8) B(m n , w n ) + F(m n ) ≤ C 1 for all n ∈ N.
As a consequence, by definition of b, (m n ) k i,j ≥ 0 for all i, j and k and (w n ) k ∈ K for all k. Since Am n + Bw n = 0, relation (3.6) implies that h 2 i,j (m n ) k i,j = 1. In particular, there exists

C 2 > 0 (independent of n) such that sup i,j,k (m n ) k i,j ≤ C 2 . Using that, if (m n ) k i,j > 0, b((m n ) k i,j , (w n ) k i,j ) ≥ |(w n ) k i,j | q qC q-1 2 ,
and that F(m n ) is uniformly bounded (because F and G are continuous and m n is bounded), relation (3.8) yields the existence of

C 3 > 0 (independent of n) such that sup i,j,k |(w n ) k i,j | ≤ C 3 .
Thus, there exists (m h,∆t , w h,∆t ) ∈ M × W such that, up to some subsequence, m n → m h,∆t and w n → w h,∆t as n → ∞. Since G(m n , w n ) = (0, m) we obtain that G(m h,∆t , w h,∆t ) = (0, m), The lower semicontinuity of B + F implies that

B(m h,∆t , w h,∆t ) + F(m h,∆t ) ≤ lim n→∞ B(m n , w n ) + F(m n ),
which implies that (m h,∆t , w h,∆t ) solves (P h,∆t ). Finally, if (m, w) ∈ M × W solves (P h,∆t ) and m k i,j = 0 for some i, j and k = 1, . . . , N T , then, by the definition of B, we must have that w k-1 i,j = 0. Thus, the constraint (Am + Bw) k-1 i,j = 0 can be written as

- m k-1 i,j ∆t -ν h 2 (m k i+1,j + m k i-1,j + m k i,j+1 + m k i,j-1 ) = w k-1,(1) i-1,j h - w k-1,(2) i+1,j h + w k-1,(3) i,j-1 h - w k-1,(4) i,j+1 h .
Since the left hand side above is non-positive and the right hand side is non-negative (by definition of K), we deduce that all the terms above are zero. By repeating the argument at the indexes neighboring (i, j), we deduce that m k ≡ 0 and so h 2 i,j m k i,j = 0 which, by (3.6), contradicts G(m, w) = (0, m). The result follows.

Proof of Theorem 3.1. By Lemma 3.2 we know that there exists a solution (m h,∆t , w h,∆t ) to (P h,∆t ) and m h,∆t i,j > 0 for all i, j. Thus, in order to conclude it suffices to show the existence of u h,∆t such that (MFG h,∆t ) holds true. For notational convenience we will omit the superindexes h and ∆t. Define the

Lagrangian L := M × W × U × R N h ×N h → R ∪ {+∞}, associated to (P h,∆t ), as (3.9) L(m, w, u, λ) := B(m, w) + F(m) -u, Am + Bw -λ, m 0 -m = B(m, w) + F(m) -A * u, m -B * u, w -λ, m 0 -m .
Note that the linear mapping M m → (Am, m) ∈ U × R N h ×N h is invertible as it is shown by its matrix representation (see (4.7) in the next section). As a consequence G is surjective and, hence, by standard arguments, there exists (u, λ)

∈ U × R N h ×N h such that (3.10) 0 = ∂ m k i,j L(m, w, u, λ) = -1 q |w k-1 i,j | q (m k i,j ) q + f (x i,j , m k i,j ) -[A * u] k i,j ∀ k = 1, . . . , N T -1, ∀ i, j, 0 = ∂ m 0 i,j L(m, w, u, λ) = -λ i,j -[A * u] 0 i,j ∀ i, j, 0 = ∂ m N T i,j L(m, w, u, λ) = -1 q |w N T -1 i,j | q (m N T i,j ) q + f (x i,j , m N T i,j ) + 1 ∆t g(x i,j , m N T i,j ) -[A * u] N T i,j ∀ i, j, 0 ∈ ∂ w k-1 i,j L(m, w, u, λ) = |w k-1 i,j | q-2 w k-1 i,j (m k i,j ) q-1 -[B * u] k-1 i,j + N K (w k-1 i,j ) ∀ k = 1, . . . , N T , ∀ i, j,
where we have used definition (3.4) and that m k i,j > 0 for all k = 1, . . . , N T and all i, j. Defining u N T i,j := g(x i,j , m N T i,j ), by the last relation in (3.2), the third relation in (3.10) can be rewritten as

-D t u N T -1 i,j -ν(∆ h u N T -1 ) i,j + 1 q |w N T -1 i,j | q (m N T i,j ) q = f (x i,j , m N T i,j ),
and hence, by the second relation in (3.2) and the first relation in (3.10), we have that

(3.11) -D t u k i,j -ν(∆ h u k ) i,j + 1 q |w k i,j | q (m k+1 i,j ) q = f (x i,j , m k+1 i,j ) ∀ k = 0, . . . , N T -1, ∀ i, j.
The last relation in (3.10) yields that for all k = 1, . . . , N T and all i, j

   (m k i,j ) q-1 |w k-1 i,j | q-2 [B * u] k-1 i,j ∈ w k-1 i,j + N K (w k-1 i,j ) if w k-1 i,j = 0, [B * u] k-1 i,j ∈ N K (0) if w k-1
i,j = 0, which, by (3.2) and under the convention (2.3), is equivalent to

(3.12) w k-1 i,j = m k i,j |P K (-[D h u] k-1 i,j )| 2-q q-1 P K (-[D h u] k-1 i,j ) = m k i,j | [D h u k-1 ] i,j | 2-q q-1 [D h u k-1 ] i,j .
Shifting the index k, the expression above yields 1 q

|w k i,j | q (m k+1 i,j ) q = 1 q | [D h u k ] i,j | q ∀ k = 0, . . . , N T -1, ∀ i, j,
which, combined with (3.11), implies the first equation in (MFG h,∆t ). The second equation in (MFG h,∆t ) is a consequence of Am + Bw = 0 and the fact that (3.12) provides the identity

(Bw) k i,j = -T i,j (u k , m k+1 ) ∀ k = 0, . . . , N T -1, ∀ i, j.
The result follows.

Remark 3.1. The proof of the existence of solutions to (MFG h,∆t ) in Theorem 3.1 provides an alternative argument to the one in [START_REF] Achdou | Mean field games: numerical methods[END_REF], based on Brouwer fixed-point theorem. The uniqueness of solutions under the assumption that f (x, •) and g(x, •) are increasing, with one of them being strictly increasing, is straightforward because, in this case, the cost functional in (P h,∆t ) is strictly convex.

A primal-dual algorithm to solve (P h,∆t )

As discussed in [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF], for solving the optimization problem where ϕ : R N → R ∪ {+∞} and ψ : R N → R ∪ {+∞} are convex l.s.c. proper functions, methods in [START_REF] Briceño-Arias | A monotone+ skew splitting model for composite monotone inclusions in duality[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF][START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF][START_REF] Glowinski | Sur l'approximation, par éléments finis d'ordre un, et la résolution, par pénalisationdualité, d'une classe de problèmes de Dirichlet non linéaires[END_REF] can be applied with guaranteed convergence under mild assumptions. In [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF], devoted to the stationary case, the method proposed in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] has the best performance when the viscosity parameter is small. This method is inspired by the first-order optimality conditions satisfied by a solution (ŷ, σ) to (4.1)-(4.2) under standard qualification conditions, which reads [34, Theorem 8]

(4.3) -σ ∈ ∂ϕ(ŷ) ŷ ∈ ∂ψ * (σ) ⇔ ŷ -τ σ ∈ τ ∂ϕ(ŷ) + ŷ σ + γ ŷ ∈ γ∂ψ * (σ) + σ ⇔ prox τ ϕ (ŷ -τ σ) = ŷ prox γψ * (σ + γ ŷ) = σ,
where γ > 0 and τ > 0 are arbitrary and, given a l.s.c. convex proper function φ :

R N → ]-∞, +∞], prox γφ x := argmin y∈R N φ(y) + |y -x| 2 2γ = (I + ∂(γφ)) -1 (x) ∀ x ∈ R N .
Given θ ∈ [0, 1], τ and γ satisfying τ γ < 1, and starting points (y 0 , ỹ0 , σ 0 ) ∈ R N × R N × R M , the iterates {(y k , σ k )} k∈N generated by (4.4)

σ k+1 = prox γψ * (σ k + γ ỹk ),
y k+1 = prox τ ϕ (y k -τ σ k+1 ), ỹk+1 = y k+1 + θ(y k+1 -y k )
converge to a primal-dual solution (ŷ, σ) to (4.1)-(4.2) (see, e.g., [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]).

In the case under study, the equations of the time-dependent discretization are very similar to their stationary counterparts (see [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF]). Specifically, the discrete linear operators A and B defined in (3.1), by an abuse of notation, are represented by real matrices A and B, of dimensions (

N T × N 2 h )×((N T + 1) × N 2 h ) and (N T × N 2 h ) × (N T × 4N 2 h )
, respectively, given by (4.5)

A =       -1 ∆t Id N 2 h νL + 1 ∆t Id N 2 h 0 • • • 0 0 -1 ∆t Id N 2 h νL + 1 ∆t Id N 2 h . . . . . . . . . . . . . . . . . . 0 0 • • • 0 -1 ∆t Id N 2 h νL + 1 ∆t Id N 2 h       and (4.6) B =      M 0 • • • 0 0 M • • • 0 . . . . . . . . . 0 • • • 0 M      , where L ∈ M N 2 h , N 2 h (R) is the matrix that represents -∆ h on the torus T 2 h and M ∈ M N 2 h , 4N 2
h (R) is the matrix representing the discrete divergence. Denoting by à and B the ( (i) The matrix à is block lower triangular with invertible diagonal blocks and, hence, it is invertible. Indeed, the first diagonal block Id N 2 h is obviously invertible and the other blocks, given by νL + 1 ∆t Id N 2 h , are also invertible because they are strictly diagonally dominant. (ii) Since à is invertible, the matrix

(N T + 1) × N 2 h )×((N T + 1) × N 2 h ) and ((N T + 1) × N 2 h ) × (N T × 4N 2 h ) real matrices (4.7) Ã = Id N 2 h 0 • • • 0 A and B = 0 • • • 0 B , the constraint G(m, w) = (0, m) in (P h,
prox ψ : (m, w) → (m, w) -C * Q -1 (C(m, w) -( m, 0)),
we have

prox γψ * : (m, w) → C * Q -1 (C(m, w) -γ( m, 0)),
where Q is defined in (4.8). By setting y 0 = (m 0 , w 0 ), ỹ0 = ( m 0 , w 0 ),

σ 0 = (n 0 , v 0 ) ∈ R N T ×N 2 h ×R N T (4N 2 h ) , (4.4) becomes (4.9)                          z [l+1] = -Q -1 Ã(n [l] + γ m[l] ) + B(v [l] + γ w[l] ) -γ( m, 0) , n [l+1] v [l+1] = Ã * z [l+1] B * z [l+1] , m [l+1] w [l+1] = prox τ ϕ m [l] + τ n [l+1] w [l] + τ v [l+1] , m[l+1] w[l+1] = m [l+1] + θ(m [l+1] -m [l] ) w [l+1] + θ(w [l+1] -w [l] ) ,
and, if γτ < 1, the convergence of (m [l] , w [l] ) to a solution ( m, ŵ) to (P h,∆t ) is guaranteed together with the convergence of (n [l] , v [l] ) to some (n, v) as l → ∞. In order to compute the Lagrange multiplier û ∈ U, which solves the first equation in (MFG h,∆t ), note that (3.9) can be written equivalently as where ( m, ŵ) is the primal solution and ẑ = ( λ, û). Therefore, in order to approximate ẑ, note that from (4.9) we have (4.12)

m [l] -m [l+1] τ + Ã * z [l+1] w [l] -w [l+1] τ + B * z [l+1] ∈ ∂ϕ(m [l+1] , w [l+1] )
and, hence, since the algorithm generates converging sequences m [l] → m and w [l] → ŵ and z

[l] → ẑ := -Q -1 ( Ã(n + γ m) + B(v + γ ŵ) -γ( m, 0)
), the closedness of the graph of ∂ϕ [8, Proposition 20.38] yields (4.11) and, hence, a good approximation of ẑ is z [l] for l large enough. For obtaining

[u * ] N T i,j , a good approximation is [u [l] ] N T i,j = g(x i,j , [m [l] ] N T i,j
). Remark 4.2. (i) In order to obtain an efficient algorithm, the computation of prox τ ϕ in (4.9) should be fast. A complete study of prox τ ϕ is presented in [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF]Section 3.2] showing that its computation depends on the resolution of a real equation, which can be efficiently solved. (ii) An important step in (4.9) is the efficient computation of the inverse of Q. Different preconditioning strategies to tackle this issue will be presented in the following section.

Preconditioning strategies

At the beginning of each iteration of the primal-dual algorithm (4.9), we require the solution of a linear system (5.1)

Qz = b.
The purpose of this section is to discuss preconditioning strategies for the solution of this linear system. For the stationary setting discussed in [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF], the solution of such a system via direct methods as the backlash (mldivide) command in MATLAB 1 was feasible for relatively fine meshes (up to the order of 100 nodes per space dimension). However, as shown in Table 1, introducing a temporal dimension and thus increasing the degrees of freedom to N 2 h × N T significantly increases the computation time. Indeed, the use of backlash on fine space and time grids -e.g. 128 2 space grid points and 40 time steps -requires an amount of RAM that is prohibitive on the machine used for our performance tests 2 , leading to "out of memory" errors. We mitigate this problem by exploring the solution of (5.1) via preconditioned iterative methods, which perform efficiently for finer space and time subdivisions and different viscosities. We begin by illustrating the difficulties associated to the conditioning of the system in (4.9). Table 2 shows the condition number of the system for different space-time discretizations and viscositiy values. Without any precoditioner, the condition numbers of different discretizations scale up to 10 8 . The same Table shows that by selecting a suitable preconditioner, such as the modified incomplete Cholesky factorization [START_REF] Benzi | Preconditioning techniques for large linear systems: A survey[END_REF] (michol in MATLAB), the conditioning of the system is improved by 4 orders of magnitude.

We have tested different choices of preconditioners and iterative methods for our problem. Since the matrix Q in our setting is sparse, symmetric, and positive-definite, we have implemented an incomplete Cholesky factorization with diagonal scaling, a modified incomplete Cholesky factorization, and multigrid preconditioning. As for the choice of the iterative method, our tests included both conjugate gradient (pcg), and the biconjugate gradient stabilized method (BiCGStab). The interested reader will find in [START_REF] Wendland | Numerical Linear Algebra[END_REF]Chapters 6 and 8] a thorough description of the aforementioned methods, and in the Appendix of this article tables with the performance of the different methods. Our findings suggest that the use of a modified incomplete Cholesky factorization in the preconditioned conjugate gradient method is satisfactory for small viscosities ( ν ≤ 0.05). However, among the tested methods, the multigrid preconditioner with BiCGStab iteration (as in Algorithm 1) is the only method which performs consistently for different viscosities and space-time discretizations. In the following, we discuss its implementation and assess its performance.

Table 2. Condition numbers for Q without preconditioning (a), and with modified incomplete Cholesky factorization preconditioning (b). 0.0004 0.0022 0.0048 5 × 10 -1 0.0004 0.0227 0.0465 0.5 0.4398 1.9248 3.6541 5.1. Multigrid preconditioner. We implement multigrid preconditioned algorithm for solving (5.1). We refer the reader to [START_REF] Trottenberg | Multigrid[END_REF] for an introduction and an overview of multigrid methods. We briefly review the main concepts behind the method. Consider two linear systems A 1 x1 = b 1 and A 0 x0 = b 0 , stemming from two discretizations of a linear PDE over the grids G 1 and G 0 , respectively. Assume also that G 1 is a refinement of G 0 . Loosely speaking, the main idea of the method is that in order to find a good approximation of the solution x1 on the finer grid, we first consider what is known as a smoothing step. This step consists on computing a few iterates x 1 1 , . . . , x η1 1 with a standard indirect method, such as Jacobi or Gauss-Seidel, and to define the residual r

1 := b 1 -A 1 x η1
1 , which is shown to be smoother (less oscillatory) than the first residual b 1 -A 1 x 1 . Then, we consider in the coarser grid G 0 the second system A 0 x0 = b 0 with b 0 = r1 , where r1 is the restriction of r 1 to G 0 . Assuming that we can compute a good approximation of x0 , which we still denote by x0 , we then extend this solution to G 1 by using a linear interpolation. Calling e 1 the resulting vector, we update x η1

1 by redefining it as x η1 1 + e 1 and we end the procedure by applying again a few iterations, say η 2 , of a smoothing method initialized at x η1 1 . This last step is called post smoothing.

The previous paragraph introduced what is known as a two grid iteration. If we consider more grids G 0 , G 1 ,. . ., G , where for each k = 0, . . . , -1, G k ⊆ G k+1 , we can proceed similarly and obtain a better approximation of the solution to A x = b . As in the previous case, we begin with the finest grid G and we perform η 1 smoothing steps to obtain the residual r := b -A x η1 whose restriction to G -1 is denoted by r . In this grid we consider the system A -1 x -1 = r and we perform again a smoothing step and a restriction of the residual to G -2 . The procedure continues until we get to the coarsest grid G 0 , where the solution e 0 to the corresponding linear system can be found easily (typically using a direct method). Next, the solution e 1 on the grid G 1 is corrected with the interpolation of e 0 . Another post smoothing is performed to the corrected solution on G 1 and using its interpolation in the grid G 2 we correct the previous solution on this grid. The smoothing, interpolation and correction iterations end once we arrive to the finest grid G to obtain the final approximation of x . The previous procedure is called a multigrid method with a V -cycle. An alternative, to obtain a more accurate solution, is to proceed as before going from G to G -1 and then for k = -1, . . . , 1 to perform two consecutive coarsegrid corrections, instead of one as in the V -cycle. The resulting procedure is known as multigrid with a W -cycle. Finally, in between the V -cycle and the W -cycle, we have the F-cycle, where in the process of going from the coarsest grid to the finest one, if a grid has been reached for the first time, another correction with the coarser grids using a V -cycle is performed.

In our context, we use one cycle of the multigrid algorithm, which is a linear operator as a function of the residual on the finest grid, as a preconditioner for solving (5.1) with the BiCGStab method. Since Q is related to the finite difference discretization of the operator -∂ 2 tt + ν 2 ∆ 2 -∆ and ν is not necessarily small, as in [START_REF] Achdou | Iterative strategies for solving linearized discrete mean field games systems[END_REF], it is natural to consider the refinements of the grid only in the space variable (we refer the reader to [START_REF] Trottenberg | Multigrid[END_REF] for semi-coarsening multigrid methods in the context of anisotropic operators). We suppose that the spatial mesh is such that N h = H2 , with H > 1 and is a positive integer (in the numerical example in the next section H will be equal to 2 or 3, H 2 being the number of spatial points in the coarsest grid).

Let us specify the main steps of the multigrid method we use as a preconditioner.

Hierachy of Grids: Semi-coarsened grids G k with size (N T + 1)H 2 2 2k for all k = 0 . . . .

Cycle:

We use the F-cycle.

Restriction operator: As in [START_REF] Achdou | Iterative strategies for solving linearized discrete mean field games systems[END_REF], in order to restrict the residual on the grid G k to the grid G k-1 , we use the second-order operator R k : R

(2 k H) 2 (N T +1) → R (2 k-1 H) 2 (N T +1) defined by (R k X) n i,j = 1 16 4X n 2i,2j + 2(X n 2i+1,2j + X n 2i-1,2j + X n 2i,2j+1 + X n 2i,2j-1 ) X n 2i-1,2j-1 + X n 2i-1,2j+1 + X n 2i+1,2j-1 + X n 2i+1,2j+1
, for n = 0, . . . , N T , i, j = 1, . . . , 2 k-1 H.

Interpolation operator: We denote by

I k : R (2 k-1 H) 2 (N T +1) → R (2 k H) 2 (N T +1
) the interpolation operator from the grid G k-1 to the grid G k . We have chosen a standard bilinear interpolation operator in the space variable, which is also a second-order operator and dual to the restriction operator (I k = 4R * k ). According to [START_REF] Brandt | Rigorous quantitative analysis of multigrid. I. Constant coefficients two-level cycle with L 2 -norm[END_REF], the sum of the orders of R k and I k has to be at least equal to the degree of the differential operator. In our context, both are equal to 4.

Linear systems on the different grids: The linear systems are defined by the matrices

Q k := A k A * k + B k B * k , k = 0, . . . , ,
where we recall that A k and B k are the finite difference discretizations of ∂ t -ν∆ and div(•), respectively, on the grid G k (see (3.1)).

Smoother: Here we have used Gauss-Seidel iterations in the lexicographic order. There is no reason for choosing the lexicographic order, other than its simplicity.

Solving the system on the coarsest grid G 0 : We can use an exact solver such as backlash in MATLAB. Indeed, in G 0 the size of the system is really small with respect to the size of the system on the grid G (in G 0 , we can even store the inverse of Q 0 and inversion at this level just becomes a matrix multiplication).

The multigrid preconditoning procedure is summarized in Algorithm 2.

Algorithm 2 Multigrid Preconditioner for Q x = b P L : y → MultigridSolver( , 0, y,cycle) x l ←BiCGStab(Q , b , P L , I d , x 0 , tol ) procedure MultigridSolver(k, x k , b k ,cycle) if k = 0 then x k ← Q -1 0 b k else x k ←Perform η 1 steps of Gauss-Seidel from x k with b k as second member. x k-1 ← 0 x k-1 ←MultigridSolver(k -1, x k-1 , R k (b k -Q k x k ),cycle) if cycle is W then x k-1 ←MultigridSolver(k -1, x k-1 , R k (b k -Q k x k ),cycle) if cycle is F then x k-1 ←MultigridSolver(k -1, x k-1 , R k (b k -Q k x k ),V) x k ← x k + I k x k-1
x k ←Perform η 2 steps of Gauss-Seidel from x k with b k as second member. return x k 5.2. Numerical Tests. In this section we present a test case considered in [START_REF] Achdou | Mean field games: numerical methods[END_REF], for which the stationary solution has been computed numerically in [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF] using the primal-dual algorithm presented above. The setting is as follows: for (x, y) ∈ R 2 and m ∈ R + , f (x, y, m) = m 2 -H(x, y), H(x, y) = sin(2πy) + sin(2πx) + cos(2πx).

We first validate the dynamic behavior of our solution. Figure 1 shows the evolution of the mass at four different time steps. Starting from a constant initial density, the mass converges to a steady state, and then, when t gets close to the final time T , the mass is influenced by the final cost and converges to a final state. This behavior is referred to as turnpike phenomenon in the literature [START_REF] Porretta | Remarks on long time versus steady state optimal control[END_REF]. It is illustrated by Figure 2, which displays as a function of time t the distance of the mass at time t to the stationary state computed as in [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF].

For the multi-grid preconditioner, Table 3 shows the computation times for different discretizations. It can be observed that finer meshes with 128 3 degrees of freedom are solvable within CPU times which outperfom others methods shown in the Appendix and in [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF]. Furthermore, the method is robust with respect to different viscosities.

From Table 3 we observe that most of the computational time is used for solving the second proximal operator (the third equality of (4.9)), which does not use a multigrid strategy but which is a pointwise operator (see Proposition 3.1 of [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF]) and thus could be fully paralellizable. Unlike the stationary case, low viscosities seem to make the algorithm be slightly slower. However, Table 4 shows that the average number of iterations of BiCGStab stays low regardless of the viscosity. Indeed Table 3 shows that more Chambolle-Pock iterations are needed to converge. The same behavior happens when we use a direct exact solver instead of the multi-grid preconditioned BiCGStab algorithm. Concluding Remarks. In this work we have developed a first-order primal-dual algorithm for the solution of second-order, time-dependent mean field games. The procedure consists of: a variational formulation for the MFG, its discretization via finite differences, the application Chambolle-Pock algorithm to the resulting minimization. While this method has been studied for stationary MFG in [START_REF] Briceño-Arias | Proximal methods for stationary mean field games with local couplings[END_REF], its numerical realization for time-dependent MFGs was prohibitive in terms of computing time, as the Chambolle-Pock iteration requires the solution of a large-scale linear system at each iteration. We have overcome this difficulty by studying different preconditioning strategies for the associated linear system. Overall, the multigrid preconditioner with a BiCGStab iteration performs satisfactorily for different discretizations and viscosity values. partially supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 639638). His participation to this project has also been partially supported by a CEMRACS 2017 scholarship.

The sixth author thanks the support from the PGMO project VarPDEMFG and from the ANR project MFG ANR-16-CE40-0015-01.

T 1 F

 1 d [b(x, m(x, t), w(x, t)) + F (x, m(x, t))] dx + T d G(x, m(x, T ))dx subject to ∂ t m -ν∆m + div(w) = 0 in T d × (0, T ), m(•, 0) = m 0 (•) in T d ,(provided that they exist). In (P), the functions b :T d × R × R d → R ∪ {+∞} and F , G : T d × R → R ∪ {+∞} are defined as follows ((x, m) := m 0 f (x, m )dm if m ≥ 0, +∞ otherwise, G(x, m) := m 0 g(x, m )dm if m ≥ 0, +∞ otherwise,

ϕ

  * (-σ) + ψ * (σ),

  ∆t ) can be rewritten as C(m, w) = ( m, 0), where C := [ Ã | B]. Remark 4.1.

( 4 . 8 )

 48 Q := CC * = Ã Ã * + B B * is positive definite and, hence, invertible. Therefore, (P h,∆t ) is a particular instance of (4.1) with ϕ(m, w) := B(m, w) + F(m), ψ(m, w) := ι ker C+{(m f ,w f )} (m, w), where (m f , w f ) is a feasible vector (provided for instance by Lemma 3.1), and ι ker C+{(m f ,w f )} is the function defined as 0 for all (m, w) ∈ ker C + {(m f , w f )} and +∞, otherwise. Since prox γψ * = Id -γ prox ψ/γ •(Id/γ) = Id -γ prox ψ •(Id/γ) (see e.g. [8, Section 24.2]) and

(4. 10 )

 10 L(m, w, u, λ) := ϕ(m, w) -(λ, u), Ãm + Bw + λ, m = ϕ(m, w) -∈ ∂ϕ( m, ŵ),

Algorithm 1

 1 Preconditioned BiCGStab x l ←BiCGStab(Q l , b l , P L , P R , x 0 , tol) procedure BiCGStab(Q l , b l , P L , P R , x 0 , tol) r 0 := p 0 := Q l x 0 -b l ; r0 := p0 := P l r 0 ; ρ0 := r 0 , r 0 ; k := 0 while r k ≥ tol do v k := QP R pk vk := P L v k αk := ρk / vk , r0 s k+1 := r k -αk v k ŝk+1 := P L s k+1 t k+1 := AP R ŝk+1 tk+1 := P L t k+1 ωk+1 := ŝk+1 , tk+1 / tk+1 , tk+1 xk+1 := xk + αk pk + ωk+1 ŝk+1 r k+1 := s k+1 -ωk+1 t k+1 rk+1 := ŝk+1 -ωk+1 tk+1 ρk+1 := rk+1 , r0 βk+1 := (α k /ω k+1 )(ρ k+1 /ρ k ) pk+1 := rk+1 + βk+1 (p k -ωk+1 vk ) k := k + 1 return x k := P R xk

Figure 1 .Figure 2 .

 12 Figure 1. Evolution of the density m obtained with the multi-grid preconditioner for ν = 0.5, T = 1, N T = 200 and N h = 128. At t = 0.12 the solution is close to the solution of the associated stationary MFG.

(

  (a) iterations to decrease the residual by a factor 10 -3 .

  ν

Table 1 .

 1 MATLAB's backslash computation times (seconds) for a single linear system solved in (4.9) within the Chambolle-Pock algorithm under a tolerance equal to 10 -4 in in normalized 2 -norm. For fine meshes [OOM] indicates an out of memory error for the tested architecture.

	2 Intel Core i7-4600U @ 2.7GHz, 16GB RAM

1 http://uk.mathworks.com/help/matlab/ref/mldivide.html

Table 3 .

 3 a) Grid with 64 × 64 × 64 points. Time (in seconds) for the convergence of the Chambolle-Pock algorithm, cumulative time of the first proximal operator with the multigrid preconditioner, and number of iterations, for different viscoty values ν and two types of grids. Here we used η 1 = η 2 = 2, T = 1 and a tolerance between two iterations of the Chambolle-Pock algorithm equal to 10 -6 in normalized 2 -norm.

						(b) Grid with 128 × 128 × 128 points.
	ν	Total time Time first prox Iterations	ν	Total time Time first prox Iterations
	0.6	116.3 [s]	11.50 [s]	20	0.6	921.1 [s]	107.2 [s]	20
	0.36	120.4 [s]	11.40 [s]	21	0.36	952.3 [s]	118.0 [s]	21
	0.2	119.0 [s]	11.26 [s]	22	0.2	1028.8 [s]	127.6 [s]	22
	0.12	129.1 [s]	14.11 [s]	22	0.12	1036.4 [s]	135.5 [s]	23
	0.046	225.0 [s]	23.28 [s]	39	0.046 1982.2 [s]	260.0 [s]	42

Table 4 .

 4 Average number of iterations of the preconditioned BiCGStab with η 1 = η 2 = 2, T = 1 and a tolerance between two iterations of the Chambolle-Pock algorithm equal to 10 -6 in normalized 2 -norm.

					(b) iterations to solve the system
					with an error of 10 -8 .
		32 × 32 × 32	64 × 64 × 64	128 × 128 × 128	ν	32 × 32 × 32	64 × 64 × 64	128 × 128 × 128
	0.6	1.65	1.86	2.33	0.6	3.33	3.40	3.38
	0.36	1.62	1.90	2.43	0.36	3.10	3.21	3.83
	0.2	1.68	1.93	2.59	0.2	3.07	3.31	4.20
	0.12	1.84	2.25	2.65	0.12	3.25	3.73	4.64
	0.046	1.68	2.05	2.63	0.046	2.88	3.59	4.67

Acknowledgments. The third author wants to acknowledge funding within the ANR project MFG ANR-16-CE40-0015-01 operated by the French National Research Agency (ANR).

Most of this work was realized while the fourth author was a postdoctoral fellow at NYU Shanghai and supported by a discretionary research fund.

During the first phase of the project, the fifth author was affiliated to the Unité de Mathématiques Pures et Appliquées (UMPA) UMR CNRS 5669, of the École Normale Supérieure de Lyon and to the Project-Team Beagle of the Inria Rhône-Alpes. He wishes to acknowledge funding within the framework of the LABEX MILYON (ANR-10-LABX-0070) of the Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007). In addition, his participation to this project has been