Enrico Facca
email: enrico.facca@uib.no

Gabriele Todeschi
email: gabriele.todeschi@univ-grenoble-alpes.fr

Andrea Natale
email: andrea.natalea@inria.fr

Michele Benzi
email: michele.benzi@sns.it

EFFICIENT PRECONDITIONERS FOR SOLVING DYNAMICAL OPTIMAL TRANSPORT VIA INTERIOR POINT METHODS

Keywords: Optimal transport, Benamou-Brenier formulation, Saddle point problem, Algebraic multigrid methods, Preconditioners

In this paper we address the numerical solution of the quadratic optimal transport problem in its dynamical form, the so-called Benamou-Brenier formulation. When solved using interior point methods, the main computational bottleneck is the solution of large saddle point linear systems arising from the associated Newton-Raphson scheme. The main purpose of this paper is to design efficient preconditioners to solve these linear systems via iterative methods. Among the proposed preconditioners, we introduce one based on the partial commutation of the operators that compose the dual Schur complement of these saddle point linear systems, which we refer as BB-preconditioner. A series of numerical tests show that the BB-preconditioner is the most efficient among those presented, despite a performance deterioration in the last steps of the interior point method. It is in fact the only one having a CPU-time that scales only slightly worse than linearly with respect to the number of unknowns used to discretize the problem.

Introduction

Optimal transport deals with the problem of finding the optimal way to reallocate one nonnegative density into another by minimizing the total cost of displacement in space. In recent years, numerous contributions have been made to the study of this problem, both on the theoretical and computational level. We suggest, for example, the monographs [START_REF] Villani | Topics in Optimal Transportation[END_REF][START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF][START_REF] Ambrosio | Lectures on Optimal Transport[END_REF][START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF] for a detailed presentation of the subject. Due to these advances, optimal transport is nowadays an established tool for many applications including, for example, the analysis of partial differential equations (PDE) [START_REF] Ambrosio | Gradient flows: in Metric Spaces and in the Space of Probability Measures[END_REF], physical modeling [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF], data science and machine learning [START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF], economics [START_REF] Galichon | Optimal transport methods in economics[END_REF], and inverse problems [START_REF] Métivier | A review of the use of optimal transport distances for high resolution seismic imaging based on the full waveform[END_REF].

When the cost of displacement per unit mass is given by the square of the Euclidean distance, the problem can be recast dynamically, as shown by Benamou and Brenier [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF]. Consider a compact and convex domain Ω ⊂ R d and two nonnegative densities ρ in and ρ f in L 1 (Ω), with Ω ρ in dx = Ω ρ f dx. To transport the former to the latter, we aim to find a time-dependent density ρ : [0, 1] × Ω → R ≥0 and a velocity field v : [0, 1] × Ω → R d that solve the following minimization problem:

min ρ,v 1 0 Ω ρ|v| 2 2 dt dx :      ∂ t ρ + div(ρv) = 0 in [0, 1] × Ω ρv • n = 0 on [0, 1] × ∂Ω ρ(0, •) = ρ in , ρ(1, •) = ρ f . (1)
The total kinetic energy represents the cost of displacement. Thanks to a change of variables (ρ, v) → (ρ, m = ρv),

(1) can be rewritten as a convex optimization problem. From the optimality conditions, one can deduce that the optimal velocity field is the gradient of a potential ϕ : [0, 1] × Ω → R. The potential ϕ and the density ρ are given as the solution of the following system of PDEs:

-∂ t ρ -div (ρ ∇ ϕ) = 0, (2a)
∂ t ϕ + ∥ ∇ ϕ∥ 2 2 + s = 0, (2b)
ρ ≥ 0, s ≥ 0, ρs = 0,

with boundary conditions ρ(0,

•) = ρ in , ρ(1, •) = ρ f , ρ ∇ ϕ • n = 0 on [0, 1] × ∂Ω.
The auxiliary variable s : [0, 1] × Ω → R ≥0 is related to the positivity constraint on ρ.

The Benamou-Brenier formulation offers several advantages. In addition to the total displacement cost, it also directly provides how to continuously reallocate the mass. This also constitutes a natural way to define interpolations between densities. Furthermore, it draws a clear link between optimal transport and continuum mechanics, and it is naturally suited for Eulerian discretizations. Finally, it can be easily generalized to other problems by penalizing/constraining the evolution ρ (such as variational mean field games and planning problems [START_REF] Achdou | Mean field games: numerical methods for the planning problem[END_REF], or unbalanced optimal transport [START_REF] Chizat | An interpolating distance between optimal transport and fisher-rao metrics[END_REF][START_REF] Baradat | Regularized unbalanced optimal transport as entropy minimization with respect to branching brownian motion[END_REF]). On the other hand, the numerical solution of (1), or its system of optimality conditions [START_REF] Achdou | Mean field games: numerical methods for the planning problem[END_REF], poses significant challenges. Although it is a convex optimization problem, it is nonlinear and, in general, non-smooth for vanishing densities. Moreover, it is a time-space boundary value problem and there is a positivity constraint on ρ to take into account.

While different strategies have been proposed to discretize the Benamou-Brenier formulation, most of these rely on staggered time discretization and a space discretization which may be based on finite differences [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF][START_REF] Papadakis | Optimal transport with proximal splitting[END_REF], finite volumes [START_REF] Gladbach | Scaling limits of discrete optimal transport[END_REF][START_REF] Natale | Computation of optimal transport with finite volumes[END_REF][START_REF] Lavenant | Unconditional convergence for discretizations of dynamical optimal transport[END_REF][START_REF] Haber | A multilevel method for the solution of time dependent optimal transport[END_REF] or finite elements [START_REF] Lavenant | Dynamical optimal transport on discrete surfaces[END_REF][START_REF]A mixed finite element discretization of dynamical optimal transport[END_REF][START_REF] Lavenant | Unconditional convergence for discretizations of dynamical optimal transport[END_REF]. Here, we focus on the framework considered in [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF], where one uses finite volumes in space with a two-level discretization of the domain Ω, in order to discretize the density ρ and the potential ϕ separately, which alleviates some checkerboard instabilities that may appear when the same grid is used to discretize both variables (a strategy first used in [START_REF] Facca | Towards a stationary Monge-Kantorovich dynamics: The Physarum Polycephalum experience[END_REF][START_REF] Facca | Numerical solution of Monge-Kantorovich equations via a dynamic formulation[END_REF] when dealing with optimal transport with unitary displacement cost given by the Euclidean distance). This choice is not restrictive since such a framework constitutes a generalization of the finite volume scheme studied in [START_REF] Gladbach | Scaling limits of discrete optimal transport[END_REF][START_REF] Lavenant | Unconditional convergence for discretizations of dynamical optimal transport[END_REF] (in which a single grid is used for density and potential), and it also contains as special cases the finite difference scheme in [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF] (when a single Cartesian grid is used in space), or the discrete transport models on networks studied in [START_REF] Erbar | Computation of optimal transport on discrete metric measure spaces[END_REF] (by an appropriate reinterpretation of the space tessellation).

From a computational point of view, the numerical solution of the resulting discrete optimization problem is generally tackled by iterative first-order primal-dual optimization schemes (see, e.g., [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF][START_REF] Papadakis | Optimal transport with proximal splitting[END_REF]). The computational cost per iteration of these methods is typically low but the total number of iterations strongly depends on the data and the tuning of the parameters, generally increasing with the problem size. In this paper we will consider instead a secondorder approach, an Interior Point (IP) method proposed in [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF], where the optimization problem eq. ( 1) is relaxed by adding a logarithmic barrier function (we refer the reader to [START_REF] Wright | Primal-Dual Interior-Point Methods[END_REF][START_REF] Boyd | Convex Optimization[END_REF][START_REF] Gondzio | Interior point methods 25 years later[END_REF] for a broad introduction of IP methods). This perturbation provides smoothness by enforcing the strict positivity of the density and uniqueness of the solution.

The problem can then be effectively solved by solving the relaxed optimality conditions via the Newton method, and the original unperturbed solution is retrieved by repeating this procedure while reducing the relaxation term to zero. Numerical experiments in [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF] show that the total number of Newton iterations remains practically constant, independently of the number of unknowns.

The most demanding task in IP methods is the solution of a sequence of saddle point linear systems in the following form

A B T B -C x y = f g ,
generated by the Newton method. In our problem all matrices A, B, and C change at each Newton iteration but they are sparse, hence iterative solvers are the natural candidate for solving these linear systems. Iterative methods were not considered in [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF], and thus the goal of this paper is to determine whether one can devise a preconditioning strategy such that, ideally, the total CPU time scales linearly with the number of unknowns.

We consider three preconditioners: a preconditioner based on the approximation of the inverse of the primal Schur complement S p = A + BC -1 B T , the SIMPLE preconditioner [START_REF] Patankar | A calculation procedure for heat, mass and momentum transfer in threedimensional parabolic flows[END_REF], and one preconditioner inspired by the ideas in [START_REF] Elman | Block preconditioners based on approximate commutators[END_REF], where we approximate the inverse of the dual Schur complement S = -C -BA -1 B T looking at the differential nature of the operators that compose S (note that, formally, we cannot write A -1 since in our problem the matrix A will be singular). A series of numerical experiments suggest that the latter, named BB-preconditioner, is the most efficient among those presented in this paper, despite some loss of robustness in the last IP iterations.

Let us stress that, while we conducted our experiments using the scheme described in [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF], we do expect that the BBpreconditioner will provide good results when adopting different discretization schemes even outside the framework considered here, e.g. when finite elements are used. This is because it is entirely derived from the differential operators involved in the continuous problem.

Related works

The literature on solving eq. ( 2) using second-order and/or Newton-based methods is relatively scarce. The most closely related works are [START_REF] Perez | Iterative strategies for solving linearized discrete mean field games systems[END_REF][START_REF] Haber | A multilevel method for the solution of time dependent optimal transport[END_REF]. In the first paper, the authors study linear algebra approaches involved in the solution of mean-field games systems via the Newton method. Problem (1) may in fact be seen as the vanishing viscosity limit of a particular mean-field game. However, their linear algebra approaches do not fit into the staggered temporal grids, which is generally required to discretize the Benamou-Brenier formulation. Moreover, they lose efficiency with vanishing viscosity. In the second paper, the authors solved precisely the Benamou-Brenier problem using an inexact Newton method combined with finite volumes on Cartesian grids. However, the discretization scheme is not designed to handle densities ρ in and ρ f with compact support. The solution of the sequence of saddle-point linear systems associated to the Newton method was performed using a preconditioning strategy that resembles the approach based on the approximation of the primal Schur complement, dropping some of the mixed time-space terms. However, this approach did not give satisfactory results in our experiments (more details in section 4.2).

Paper structure

In section 2 we summarize the discretization method proposed in [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF], the IP approach used to solve the optimization problem, and the nonlinear problems associated with it. Then, in section 3, we present the linear system studied in this paper. Finally, in section 4 we present the preconditioners described in this paper, together with some numerical experiments where we solve a specific test case for different time and space refinements. In particular, in section 4.5 we compare the CPU time required by the preconditioners in solving different test cases.

Discrete problem and Interior Point method

In this section we present the discrete counterpart of system eq. ( 2) considered in [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF]. This combines the use of finite volumes with staggered temporal grids. We further present the IP strategy adopted for solving the problem.

Spatial discretization

In [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF], the authors considered two different discretizations of the domain Ω, which is assumed to be polygonal, both admissible for Two Point Flux Approximation (TPFA) finite volumes, according to [24, Definition 9.1]. At each time t ∈ [0, 1], the variable ρ(t) is discretized on a Delaunay triangulation, with the further hypothesis that only acute angles appear. The variable ϕ(t) is discretized on a finer grid obtained from the previous one by dividing each triangular cell into three quadrilateral cells, joining the edges midpoints to the triangle's circumcenter [38, Figure 1]. Both meshes consist of two sets (T , E), the set of cells c and edges e, respectively. To distinguish the coarser mesh from the finer one, we denote the former by (T ′ , E ′ ), with the same notation for all its elements. Due to the no-flux boundary condition, boundary edges are not relevant to the discrete model. We will then consider, by convention, the sets E and E ′ without boundary edges. Let us denote by N T ′ and N E ′ the total number of cells and (internal) edges of the coarser mesh. The total number of cells and edges of the finer mesh are then

N T = 3N T ′ and N E = 2N E ′ + 3N T ′ .
We define

M := Diag (|c|) , |c| := (|c i |) N T i=1 , M ′ := Diag (|c ′ |) , |c ′ | := (|c i ′ |) N T ′ i=1 ,
where |c i | denotes the area of the triangle c i . Since the discrete model involves two different spatial discretizations, we also introduce the injection operator J ∈ R N T ,N T ′ ,

J [i, j] = 1 if c i ⊂ c ′ j , 0 else . (3)
The operator J injects piecewise constant functions from the coarse grid to the fine one. The operator (M ′ ) -1 J T M is an averaging operator that acts in the opposite direction. Both operators are the identity when the same grid is used for ϕ and ρ.

In order to define the discrete differential operators associated to the finite volume discretization, we fix an arbitrary orientation on the set of edges E and define the matrix E ∈ R N T ,N E given by

E[i, k] = 1 if cell c i is the left side of edge e k , -1 if cell c i is the right side of edge e k .
The discrete gradient and divergence, ∇ ∈ R N E ,N T and div ∈ R N T ,N E , are given by

∇ = ∇ E,T = Diag (|w|) -1 E T , div = div T ,E = -∇ T Diag (|w| ⊙ |e|) = -EDiag (|e|) ,
where |w| ∈ R N E and |e| ∈ R N E are the vectors of distances between the circumcenters of adjacent cells and the lengths of the edges, respectively (we used the symbol ⊙ to denote the pointwise multiplication of vectors).

The TPFA finite volume discretization and the two-level grids require the introduction of two additional operators to match the dimension between the different spaces. The first is a reconstruction operator R E : R N T ′ → R N E , mapping a positive density ρ defined on the cells of the coarse grid to a variable defined on the edges of the finer one. In [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF] this reconstruction first lifts the density ρ into the finer space via the operator J and then averages the values of adjacent cells. The authors considered two types of averages, a weighted arithmetic mean and a weighted harmonic mean (the latter is not considered in this paper to avoid complicating the exposition). Using the former, the linear reconstruction operator explicitly writes (R E ρ) e := λ e (J ρ) ci + (1λ e )(J ρ) cj for the two cells c i , c j sharing the edge e, where λ e is a weight that depends on the mesh geometry. Moreover, they introduced a further operator R T : R N E → R N T ′ that maps the variables defined at the edges of the finer grid to the variables defined at the cells of the coarser one. To preserve the variational structure of the discrete problem, it is defined as

R T := (R E ) T Diag (|w| ⊙ |e|) .

Temporal discretization

In [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF], the temporal discretization is based on staggered temporal grids. The time interval [0, 1] is divided into K + 1 sub-intervals [t k , t k+1 ] of equal length ∆t = 1/(K + 1), for k = 0, . . . , K ≥ 1, with t 0 = 0 and t K+1 = 1. The variables ρ and s are defined at time t k for k = 1, . . . , K, while ϕ is discretized at each instant (t k + t k+1 )/2 for k = 0, . . . , K.

Combining this temporal and spatial discretization, the discrete counterpart of the potential ϕ, the density ρ, and the slack variable s in equation eq. ( 2) are the vectors

ϕ ∈ R n , ρ ∈ R m ≥0 , and s ∈ R m ≥0 with n = N T (K + 1) and m = N T ′ K given by ϕ = ϕ 1 ; . . . ; ϕ K+1 , ϕ k ∈ R N T , ρ = ρ 1 ; . . . ; ρ K , ρ k ∈ R N T ′ ≥0 , s = s 1 ; . . . ; s K , s k ∈ R N T ′ ≥0 ,
where we use the symbol ; to denote the concatenation of vectors. Moreover, we will denote by x k the k-slice of a concatenated vector x, which corresponds to its k-th time portion. The discrete counterparts of the initial and final densities are the two vectors ρ 0 and ρ K+1 in R N T ′ given by

ρ 0 i = 1 |c i ′ | ci ′ ρ in dx , ρ K+1 i = 1 |c i ′ | ci ′ ρ f dx, i = 1, . . . , N T ′ .

Discrete nonlinear system of equations

With the aforementioned discretization, the discrete counterpart of the nonlinear system eq. ( 2) is given by

F ϕ (ϕ, ρ) := F 1 ϕ ; . . . ; F K+1 ϕ =0 ∈ R N T (K+1) , F ρ (ϕ, ρ, s) := F 1 ρ ; . . . ; F K ρ =0 ∈ R N T ′ K , F s (ρ, s) := F 1 s ; . . . ; F K s =0 ∈ R N T ′ K , (4)
where for k = 1, . . . , K + 1 the functions F k ϕ are given by

F k ϕ (ϕ, ρ) = -M J ρ k -ρ k-1 ∆t -div R E ρ k +ρ k-1 2 ⊙ ∇ϕ k , (5)
while for k = 1, . . . , K the functions F k ρ and F k s are given by

F k ρ (ϕ, ρ, s) = M ′ J T ϕ k+1 -ϕ k ∆t + 1 4 R T ∇ϕ k 2 + ∇ϕ k+1 2 + M ′ s k , (6)
F k s (ρ, s) = ρ k ⊙ s k , (7)
with the additional (component-wise) requirement ρ, s ≥ 0 . System eq. ( 4) is the system of optimality conditions for the discrete counterpart of the problem in eq. ( 1). Remark 1. Thanks to the finite volume discretization, the conservative structure of the continuity equation is preserved. Then, due to the no-flux boundary conditions, which are encoded explicitly in the divergence operator, the equation

F k ϕ (ϕ, ρ) = 0 implies 1 T F k ϕ = |c| T J ρ k -ρ k-1 = |c ′ | T ρ k -ρ k-1 = 0 ∀k = 1, . . . , K .
Hence, at each intermediate time step, the discrete mass |c ′ | T ρ k is preserved and is equal to the mass of the discrete initial and final densities.

Interior Point method

Due to the IP strategy, the discrete optimization problem is perturbed by adding a logarithmic barrier function,

-µ K k=1 ∆t N T ′ i=1 log(ρ k i )|c i |
tuned by a parameter µ > 0. This turns the system of optimality conditions eq. ( 4) into

F ϕ (ϕ, ρ) = 0 ∈ R N T (K+1) , F ρ (ϕ, ρ, s) = 0 ∈ R N T ′ K , F s (ρ, s) -µ1 = 0 ∈ R N T ′ K , (8)
where the complementarity constraint between ρ and s in eq. ( 7) is now relaxed. In this way, the two variables are forced to be strictly positive, making the problem easier to solve by using Newton methods.

The solution of the original problem is recovered in the limit µ → 0 and the parameter µ directly provides a bound on the sub-optimality of the discrete solution of eq. ( 8) for the original unperturbed problem [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF]Section 4]. From the practical point of view, it is computed via a continuation method: for each value of µ l of a sequence of relaxation parameters (µ l ) l≥0 , the problem is solved with a Newton method initialized with the previously computed solution.

We refer to [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF] for the precise set up of the algorithm. We stress only that for each value of µ, problem eq. ( 8) is solved up to the relatively low tolerance 1e -6 in order to obtain a more robust implementation. Furthermore, within each Newton cycle, we ensure that any Newton update does not introduce negative entries in the vectors ρ and s with a simple line-search procedure.

3 The linear algebra problem and the KKT system

The nonlinear system of equations eq. ( 8) is solved using an inexact Newton method. Each Newton iteration requires the solution of a linear system (referred to as KKT system) in the form

  A B T B M ′ Diag (s) Diag (ρ)   δϕ δρ δs = f g h = - F ϕ F ρ F s -µ1 , (9)
where the block matrix in equation eq. ( 9) is the Jacobian matrix of (F ϕ ; F ρ ; F s -µ1). We denote by M ′ and M the matrices given by

M ′ = Block Diag (M ′ ) K k=1 , M = Block Diag (M ) K+1 k=1 .
The matrices A, B, B T in eq. ( 9) are the finite-dimensional version of the following differential operators (up to a multiplication by a mass matrix)

A ≈ -div(ρ ∇) , B ≈ ∂ t + ∇ ϕ • ∇ , B T ≈ -∂ t -div(• ∇ ϕ) .
According to the discretization described in section 2, matrix A ∈ R n,n is a symmetric block diagonal matrix with K + 1 blocks. Each block is a weighted Laplacian matrix given by

A k = -divDiag ρk ∇ ∈ R N T ,N T , ρk := R E ρ k + ρ k-1 2 ,
which we can write equivalently as A = -Div x Diag ( ρ1 ; . . . ; ρk ; . . . ; ρK+1 ) D x where

Div x = Block Diag (div) K+1 k=1 , D x = Block Diag (∇) K+1 k=1 .
Matrix B is a block bi-diagonal matrix in R m,n given by

B = J T D t M + HGD x ,
where the matrices

D t ∈ R KN T ,(K+1)N T , H ∈ R KN T ′ ,(K+1)N T ′ , J ∈ R (K+1)N T ,(K+1)N T ′ , and G ∈ R (K+1)N T ′ ,(K+1)N E ′ are D t = 1 ∆t    -I N T I N T . . . . . . -I N T I N T    , J = Block Diag (J ) K+1 k=1 , (10)
H = 1 2    I N T ′ I N T ′ . . . . . . I N T ′ I N T ′    , G = Block Diag (G k ) K+1 k=1 , (11)
with (G k ) k=1,...,K+1 ∈ R N T ′ ,N E given by

G k := R T Diag ∇ϕ k .
Remark 2. The kernel of the block diagonal matrix A has dimension K + 1, since each block A k has the constant vectors as kernel. On the other hand, Ker(A) is included in Ker(HGD x ) and therefore

Ker(A) ∩ Ker(B) = Ker(A) ∩ Ker(J T D t M) = ⟨1⟩ ∈ R n ,
thus the Jacobian matrix is singular. This singularity can be removed grounding one entry of the solution δϕ or via regularization techniques [START_REF] Lehoucq | On the finite element solution of the pure neumann problem[END_REF]. However, we prefer avoiding these approaches in this paper, since iterative methods work also when dealing with singular matrices [START_REF] Van Der | Iterative Krylov Methods for Large Linear Systems[END_REF][START_REF] Elman | Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics[END_REF].

Reduction to a saddle point linear system

In order to solve the linear system in eq. ( 9) it is convenient to eliminate the variable δs writing

δs = (Diag (ρ)) -1 (h -Diag (s) δρ)
and reduce it into the following saddle point system

J δϕ δρ = A B T B -C δϕ δρ = f g , (12)
where the matrix C ∈ R m,m and the vector g ∈ R m are given by

C := M ′ Diag (ρ) -1 Diag (s) , g := g -M ′ Diag (ρ) -1 h .
The matrix J in eq. ( 12) has typically higher conditioning number than the Jacobian matrix in eq. ( 9), in particular when approaching the optimal solution, i.e. as s ⊙ ρ ≈ µ1 → 0. However, designing an efficient preconditioner for the fully coupled system in eq. ( 9) can be harder than for the standard saddle point linear system in eq. ( 12), for which different preconditioning approaches are described in [START_REF] Benzi | Numerical solution of saddle point problems[END_REF][START_REF] Bergamaschi | Preconditioning indefinite systems in interior point methods for optimization[END_REF][START_REF] Morini | A comparison of reduced and unreduced KKT systems arising from interior point methods[END_REF].

We adopt an Inexact Newton approach, which means that we seek for a solution (δϕ, δρ) such that

∥J(δϕ; δρ) -(f ; g)∥ ≤ ε out ∥(f ; g; h)∥ .
The linear system residual is scaled by the right-hand side of the original system to avoid over-solving issues we faced using ∥(f ; g)∥.

In our experiments we use a fixed tolerance ε out = 1e -5. This value may be quite small compared to the literature for an inexact Newton approach (for example [START_REF] Bellavia | Inexact interior-point method[END_REF][START_REF]Convergence analysis of an inexact feasible interior point method for convex quadratic programming[END_REF]), but we want to avoid more aggressive inexact strategies which could undermine the effectiveness of the nonlinear solver and focus on the linear algebra only.

Remark 3. From remark 1, we know the solution ρ to equation eq. ( 8)

must satisfy |c ′ | T ρ k = |c ′ | T ρ 0 for all k = 1, . . . K.
If this condition is satisfied for the starting point of the Newton scheme, then the increment δρ solution to equations [START_REF] Benzi | Numerical solution of saddle point problems[END_REF]

satisfies |c ′ | T δρ k = 0 ∀k = 1, . . . K ,
at each Newton iteration. To see this, sum both sides of the first block of equations of system [START_REF] Bergamaschi | Preconditioning indefinite systems in interior point methods for optimization[END_REF], for each time step k, and use the no-flux boundary conditions in the divergence operator. Due to our inexact Newton steps, this condition is not verified exactly, but only up to the tolerance chosen. We impose it by renormalizing ρ after each Newton update.

Preconditioning approaches and numerical results

In this section we present three preconditioners to solve the linear system in eq. ( 12) and a series of numerical experiments used to measure their performance. All of these preconditioners involve the usage of inexact inner solvers, therefore, they may not be linear and must be applied as right preconditioners within a FGMRES cycle [START_REF] Saad | A flexible inner-outer preconditioned GMRES algorithm[END_REF]. In all cases, the inner solver we adopted is the AGgregation-based Algebraic MultiGrid (AGMG) method described in [START_REF] Notay | An aggregation-based algebraic multigrid method[END_REF].

We adopt three measures to evaluate the preconditioner performance for each Newton cycle associated with each IP iteration:

1. Outer/Lin.sys.: the average number of outer FGMRES iterations. It is calculated as the number of outer iterations divided by the number of linear systems solved for each Newton cycle.

2. CPU/Lin.sys.: the average CPU time (measured in seconds) required for the solution of the linear systems.

It is computed as the total CPU time required to solve all linear systems within a Newton cycle divided by the number of linear systems solved (all experiments were conducted single core, on a machine equipped with an Intel ® Xeon ® 2.8 GHz CPU and 128 GiB of RAM memory).

3. Inner/Outer: the average number of inner iterations per preconditioner application. It is computed by dividing the cumulative number of AGMG iterations required to solve the linear systems involved in the preconditioner application (for each preconditioner, we will specify the linear systems to which we will refer) by the total number of outer iterations for each Newton cycle.

Test cases

The IP algorithm and the linear solver strategies described in the next sections are tested on four numerical experiments adapted from [38, Section 5.2]. The first three are set in Ω = [0, 1] 2 whereas the fourth one, three dimensional, in Ω = [0, 1] 3 . They are:

1. Gaussian densities ρ in and ρ f centered at (0.3, 0.3) and (0.7, 0.7) respectively with variance 0.1. This ensures that ρ in and ρ f are lower bounded by 1e -3.

2. Translation of a compactly supported smooth sinusoidal density ρ in .

3. Compression of a compactly supported smooth sinusoidal density ρ in . This test case has been shown experimentally in [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF] to exhibit severe instabilities in the discrete solution, and this, in turn, has a negative effect on the discrete solver.

4. Same as test 1 but in three dimensions. The Gaussian densities ρ in and ρ f are centered at (0.3, 0.3, 0.3) and (0.7, 0.7, 0.7).

Test cases 1, 2 and 3 are discretized using four different meshes, taken from [1], having N T ′ = 224, 896, 3584, 14336 cells, where each mesh is a refinement of the other. These meshes will be denoted by T 0 , T 1 , T 2 , T 3 . The corresponding number of subcells is N T = 672, 2688, 10752, 43008. For each grid, we consider ∆t = 1/(K + 1) with K + 1 = 16, 32, 64, 128. The typical mesh size h of the mesh T 0 is approximately 1/16, so that the diagonal pairing of the space and time steps provides a uniform discretization of the time-space domain [0, 1] × Ω. The total number of ϕ and ρ degrees of freedom for these sequence of problems is n + m ≈ (1.4e4, 1.1e5, 9.1e5, 7.3e6), which is the size of the linear system in eq. ( 12). The last test case is solved on Cartesian grids with cubic cells. In this case the operator J ineq. ( 3) is simply the identity. We consider 2 refinements (both in time and in space) of an initial experiment set on a 8 × 8 × 8 grid and using ∆t = 1/(K + 1) with K + 1 = 8. In this case, the total number of ϕ and ρ degrees of freedom for a uniform time-space refinement is n + m ≈ (7.7e3, 1.3e5, 2.0e6).

Convergence is achieved when the relaxation parameter µ is below 1e -6, which corresponds to ten iterations of IP, with the choice of parameters described in [START_REF] Natale | Computation of optimal transport with finite volumes[END_REF]. The nonlinear system eq. ( 8) is solved with rather tight tolerance 1e -6 in order to have a robust implementation for the IP solver. Each IP step requires between 3 and 7 inexact Newton iterations. These numbers are practically insensitive to the mesh and the time step adopted. The total number of Newton steps (which corresponds to the number of linear systems to be solved) ranges between 40 and 50. The preconditioning approaches used to solve these linear systems are described in the following sections.

For each preconditioner, we will present a table reporting the metrics mentioned above, while the relaxation parameter µ is reduced for all combinations of spatial and temporal discretizations. We will show these detailed results only for the second test case, since it captures the main challenges involved in the linear algebra problems. The numerical results for the other three cases will be summarized in section 4.5, where we will compare the performance of the preconditioners in terms of total CPU time.

The code is written in MATLAB ® and its source code is available at the link https://github.com/gptod/OT-FV, where it is possible to reproduce the experiments presented in this paper. The AGMG solver is interfaced via MEX ® .

Preconditioner based on the primal Schur complement

In this section we consider a preconditioner for the linear system in eq. ( 12) based on the primal Schur complement S p = A+B T C -1 B, similarly to the approach used in [START_REF] Facca | Fast iterative solution of the optimal transport problem on graphs[END_REF] for the solution of the optimal transport problem on graphs with cost equal to the shortest path distance. It is given by

P -1 = I C -1 B I S -1 p -C -1 I B T C -1 I . ( 13
)
The application of the preconditioner in eq. ( 13) requires the inversion of C, which is diagonal, and the (approximate) solution of the linear system

S p x = b , (14)
with x, b ∈ R n . Using this preconditioner is essentially equivalent to reducing the linear system in eq. ( 12) to the variable δϕ only, an approach referred to in the literature as fully reduced [START_REF] Benzi | Multilevel algorithms for large-scale interior point methods[END_REF] or condensed [START_REF] Apuzzo | On mutual impact of numerical linear algebra and large-scale optimization with focus on interior point methods[END_REF]. The matrix S p is a block tridiagonal matrix. It can be written as

S p = A -Div x Diag ( ρ) D x =:Sxx -Div x G T H T C -1 HGD x + =:Stt MD t T J C -1 J T D t M + MD t T J C -1 HGD x + (MD t T J C -1 HGD x ) T =:Stx+S T tx ,
and it may be seen as the discretization of an elliptic time-space operator. In fact, it is the sum of a ρ-weighted spatial Laplacian matrix A, an anisotropic spatial Laplacian S xx ≈div(ρ/s ∇ ϕ ⊗ ∇ ϕ ∇), a weighted temporal Laplacian with Neumann boundary condition S tt , and a time-space operator S tx + S T tx . We adopt the AGMG solver to solve the linear system in eq. ( 14), with relative accuracy ε in = 1e -1 > ε out = 1e -5. The value of the inner tolerance ε in has been determined experimentally to distribute the work load between the inner and outer solver in a balanced way.

In table 1 we summarize the results obtained using the preconditioner presented in this section for different combinations of time and space discretizations. We use the metrics described in section 4, where the average number of inner iterations, denoted by Inner/Outer, refers to the linear system in eq. ( 14).

The number of average outer iterations, denoted by Outer/Lin.sys., remains practically constant only for the coarsest and most uniform combinations of grids and time steps. However, in some cases, also Outer/Lin.sys. increases when ). Each sub-table reports, while µ is reduced (leftmost column), the averaged outer iteration and CPU time per linear system (Outer/Lin.sys. and CPU/Lin.sys.) and the averaged inner iterations per outer iteration Inner/Outer for solving the linear system in eq. ( 14) (metrics defined in section 4). A final row summarizes the averages on the whole simulation. We highlighted in gray the time-space combination providing the uniform discretization. The † symbol denotes those IP steps where the linear solver failed.

the average number of inner iterations Inner/Outer get closer to 200 (the limit we fixed for the inner solver), which means that, in some cases, the inner solver did not reach the accuracy ε in = 1e -1. For the cases with more degrees of freedom, this led to a huge CPU time (almost 1e5 seconds per linear system) and to failures occurring as µ is reduced.

The average number of inner iterations Inner/Outer is strongly influenced by the relationship between the mesh and the time step ∆t = 1/(K + 1). The phenomenon is particularly evident for mesh T 2 . For K + 1 = 32, the average number of inner iterations increases progressively while µ is reduced, reaching ≈ 179 at the last IP iteration. For K + 1 = 128, Inner/Outer increases reaching the 200 iterations limit. For K + 1 = 64 instead, which corresponds to the most balanced time-space scaling, this number remains between 2.0 and 2.7. The same phenomenon occurs using T 1 , passing from K + 1 = 16 to K + 1 = 32. We attribute this to an improper construction of the coarse matrix sequence used by the multigrid solver due to the poor scaling of the temporal and spatial components of S p . The AGMG is based on an aggregation-based coarsening approach, hence if these components are unbalanced the coarse matrices may not reflect the infinite-dimensional operators discretized by our problem, compromising the performance of the multigrid solver.

This phenomenon can only get worse in the last IP iterations due to presence of the term

C -1 = M ′ -1 Diag (ρ) Diag (s) -1 in S p = A + BC -1 B T .
In fact, the term C -1 contains entries varying by several order of magnitude due to the relaxed complementarity condition s ⊙ ρ ≈ µ1. This also means that, on those entries where ρ remain strictly positive, the anisotropic Laplacian S xx scales like 1/µ, making the linear system eq. ( 14) more difficult to solve since this anisotropic term becomes the dominant part of S p . Similar matrices and analogous considerations can be found in [START_REF] Benzi | Multilevel algorithms for large-scale interior point methods[END_REF][START_REF]A preconditioning technique for a class of pde-constrained optimization problems[END_REF][START_REF] Haber | A multilevel method for the solution of time dependent optimal transport[END_REF][START_REF] Morini | A comparison of reduced and unreduced KKT systems arising from interior point methods[END_REF]. Furthermore, even in best-case scenarios, we experimented that the preprocessing phase is the most demanding in terms of CPU time. Typically, more than 50% of the time is spent building the sequence of coarse matrices, reaching peaks of 90% during the last IP steps using the finest discretization.

In order to cope with all the issues presented so far, we tried to neglect some components of S p . We tried to remove the extra-diagonal or lower-diagonal part of S p , or the time-space operator S tx + S T tx so that the remaining terms form a weighted time-space Laplacian (similarly to what done in [START_REF] Haber | A multilevel method for the solution of time dependent optimal transport[END_REF]). However, none of these approaches worked. All numerical experiments suggest that none of these terms can be neglected without altering the spectral property of the matrix, and consequently affecting the effectiveness of the preconditioner.

Summarizing the results obtained in this section, the preconditioner based on the primal Schur complement S p = A + B T C -1 B lacks both in efficiency and in robustness, due to the complexity of the PDE behind the linear system in eq. ( 14), in particular when the relaxation parameter µ goes to zero.

The SIMPLE preconditioner

In this section, we recall the classical SIMPLE preconditioner described in [START_REF] Patankar | A calculation procedure for heat, mass and momentum transfer in threedimensional parabolic flows[END_REF][START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF] and based on the approximation of the inverse of the dual Schur complement S = -(C + BA -1 B T ) (we recall that, formally, we cannot write A -1 since the matrix A is singular). The SIMPLE preconditioner is given by

P -1 = I - Â-1 B T I Â-1 Ŝ-1 I -B Â-1 I where  = Diag (A) , Ŝ = -(C + B Â-1 B T ) .
This preconditioner is referred to as constraint preconditioner in the literature studying linear solvers for IP methods [START_REF] Morini | A comparison of reduced and unreduced KKT systems arising from interior point methods[END_REF]. The main computational cost of applying this preconditioner is to solve linear systems in the form

Ŝy = c . (15)
The matrix Ŝ can be formed explicitly and is block tridiagonal.

During the latest steps of the IP method, when µ → 0, linear system eq. ( 15) may become ill-conditioned since the diagonal matrix C = M ′ Diag (ρ) -1 Diag (s) can contain terms that vary by several orders of magnitude, since ρ ⊙ s ≈ µ1 → 0. To cope with this issue, we scale both sides of eq. ( 15) by Diag (ρ) and we solve the linear system

-(M ′ Diag (s) + Diag (ρ) B Â-1 B T )y = Diag (ρ) c . ( 16
)
This scaling may seem a dangerous procedure if small entries appears in ρ, however it is used within the preconditioner application, where the potential inaccuracies do not affect the convergence of FGMRES to the solution of the linear system. We use the AGMG solver with tolerance ε in = 1e -1 to solve the linear system in eq. ( 16) and denote by Ŝ-1 εin the resulting (nonlinear) operator. We found experimentally that this value provided the best performance. Further reductions of ε in led only to an increase of the number of inner iterations but no reduction in the outer loop iterations.

The results for the SIMPLE preconditioner are summarized in table 2. Despite its simplicity, the proposed preconditioner turned out to be robust (in particular for small values of µ) but not particularly efficient. Some failures occurred only at the initial IP steps using the finest grid. In this cases, we restarted the solver using the data obtained with the BB-preconditioner in order to study the behavior of the SIMPLE preconditioner close to the optimal solution. The average number of inner iterations Inner/Outer increases only sightly as µ → 0, independently of the time step and the mesh size used or the IP step considered. The preprocessing time required by the AGMG solver is limited, approximately accounting for 1% of the CPU time spent in the linear system solution. The average number of outer iterations Outer/Lin.sys. is affected only mildly by the number of time steps K + 1.

Remarkably, the preconditioner becomes more efficient as the IP relaxation term µ goes to zero. However, Outer/Lin.sys. more than doubles at each mesh refinement. Both phenomena can also be explained by looking at the structure of the preconditioned matrix, which is given by

A B T B -C P -1 = I + (A Â-1 -I )(I + B T Ŝ-1 εin B Â-1 ) (I-A Â-1 )B T Ŝ-1 εin (I -Ŝ Ŝ-1 εin )B Â-1 Ŝ Ŝ-1 εin ,
and estimating its eigenvalues. In order to do this, let us first observe that the preconditioned matrix is approximately a block upper triangular matrix, whose bottom-right block becomes close to the identity plus a perturbation of order ε in = 1e -1, hence there are m eigenvalues close to one.

The upper left block is the identity plus a perturbation matrix (which should be as small as possible to have an ideal preconditioner) given by product of two factors. The first factor is

A Â-1 -I = Block Diag A k Diag A k -1 -I N T .
Since the matrices (A k ) k=1,...,K+1 are weighted Laplacian matrices, a diagonal preconditioner becomes a poor preconditioner for large A k (see [START_REF] Kamenski | Conditioning of finite element equations with arbitrary anisotropic meshes[END_REF] for precise estimates). Hence the largest eigenvalue of the first factor will increase while refining in space, but not in time, by the block diagonal structure of matrix A. The second factor, which we would like to keep as small as possible to compensate the first one, is

I + Q with Q = B T Ŝ-1 εin B Â-1
. Using the formula in eq. ( 16) and the complementarity condition ρ ⊙ s ≈ µ1, we can approximate it as follows:

I + Q = I -B T (M ′ Diag (s) + Diag (ρ) B Â-1 B T ) -1 Diag (ρ) B Â-1 ≈ I -B T (M ′ µ + Diag (ρ) 2 B Â-1 B T ) -1 Diag (ρ) 2 B Â-1 .
If we set µ = 0 in this formula, the resulting matrix has 0 or 1 eigenvalues since it is idempotent. This gives a qualitative explanation of why the SIMPLE preconditioner performs better as µ → 0. Nevertheless, this approach is not viable to tackle large problems due to the quadratic complexity with respect to the number of spatial degrees of freedom.

The BB-preconditioner

We present now the third preconditioner considered in this paper. Its design started from the idea of using the following block triangular preconditioner (see for example [START_REF] Benzi | Numerical solution of saddle point problems[END_REF]),

P t = A B T -Ŝ ,
where Ŝ should ideally equal the dual Schur complement S = -(C + BA -1 B T ) or some approximation of it (we again recall that, formally, we cannot write A -1 since the matrix A is singular). Given a vector (c; d) ∈ R n+m , one application of this preconditioner requires to compute the vector (x; y) ∈ R n+m solving

-Ŝy = d, (17)
Ax = b = c -B T y. ( 18
)
This approach is particularly attractive from the computational point of view, since matrix A in eq. ( 28) is block diagonal and each block A k is a weighted Laplacian. Thus, we can efficiently get the solution x by solving separately K + 1 linear systems in the form A k x k = b k using multigrid solvers. However, two issues arise. 1 108 110 106 107 3.8e-1 7.3e-1 1.6e0 4.1e0 1.0 1.0 1.8 2.1 2e-1 72 72 73 74 3.1e-1 6.3e-1 1.4e0 3.4e0 1.1 1.6 2.4 2.9 4e-2 65 65 66 66 3.1e-1 6.4e-1 1.5e0 3.5e0 1.7 2.7 3.5 4.2 8e-3 56 58 60 60 2.9e-1 6.4e-1 1.5e0 3.8e0 2.2 3.1 4.0 4.8 2e-3 52 54 48 52 2.9e-1 6.4e-1 1.4e0 3.7e0 2.3 3.1 4.2 5.2 3e-4 52 43 47 48 3.1e-1 5.8e-1 1.4e0 3.8e0 2.5 3.4 4.4 6.5 6e-5 42 41 41 39 3.0e-1 5.8e-1 1.3e0 3.1e0 2.9 4.4 5.6 7.3 1e-5 41 42 40 43 2.8e-1 6.3e-1 1.3e0 3.8e0 3.8 5.5 6.6 8.5 3e-6 43 55 50 68 3.2e-1 2.0e0 1.8e0 6.5e0 3.8 19.5 7.9 9.7 5e-7 54 72 67 74 3.8e-1 2.0e0 2.7e0 9.7e0 3.8 14.0 8.

.0 2e-1 1124 † † † 2.7e2 † † † 1.0 † † † 4e-2 † † 540 503 † † 6.9e2 2.0e3 † † 1.
† † † † † † † † † † †
Table 2: Numerical results using the SIMPLE preconditioner. Each sub-table reports the results for a given mesh, from the coarsest (top) to the finest (bottom), and for different time discretizations (∆t = 1/(K + 1)). Each subtable reports, while µ is reduced (leftmost column), the averaged outer iteration and CPU time per linear system (Outer/Lin.sys. and CPU/Lin.sys.) and the averaged inner iterations per outer iteration Inner/Outer for solving the linear system in eq. ( 16) (metrics defined in section 4). A final row summarizes the averages on the whole simulation. We highlighted in gray the time-space combination providing the uniform discretization. The † symbol denotes those IP steps where the linear solver failed.

The first one is due to the fact that each block A k is singular, thus we need to ensure that

b k ∈ Im(A k ) = Ker(A k ) ⊥ = 1 ∈ R N T ⊥ , k = 1, . . . K + 1 , (19)
for the linear systems to be admissible. This condition may not hold for the right-hand side b = c -B T y in eq. ( 18) for any application of the preconditioner, because Im(B T ) is not orthogonal to Ker(A) (see remark 2). Even if the orthogonality condition in eq. ( 19) is satisfied, or imposed by orthogonalization, there are multiple solutions, which differ by an additive constant, and thus we also need a criteria to select one. This is not trivial since the solution δϕ of eq. ( 12) must also solve the second block of equations, Bδϕ -Cδρ = g. In the next section, we describe how to deal with this compatibility issue with a reformulation of the continuous nonlinear system in eq. ( 2).

A second, harder, issue is the design of a good approximation of the inverse of the matrix S = -(C + BA -1 B T ), used in the solution of the linear system eq. ( 17). First, S is not well defined since A is singular and Im(B T ) is not orthogonal to Ker(A). Moreover, even overcoming this problem, the matrix S is block tridiagonal but with dense blocks, due to the presence of the pseudo-inverse of the Laplacian matrices A k , thus too costly to form and invert. This second issue is addressed in section 4.4.3.

Splitting formulation

In order to cope with the first issue, we split the solution ϕ of eq. ( 2) into two components, a new potential φ : Ω × [0, 1] → R satisfying the constraint Ω φ(t) dx = 0 for all t ∈ [0, 1] and a correction term:

ϕ(t, x) = φ(t, x) + t 0 λ(s)ds (20)
where λ : [0, 1] → R depends only on the time variable. In the new unknowns ( φ, ρ, s, λ) the system of PDEs (2) becomes

-∂ t ρ -div x ρ ∇ x φ = 0, ∂ t φ + | ∇ x φ| 2 2 + s + λ = 0, ρ ≥ 0, s ≥ 0, ρs = 0, Ω φ dx = 0. (21)
Following the same steps described in section 2, the (relaxed) discrete counterpart of equation in eq. ( 21) is a nonlinear system of equations with unknowns ( φ, ρ, s, λ)

∈ (R n × R m × R m × R K ) given by F ϕ ( φ, ρ) = F 1 ϕ ; . . . ; F K+1 ϕ =0 ∈ R n , Fρ ( φ, ρ, s, λ) = F 1 ρ ; . . . ; F K ρ =0 ∈ R m , F s (ρ, s) = F 1 s ; . . . ; F K s -µ1=0 ∈ R m , F λ ( φ) = F 1 λ ; . . . ; F K λ =0 ∈ R K+1 ,
where F k ϕ and F k s are defined in eqs. ( 5) and ( 7), while F k ρ and F k λ are given by

F k ρ ( φ, ρ, s, λ) := F k ρ (ρ, φ, s) + |c ′ |λ k , k = 1, . . . , K, F k λ ( φk ) := |c| T φk , k = 1, . . . , K + 1,
with F k ρ as given in eq. ( 6). Note that there are n + 2m + K + 1 equations and n + 2m + K unknowns. The additional equation fixes the global constant for the potential φ. The linear system arising from the Newton method becomes

    A B T B M ′ M ′ E ′ T Diag (s) Diag (ρ) EM        δ φ δρ δs δλ    =    f ḡ h i    = -    F ϕ F ρ F s -µ1 F λ    , (22)
where matrices E ∈ R K+1,(K+1)N T and E ′ ∈ R K,KN T ′ are given by

E =    1 T N T . . . 1 T N T    , E ′ =    1 T N T ′ . . . 1 T N T ′    .
Once the vectors δ φ and δλ in eq. ( 22) are found, δϕ can be retrieved in the original system eq. ( 12) using the formula

δϕ k = δ φk + k i=1 ∆t i δλ i ,
which is the finite-dimensional counterpart of eq. ( 20).

New saddle point linear system

After eliminating the unknown δs as in section 3.1, the linear system in eq. ( 22) becomes

  A B T B -C M ′ E ′ T EM     δ φ δρ δλ   = f g i . ( 23
)
The variable δλ can be eliminated as well by multiplying by E ′ the second block of equations in eq. ( 22), yielding

E ′ (Bδ φ -Cδρ) + E ′ M ′ E ′ T |Ω|I K δλ = E ′ g .
Plugging this expression in eq. ( 23), the linear system in the unknowns (δ φ; δρ) becomes

  A B T PB -PC EM   δ φ δρ = f P g i . ( 24
)
where the matrix P, which is given by

P := I - 1 |Ω| M ′ E ′ T E ′ ,
is a projector (indeed P 2 = P). Moreover, its transpose P T is the discretization of the zero mean-projector, which maps a function g : [0, 1] × Ω → R to g -Ω g/|Ω|.

The linear system in [START_REF] Eymard | Finite volume methods[END_REF] is not over-determined since Ker(A) ⊂ Ker(PB) and the rows of EM are K + 1 vectors linearly independent from the rows of A or PB. Recalling remark 3, δρ ∈ Ker(P) or equivalently δρ ∈ Im(P T ) and, since P T is also a projector, we can write δρ = P T δρ. Thus, we can rewrite system [START_REF] Eymard | Finite volume methods[END_REF] as

A B T P T PB -PCP T δ φ δρ = f P g , (25)
where we neglect the constraints EMδ φ = i. For any solution (δ φ; δρ), the solution (δ φ; δρ) of eq. ( 24) can be recovered by simply correcting the mean of δ φ.

The linear system in eq. ( 25) does not suffer from the first issue mentioned above associated to the singularity of matrix A. In fact, note that Im(B T P T ) ⊥ Ker(A), since Ker(A) ⊂ Ker(PB). This ensures that the dual Schur complement S = -P(C + BA + B T )P T (where A + denotes the Moore-Penrose pseudo-inverse of A) is well defined. Moreover the following (ideal) block triangular preconditioner

P t = A B T P T -S (26)
is well defined as well. In fact, given a vector (c; d) with c ∈ Ker(A) ⊥ and d ∈ Im(P), the application of this preconditioner requires to compute the vector (x; y) solving the following two linear systems

-Sy = d , (27)
Ax = b = c -B T P T y . (28)
The linear system in eq. ( 28) is now well defined, since b ∈ Ker(A) ⊥ at any preconditioner application. In fact, both vectors B T P T y and c belong to Ker(A) ⊥ , the first vector because Im(B T P T ) ⊂ Ker(A) ⊥ , the second vector since it is a linear combination of the vector f in eq. ( 25) and vectors in Im(A) ∪ Im(B T P T ), all belonging to Ker(A) ⊥ .

Approximating the inverse of the dual Schur complement

We turn now to the issue of approximating the Schur complement. Our first attempt was to approximate it with P(C + B(Diag (A)) -1 B T )P T , but it led to poor performance. We tried to apply the Algebraically stabilized Least Square Commutator proposed in [START_REF] Elman | Least squares preconditioners for stabilized discretizations of the Navier-Stokes equations[END_REF]Section 4.2], but this leads to poor results, for different combinations of the relaxation parameters that appear in such approach. We attribute this phenomenon to the fact that the matrix C is not a stabilization operator for the saddle point system, as assumed in [START_REF] Elman | Least squares preconditioners for stabilized discretizations of the Navier-Stokes equations[END_REF].

In order to devise a better approximation of the Schur complement, here we follow the idea in [START_REF] Elman | Block preconditioners based on approximate commutators[END_REF][START_REF] Elman | Least squares preconditioners for stabilized discretizations of the Navier-Stokes equations[END_REF][START_REF] Elman | Preconditioning strategies for models of incompressible flow[END_REF] of looking at the infinite-dimensional differential operators associated to the matrices that compose it and try to deduce a possible approximation of its inverse. The following proposition shows a differential identity that goes in this direction. Proposition 1. Let ρ : [0, 1] × Ω → R >0 , ϕ : [0, 1] × Ω → R, smooth enough. Denoting by -∆ ρ =div(ρ ∇), the following operator identity holds:

(∂ t + div(• ∇ ϕ)) (-∆ ρ ) = -∆ ρ (∂ t + ∇ ϕ • ∇) -div ((∂ t ρ + div(ρ ∇ ϕ)) ∇) +2 div(ρ ∇ 2 ϕ ∇) . (29)
Proof. Consider any function g 0 : [0, 1] × Ω → R, such that Ω g 0 (t, •) = 0 for all t ∈ [0, 1] and let us define

u : [0, 1] × Ω → R given by -div(ρ ∇ u) = g 0 (30)
with zero Neumann boundary conditions. Multiplying both sides of equation eq. ( 30) by ∂ i ϕ (where ∂ i = ∂ xi ) and taking the divergence, we obtain i,j

∂ i (∂ i ϕ∂ j (ρ∂ j u)) = div(g 0 ∇ ϕ).
Using the identity

∂ i ϕ∂ j (ρ∂ j u) = ∂ j (∂ i ϕρ∂ j u) -∂ j (∂ i ϕ)ρ∂ j u, we get div(g 0 ∇ ϕ) = - i,j ∂ i (∂ j (∂ i ϕρ∂ j u) -∂ j,i ϕρ∂ j u) = - j,i ∂ j (∂ i,j uρ∂ i ϕ) - j,i ∂ j (∂ i (ρ∂ i ϕ)∂ j u) + div(ρ ∇ 2 ϕ ∇ u) = - j ∂ j i ∂ i,j uρ∂ i ϕ - j ∂ j i ∂ i (ρ∂ i ϕ)∂ j u + div(ρ ∇ 2 ϕ ∇ u) = -div(ρ ∇ 2 u ∇ ϕ) -div(div(ρ ∇ ϕ) ∇ u) + div(ρ ∇ 2 ϕ ∇ u) . Using the identity ∇(∇ ϕ • ∇ u) = ∇ 2 ϕ ∇ u + ∇ 2 u ∇ ϕ, we obtain div(g 0 ∇ ϕ) = -div(ρ ∇(∇ ϕ • ∇ u)) -div(div(ρ ∇ ϕ) ∇ u) + 2 div(ρ ∇ 2 ϕ ∇ u) . (31)
Now, differentiating with respect to time the expression in eq. ( 30), we obtain

-div (∂ t ρ ∇ u + ρ ∇ ∂ t u) = ∂ t g 0 . (32)
Combining eq. ( 31), eq. ( 32) we obtain

∂ t g 0 + div(g 0 ∇ ϕ) = -div(ρ ∇(∂ t u + ∇ ϕ • ∇ u)) -div ((∂ t ρ + div(ρ ∇ ϕ)) ∇ u) +2 div(ρ ∇ 2 ϕ ∇ u) .
Since g 0 = -∆ ρ u by eq. ( 30), we proved eq. ( 29).

Our approximation of the Schur complement is based on neglecting the termsdiv ((∂ t ρ + div(ρ ∇ ϕ)) ∇ and 2 div(ρ ∇ 2 ϕ ∇) on the right-hand side of eq. ( 29). The first term may be assumed to be small since it contains the continuity equation, which is the first nonlinear equation in the system eq. (2a). This term is small at the initial and final Newton steps within each IP iteration, and experimental observations showed it remains small during the intermediate Newton iterations. The second term can be neglected under certain assumptions. If the optimal transport is a translation this term is null. When the transported measures ρ in and ρ f are Gaussian [START_REF] Valdimarsson | On the hessian of the optimal transport potential[END_REF], or more generally logconcave probability densities [START_REF] Caffarelli | Monotonicity properties of optimal transportation and the FKG and related inequalities[END_REF][START_REF] Colombo | Lipschitz changes of variables between perturbations of logconcave measures[END_REF], the largest eigenvalue of ∇ 2 ϕ(0, x) is bounded. However, note that in general

Numerical results

In this section we present the results obtained on test case 2 using the preconditioner described in section 4.4.3. Its application requires the approximate solution of the linear systems in eqs. ( 28) and [START_REF] Métivier | A review of the use of optimal transport distances for high resolution seismic imaging based on the full waveform[END_REF]. The first is solved using the AGMG solver with inner tolerance ε S in = 1e -1 (determined experimentally). The second involves the solution of K + 1 weighted Laplacian systems. We solved them with the AGMG solver with tolerance ε A in = 5e -2. This value, slightly smaller than those adopted in previous cases, improves the robustness of the preconditioner. The number of multigrid iterations for these systems ranges between 2 and 3 in all cases.

In table 3 we summarized the results for the metrics defined in section 4 that we used to measure the preconditioner performance. The number of averaged inner iterations Inner/Outer refers to the solution of the linear system eq. ( 36). This number remains bounded between 2.1 and 8.5 iterations per preconditioner applications, it tends to increase slightly with µ → 0 but it is not affected by the size of the problem.

Unfortunately, the number of outer iterations Outer/Lin.sys. tends to increase as µ → 0, exceeding in some cases the 400 limit we fixed for the FGMRES iterations. This phenomenon is more pronounced when the temporal scale is finer than the spatial one to. This loss in efficiency for µ → 0 is more evident in the test cases 2 and 3 whose solutions have compact support, while for the test case 1 it is less pronounced. This suggests that this phenomenon can be related to the degeneracy of the linear system eq. ( 36), which becomes underdetermined due to the presence of terms that go to zero on both sides of the equation.

Nevertheless, for µ ≈ 1e -5, this preconditioner is the only one scaling well with respect to the temporal and spatial discretization. The number of averaged outer iterations Outer/Lin.sys. increases only slightly when the mesh is refined or the time steps is halved. This results in a CPU time that scales sightly worse than linearly with respect to the number of degrees of freedom used to discretize the problem.

Summary of numerical results

In this section, we summarize the numerical results to compare the pros and cons of the preconditioning approaches proposed in this paper. In fig. 1 we compare the total CPU time (y-axis) required to achieve IP relaxation µ ≈ 1e -5 (eight IP iterations) with respect to the total number of degrees of freedom used, while halving the time step and the mesh size used. This summarizes the data in tables 1 to 3 with gray background color. We included in this comparison the results obtain for the test cases 1, 3 and 4.

The behavior of the preconditioner based on the primal Schur complement described in section 4.2 is determined by the inner solver used to solve the block linear system 14, which is the finite-dimensional counterpart of a time-space weighted and anisotropic Laplacian. This linear system was solved efficiently with the AGMG solver only for the smallest test cases during the initial IP steps, when the relaxation µ is relatively high. In fact, it did not reach the threshold µ ≈ 1e -5 in the finest discretizations for all test cases. The SIMPLE preconditioner is rather robust with respect to different time steps and relaxation parameters µ. However, it becomes inefficient when large grids are used. In fact, the number of outer iterations and the CPU time approximately scale linearly with respect to the number of time steps and the number of interior point steps, but quadratically with respect to the number of cells. Moreover, it had some failures in the first IP steps for all test cases in the finest discretization. However, the SIMPLE preconditioner is the only one that becomes faster in the last IP iterations. Thus, it may be the only one providing a viable option when accurate solutions of the optimal transport problem are required, possibly combined with more efficient approaches for the initial IP steps.

The BB-preconditioner described in section 4.4 is the only one able to tackle efficiently large scale problems, since the number of inner iterations required per each Krylov step remains rather constant, while the number of outer iterations increases only slightly using smaller time steps and finer grids. This holds as long as the IP tolerance is approximately 1e -5. Using smaller values, the number of inner and outer iterations tends to increase, leading in some cases to the solver's failure. The reason causing this performance degradation is not clear and needs further investigations.

Conclusions

We presented different preconditioners for solving via iterative methods the saddle point linear systems arising when computing solutions of the Benamou-Brenier formulation using an IP method. The first preconditioner presented, based on the primal Schur complement, is neither efficient nor robust since its application requires the solution of a complicated linear system involving an anisotropic weighted Laplacian in the time-space domain. The second preconditioner, called SIMPLE, was shown to be remarkably robust but with a CPU time that scales quadratically with respect to the number of spatial degrees of freedom. The most efficient approach turned out to be a block triangular 3) and BB (the BBpreconditioner). The results refer to the test cases 1,2,3 and 4 (from top to bottom panel) and using different time-space discretization (left to right). The columns height is normalized by the maximum among the three preconditioners. The gray portion of the column denotes the part of the preprocessing time of each preconditioner. If µ ≈ 1e -5 was not achieved, we report the value µ reached. For the SIMPLE preconditioner some intermediate IP steps failed. We report the failures number and we do not sum any CPU time to the total cost, since the BB-preconditioner appears to be the best among the three presented anyway.

preconditioner where the inverse of the dual Schur complement is approximated exploiting a partial commutation of its components, which we named BB-preconditioner. Although this approach loses efficiency in the latest IP steps, it is the only one scaling well with respect to the time and space discretization size, both in 2d and 3d problems. Up to a reasonable value µ ≈ 1e -5 for the IP relaxation parameter, our numerical experiments showed that the problem can be solved with a good scaling of the CPU time with respect to the number of degrees of freedom. Combining it with further tuning strategies of the IP method (e.g. predictor-corrector methods, multilevel approach, parallelization) we believe it has the potential to provide highly efficient solvers for the dynamical optimal transport problem. Moreover, we remark that the same strategy can be applied to variations of the problem which consider further penalization of the density in eq. ( 1) (such as those considered in [START_REF] Porretta | Regularizing effects of the entropy functional in optimal transport and planning problems[END_REF]) without major modifications. The deterioration of the preconditioner for smaller values of µ requires further investigation and it is left for future works.

7 T 1 ( 7 T 2 ( 1 T 3 (

 717213 NT = 896, NT ′ = 2688) #DOF=(5.6e4, 1.1e5, 2.3e5, 4.6e5) 1 250 230 229 230 3.1e0 5.9e0 1.2e1 2.8e1 1.0 1.0 1.1 1.9 2e-1 136 144 147 148 1.9e04.3e0 9.8e0 2.3e1 1.0 1.2 1.6 2.4 4e-2 114 116 119 119 1.7e0 4.0e0 9.7e0 2.3e1 1.4 1.8 3.0 3.9 8e-3 90 107 107 106 1.5e0 4.0e0 9.9e0 2.4e1 1.8 2.2 3.8 4.7 2e-3 89 91 89 91 1.6e0 3.8e0 9.5e0 2.4e1 2.0 2.9 4.3 5.3 3e-4 102 85 84 79 1.9e0 3.8e0 9.7e0 2.3e1 2.0 3.0 4.5 5.4 6e-5 105 97 82 74 2.0e0 4.7e0 1.0e1 2.3e1 2.0 3.1 4.8 5.9 1e-5 100 86 79 78 2.1e0 4.2e0 1.1e1 2.6e1 2.2 3.2 6.1 7.1 3e-6 91 74 78 85 2.1e0 4.2e0 1.2e1 3.4e1 2.7 4.7 7.9 10.2 5e-7 88 74 96 115 2.1e0 4.7e0 1.6e1 5.2e1 3.1 5.3 8.3 12.3 120 118 118 119 2.0e0 4.4e0 1.1e1 2.7e1 1.7 2.2 3.4 4.NT = 3584, NT ′ = 10752) #DOF=(2.3e5, 4.6e5, 9.1e5, 1.8e6) 1 516 459 497 475 2.6e1 4.8e1 1.2e2 2.7e2 1.0 1.0 1.0 1.0 2e-1 300 298 319 317 1.7e1 3.5e1 8.4e1 2.1e2 1.0 1.0 1.1 1.6 4e-2 206 222 234 221 1.3e1 2.9e1 7.2e1 1.8e2 1.5 1.4 1.9 3.0 8e-3 190 164 167 169 1.3e1 2.4e1 5.8e1 1.5e2 1.8 1.8 2.4 4.0 2e-3 185 161 150 154 1.4e1 2.5e1 5.6e1 1.6e2 2.6 2.0 3.2 4.8 3e-4 223 170 168 139 1.8e1 2.8e1 6.9e1 1.6e2 3.0 2.5 3.5 5.2 6e-5 228 222 157 141 2.0e1 3.8e1 7.1e1 1.8e2 3.0 2.6 3.7 5.5 1e-5 255 183 156 145 2.3e1 3.4e1 7.4e1 1.9e2 3.0 2.7 3.9 5.7 3e-6 250 189 162 136 2.4e1 3.8e1 8.1e1 1.9e2 3.1 2.9 4.1 6.2 5e-7 252 160 140 140 2.5e1 3.3e1 7.4e1 2.2e2 3.1 3.0 4.2 6.9 272 236 237 224 1.9e1 3.4e1 7.8e1 2.0e2 1.9 1.7 2.1 3.NT = 14336, NT ′ = 43008) #DOF=(9.0e5, 1.8e6, 3.7e6, 7.3e6) 1 1003 919 959 1045 2.2e2 4.9e2 1.1e3 3.8e3 1.0 1.0 1.0 1

Figure 1 :

 1 Figure 1: Comparison of CPU time spent in solving linear algebra problems to reach an IP relaxation µ ≈ 1e -5 using the preconditioners approaches presented in this paper. They are denoted by Prim. (the preconditioner based on the primal Schur complement in section 4.2, Simp. (the SIMPLE preconditioner in section 4.3) and BB (the BBpreconditioner). The results refer to the test cases 1,2,3 and 4 (from top to bottom panel) and using different time-space discretization (left to right). The columns height is normalized by the maximum among the three preconditioners. The gray portion of the column denotes the part of the preprocessing time of each preconditioner. If µ ≈ 1e -5 was not achieved, we report the value µ reached. For the SIMPLE preconditioner some intermediate IP steps failed. We report the failures number and we do not sum any CPU time to the total cost, since the BB-preconditioner appears to be the best among the three presented anyway.

Table 1 :

 1 Numerical results using the preconditioner based on the primal Schur complement. Each sub-table reports

	K +1 16 32 64 128	16	32	64 128	16	32	64 128
	µ Outer/Lin.sys.			CPU/Lin.sys.			Inner/Outer
		T 0 (NT = 224, NT ′ = 672) #DOF=(1.4e4, 2.8e4, 5.7e4, 1.1e5)
	1	4 4	2 1.8e-1 3.1e-1 9.7e-1 2.3e1	2.2	2.1	2.7	1.5
	2e-1	4 3	1 2.1e-1 6.1e-1 5.8e0 1.9e1	2.7	2.3	2.6	1.0
	4e-2	5 4	1 3.3e-1 1.1e0 7.4e0 2.0e1	1.9	2.3	1.9	1.0
	8e-3	4 4	1 4.8e-1 1.7e0 9.0e0 1.8e1	1.6	2.0	1.8	1.0
	2e-3	5 4	2 5.9e-1 1.9e0 8.8e0 2.0e1	1.6	1.8	1.8	1.1
	3e-4	4 4	2 6.4e-1 2.1e0 8.4e0 1.9e1	1.7	2.0	1.7	1.2
	6e-5	4 4	2 6.8e-1 2.1e0 8.8e0 2.2e1	1.7	1.8	2.0	1.2
	1e-5	4 4	2 7.0e-1 2.2e0 9.1e0 2.4e1	1.7	2.1	1.9	1.2
	3e-6	4 4	2 7.1e-1 2.4e0 7.9e0 2.5e1	1.8	1.8	1.7	1.0
	5e-7	4 4	2 7.8e-1 2.5e0 9.2e0 2.3e1	1.7	1.9	1.9	1.1
		4.5 4.2 4.0 1.7 4.9e-1 1.6e0 7.4e0 2.1e1	1.9	2.0	2.0	1.1
		T (NT = 896, NT ′ = 2688) #DOF=(5.6e4, 1.1e5, 2.3e5, 4.6e5)
	1	4 4	4 5.7e-1 1.4e0 2.6e0 1.4e1	2.7	2.3	2.1	3.5
	2e-1	3 3	3 6.5e-1 1.7e0 1.2e1 5.8e2	2.7	2.7	2.9	2.9
	4e-2	5 5	4 1.1e0 3.8e0 5.1e1 9.0e2	2.6	2.1	2.3	2.2
	8e-3	6 5	4 1.9e0 8.6e0 1.0e2 1.1e3	3.3	2.1	1.9	1.7
	2e-3	6 5	4 3.4e0 1.3e1 1.4e2 1.0e3	6.0	2.0	1.8	1.6
	3e-4	6 5	4 6.3e0 1.5e1 1.5e2 1.0e3 14.6	2.1	2.0	1.9
	6e-5	5 4	4 1.1e1 1.8e1 1.6e2 1.0e3 32.2	2.5	2.1	1.9
	1e-5	5 4	4 3.0e1 1.9e1 1.7e2 1.1e3 95.0	2.5	2.1	1.8
	3e-6	5 4	4 4.2e1 2.2e1 1.7e2 1.0e3 122.8 2.3	1.9	1.9
	5e-7	5 4	4 7.9e1 2.7e1 1.9e2 1.1e3 151.4 2.5	2.1	2.7
		5.5 4.8 4.3 3.7 1.5e1 1.2e1 1.0e2 8.6e2 49.0	2.3	2.1	2.2
		T (NT = 3584, NT ′ = 10752) #DOF=(2.3e5, 4.6e5, 9.1e5, 1.8e6)
	1	3 4	4 2.3e0 5.0e0 1.0e1 2.2e1	3.2	2.9	2.4	2.3
	2e-1	3 3	3 2.9e0 6.3e0 1.5e1 2.1e3	3.4	3.2	3.1	3.3
	4e-2	4 5	9 4.9e0 2.6e1 2.1e2 1.3e3	4.0	4.2	2.5 128.5
	8e-3	5 6 20 1.1e1 9.6e1 1.2e3 2.0e3	4.2	5.0	2.4 160.1
	2e-3	6 6 55 1.9e1 1.3e2 1.5e3 6.2e3	3.1	7.7	2.0 198.8
	3e-4	8 6	† 2.5e1 2.0e2 2.3e3	†	3.3 18.1	2.2	†
	6e-5	6 4	† 3.1e1 3.8e2 2.2e3	†	6.4 65.2	2.6	†
	1e-5	7 4	† 4.5e1 7.1e2 1.7e3	† 13.1 133.8 2.7	†
	3e-6	9 4	† 7.3e1 1.1e3 1.7e3	† 28.9 161.5 2.9	†
	5e-7	16 4	† 1.3e2 2.1e3 1.8e3	† 70.9 179.3 2.6	†
		4.7 6.3 4.6	† 2.7e1 4.1e2 1.1e3	† 12.2 79.4	2.5	†
	T 3 (NT = 14336, NT ′ = 43008) #DOF=(9.0e5, 1.8e6, 3.7e6, 7.3e6)
	1	3 3	3 9.3e0 2.0e1 5.0e1 9.8e1	3.5	3.4	2.9	2.8
	2e-1	3 3	3 1.1e1 2.5e1 6.0e1 4.7e2	4.0	3.9	3.6	3.3
	4e-2	4 4	7 2.7e1 1.5e2 3.4e3 1.8e3	5.6	4.6 22.9 84.2
	8e-3	5 9 12 7.4e1 5.0e2 2.6e3 3.7e3	5.7	5.2 154.9 105.4
	2e-3	7 18 26 1.4e2 2.4e3 4.9e3 1.1e4	4.5	3.3 182.5 152.1
	3e-4	6 58 59 2.1e2 2.7e3 1.5e4 3.6e4	4.0	3.2 199.0 198.9
	6e-5	6 † 162 2.6e2 2.1e3	† 9.5e4	4.3	4.1	† 200.0
	1e-5	5 †	† 3.2e2 2.4e3	†	†	4.5	7.1	†	†
	3e-6	5 †	† 3.5e2 2.7e3	†	†	4.4 14.9	†	†
	5e-7	5 †	† 3.6e2 3.7e3	†	†	4.6 41.8	†	†
		4.3 4.7 †	† 1.5e2 1.4e3	†	†	4.6	7.9	†	†

the results for a given mesh, from the coarsest (top) to the finest (bottom), and for different time discretizations (∆t = 1/(K + 1)

Acknowledgments

The bulk of this work was completed while E.F was affiliated with Univ. Lille, Inria, CNRS, UMR 8524 -Laboratoire Paul Painlevé, as member of the RAPSODI team. E.F. was partially funded by the Marie Skłodowska-Curie Action (HE MSCA PF 101103631). A.N. was partially funded by the Labex CEMPI (ANR-11-LABX-0007-01). The work of M. B. was supported in part by PNRR MUR Project PE0000013-FAIR.

this term may be non null, as for example in the test case 3. Supposing that both terms can be neglected, we get the following approximate identity:

In our finite-dimensional problem, this expression translates into the approximate matrix identity

where à ∈ R m,m and B ∈ R n,m discretize the operators -∆ ρ and ∂ t + ∇ ϕ • ∇, respectively, similarly to the matrices A and B but with different dimensions. The scaling factors M ′ -1 and M in eq. ( 33) are introduced to match the scaling dimensions of the matrices A, B, and B T . The most natural definition of à and B is the following:

where the matrices J , H, and G are defined in eqs. [START_REF] Benzi | Multilevel algorithms for large-scale interior point methods[END_REF] and [START_REF]A preconditioning technique for a class of pde-constrained optimization problems[END_REF], while the matrices Dt ∈ R KN T ′ ,(K+1)N T ′ and

Using the approximate identity in eq. ( 33), the linear system in eq. ( 27) becomes

Since Ker(PM ′ ) = Im(E T ) are discrete functions constant in space, we have that Im(M ′ P T ) ⊂ Im( Ã) and the last expression is equal to

Hence, we only have to solve the linear system

which is compatible since d ∈ Im(P). Note that the dense matrix (C + BA + B T ) in eq. ( 34) is replaced by the operator (CM ′ -1 Ã -BM -1 B) Ã+ M ′ which is dense by the presence of the pseudo-inverse of Ã. But now this operator is composed by two factors. This means that we can approximate the action of their respective inverse by solving a sparse linear system followed by a matrix-vector product. In fact, since the solution δρ in eq. ( 25) belongs to the image of P T , we may look for a solution y ∈ Im(P T ). By setting

we can compute y by first solving the linear system

and then computing y = M ′ -1 Ãz.

Moreover, in order to cope with the ill-conditioning derived from the presence of the matrix C = M ′ Diag (s/ρ), we scale eq. ( 35) by Diag (ρ) and we solve the linear system

similarly to what done in eq. ( 16) for the SIMPLE preconditioner.

We refer to the resulting preconditioner as "BB-preconditioner", due to presence of this composed operator in the approximate factorization of the Schur complement. This choice was done in analogy with the "BF Bt-preconditioner" introduced in [START_REF] Elman | Preconditioning strategies for models of incompressible flow[END_REF], whose ideas inspired the preconditioner described in this section.


T Le Bourgeois
email: thomas.le_bourgeois@cirad.fr

S Bazan

VÉGÉPHYL -25 ème CONFÉRENCE DU COLUMA JOURNÉES INTERNATIONALES SUR LA LUTTE CONTRE LES MAUVAISES HERBES ORLÉANS -5, 6 et 7 DÉCEMBRE 2023 LA FLORE ADVENTICE DES PRINCIPALES CULTURES EN GUYANE FRANCAISE

Keywords: flore adventice, Guyane Française, systèmes de culture, diversité, indicateur écologique weed flora, French Guiana, cropping systems, diversity, ecological indicator

A partir de 9 études préalables, nous avons obtenu une vision synoptique de la flore adventice des principaux systèmes de culture de Guyane Française (vivrier, maraîchage, arboriculture, pâturage, riz irrigué). La flore globale est riche de 334 espèces, mais au niveau de chaque culture, la diversité floristique est assez semblable avec près de 100 espèces. Une vingtaine d'entre elles apparaissent majeures dans différentes cultures, tandis que 30 semblent spécifiques à un contexte agricole particulier. Euphorbia heterophylla est majeure en vivrier, maraîchage et arboriculture, Cyperus aromaticus est commune au maraîchage, à l'arboriculture et aux pâturages, tandis que les Ludwigia spp. se retrouvent en maraîchage, riziculture et pâturages. Spermacoce verticillata et Mimosa pudica sont majeures en pâturages et en arboriculture, Asystasia gangetica, Digitaria horizontalis, Cyanthillium cinereum, Paspalum conjugatum sont majeures en maraîchage et arboriculture tandis que Fimbristylis littoralis, Eleocharis filiculmis, et Fuirena umbellata sont majeures en pâturage et riziculture. Certaines espèces sont des indicateurs notamment de l'humidité du sol.

Question et cadrage

En Guyane Française comme partout ailleurs dans le monde, la maîtrise des adventices représente une problématique majeure en agriculture quel que soit les systèmes de culture. Que la gestion des plantes indésirables soit manuelle, mécanique, chimique ou agronomique, elle constitue une charge très importante en temps de travail, équipements spécifiques et intrants, nouvelle organisation de l'exploitation et coût financier. Depuis la fin des années 1990, plusieurs études de caractérisation de la flore adventice ont été réalisées dans les principaux systèmes de culture du département. A partir de ces travaux, l'objectif de cette communication est de présenter une synthèse des connaissances sur la flore adventice des cultures de Guyane : flore globale, diversité par cultures, espèces généralistes ou au contraire spécifiques de chaque culture, en fonction des contextes écologiques et des pratiques de gestion qui leur sont propres. Une bonne connaissance de ces plantes est primordiale pour adapter efficacement les méthodes de gestion des adventices au sein des systèmes de culture.

Approche méthodologique

La Guyane Française : contexte et production agricole La Guyane se situe au nord-est du continent sud-américain. Elle s'intègre dans le plateau des Guyanes qui s'étend du sud du Venezuela au nord-est du Brésil. La Guyane occupe seulement 0,5 % de la surface du continent alors qu'avec une superficie de 83 846 km² elle est de loin le plus vaste des départements français (16 % du territoire de la métropole). Le relief y est peu marqué, avec une zone de collines littorales et quelques points dépassant 600 m d'altitude dans le sud, le plus haut culminant à 830 m. La Guyane est recouverte à 96 % par la forêt amazonienne qui est sillonnée de rivières et de fleuves. Les 378 km de littoral, très vaseux (sédiments charriés par les fleuves), sont en grande partie couverts de mangroves. La bande côtière, qui a subi une forte déforestation, se présente sous forme de savanes.

La Guyane bénéficie d'un climat de type équatorial humide. Sa position proche de l'équateur ainsi que sa façade océanique lui confèrent une grande stabilité climatique, marquée notamment par la faiblesse des vents et la faible amplitude des températures. En revanche, les précipitations connaissent de fortes variations annuelles et déterminent le rythme des saisons. L'année est marquée par une saison humide (ou saison des pluies) qui s'étale de décembre à juillet, entrecoupée par une petite saison sèche aux alentours de mars (communément appelée « petit été de mars ») et une plus longue saison sèche d'août à novembre. La pluviométrie annuelle oscille de près de 2 000 mm dans le nordouest à plus de 4 000 mm dans la région située entre Roura et Régina. Elle est de 3 000 mm en moyenne sur la bande côtière de Kourou à Cayenne et un peu moindre dans les régions de l'intérieur (2 500 mm) [START_REF] Minagri | Ministère de l'Agriculture et de l'Alimentation[END_REF]. L'humidité relative moyenne est élevée, entre 80 % et 90 % selon la saison. Néanmoins, la Guyane dispose d'un ensoleillement important, avec en moyenne 2 200 heures d'insolation annuelle, les maxima étant situés sur la bande côtière (IEDOM, 2017).

Les terres agricoles sont détenues à 90% par l'État et allouées aux agriculteurs sous la forme de concessions. Le secteur primaire ne représente que 4% de la valeur ajoutée dans l'économie guyanaise en 2014 (IEDOM, 2017), et les emplois salariés en agriculture ne représentent que 1% des emplois du département [START_REF] Minagri | Ministère de l'Agriculture et de l'Alimentation[END_REF]. Les exploitations traditionnelles sur abattis, essentiellement localisées le long des fleuves Maroni à l'ouest et Oyapock à l'est, sont dédiées à une production vivrière, notamment de manioc. La production maraîchère et fruitière est localisée à Régina, Cacao, Matiti, Javouhey et Saint-Laurent-du-Maroni, tandis que la production rizicole (maintenant arrêtée) était uniquement localisée au polder de Mana. Enfin, les exploitations d'élevage bovine sont réparties le long de la côte de Régina à Mana et un peu dans l'ouest à Saint-Laurent-du-Maroni. En 2022, 6 140 exploitations agricoles étaient recensées en Guyane, toutes cultures confondues, la grande majorité d'entre elles (80%) ayant une SAU comprise entre 1 et 5 ha. La surface agricole totale de la Guyane consacrée à la production végétale représente 36 430 ha, dont 500 ha sont consacrés au maraichage, 8 190 ha à l'arboriculture fruitière, et près de 18 190 ha aux pâturages. Les cultures vivrières, notamment les tubercules comme le manioc représentent 4 070 ha. La riziculture représentait 5 500 ha mais elle est arrêtée depuis 2010. A ces cultures, s'ajoutent quelques 470 ha de canne à sucre, tandis que des cultures comme café, cacao et hévéa sont extrêmement limitées [START_REF] Minagri | Ministère de l'Agriculture et de l'Alimentation[END_REF].

Méthodologie

Nous avons exploité 9 études réalisées depuis la fin des années 1990 sur les adventices des 5 principales cultures de Guyane: le manioc (en tant que culture vivrière traditionnelle), le riz irrigué, le maraîchage, l'arboriculture fruitière, et les pâturages (tableau I). Ces sources d'informations comprennent à 1 enquête, 1 expérimentation, 1 guide de terrain, 1 site Web et 4 études par relevés phytoécologiques, dont 3 sont valorisables en méta-analyse (en grisé dans le tableau I) grâce à la disponibilité des données élémentaires. Enfin, 51 espèces majeures dans une culture sont mineures dans les autres cultures, ce qui représente la moitié de la liste des majeures de l'agriculture guyanaise dans notre synthèse.

Discussion -Conclusion

Il ressort de cette synthèse que la flore adventice globale des principales cultures de la Guyane française est très diversifiée avec 334 espèces inventoriées. Ce niveau de diversité est assez proche de celui de la flore adventice de La Réunion qui sans riziculture compte 273 espèces (Dodet et al., unpublished). Chaque culture en Guyane a une flore spontanée de 100 à 166 espèces. En comparaison, à La Réunion, on recense 220 espèces en canne à sucre (Le [START_REF] Bourgeois | Caractérisation des enherbements en culture de canne à sucre à la Réunion. 19ème Conférence Internationale du Columa[END_REF], Le Bourgeois et al., 2005) dont 165 considérées comme majeures [START_REF] Marnotte | La flore des mauvaises herbes de la Canne à Sucre à La Réunion -Caractérisation à partir des témoins des essais d'herbicides de[END_REF], 95 dans les pâturages (Le [START_REF] Bourgeois | Adventices des cultures des pâturages et des milieux naturels, situation et programmes de recherche à la Réunion[END_REF], 137 en ananas (Lebreton et Le Bourgeois, 2005) et 84 en lentille, culture très particulière en altitude sur une faible surface (Lebreton et Le Bourgeois, 2005). Par ailleurs, en Guadeloupe, la flore adventice de la canne à sucre comprend 107 espèces [START_REF] Hatil | Caractérisation de l'enherbement de la canne à sucre en Guadeloupe[END_REF][START_REF] Hatil | Characterization of sugarcane weed communities in Guadeloupe[END_REF], la flore des cultures maraîchères à l'île Maurice compte 131 espèces [START_REF] Minagri | Ministère de l'Agriculture et de l'Alimentation[END_REF], tandis qu'au Nord Cameroun, dans une grande région à fort gradient climatique et pédologique, la flore adventice du cotonnier compte 280 espèces (Le [START_REF] Bourgeois | Les mauvaises herbes dans la rotation cotonnière au Nord-Cameroun (Afrique). Amplitude d'habitat -Degré d'infestation[END_REF]. Ces quelques éléments de comparaison, montrent qu'à l'échelle d'un bassin de production ou d'une région tropicale au contexte écologique peu diversifié, comme peuvent l'être les DOM français, ou les îles tropicales comme Maurice, la flore adventice d'un système de culture donné correspond environ à une centaine d'espèces dont 20 à 40 espèces sont considérées comme des espèces d'importance majeure d'un point de vue agronomique. Il apparaît donc important de se focaliser sur ces quelques espèces majeures et de bien connaître leur comportement et leur comportement biologique pour élaborer une stratégie de gestion de l'enherbement la plus pertinente possible pour une culture donnée en fonction des conditions environnementales de la région concernée.

  Tableau I : Liste des études utilisées, cultures concernées, contenu et références (List of studies used, crop studied, content and reference) Liliopsida et Ptérytophytes) et de familles telles que les Poaceae et les Cyperaceae au sein des flores par cultures et de la flore globale. Dans un deuxième temps, nous avons identifié les adventices majeures de chaque culture ou système de culture (les plus fréquentes et/ou abondantes) puis nous avons consolidé et rapproché ces 5 listes d'espèces, majeures pour au moins une culture, et mineures ou absentes pour d'autre cultures. Des similitudes et des singularités ont été repérées, vraisemblablement caractéristiques de conditions agroécologiques partagées ou particulières ((valant indicateurs écologiques) Tous ces éléments de floristique ont été synthétisés dans un jeu de données publié sur le Dataverse Amatrop -Tropical weed studies du Cirad (Le Bourgeois et al., 2023) et sont libres d'accès. Diversité de la flore adventice des cultures en Guyane française La diversité de la flore adventice des principaux systèmes de culture de Guyane et la représentativité des grands groupes taxonomiques (Dicotylédones, Monocotylédones et Ptérytophytes) et de familles telles que les Poaceae et les Cyperaceae sont présentées dans le tableau II.La flore adventice de la Guyane ainsi consolidée, dite globale dans la suite de l'article, comprend 334 espèces, dont 1/3 de monocotylédones constituées essentiellement de 58 espèces de Poaceae et 42 espèces de Cyperaceae et 2/3 de dicotylédones. Les fougères ne représentent que 5 espèces.[START_REF] Hauswirth | L'enherbement du polder de Mana -Etat des lieux et perspectives[END_REF], la liste des plantes des rizières inclut, 45 adventices du riz dans les casiers, et 119 rudérales des canaux, des diguettes et des parcours entre les casiers (grande diversité floristique en lien avec un important gradient écologique, depuis les milieux toujours exondés jusqu'aux milieux inondés). A l'inverse, avec 108 espèces, la diversité de la flore adventice des pâturages apparaît assez faible, car 30 Poaceae et 12 Fabaceae spontanées ou implantées, y sont considérées comme des espèces fourragères et non des adventices. Cette situation illustre bien la versatilité du statut des plantes en agriculture (non désirables vs désirables, adventices vs cultivées) en fonction du contexte dans lequel elles se développent et des objectifs de production de l'agriculteur. , depuis des pentes toujours exondées jusqu'à des bas-fonds régulièrement inondés. Les fougères (Ptérydophytes) sont très peu nombreuses, que ce soit en conditions terrestres ou aquatiques. Euphorbia heterophylla, Digitaria horizontalis et Paspalum conjugatum sont communes au vivrier, au maraîchage et à l'arboriculture ; Mimosa pudica et Scleria gaertneri sont communes au vivrier, à l'arboriculture et aux pâturages. Commelina erecta est commune au vivrier et à l'arboriculture. Amaranthus dubius, Asystasia gangetica, Bidens pilosa, Cyanthillium cinereum, Cynodon dactylon, Digitaria ciliaris, Emilia sonchifolia et Panicum pilosum sont communes au maraîchage et à l'arboriculture. Cyperus aromaticus (espèce exotique introduite récemment) est rapidement devenue une adventice majeure en maraîchage, arboriculture et pâturages. Fimbristylis littoralis est commune au maraîchage, aux pâturages et au riz, notamment dans les situations de bas-fonds humides à inondés, de même que Ludwigia octovalvis en maraîchage et riz. Spermacoce verticillata est une adventice majeure en arboriculture enherbée et en pâturage, c'est une espèce typique des enherbements permanents où elle peut former de grosses populations et représente une très grosse contrainte. Fuirena umbellata, Cyperaceae des milieux humides à inondés est commune aux pâturages temporairement inondés et au riz. A l'opposé, 30 espèces majeures sont spécifiques d'une seule culture et sont donc absentes des autres cultures, ce qui en fait de bonnes espèces indicatrices. Cela représente 1 espèce en arboriculture (Rottboellia cochinchinensis), 3 espèces en maraîchage (Leptochloa mucronata, Teranaya aculeata et Vandellia diffusa), 6 en culture vivrière (Cyathula prostrata, Digitaria fuscescens, Ipomoea quamoclit et 3 taxons identifiés au genre Lindernia sp., Mikania sp. et Solanum sp.). Il y en a 8 en pâturage qui sont essentiellement des Cyperaceae et des Melastomataceae, et jusqu'à 12 espèces en riziculture, qui sont des espèces typiquement aquatiques, sélectionnées par les conditions d'inondation longue des parcelles rizicoles.Le cas de l'arboriculture fruitière à enherbement permanent est particulier car cet enherbement est une volonté de l'agriculteur pour maintenir la couverture du sol. Aussi, cette communauté végétale, la plupart du temps spontanée et généralement entretenue par fauchage, joue également un rôle de service.

			Vivrier Vivrier
	Espèce Espèce		(Manioc) Maraîchage Arboriculture Pâturages (Manioc) Maraîchage Arboriculture Pâturages	Riziculture Riziculture
	Eleusine indica Phyllanthus urinaria Principaux résultats et interprétation Emilia fosbergii Pleroma tibouchinum Espèces majeures des différentes cultures Emilia sonchifolia Pontederia crassipes	* * * *	M * M M	M M		* * M *	* * * * M
	Eragrostis ciliaris Portulaca oleracea Les espèces considérées comme majeures dans au moins une des 5 cultures principales de l'agriculture M * M guyanaise d'un point de vue agronomique, du fait de la combinaison de leur fréquence et de leur Erechtites hieracifolia M * * Pterolepis glomerata * * M	*
	abondance, sont présentées en fonction des différentes cultures dans le tableau III. Euphorbia heterophylla M M M Rhynchanthera grandiflora	M	*
	Euphorbia hirta Rhynchospora holoschoenoides		*	*	M *		M	* *
	Tableau II Tableau III Euphorbia hyssopifolia Diversité floristique de la flore adventice et représentativité de certains groupes taxonomiques à l'échelle de la Guyane et par culture (floristical diversity of weed flora and representativeness of some taxonomic groups at French Guiana level and by crop) Statut des adventices par culture : M = Majeure, * = mineure, cellule vide = absente * * M * Rhynchospora pubera * M (non observée) ; le fond coloré indique une espèce majeure dans la culture concernée Fimbristylis dichotoma * M Rhynchospora trispicata M et absente des autres cultures Fimbristylis littoralis M M Rolandra fruticosa M * (Status of weeds according to crops, M = Major, * = present, empty cell = absent (not recorded), color = Major weed for a crop but never recorded elsewhere) Fuirena umbellata M Rottboellia cochinchinensis M	* M M
	Grona heterocarpos subsp. ovalifolia Sagittaria guayanensis		*	M		*	M
	Vivrier Hymenachne amplexicaulis Salvinia auriculata	Vivrier	Arboricul-			Flore M M
	ETUDE Espèce Hyptis atrorubens Sauvagesia erecta	CULTURE (Manioc) Maraichage TYPE (Manioc) Maraîchage Arboriculture Pâturages REFERENCE ture Pâturage Riziculture M * M * M	totale Riziculture *
	Mission malherbologie. Projet recherches rurales Dicotylédones 79 Amaranthus dubius Hyptis lanceolata Scleria gaertneri L'enherbement du polder de 80% (Vivrier) 88 Manioc Mana Riz Monocotylédones 20 20% 30 Amaranthus viridis Imperata brasiliensis Scoparia dulcis Ptérydophytes 0 0% 3 Asystasia gangetica Ipomoea quamoclit Sieruela rutidosperma Nuisibilité de l'enherbement en riziculture Riz Plantes des rizières de Riz Poaceae 10 10% 14 Cyperaceae 6 6% 13 Axonopus compressus Ipomoea setifera Solanum sp. Nombre total 99 101 Bidens cynapiifolia Ischaemum rugosum Spermacoce ocymifolia Brachiaria sp. Lantana camara Spermacoce verticillata	* M * M * M * M M * M *	Enquête 87% 78 M 215 Relevés 30% 27 M 3% 1 M M phytoècologiques 74% 70 65% 100 61% 212 63% (Marnotte, 1999) M M * * M M * 25% 38 35% 61 37% 117 35% * M * (Hauswirth, 2002) 1% 0 0% 3 2% 5 1% M * Expérimentation (Marnotte et Téoulet, 2004) Site web (Marnotte et Carrara, 14% 18 17% 0 0% 30 18% 58 17% 13% 8 8% 35 32% 18 * * M * * 11% 42 13% 106 108 164 M M * M * M * 334 M * * * M M *
	Guyane Cecropia peltata Leptochloa mucronata Sphagneticola trilobata		M *	* M M	* *	2007)	* *	*
	Gestion des adventices des pâturages Ceratopteris thalictroides Leptochloa scabra Sphenoclea zeylanica Caractérisation des enherbements en Chromolaena odorata Lindernia sp. Tarenaya aculeata Cissus verticillata Ludwigia erecta Tonina fluviatilis Commelina diffusa Ludwigia hyssopifolia Torenia crustacea maraîchage et arboriculture Guide des adventices des cultures maraichères et Commelina erecta Ludwigia octovalvis Trema micrantha Cyanthillium cinereum Luziola subintegra Vandellia diffusa	Pâturages Maraîchage Fruitiers M M M Maraîchage Fruitiers M * M *	90 Relevés * * phytoécologiques 60 Relevés * M phytoécologiques * * * Guide de terrain * M M M	* M M M	(Huguenin, 2008) (Le Bourgeois, 2018) M M M M * (Berton, 2020) *	M M M * * * * * M * * M
	fruitières Cyathula prostrata Marsypianthes chamaedrys Xyris jupicai		M *	*	M		* M	*
	Flore des pâturages Cynodon dactylon Miconia dependens Xyris savanensis Cyperus aromaticus Miconia hirta Nombre d'espèces majeures par Pâturages		90 Relevés M phytoécologiques M	M M M M	(Favale, 2022) * M *	* M *
	Cyperus brevifolius Mikania sp. culture			M 23	26	35		M 26	18
	Cyperus haspan Cyperus luzulae Mimosa pudica Nombre d'espèces majeures Momordica charantia spécifiques d'une seule culture		* M M 6	* * 3	* M * 1		M M M * 8	* * * 12
	Les Dicotylédones représentent 80 % des espèces adventices des cultures vivrières et maraîchères vs 60 % en arboriculture, dans les pâturages et en riziculture. Les Monocotylédones sont beaucoup moins représentées avec 20 % des espèces sauf dans les pâturages et en rizicultures où elles représentent environ 30% de la flore. Dans ce groupe, les Poaceae, absentes en tant qu'adventices des Cyperus polystachyos M Cyperus rotundus M * Cyperus sphacelatus * * M Grona adscendens * M Murdannia nudiflora * M Nymphaea rudgeana Les 101 adventices majeures du fait de la combinaison de leur fréquence et de leur abondance, sont Oldenlandia corymbosa * M considérées comme problématiques pour les agriculteurs. Elles représentent 29% de la flore globale de la Guyane. Par culture (tab.III, avant-dernière ligne), on dénombre 18 espèces pour le riz, près de * Oryza sativa 25 espèces pour le vivrier, le maraîchage et les pâturages et 35 espèces en arboriculture. A noter que Digitaria ciliaris M M * Panicum maximum * * M * certaines espèces sont majeures dans plusieurs cultures.	* * * M * M *
	pâturages,sont peu nombreuses en culture vivrière et en maraîchage (10 à 14%), et plus diversifiées en arboriculture du fait de l'enherbement permanent (17 %), et en riziculture (18 %) ou la lutte contre Digitaria fuscescens M Digitaria horizontalis M M M * Panicum pilosum M M * Paspalum conjugatum M M M * 20 espèces sont majeures dans 2 cultures et 7 espèces dans 3 cultures. Ainsi,	*
	Echinochloa colonum Passiflora foetida			M	*	*		M *
	Echinochloa crus-pavonis Passovia pyrifolia			*	M		M
	Eleocharis filiculmis Phyllanthus amarus			*	M			*	M *

Pour établir la flore adventice globale de la Guyane en agrégeant les flores des différentes cultures , il a fallu au préalable homogénéiser la nomenclature botanique en l'actualisant à partir du référentiel taxonomique de GBIF

(GBIF, 2023)

. Une fois établie, nous avons comparé la diversité floristique par culture et comparé la représentativité des grands groupes taxonomiques (Dicotylédones -Magnolopsida, Monocotylédones -A l'échelle de chacune des cultures, la diversité de flore apparaît assez semblable avec environ une centaine d'espèces soit 1/3 de la flore adventice globale. Le nombre plus élevé d'espèces recensé en riziculture (165 espèces) vient du fait que, d'après (ce groupe comprenant notamment les riz sauvages est particulièrement difficile. Les Cyperaceae sont peu diversifiées (environ 6 à 13 % des espèces) dans les cultures toujours exondées (vivrier, maraîchage, arboriculture) ou dans les cultures toujours inondées comme le riz. En revanche, elles sont très nombreuses (32 %) dans les pâturages guyanais situés sur tout un gradient de situations d'humidité