Jean-Marie Mirebeau

Minimal Stencils for Discretizations of Anisotropic PDEs Preserving Causality or the Maximum Principle

We consider discretizations of anisotropic diffusion and of the anisotropic eikonal equation, on two dimensional cartesian grids, which preserve their structural properties: the maximum principle for diffusion, and causality for the eikonal equation. These two PDEs embed geometric information, in the form of a field of diffusion tensors and of a Riemannian metric respectively. Common knowledge is that, when these tensors are strongly anisotropic, monotonous or causal discretizations of these PDEs cannot be strictly local: numerical schemes need to involve interactions between each point and the elements of a stencil, which is not limited to its immediate neighbors on the discretization grid. Using tools from discrete geometry we identify the smallest valid stencils, in the sense of convex hull inclusion. We also estimate, for a fixed condition number but a random tensor orientation, the worst case and average case radius of these minimal stencils, which is relevant for numerical error analysis.

Introduction

The Partial Differential Equation (PDE) of diffusion obeys an important structural property, named the maximum principle. Numerical schemes which preserve this structure benefit from the discrete maximum principle [START_REF] Ciarlet | Discrete maximum principle for finite-difference operators[END_REF], a strong stability guarantee. The eikonal equation is the PDE formulation of an optimal control problem: to find the shortest exit path from a given domain. Discretization schemes which preserve its causal structure can be solved in a single pass using the Fast Marching algorithm [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF][START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations[END_REF], which has a quasi-linear complexity. Motivation for structure preservation in the discretization of PDEs is therefore plentiful, and stems from theoretical as much as practical considerations. In this intention, a variety of numerical schemes have been developed; for instance and without exhaustivity [MW53, BOZ04, Obe06, Wei98, FM13] for anisotropic diffusion, and [SV01, [START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Convex Hamilton-Jacobi Equations[END_REF][START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF] for anisotropic eikonal equations.

Non-isotropic PDEs have numerous applications, of which we can only give a glimpse. Anisotropic diffusion is required in porous media simulation [START_REF] Droniou | Finite volume schemes for diffusion equations: Introduction to and review of modern methods[END_REF], or stochastic control [START_REF] Bonnans | A fast algorithm for the two dimensional HJB equation of stochastic control[END_REF], but is also fundamental in image processing [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF]. The Anisotropic eikonal equation is naturally required for trajectory planning, but is also relevant for seismic imaging [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: theory and algorithms[END_REF], and extensively used in medical image segmentation [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF]. When discretizing PDEs, anisotropy usually comes with technical difficulties. Indeed, monotone or causal numerical schemes cannot be strictly local, but need to introduce interactions between each point and a stencil of potentially distant neighbors [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: theory and algorithms[END_REF][START_REF] Kocan | Approximation of viscosity solutions of elliptic partial differential equations on minimal grids[END_REF]. The objective of this paper is to provide a complete qualitative and quantitative understanding of this non-locality, in the limited setting of two dimensional discretizations on cartesian grids: we identify the smallest stencils, in the sense of convex hull inclusion, and we sharply estimate their worst case and average radius. Note that minimal stencils were similarly identified in [START_REF] Benamou | Monotone and Consistent discretization of the Monge-Ampere operator[END_REF] for the monotone discretization of the two dimensional Monge-Ampere equation.

We introduce in §1.1 the concept of D-diffusion stencil, where D is a positive definite diffusion tensor, which allows to discretize anisotropic diffusion PDEs while obeying the maximum principle. Likewise, the concept of M -eikonal stencil, where M is a positive definite matrix encoding the metric structure, allows to discretize anisotropic eikonal PDEs while preserving their causal structure.

We identify in §1.2 the smallest D-diffusion stencil, and the smallest M -eikonal stencil, in the sense of convex hull inclusion, see Theorems 1.12 and 1.13. These discretizations turn out to coincide with numerical schemes previously introduced by the author in [START_REF] Fehrenbach | Sparse Non-negative Stencils for Anisotropic Diffusion[END_REF][START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF], see Proposition 1.17. We refer to these earlier works for implementation details, extensive numerical studies, and comparisons with several competing methods.

We present in §1.3 quantitative estimates of the size of the optimal D-diffusion and Meikonal stencils, which are identified in the previous subsection §1.2. For a given tensor condition number κ, the stencil radius strongly depends on the angle θ between the grid axes and the tensor eigenvectors. Our estimates provide, for a given κ, both the worst and the average stencil radius w.r.t. θ, see Theorem 1.19.

We briefly discuss in §1.4 the case of PDEs with non-constant coefficients. Indeed the results of this paper, being strictly local, are for simplicity stated for PDEs with constant coefficients. We describe problem instances to which our results apply particularly well, and show how Theorem 1.19 could be the starting point of a numerical error analysis.

Notation: Let f, g and det(f, g) denote respectively the scalar product and the determinant of two vectors f, g ∈ R 2 . The perpendicular vector f ⊥ to f obeys f ⊥ , g = det(f, g). Let S 2 denote the collection of symmetric 2 × 2 matrices, and let S + 2 denote the subset of those which are positive definite. To each D ∈ S + 2 we associate the norm defined for all e ∈ R 2 by e D := e, De .

Let e T (resp. A T ) denote the transpose of a vector e (resp. matrix A), and let e⊗e := ee T denote the tensor product of the vector e ∈ R 2 with itself, which is a positive semi-definite matrix. We denote by Hull(V ) the convex hull of a set V ⊆ R 2 , and use the shorthand

[x 1 , • • • , x r ] := Hull({x 1 , • • • , x r }),
for any x 1 , • • • , x r ∈ R 2 . We say that some vectors e 1 , • • • , e r ∈ R 2 are ordered trigonometrically if they are non-zero, pairwise non-positively collinear, and if the normalized vectors e 1 / e 1 , • • • , e r / e r are indexed according to the cyclic order on the unit circle. In that case we say that e i and e i+1 are trigonometrically consecutive within this family of vectors, for any 0 ≤ i < r and with the convention e 0 := e r . Finally, we denote by Z the collection of elements of Z 2 which have co-prime coordinates Z := {e ∈ Z 2 ; gcd(e) = 1}, where gcd(e) := gcd(a, b) if e = (a, b) ∈ Z 2 . Note that (0, 0) / ∈ Z.

Structure preserving schemes for the Diffusion and Eikonal equations

We introduce in this section discretization schemes for anisotropic diffusion and eikonal PDEs, based on adequate notions of stencils. We focus on constant PDE coefficients, see §1.4 for a discussion on non-constant coefficients. This restriction is not limiting since the definitions and results presented in §1.1, §1.2 and §1.3 are purely local. As mentioned in the introduction, our main results are limited to discretizations on two dimensional cartesian grids, which take the form

x 0 + hRZ 2 ,
where x 0 is an offset, h > 0 a scale, and R a rotation. Without loss of generality, we assume that these parameters take their canonical values, so that the discretization grid is simply Z 2 . The results presented in this subsection are fairly classical, but are nevertheless required to justify the axioms chosen for the definition of diffusion and eikonal stencils. Classical references on these topics include [START_REF] Ciarlet | Discrete maximum principle for finite-difference operators[END_REF] on the discrete maximum principle, [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF] on anisotropic diffusion in image processing, and [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: theory and algorithms[END_REF] on anisotropic eikonal equations. The results of this subsection are established in Appendix A.

Discretization of anisotropic Diffusion. We introduce the concept of D-Diffusion stencil in Proposition 1.1, show how it leads to a natural discretization of diffusion PDEs in Proposition 1.2, and justify the so-called Non-Negativity axiom by a stability property in Proposition 1.4. Throughout this section, the diffusion tensor D is given and assumed to be positive definite, but it could be severely ill conditioned.

Definition 1.1. A D-diffusion stencil, D ∈ S +
2 , is a finitely supported map γ : Z → R obeying:

• Consistency: 1 2 e∈Z γ(e) e ⊗ e = D.

• Symmetry: γ(-e) = γ(e) for all e ∈ Z.

• Non-negativity: γ(e) ≥ 0 for all e ∈ Z.

A D-diffusion stencil γ may also be regarded as the finite collection of points supp(γ), together with the weights γ(e), e ∈ supp(γ). Our next result justifies the Consistency axiom.

Proposition 1.2. Consider a finitely supported γ : Z → R obeying property Symmetry of Definition 1.1. Then property Consistency is equivalent to each of the following properties:

• Consistency':

e∈Z γ(e)(u(x + e) -u(x)) = Tr(D∇ 2 u(x)) for any quadratic u : R 2 → R.

• Consistency":

1 2 e∈Z γ(e)(u(x + e) -u(x)) 2 = ∇u(x) 2 D for any linear u : R 2 → R.
To each finitely supported γ : Z → R we attach the linear operator L γ defined for all u : Z 2 → R and x ∈ Z 2 by

L γ u(x) := e∈Z γ(e)(u(x + e) -u(x)).
(2)

If γ obeys properties Symmetry and Consistency of Definition 1.1, then L γ is by Proposition 1.2 a natural discretization of Tr(D∇ 2 u) and of div(D∇u). Note that these two differential operators are identical, as well as their discretizations, since the diffusion tensor D is constant over the domain. Divergence form and non-divergence form diffusion PDEs differ however, as well as their discretizations, in the case of a space varying tensor field, see §1.4.

Observing that Id = e 1 ⊗ e 1 + e 2 ⊗ e 2 , where (e 1 , e 2 ) denotes the canonical basis of R 2 , we obtain the Id-diffusion stencil defined by γ(±e 1 ) = γ(±e 2 ) = 1 and γ = 0 elsewhere. The operator (2) associated to this specific stencil is the standard five points discretization of the Laplace operator ∆.

The maximum principle is a fundamental property of the diffusion PDE:

if u ∈ C 2 (R + × R 2 ) is bounded and obeys ∂ t u = div(D∇u), then for any t ≥ 0 one has inf R 2 u 0 ≤ u t ≤ sup R 2 u 0 ,
where u t : x → u(t, x). Here and after, if u : X → R is a map and c ∈ R is a constant, then "u ≥ c" stands for "u(x) ≥ c for all x ∈ X". The next proposition shows that the Non-Negativity axiom of D-diffusion stencils is equivalent (under a CFL condition) to the discrete maximum principle for the operator Id +δL γ which is associated to an explicit time step of discretized diffusion.

Definition 1.3. In the following, an operator is a continuous map from L ∞ (Z 2 ) to itself. An operator A obeys the discrete maximum principle iff for all u ∈ L ∞ (Z 2 ) one has on

Z 2 inf Z 2 u ≤ Au ≤ sup Z 2 u.
Proposition 1.4. Let γ : Z → R be finitely supported. Then the following are equivalent:

• The weights γ are non-negative.

• There exists δ > 0 such that Id +δL γ obeys the maximum principle.

Assume now that the weights γ are non-negative. Then maximum principle for Id +δL γ is equivalent to the CFL condition 0 ≤ δ e∈Z γ(e) ≤ 1. In addition the operator (Id -δL γ ) -1 , describing an implicit time step of discretized diffusion, obeys the maximum principle for any δ > 0.

The discrete maximum principle also makes sense for static boundary value problems, without time evolution. See the discussion in [START_REF] Ciarlet | Discrete maximum principle for finite-difference operators[END_REF] which is based on the concept of M -matrix.

Remark 1.5 (Stencils with support outside of Z). It may seem relevant to consider D-diffusion stencils γ : Z 2 \{0} → R, obeying the axioms of Definition 1.1 up to the replacement of Z with the larger set Z 2 \ {0}. Note however that a standard D-diffusion stencil γ : Z → R may then be defined by γ(e) := λ>0 λ 2 γ(λe). In addition γ is preferable to the original γ, since it is supported on a smaller set in the sense of convex hull inclusion. Indeed e/ gcd(e) ∈ [-e, e] ⊆ Hull(supp(γ)) for any e ∈ supp(γ).

Discretization of Anisotropic Eikonal equations. Given a bounded domain Ω ⊆ R 2 , and a positive definite tensor D ∈ S + 2 , one may consider the eikonal equation

∀x ∈ Ω, ∇u(x) D = 1, ∀x ∈ R 2 \ Ω, u(x) = 0. ( 3
)
This PDE is known to admit a unique viscosity solution [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF], which is the escape time u : R 2 → R + from the domain Ω, where path length is measured using the metric • M , with M := D -1 . The reason for the matrix inversion D = M -1 is that the norm • M is intended for measuring vectors, whereas the dual norm • D is intended for measuring gradients (co-vectors).

Bellman's optimality principle, illustrated in Figure 1, expresses that in order to escape the domain Ω one must first escape any sub-domain

V . Hence if x ∈ V ⊆ Ω, then u(x) = min y∈∂V x -y M + u(y). (4)
The next definition introduces the notion of eikonal stencil, used in Definition 1.7 to mimic Bellman's optimality principle at the discrete level.

Definition 1.6. An M -eikonal stencil is a finite collection of points e 1 , • • • , e r ∈ Z obeying:

• Ordering: the points are indexed according to the cyclic trigonometric order. Let e 0 := e r .

• Orientation: det(e i , e i+1 ) > 0 for all 0 ≤ i < r.

• Acuteness: e i , M e i+1 ≥ 0 for all 0 ≤ i < r.

The properties Ordering and Orientation of this definition guarantee that the union of segments [e 0 , e 1 ] ∪ • • • ∪ [e r-1 , e r ] delimits a neighborhood V 0 of the origin. The Hopf-Lax operator, introduced below, mimics the r.h.s. of (4) at the discrete level, on the neighborhood V 0 translated around a point x ∈ Z 2 of interest, and denoted by V (x) in Figure 1. Hopf formulas were introduced by Evans [?] as explicit representations of some Hamilton-Jacobi equations with constant coefficients. x -y i,t M + (1 -t)u(y i ) + tu(y i+1 ),

where y i := x + e i and y i,t := (1 -t)y i + ty i+1 , for all 0 ≤ i < r and all t ∈ [0, 1].

Semi-Lagrangian schemes [START_REF] Falcone | Semi-Lagrangian schemes for Hamilton-Jacobi equations, discrete representation formulae and Godunov methods[END_REF][START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF] for the eikonal equation take the form of a fixed point problem: find u :

Z 2 → R such that ∀x ∈ X, u(x) = Λu(x), ∀x ∈ Z 2 \ X, u(x) = 0, (6)
where the discrete domain X ⊆ Z 2 approximates the PDE domain Ω in (3). The eikonal PDE has a causal structure: denoting by y * a point at which the minimum in Bellman's optimality principle (4) is attained, the PDE solution satisfies u(x) > u(y * ). The discrete counterpart of this property, which may or may not hold, is introduced in the next definition.

Definition 1.8. The Hopf-Lax operator (5) has the causality property iff, for any u : Z 2 → R and any x ∈ Z 2 , one has denoting by 0 ≤ i < r and t ∈ [0, 1] the minimizers in (5)

Λu(x) > u(y i ) if t < 1, and Λu(x) > u(y i+1 ) if t > 0.
Numerical solvers of (6) involve a mutable map u : Z 2 → R ∪ {+∞}, initialized by u = 0 on Z 2 \ X and u = +∞ on X. They iteratively substitute u(x i ) with Λu(x i ), an operation referred to as an Hopf-Lax update [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF], for some sequence of points (x i ) i≥0 in the discrete domain X, until a convergence criterion is met. The point sequence is fixed a-priori or determined at run time.

If the causality property holds, then (6) can be efficiently solved using the Fast Marching algorithm [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF], which is a variant of Dijkstra's algorithm on graphs, see Figure 1. This method

x y Γ V x V x X x y1 y2 V x
Figure 1: From left to right. I: Bellman's optimality principle (4) expresses that any path starting from x, and escaping Ω, must intersect ∂V at some point y. II: The discretization of the eikonal PDE mimics this principle, see ( 5) and (6). III: The minimum defining Λu(x) is attained on a facet of ∂V (x), here a segment of endpoints y 1 , y 2 . The value of the solution u to (6) at x is thus sensitive to the values at y 1 and y 2 . IV: The discrete Causality property guarantees that the solution to (6) obeys u(x) > u(y), whenever the dependency graph features an arrow from x to y. This is an essential requirement for the fast marching algorithm.

determines the sequence of points (x i ) i≥0 at run-time using a priority queue, and guarantees that the number of Hopf-Lax updates required at any point x ∈ X to reach convergence is bounded a-priori by the number of elements of its stencil. Without the causality property, alternative methods must be used [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF]?], which are inspired by the Bellman-Ford algorithm on graphs. They are typically slower and offer no a-priori bound on the number of Hopf-Lax updates required for each point. The causality property, which is thus highly desirable from the algorithmic point of view, turns out to be equivalent to a geometric property of the stencils, as shown in the next proposition. This proposition can be extended to arbitrary dimension, and to arbitrary norms, which need not be of the form (1) but may be non-smooth or even asymmetric [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: theory and algorithms[END_REF][START_REF] Vladimirsky | Label-setting methods for Multimode Stochastic Shortest Path problems on graphs[END_REF][START_REF] Mirebeau | Efficient fast marching with Finsler metrics[END_REF].

Proposition 1.9 (Causality is equivalent to Acuteness [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: theory and algorithms[END_REF][START_REF] Vladimirsky | Label-setting methods for Multimode Stochastic Shortest Path problems on graphs[END_REF]). Under the assumptions of Definition 1.7, the causality of the Hopf-Lax operator is equivalent to the Acuteness property of the stencil in Definition 1.6.

Remark 1.10 (Stencils with points outside Z). It may seem relevant to consider an M -eikonal stencil e 1 , • • • , e r ∈ Z 2 \ {0} obeying the axioms of Definition 1.6 up to the replacement of Z with the larger set Z 2 \ {0}. Note however that a standard M -eikonal stencil may then be defined by e 1 / gcd(e 1 ), • • • , e r / gcd(e r ). In addition, the new stencil is preferable since it is supported on a smaller set in the sense of convex hull inclusion: indeed

e i / gcd(e i ) ∈ [0, e i ] ⊆ [e 1 , • • • , e r ].

Minimal stencils in the sense of convex hull inclusion

We identify in this section the smallest D-diffusion and M -eikonal stencils in the sense of convex hull inclusion, see Theorems 1.12 and 1.13, for any D, M ∈ S + 2 . Empirical experience tells that the most robust and accurate PDE discretizations are typically achieved with the smallest stencils. Small stencils also limit discretization issues close to the domain boundary, and ease parallel implementations. In this light, our results characterize the best possible stencils for the discretization of anisotropic diffusion and anisotropic eikonal equations on cartesian grids -unless one is looking for discretizations with higher order consistency. The statement of our results requires the introduction of some elementary notions of lattice geometry, illustrated in Figures 2 and3. The results presented in this section are established in §2.

Vor(M ) := {p ∈ R 2 ; ∀e ∈ Z 2 , p M ≤ p -e M }, Vor(M ; e) := {p ∈ Vor(M ); p M = p -e M }. An M -Voronoi vector is an element e ∈ Z 2 \ {0} such that Vor(M ; e) = ∅. It is said strict iff Vor(M ; e) is not reduced to a point, in which case it is a 1-dimensional segment.
Note that p M ≤ p -e M if and only if 2 p, M e ≤ e 2 M , for any p, e ∈ R 2 . Hence Vor(M ) is an intersection of half-spaces, thus a convex polytope, which has one facet Vor(M ; e) with exterior normal is M e for each strict M -Voronoi vector e. The next two theorems show that the smallest D-diffusion and M -eikonal stencils are supported precisely on the set of strict M -Voronoi vectors, with D := M -1 . Theorem 1.12. For any D ∈ S + 2 , there exists a unique D-diffusion stencil supported on the set V of strict M -Voronoi vectors, where M := D -1 . Furthermore, this stencil is minimal in the following sense: any other D-diffusion stencil γ satisfies Hull(V ) ⊆ Hull(supp(γ)).

Theorem 1.13. For any M ∈ S + 2 , the set V of strict M -Voronoi vectors, ordered trigonometrically, is an M -eikonal stencil. Furthermore this stencil is minimal in the following sense: any other

M -eikonal stencil e 1 , • • • , e r satisfies Hull(V ) ⊆ Hull({e 1 , • • • , e r }).
The next two definitions distinguish, within the set S + 2 of symmetric positive definite matrices, an open dense subset of generic matrices, and a complementary subset of exceptional matrices obeying a geometric property. This allows in Proposition 1.16 to characterize strict M -Voronoi vectors as the smallest elements of Z w.r.t. the metric • M , and in Proposition 1.17 to relate our stencil constructions with earlier works.

Definition 1.14. A basis of Z 2 is a pair (f, g) ∈ (Z 2 ) 2 such that | det(f, g)| = 1. The basis is said M -orthogonal, where M ∈ S + 2 , iff f, M g = 0.
Definition 1.15. A matrix M ∈ S + 2 is said generic iff there exists no M -orthogonal basis of Z 2 .

Proposition 1.16. The set of strict M -Voronoi vectors, where M ∈ S + 2 , consists of

κ = 5, θ ∈ [0, π/2]. Right: κ ∈ [1, 16], θ = π/3). • if M is non-generic: the four smallest elements of Z w.r.t. • M . • if M is generic: the six smallest elements of Z w.r.t. • M .
As far as generic matrices are concerned, the stencil constructions of Theorems 1.12 and 1.13 are not new, but coincide with those of [START_REF] Fehrenbach | Sparse Non-negative Stencils for Anisotropic Diffusion[END_REF] for diffusion and [START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF] for eikonal PDEs. We refer to these works for convergence results, numerical experiments and comparisons with other approaches. The stencil constructions of [START_REF] Fehrenbach | Sparse Non-negative Stencils for Anisotropic Diffusion[END_REF][START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF] are based on another tool of lattice geometry, called obtuse superbases see §2.2, and also make sense in three dimensions. However no optimality result is known in this latter case.

Proposition 1.17. If M ∈ S + 2 is generic, then the D-diffusion stencil of Theorem 1.12 is the one presented in [START_REF] Fehrenbach | Sparse Non-negative Stencils for Anisotropic Diffusion[END_REF], and the M -eikonal stencil of Theorem 1.13 is the one presented in [START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF].

Our D-diffusion stencils are also equivalent to the (two dimensional only) construction of [START_REF] Bonnans | A fast algorithm for the two dimensional HJB equation of stochastic control[END_REF]. Non-equivalent constructions can be found in [START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Convex Hamilton-Jacobi Equations[END_REF] for M -eikonal stencils, and [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF] for D-diffusion stencils. They are, inevitably, less local and accurate, as exposed by the comparisons in [START_REF] Fehrenbach | Sparse Non-negative Stencils for Anisotropic Diffusion[END_REF][START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF].

Average size of minimal stencils

We introduced in §1.1 the concepts of D-diffusion stencil and M -eikonal stencil, and identified in §1.2 the smallest such stencils in the sense of convex hull inclusion. They turn out to be both supported on the set of strict M -Voronoi vectors, with M := D -1 . Our next result, Theorem 1.19 is a quantitative estimate of the size these stencils, which supplements their qualitative optimality property. The result presented in this section is established in §3.

Definition 1.18. For each κ ∈ [1, ∞[, θ ∈ R, we introduce the positive definite 2 × 2 matrix M κ (θ) := κ -1 e θ ⊗ e θ + κ e ⊥ θ ⊗ e ⊥ θ ,
with e θ := (cos θ, sin θ).

We denote by V κ (θ) the collection of all strict M κ (θ)-Voronoi vectors, and define

R κ (θ) := max e∈Vκ(θ) e , S κ (θ) := max e∈Vκ(θ) e Mκ(θ) . (8)
We allow ourselves a slight abuse of notation: "e θ " always refers to the unit vector (7) of direction θ, although at times we introduce families e 1 , • • • , e r of vectors in Z 2 .

The radius of M -diffusion stencils is measured in (8) using two different norms: the "extrinsic" euclidean norm (8, left), and the "intrinsic" • Mκ(θ) -norm (8, right) which is tied to the on the variable θ is highly irregular, and reflects how well the slope tan θ is approximated by rationals with small denominator. Both radii measures R κ and S κ have their relevance and merits, as discussed in Appendix B of [START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF] with an heuristic analysis of the accuracy of several discretizations of the eikonal PDE.

0 Π 4 Π 2 5. 2. 3. R Κ Θ , Κ 10 0 Π 4 Π 2 10. 5. 20. 3. 15. 7. R Κ Θ , Κ 50 0 Π 4 Π 2 2. 3. 1.5 S Κ Θ , Κ 10 0 Π 4 Π 2 5. 2. 3. 1.5 7. S Κ Θ , Κ 50
The maximum condition number κ, of the tensor field associated to an anisotropic PDE, is generally well known: it is problem data, reflecting a continuous model. In contrast we regard the angle θ, between the cartesian grid axes and the tensor eigenvectors, as a random quantity uniformly distributed in [0, π[. In other words the preferred directions of the grid and of the PDE are viewed as independent of each other. Our main result, Theorem 1.19 illustrated in Figure 5, is thus devoted to the estimation of averaged quantities w.r.t. the angular variable θ, namely the L p ([0, π[)-norms of R κ and S κ . The uniform equivalence of two expressions, w.r.t. the parameter κ, is denoted as follows:

A(κ) ≈ B(κ) ⇔ ∃C, c > 0, ∀κ ≥ 1, cB(κ) ≤ A(κ) ≤ CB(κ). (9)
Theorem 1.19. For any p ∈ [1, ∞] one has uniformly in κ (L p norms are on the interval [0, π[)

R κ L p ≈ κ 1 2 S κ L p , S κ L p ≈      κ 1 2 -1 p if p > 2, (ln κ) 1 2 if p = 2, 1 if p < 2.
In applications, κ = 10 is already a pronounced anisotropy, and κ = 100 is presumably the most degenerate anisotropy that can conceivably be handled with the proposed discretizations of the diffusion and eikonal PDEs. More specialized approaches are recommended for stronger anisotropies, such as asymptotic preserving formulations [?]. Nevertheless Theorem 1.19 is a good indicator of the effective spread of our diffusion stencils, since the asymptotic behavior is quickly attained as shown in Figure 5.

PDEs with non-constant coefficients

We shortly discuss some PDE models interest, which are intrinsically anisotropic, often discretized on cartesian grids, and thus directly benefit from the results obtained in this paper.

Divergence form anisotropic diffusion. Let Ω ⊆ R 2 be a bounded domain, and let D : Ω → S + 2 be a continuous field of diffusion tensors. Let h > 0 be a discretization scale, and let Ω h := Ω∩hZ 2 be the discretization grid. The elliptic energy E D associated to D has a natural discretization

E h . For u : Ω → R E D (u) := ˆΩ ∇u(x) 2 D(x) dx, E h (u) := 1 2 x∈Ω h e∈Z γ D(x) (u(x + he) -u(x)) 2 .
We denoted by γ D : Z → R, where D ∈ S + 2 , the minimal D-diffusion stencil of Theorem 1.12. The contributions to E h involving the evaluation of u outside Ω should be omitted for Neumann boundary conditions, or included by extending u as 0 on R 2 \ Ω for Dirichlet boundary conditions. Taking the gradient flow of E D (resp. E h ) w.r.t. the L 2 metric, one obtains the divergence form diffusion ∂ t u(t, x) = div(D(x)∇u(t, x)) (resp. a natural discretization of this PDE obeying the maximum principle). Before discussing specific applications, we illustrate the relevance of Theorem 1.19 through an heuristic accuracy analysis. By a fourth order Taylor expansion of u close to x, we obtain

1 2h 2 e∈V γ D (u(x + he) -u(x)) 2 = ∇u(x) 2 D + h 2
|E D (u) -E h (u)| ≤ ˆΩ ∇u(x) 2 D(x) dx -h 2 x∈Ω h ∇u(x) 2 D(x) + h 2 max x∈Ω Λ(x) D(x)   h 2 x∈Ω h r(D(x)) 2   + o(h 2 ).
The first contribution, a quadrature error, is scheme independent and expected to be O(h 2 ) (at least on square domains, for which the boundary representation is exact). The second contribution, an O(h 2 ) discretization error, is proportional to the average squared stencil radius. This shows the relevance of the estimate R κ L 2 ≈ √ κ ln κ obtained in Theorem 1.19. Note that the worst case stencil radius is significantly larger: R κ L ∞ ≈ κ .

Non-linear anisotropic diffusion PDEs take the form ∂ t u(t, x) = div(D u (t, x)∇u(t, x)), with a solution dependent tensor field D u . They play an important role in image processing, being used for image deblurring, sharpening, compression, . . . see [START_REF] Weickert | Theoretical foundations of anisotropic diffusion in image processing[END_REF] for an overview. From a practical point of view, the PDE non-linearity merely requires to update the diffusion tensor field every few time steps. Our D-diffusion stencils are used to solve numerically non-linear anisotropic diffusion PDEs in [START_REF] Fehrenbach | Sparse Non-negative Stencils for Anisotropic Diffusion[END_REF]. Note that our results are irrelevant for non-linear isotropic diffusion, such as the Perona-Malik model1 with diffusion tensors D u (t, x) := Id /(1 + ∇u(t, x) ).

We end this paragraph by mentioning a celebrated model of the first layer of the visual cortex [START_REF] Petitot | The neurogeometry of pinwheels as a sub-Riemannian contact structure[END_REF] as the three dimensional manifold V 1 := R 2 × P 1 , where P 1 denotes [0, π[ equipped with periodic boundary conditions. This manifold is equipped with a degenerate hypo-elliptic diffusion operator, which can be approximated using a strongly anisotropic diffusion tensor field: for all (x 0 , x 1 , θ) ∈ V 1 , denoting by κ 1 a large parameter

D(x 0 , x 1 , θ) =   e θ ⊗ e θ + κ -2 e ⊥ θ ⊗ e ⊥ θ 1   = (κM κ (θ)) -1 1 .
The two-dimensional discretization scheme proposed in this paper for anisotropic diffusion easily extend to this case, thanks to the block diagonal structure. This model fits Theorem 1.19 particularly well, since the condition number κ is fixed, and all directions θ are equally represented.

Non-Divergence form diffusion. Monotone second-order differential operators can, under mild assumptions, be expressed in Hamilton-Jacobi-Bellman (HJB) form as infs and sups of linear non-divergence form diffusion operators. Denoting by A and B two parameter spaces, and by D : A × B → S + 2 a family of diffusion tensors, the continuous operator and its discretization read:

inf a∈A sup b∈B Tr(D(a, b)∇ 2 u(x)), inf a∈A sup b∈B e∈Z γ D(a,b) (e) u(x + he) -u(x) h 2 .
This strategy is applied in [START_REF] Bonnans | A fast algorithm for the two dimensional HJB equation of stochastic control[END_REF] to the HJB equation of stochastic control. See also [START_REF] Benamou | Monotone and Consistent discretization of the Monge-Ampere operator[END_REF] for a monotone and consistent discretization of the Monge-Ampere operator, which fits in this setting thanks to the identity det(∇ 2 u) = inf{Tr(D∇ 2 u); D ∈ S + 2 , det D = 1}. Note that the constraint det(D) = 1 includes diffusion tensors D ∈ S + 2 of arbitrary condition number, in all orientations.

Anisotropic Eikonal equations. Anisotropic eikonal equations allow to estimate riemannian distances, and compute the associated shortest paths. Anisotropy is particularly useful in applications to image segmentation, where it is used to guide the shortest paths along one dimensional structures of interest, such as blood vessels or organ boundaries [START_REF] Chen | Vessel Extraction Using Anisotropic Minimal Paths and Path Score[END_REF].

Strongly anisotropic riemannian metrics can also be used to approximate degenerate, subriemannian geometries. For instance, the Reed-Shepp car model is posed on the manifold R 2 ×P 1 , where again P 1 denotes [0, π[ equipped with periodic boundary conditions. This model involves a degenerate metric which can be approximated as follows: for all (x 0 , x 1 , θ) ∈ R 2 × P 1 , denoting by κ 1 a large parameter The two-dimensional numerical scheme proposed in this paper for anisotropic eikonal equations easily extends to this model, thanks to the block diagonal structure of the metric. See [SBD + 15] for an application to blood vessel segmentation. This model fits Theorem 1.19 particularly well, since the condition number is fixed, and all anisotropy orientations are equally represented.

M (x 0 , x 1 , θ) =   e θ ⊗ e θ + κ 2 e ⊥ θ ⊗ e ⊥ θ 1   = κM κ (θ) 1 .

Correctness and minimality of Voronoi based stencils

We establish in this section the results announced in §1.2. Specifically, Theorem 1.13 is proved in §2.1, Proposition 1.16 in §2.2, and Theorem 1.12 in §2.3. The results of each subsection are used in the following subsection, and in particular Theorem 1.13 is used to prove Theorem 1.12. The statement of Proposition 1.17 on eikonal stencils (resp. diffusion stencils) is proved in §2.2 (resp. §2.3). We denote by Cone(f, g) the convex cone spanned by two vectors f, g ∈ R 2 .

Correctness and minimality of eikonal stencils

We establish in Corollary 2.3 the correctness part of Theorem 1.13: strict M -Voronoi vectors, when ordered trigonometrically, do define an M -eikonal stencil. The rest of Theorem 1.13 is established in Proposition 2.6. Two preliminary lemmas are required.

Lemma 2.1. All M -Voronoi vectors belong to the set Z, for all M ∈ S + 2 .

Proof. Any element of Z 2 \ {0} can be written under the form λe, where λ is a positive integer and e ∈ Z. If λ > 1, then for any p ∈ R 2 we obtain

e -p 2 = 1 λ λe -p 2 + 1 - 1 λ p 2 -(λ -1) e 2 < max{ λe -p , p 2 },
which shows that Vor(M ; λe) = ∅, hence that λe is not a Voronoi vector.

Lemma 2.2. Let f, g be M -Voronoi vectors such that Vor(M ; f )∩Vor(M ; g) = ∅, where M ∈ S + 2 . Then f, M g ≥ 0.

Proof. Consider x ∈ Vor(M ; f ) ∩ Vor(M ; g). Then x -(f + g) M ≥ x -f M = x -g M = x M . Thus, as announced 0 ≤ x -(f + g) 2 M -x -f 2 M -x -g 2 M + x 2 M = 2 f, M g .
Corollary 2.3. For any M ∈ S + 2 , the collection of strict M -Voronoi vectors (e 1 , • • • , e r ), ordered trigonometrically, defines an M -eikonal stencil.

Proof. By Lemma 2.1, M -Voronoi vectors belongs to Z as required. The sets Vor(M ; e i ), 1 ≤ i ≤ r, are the faces of the convex polytope Vor(M ), ordered trigonometrically. The exterior normal to Vor(M ; e i ) is M e i , hence 0 < det(M e i , M e i+1 ) = det(M ) det(e i , e i+1 ) as required (Orientation), for all 1 ≤ i ≤ r. Since two consecutive faces have non-empty intersection, one has e i , M e i+1 by Lemma 2.2, as required also (Acuteness), which concludes the proof.

Our next objective is to establish the minimality property of the M -eikonal stencils of Theorem 1.13. For that purpose, achieved in Proposition 2.6, two preliminary results are required.

Recall that [p 1 , • • • , p n ] := Hull({p 1 , • • • , p n }) for any p 1 , • • • , p n ∈ R 2 . Lemma 2.4. Let f, g ∈ Z 2 be such that | det(f, g)| > 1. Then the triangle [0, f, g] contains a point e ∈ Z distinct from its vertices. Proof. Since | det(f, g)| > 1, the map (α, β) ∈ Z 2 → αf + βg ∈ Z 2 is not surjective. Hence there exists (α, β) ∈ Q 2 ,
at least one of them non-integer, such that αf + βg ∈ Z 2 . Up to replacing (α, β) with (α -m, β -n), for some (m, n) ∈ Z 2 , we may assume that α, β ∈ [0, 1]. Up to replacing (α, β) with (1 -α, 1 -β), we may assume that α + β ≤ 1. The point e := αf + βg ∈ Z 2 thus belongs to T and is distinct from its vertices. This establishes the announced result, unless gcd(e) > 1, in which case we can consider e/ gcd(e). Proof. Let e be a strict M -Voronoi vector, and let f = e i and g = e i+1 be trigonometrically consecutive elements of the M -eikonal stencil (e 1 , • • • , e r ) such that e ∈ Cone(f, g). Then e = αf + βg for some α, β ∈ R + . By Corollary 2.5 we can assume without loss of generality that det(f, g) = 1, thus the vectors (f, g) form a basis of Z 2 , hence α, β ∈ Z. Since e ∈ Z, one has gcd(α, β) = 1. Assuming for contradiction that α ≥ 1 and β ≥ 1 we obtain for any p ∈ R 2

α p -f 2 + β p -g 2 = (α + β) p 2 -2 αf + βg, p + α f 2 + β g 2 = (α + β -1) p 2 + p -e 2 -2αβ f, g .
Since αβ f, g ≥ 0, and since α + β = (α + β -1) + 1, it follows that min{ p -f , p -g } ≤ max{ p , p -e }. Choosing p in the relative interior of Vor(M ; e), so that the elements of Z 2 closest to p are precisely 0 and e, we obtain a contradiction. Thus (α, β) equals (1, 0) or (0, 1), therefore e equals f or g, which shows as announced that the Voronoi vector e belongs to the stencil convex hull [e 1 , • • • , e r ].

Description of strict M -Voronoi vectors

In this section, we introduce the concept of superbase of the lattice Z 2 . This allows us to establish the correspondence between the M -eikonal stencil of Theorem 1.13 and the one constructed in [START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF], as announced in Proposition 1.17. We also characterize strict M -Voronoi vectors as the smallest elements of Z w.r.t. • M , as announced in Proposition 1.16. Definition 2.7. A superbase of Z 2 is a triplet (e 0 , e 1 , e 2 ) such that (e 1 , e 2 ) is a basis of Z 2 and e 0 + e 1 + e 2 = 0. It is said M -obtuse, where M ∈ S + 2 , iff e i , M e j ≤ 0 for all 0 ≤ i < j ≤ 2.

It is known that for each M ∈ S + 2 there exists at least one M -obtuse superbase. Interestingly, a similar property holds in dimension 3, but not in dimension 4 or higher. The proof is constructive, and follows from the analysis of an algorithm by Selling [START_REF] Selling | Ueber die binären und ternären quadratischen Formen[END_REF]. Obtuse superbases are an important tool in lattice classification [START_REF] Conway | Low-Dimensional Lattices. VI. Voronoi Reduction of Three-Dimensional Lattices[END_REF]. They are also at the foundation of the M -eikonal stencil of [START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF], described in Proposition 2.8 (ii), and of the D-diffusion stencil of [START_REF] Fehrenbach | Sparse Non-negative Stencils for Anisotropic Diffusion[END_REF], defined in Corollary 2.11. The proof of the following proposition is left as an easy exercise.

Proposition 2.8. (i) Let M ∈ S + 2 be non-generic, and let (e 1 , e 2 ) be an M -orthogonal basis, sorted so that det(e 1 , e 2 ) = 1. Then (e 1 , e 2 , -e 1 , -e 2 ) is an M -eikonal stencil.

(ii) Let M ∈ S + 2 , and let (e 0 , e 1 , e 2 ) be an M -obtuse superbase, sorted so that det(e 1 , e 2 ) = 1. Then (e 0 , -e 2 , e 1 , -e 0 , e 2 , -e 1 ) is an M -eikonal stencil.

Proposition 2.9. Let M ∈ S + 2 be non-generic and let n = 4 (resp. generic and let n = 6). Then the set of strict M -Voronoi vectors precisely consists of the n members

(f 1 , • • • , f n ) of the M -eikonal stencil constructed in Proposition 2.8.
Proof. We first observe that, if (e 1 , e 2 ) is a basis of Z 2 (resp. (e 0 , e 1 , e 2 ) is a superbase), then [e 1 , e 2 , -e 1 , -e 2 ] ∩ Z = {e 1 , e 2 , -e 1 , -e 2 } (resp. [e 0 , -e 2 , e 1 , -e 0 , e 2 , -e 1 ] ∩ Z = {e 0 , -e 2 , e 1 , -e 0 , e 2 , -e 1 }. Indeed, by linear invariance, it suffices to check the case where e 1 = (1, 0) and e 2 = (0, 1) (resp. and e 0 = (-1, -1)), which is obvious.

By Theorem 1.13, strict M -Voronoi vectors form an eikonal stencil, which is minimal for convex hull inclusion hence must be a subset of

[f 1 , • • • , f n ]. Thus the set V of M -Voronoi vectors is included in [f 1 , • • • , f n ] ∩ Z = {f 1 , • • • , f n }. Since f i-1 , M f i+1 < 0 for all 0 ≤ i < n,
and in view of the Acuteness property of the M -eikonal stencil built of the strict M -Voronoi vectors, this inclusion is an equality, which concludes the proof.

We finally characterize strict M -Voronoi vectors as announced in Proposition 1.16.

Proof of Proposition 1.16. Let f 1 , • • • , f n be the M -eikonal stencil of Proposition 2.8, which by Proposition 2.9 is also the set of strict M -Voronoi vectors. Let e ∈ Z \ {f 1 , • • • , f n }, and let 0 ≤ i < n be such that e ∈ Cone(f, g) with f := f i and g := f i+1 . Then e = αf + βf for some coefficients α, β which must be integer (since (f, g) is a basis), non-negative (since e ∈ Cone(f, g)), co-prime (since e ∈ Z), and obey |α| + |β| > 1 (since e / ∈ {f, g}). Thus α ≥ 1 and β ≥ 1, and therefore e M ≥ αf + βg M ≥ f + g M since f, M g ≥ 0.

In the non-generic case, this implies e M ≥ e 1 + e 2 M = e 1 2 M + e 2 2 M > max{ e 1 M , e 2 M }, where (e 1 , e 2 ) is an M -orthogonal basis. In the generic case this implies e M ≥ e ie j M for some 0 ≤ i < j ≤ 2, where (e 0 , e 1 , e 2 ) is an M -obtuse superbase. Since e i , M e j < 0, we obtain e M ≥ e i -e j M > max{ e i M , e j M , e i + e j M } = max{ e 0 M , e 1 M , e 2 M }. We have proved in both cases that e M > max{

f 1 M , • • • , f n M } for all e ∈ Z \ {f 1 , • • • , f n }, which concludes the proof.

Correctness and minimality of diffusion stencils

We establish in Corollary 2.11 that there exists a D-diffusion stencil supported on the set of strict M -Voronoi vectors, as announced in Theorem 1.12. Its minimality, in the sense of convex hull inclusion, is established in Corollary 2.13. Our first step is an elementary algebraic lemma.

Proposition 2.10. Let (e i ) 2 i=0 be a superbase. Then (e ⊥ i ⊗ e ⊥ i ) 2 i=0 is a basis of S 2 , and for any

D ∈ S 2 one has D = - 0≤i≤2 e i+1 , De i+2 e ⊥ i ⊗ e ⊥ i . (11)
Proof. For any i, j ∈ {0, 1, 2}, one has denoting by δ ij the Kronecker symbol e j , e ⊥ i ⊗ e ⊥ i e j = e ⊥ i , e j 2 = det(e i , e j ) 2 = 1 -δ ij .

Therefore denoting by D the r.h.s. of ( 11) we obtain for any 0 ≤ j ≤ 2 e j , De j = -e j , De j+1 -e j-1 , De j = -e j+1 + e j-1 , De j = e j , De j .

In particular 0 = e 1 , (D -D)e 1 = e 2 , (D -D)e 2 = (e 1 + e 2 ), (D -D)(e 1 + e 2 ) . Since (e 1 , e 2 ) is a basis, this implies D = D as announced. As a result the three element family (e ⊥ i ⊗ e ⊥ i ) 2 i=0 generates the three dimensional space S 2 , hence it is a basis.

The next corollary describes the D-diffusion stencil used in [START_REF] Fehrenbach | Sparse Non-negative Stencils for Anisotropic Diffusion[END_REF], and shows that it is supported on strict M -Voronoi vectors, with M := D -1 , as announced in Theorem 1.12.

Corollary 2.11. Let (e 0 , e 1 , e 2 ) be a D-obtuse superbase. Then defining γ :

Z → R by γ(±e ⊥ i ) := -e i+1 , De i+2 , 0 ≤ i ≤ 2, and γ = 0 on Z \ {±e ⊥ 0 , ±e ⊥ 1 , ±e ⊥ 2 }
, we obtain the unique D-diffusion stencil which is supported on the set of strict M -Voronoi vectors, with M := D -1 .

Proof that γ defines a D-diffusion stencil. Since (e 0 , e 1 , e 2 ) is a superbase one has {±e ⊥ 1 , ±e ⊥ 2 , ±e ⊥ 3 } ⊆ Z, hence that γ : Z → R is well defined. Since it is D-obtuse one has e i+1 , De i+2 ≤ 0, hence γ is non-negative. By Proposition 2.10, this function γ defines a D-diffusion stencil.

Proof that γ is supported on the set of strict M -Voronoi vectors. For any f, g ∈ Z 2 one has

f, Dg = f ⊥ , M g ⊥ / det(M ).
Hence {e ⊥ 0 , e ⊥ 1 , e ⊥ 2 } is an M -obtuse superbase, and therefore (e ⊥ 0 , -e ⊥ 2 , e ⊥ 1 , -e ⊥ 0 , e ⊥ 2 , -e ⊥ 1 ) is an M -eikonal stencil by Proposition 2.8. By Proposition 2.6 all strict M -Voronoi vectors belong to

[e ⊥ 0 , -e ⊥ 2 , e ⊥ 1 , -e ⊥ 0 , e ⊥ 2 , -e ⊥ 1 ] ∩ Z = {e ⊥ 0 , -e ⊥ 2 , e ⊥ 1 , -e ⊥ 0 , e ⊥ 2 , -e ⊥ 1 } ⊇ supp(γ).
If M is generic, then by Proposition 2.9 there are six strict M -Voronoi vectors, and the result is proved. If M is non-generic, then there are only four strict M -Voronoi vectors, which are the members of an M -orthogonal basis -w.l.o.g. (e ⊥ 1 , e ⊥ 2 ) -and their opposites. Observe now that γ(e

⊥ 0 ) = -e ⊥ 1 , De ⊥ 2 = -e ⊥ 1 , M e ⊥ 2 / det(M ) = 0. Thus γ(±e ⊥ 0 ) = 0 and therefore γ is supported on the set {±e ⊥ 1 , ±e ⊥ 2 } of strict M -Voronoi vectors.
Finally, uniqueness follows from the Symmetry axiom of D-Diffusion stencils, and from Proposition 2.10 which states that

(e ⊥ i ⊗ e ⊥ i ) 2 i=0 is a basis of S 2 .
We show in Corollary 2.13 that the support of a D-diffusion stencil, ordered trigonometrically, defines an M -diffusion stencil, with M = D -1 . The minimality of the collection of strict M -Voronoi vectors as a D-diffusion stencil, announced in Theorem 1.12, thus follows from its minimality as an M -eikonal stencil, established in §2.1. Note that, in contrast with the two-dimensional setting studied in this paper, the three-dimensional D-diffusion and M -eikonal stencils of [START_REF] Fehrenbach | Sparse Non-negative Stencils for Anisotropic Diffusion[END_REF][START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF] are distinct, and in particular they do not have the same number of vertices. We denote S 1 := {z ∈ R 2 ; z = 1}. Lemma 2.12. Let γ : S 1 → R be symmetric and non-negative, i.e. γ(-e) = γ(e) ≥ 0 for all e ∈ S 1 . Assume that γ is supported on a finite set V ⊆ S 1 , and that e∈V γ(e)e ⊗ e = Id. Then any trigonometrically consecutive f, g ∈ V satisfy f, g ≥ 0.

Proof. We claim that for any ξ ∈ S 1 there exists e ∈ V such that ξ, e ≥ 1/ √ 2. Indeed, since V is symmetric, we obtain otherwise ξ, e 2 < 1/2 for all e ∈ V , but this implies the contradiction Let f, g ∈ V be trigonometrically consecutive, and let ξ := (f + g)/ f + g . By construction ξ, f = ξ, g ≥ ξ, e for all e ∈ V . Thus 1/ √ 2 ≤ ξ, f = ξ, g by the above. Therefore the angle between f and g is at most 2 arccos(1/ √ 2) = π/2, which concludes the proof.

1 = ξ 2 = e∈V γ(e) e,
Corollary 2.13. Let M ∈ S + 2 and let D := M -1 . Then the support of any D-diffusion stencil, ordered trigonometrically, defines an M -eikonal stencil. Therefore its convex hull contains all strict M -Voronoi vectors.

Proof. Let γ : Z → R be the D-diffusion stencil, and let A be an invertible matrix of positive determinant such that M = A T A. Apply Lemma 2.12 to γ :

S 1 → R defined by γ(Ae/ Ae ) = 1 2 γ(e) Ae 2 for all e ∈ Z. Then 0 ≤ Af, Ag = f, M g /( f 2 M g 2 M )
for any trigonometrically consecutive f, g ∈ supp(γ), thus f, M g ≥ 0. One has det(f, g) ≥ 0 since the set supp(γ) is symmetric, hence det(f, g) > 0 since otherwise f, g ∈ Z would be opposite, in contradiction with f, M g ≥ 0. As announced, we recognize an M -eikonal stencil. Finally, the convex hull of an M -eikonal stencil contains all strict M -Voronoi vectors by Proposition 2.6.

Average norm of Voronoi vectors

We estimate the norm of M κ (θ)-Voronoi vectors, as announced in Theorem 1.19, for a fixed condition number κ 2 , but in the L p ([0, π[) sense w.r.t. the anisotropy orientation θ. This type of result belongs to the well established field of lattice geometry, dating back to Lagrange [?]. Lattice geometry is currently the object of an important research activity, motivated by the numerous applications of the LLL algorithm [?] in particular in the field of cryptography. Nevertheless, the rather elementary Theorem 1.19 is original (to the author's knowledge), for the following reasons:

(a) Research on lattice geometry is mainly focused on high dimensions.

(b) For many applications of current interest, such as cryptography, it is natural to assume that the lattice Gram matrix M has integer coefficients, in contrast with the present setting where M = M κ (θ).

(c) The proposed random lattice model, which can be reformulated as Λ = DRZ d with D diagonal of eigenvalues κ ± 1 2 and R a random rotation, differs from classical [?] studies.

(d) The proposed result intertwines two mutually "singular" norms: the euclidean norm and the anisotropic • M norm. Indeed we estimate the euclidean norm of M -Voronoi vectors. In contrast, research in lattice geometry most commonly focuses on a single norm.

Reduced basis of a lattice

We introduce the concept of Minkowski basis [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF], and use it to describe the set of M -Voronoi vectors and to derive uniform bounds on their norms. To each M ∈ S + 2 we associate an • M -shortest non-zero vector e 1 (M ) ∈ Z 2 , and an • M -shortest linearly independent vector e 2 (M ) ∈ Z 2 e 1 (M ) ∈ argmin

e∈Z 2 \{0} e M , e 2 (M ) ∈ argmin f ∈Z 2 \e 1 (M )Z f M . (12)
The vectors e 1 (M ) and e 2 (M ) are uniquely determined, up to a change of sign, for all matrices M except for a Lebesgue negligible set. For non-generic matrices, an arbitrary minimizer is selected. In applications, the vectors (13) can be computed via Lagrange's algorithm [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF], a generalization of Euclid's gcd algorithm. The next lemma shows that the two vectors (12) form a basis of Z 2 , often referred to as a Minkowski reduced basis [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF].

Proposition 3.1. For any M ∈ S + 2 , the Minkowski vectors e 1 (M ), e 2 (M ) form a basis of Z 2 .

Proof. Let e 1 := e 1 (M ) and e 2 := e 2 (M ). One has gcd(e 1 ) = gcd(e 2 ) = 1, since otherwise the vector e i / gcd(e i ) would violate the minimality property of e i , i ∈ {1, 2}, see (12). Since e 2 ∈ Z 2 \ e 1 Z and gcd(e 1 ) = 1, the vectors e 1 , e 2 are not collinear, hence det(e 1 , e 2 ) = 0. Assume for contradiction that | det(e 1 , e 2 )| > 1 and, by Lemma 2.4, consider an element e ∈ Z 2 of the triangle [0, e 1 , e 2 ] distinct from its vertices: e = αe 1 + βe 2 , α, β ≥ 0, α + β ≤ 1. Since gcd(e 1 ) = gcd(e 2 ) = 1, the vector e does not lie on an edge [0, e i ], i ∈ {1, 2}. Hence by the triangular inequality, which is strict since e 1 , e 2 are not collinear and α, β are positive, we obtain e M < α e 1 M + β e 2 M ≤ e 2 M . This contradicts the minimality of e 2 M , hence | det(e 1 , e 2 )| = 1 which concludes the proof.

The following lemma upper bounds the orthogonality defect of a Minkowski reduced basis. See [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF] for extensions of this result to higher dimension. (ii) if e 1 , M e 2 < 0 then {±e 0 , ±e 1 , ±e 2 }, with e 0 := -(e 1 + e 2 ). In addition e 0 , M e 1 < 0 and e 0 , M e 2 < 0. (In the case e 1 , M e 2 > 0, replace e 2 with its opposite.)

Proof. If M is non-generic, then strict Voronoi vectors are the four smallest elements of Z by Proposition 1.16, hence they are precisely {±e 1 , ±e 2 } by construction (12). Furthermore strict Voronoi vectors feature an M -orthogonal basis by Proposition 2.9, hence case (i) holds. If M -is generic, then e 1 , M e 2 = 0 hence e 1 , M e 2 < 0 up to replacing e 2 with -e 2 . By Lemma 3.2 we obtain e 0 , M e 1 = -e 2 , M e 1 -e 1 2 M < 0 and e 0 , M e 2 = -e 1 , M e 2e 2 2 M < 0, with e 0 := -(e 1 + e 2 ). This shows that (e 0 , e 1 , e 2 ) is an M -obtuse superbase, and case (ii) holds by Proposition 2.9.

The scalars λ i (M ) are called the Minkowski minima, and depend continuously on the matrix M ∈ S + 2 , in contrast with r i (M ) and e i (M ) which take discrete values. Minkowski minima are tied together by Minkowski's second theorem [?], a fundamental result of lattice geometry stating that:

2 ≤ π √ det M λ 1 (M )λ 2 (M ) ≤ 4. (14)
For each κ ≥ 1, θ ∈ R, i ∈ {1, 2}, we denote by e κ i (θ), λ κ i (θ), r κ i (θ), the vectors (12) and their norms (13) attached to the matrix M κ (θ), defined in (7). Inequality (14) rewrites as

2/π ≤ λ κ 1 (θ)λ κ 2 (θ) ≤ 4/π. (15)
This explains the approximate vertical symmetry in the plot of λ κ 1 , λ κ 2 , Figure 7 (left, log scale). We denote by (f, g) ∈ [0, π/2] the unoriented angle between two vectors f, g ∈ R 2 \ {0}, considered up to a change of sign:

(f, g) := arccos | f, g | f g ∈ [0, π/2]. (16)
Let e ∈ R 2 \ {0}, let θ ∈ R, and let ϕ = (e, e θ ) where e θ := (cos θ, sin θ). Then

e 2 Mκ(θ) = e 2 (κ sin 2 ϕ + κ -1 cos 2 ϕ), hence κ -1 2 e ≤ e Mκ(θ) ≤ κ 1 2 e . (17)
Lemma 3.4 (Uniform bounds on Minkowski's vectors). For any κ ∈ [1, ∞[, and any θ ∈ R (omitting the argument θ for readability)

λ κ 1 ∈ [κ -1 2 , c 1 ] λ κ 2 ∈ [c 2 , κ 1 2 ], r κ 1 ∈ [1, c 1 κ 1 2 ] r κ 2 ∈ [1, κ],
where c 1 := 2/ √ π, c 2 := 2/π.

Proof. First line. Since λ κ 1 ≤ λ κ 2 by construction, and c 2 2 ≤ λ κ 1 λ κ 2 ≤ c 2 1 by (15), we obtain λ κ 1 ≤ c 1 and λ κ 2 ≥ c 2 . By minimality (13) of e κ 2 Mκ , and since at least one of the two independent vectors (1, 0) or (0, 1) is non-collinear with e κ 1 , one has λ κ 2 = e κ 2 Mκ ≤ max{ (1, 0) Mκ , (0, 1) Mκ } ≤ κ 1 2 max{ (1, 0) , (0, 1) } = κ 1 2 as announced. The lower bound for λ κ 1 is shown below. Second line. One has 1 ≤ e κ 1 =: r κ 1 since e κ 1 ∈ Z 2 \ {0}, and likewise r κ 2 ≥ 1. From (17, right) we obtain λ 1 κ ≥ κ -1 2 r κ 1 ≥ κ -1 2 which completes the first line, and

r κ 1 ≤ κ 1 2 λ κ 1 ≤ c 1 κ 1 2 and r κ 2 ≤ κ 1 2 λ κ 2 ≤ κ which concludes the proof.
Corollary 3.5. With the notations of Lemma 3.4, one has for κ ≥ c 2 1

1 2 R κ ≤ max{r κ 1 , r κ 2 } ≤ R κ , 1 2 S κ ≤ λ κ 2 ≤ S κ , (18)
R κ ∈ [1, 2κ], S κ ∈ [c 2 , 2κ 1 2 ].
Proof. By Proposition 3.3, one has max{r κ 1 , r κ 2 } ≤ R κ ≤ max{r κ 1 , r κ 2 , r κ 1 + r κ 2 } ≤ 2 max{r κ 1 , r κ 2 }, as announced in (18, left). Proceeding likewise with S κ , λ κ 1 , λ κ 2 , and observing that λ κ 1 ≤ λ κ 2 by construction (12), we obtain (18, right). The announced uniform bounds on R κ , S k , follow from those of Lemma 3.4.

We introduce a set Z, which collects one representative vector for each rational direction in the plane Z := {e ∈ Z 2 ; gcd(e) = 1, e (0, 0) lexicographically}.

For any e ∈ Z 2 \ {0}, the line eR intersects Z at a single point. Note also that {-1, 1} × Z → Z : (ε, e) → εe is clearly a bijection, where Z is defined in (1). Up to changing their sign, we assume from this point that Minkowski's vectors e i (M ), i ∈ {1, 2}, M ∈ S + 2 , belong to Z. We end this subsection by estimating a sum, over pairs of co-prime integers, which repeatedly appears in the sequel. In particular, there exists 0 < c ≤ C < ∞ such that cr ≤ Σ(r) ≤ Cr for all r ≥ 1.

Proof. Denote by p(r) the the relative density of points with co-prime coordinates in the ball of radius r: p(r) := 1 πr 2 #{e ∈ Z 2 ; 0 < e ≤ r, gcd(e) = 1}.

It is well known [START_REF] Hardy | An Introduction to the Theory of Numbers[END_REF] that p(r) → 1/ζ(2) = 6/π 2 as r → ∞. The announced result then follows from the observations

1 e = ˆr e 1 s 2 ds + 1 r , hence 2Σ(r) πr = 1 r ˆr 0 p(s)ds + p(r) = 2 ζ(2)
+ o(1).

Probabilistic estimate of the Minkowksi minima

We estimate the cumulative distribution of λ κ 1 , seen as a random variable of θ ∈ [0, π[, for a fixed κ ≥ 1, see Figure 9. The tail distribution of λ κ 2 is then inferred using Minkowski's second theorem (14). For that purpose we introduce the following angular sectors, illustrated on Figure 8

e 2 Mκ(θ) ≤ e 2 κ c e √ κ 2 + κ -1 = c 2 + e 2 κ -1 ≤ 2c 2 .
Corollary 3.8. For any κ ≥ 1, and any λ ≥ κ -1 2 , one has with λ

:= λ/ √ 2 e∈Z e ≤λ √ κ Φ κ e (λ ) ⊆ {θ ∈ [0, π[; λ κ 1 (θ) ≤ λ} ⊆ e∈Z e ≤λ √ κ Φ κ e (λ). (20)
Proof. The right implication in Lemma 3.7, with c = λ = λ/ √ 2, yields the first inclusion, since λ κ 1 (θ) ≤ e Mκ(θ) . The left implication in Lemma 3.7, with c = λ, yields the second inclusion.

We estimate in Proposition 3.10 the distribution function P (λ κ 1 ≤ λ). The main ingredient of the proof is that, under suitable assumptions, the union appearing in (20, left) is of disjoint sets. Probabilities considered here and in the following are for a uniformly distributed angle θ ∈ [0, π[. We recall that for any ϕ ∈]0, π/2] one has 2 π ϕ ≤ sin ϕ < ϕ.

Lemma 3.9. Let f, g be distinct elements of Z such that such that max{ f , g } ≤ c √ κ, where c = 1/π and κ ≥ 1. Then Φ κ f (1) ∩ Φ κ g (1) = ∅.

Proof. Since f, g are distinct elements of Z, they are not collinear. Hence det(f, g) = 0, and therefore | det(f, g)| ≥ 1 since this determinant is a positive integer. Thus Assuming that θ ∈ Φ κ f (1) ∩ Φ κ g (1), we obtain using (21) the contradiction

(f, g) > sin (f, g) = | det(f, g)| f g ≥ 1 f g . 1. 0.
(f, g) ≤ (f, e θ ) + (e θ , g) ≤ π 2 (sin (f, e θ ) + sin (e θ , g)) ≤ π 2 1 f √ κ + 1 g √ κ ≤ πc f g .
Proposition 3.10 (Cumulative distribution of λ κ 1 ). There exists constants C, c > 0 such that: for all κ ≥ 1, and all λ ∈

[Cκ -1 2 , c], one has cλ 2 ≤ P (λ κ 1 ≤ λ) ≤ Cλ 2 .
Proof. Let c 0 be the constant of Lemma 3.9, and let λ := λ/ √ 2. For λ ≤ c 0 , the left term of (20) is a disjoint union, by Lemma 3.9. Hence we obtain successively, using (22) in the second line, and introducing Σ of Proposition 3.6 in the third

e∈Z e ≤λ √ κ |Φ κ e (λ )| ≤ πP (λ κ 1 ≤ λ) ≤ e∈Z e ≤λ √ κ |Φ κ e (λ)|, e∈Z e ≤λ √ κ 2λ e √ κ ≤ πP (λ κ 1 ≤ λ) ≤ e∈Z e ≤λ √ κ πλ e √ κ , 2 λ √ κ Σ(λ √ κ) ≤ πP (λ κ 1 ≤ λ) ≤ π λ √ κ Σ(λ √ κ).
Assume in addition that λ √ κ ≥ 1. Then denoting by c 1 , C 1 the constants of Proposition 3.6

c 1 λ 2 = 2 λ √ κ c 1 λ √ κ ≤ πP (λ κ 1 ≤ λ) ≤ π λ √ κ C 1 λ √ κ = πC 1 λ 2 .
This concludes the proof, with c = min{c 0 √ 2, c 1 /π} and C = max{ √ 2, C 1 }.

Corollary 3.11 (Tail distribution of λ κ 2 ). There exists constants C, c > 0 such that: for any κ ≥ 1, and any λ ∈ [C, cκ

1 2 ] one has cλ -2 ≤ P (λ κ 2 ≥ λ) ≤ Cλ -2 .
Proof. The result follows from Proposition 3.10, and from (15) which implies for any λ > 0:

P (λ κ 1 ≤ 2/(λπ)) ≤ P (λ κ 2 ≥ λ) ≤ P (λ κ 1 ≤ 4/(λπ)).

Euclidean norm of Minkowski's shortest vector

We estimate the distribution of the euclidean norm r κ 1 (θ) of Minkowski's vector e κ 1 (θ) ∈ Z, which is defined in (12) as the shortest vector with respect to the anisotropic norm • Mκ(θ) . For that purpose, we need to precisely identify this vector, hence we introduce for each κ ≥ 1, e ∈ Z, the set

Θ κ e := {θ ∈ [0, π[; e = e κ 1 (θ)}, (23)
Our first two lemmas compare Θ κ e with the angular sector Φ κ e (c) introduced in (19). Our next proposition estimate the distribution of r κ 1 . Strictly speaking this is not needed for the proof of Theorem 1.19, but we provide for completeness the argument.

√ κ ≥ 1 e∈Z e ≤λ √ κ 2c 0 e √ κ ≤ πP (r κ 1 ≤ λ √ κ) ≤ e∈Z e ≤λ √ κ πC 0 e √ κ . 2c 0 c 1 λ ≤ 2c 0 Σ(λ √ κ) √ κ ≤ πP (r κ 1 ≤ λ √ κ) ≤ πC 0 Σ(λ √ κ) √ κ ≤ πC 0 C 1 λ.

Euclidean norm of Minkowski's second shortest vector

We estimate the tail distribution of the euclidean norm r κ

Observe that κ 2 -1 ≥ κ 2 /2, simplify powers of κ and e , and bound above the sum in the r.h.s. by twice its largest element. Then (26) is the consequence of the equivalent conditions.

| tan -1 (I)| + | tan -1 (-I)| ≥ 2 1 κ ρ - 1 ρ × (tan -1 ) ρ κ ≥ ρ 2κ
We used fact that ρ ≥ 2, so that ρ -1 ρ ≥ ρ/2. In addition ρ κ = 2c λ e √ κ ≤ c, since λ ≥ 2, e ≥ 1 and κ ≥ 1, hence (tan -1 ) (ρ/κ) ≥ 1/(1 + c 2 ) ≥ 1/2. This concludes the proof. Proof. Let c 0 , C 0 be the constants of Lemma 3.18, let κ, λ, e, satisfy the conditions of its statement, and let θ be in the set (24). Then

µ κ (θ) ≥ λ √ κ e ≥ λ 2 c 0 ≥ 1, thus µ κ (θ) - 1 2 ≥ 1 2 µ κ (θ) ≥ λ √ κ 2 e .
Therefore, using Corollary 3.17

A Discretizations of anisotropic PDEs based on stencils

We establish in this appendix the results announced in §1.1, on the consistency of numerical schemes based on D-diffusion and M -eikonal stencils, and their structure preservation properties. These results are fairly classical, but the proofs are provided for completeness. Results of similar nature can be found in [START_REF] Weickert | Theoretical foundations of anisotropic diffusion in image processing[END_REF] for diffusion with a focus on image processing, [START_REF] Ciarlet | Discrete maximum principle for finite-difference operators[END_REF] for static diffusion problems, or [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: theory and algorithms[END_REF][START_REF] Vladimirsky | Label-setting methods for Multimode Stochastic Shortest Path problems on graphs[END_REF] for eikonal equations. Note also the restriction to two dimensional discretizations on cartesian grids is not necessary, and that in both cases extensions to arbitrary dimensional discretizations on unstructured point sets can easily be formulated.

A.1 Diffusion stencils

We first establish that the Consistency axiom of D-diffusion stencils, is indeed equivalent to the consistency of two discretizations of ∇u(x) 2 D and Tr(∇ 2 u(x)) respectively. Thus the Consistency axiom implies properties Consistency' and Consistency". Conversely, observing that the quadratic form L ∈ R 2 → L, DL (resp. the linear form Q ∈ S 2 → Tr(DQ)) uniquely determines the matrix D ∈ S 2 , we obtain that Consistency' (resp. Consistency") implies Consistency, which concludes the proof.

Our next objective is to establish Proposition 1.4 on the stability of the explicit and implicit time steps for anisotropic diffusion. A definition and two intermediate results are required. In this paper, the term operator always refers to a continuous map from L ∞ (Z 2 ) to itself, and the maximum principle is understood in the sense of Definition 1.3. Definition A.1. An operator L preserves constants iff Lu = u for any identically constant u ∈ L ∞ (Z 2 ). An operator L is non-negative iff u ≥ 0 ⇒ Lu ≥ 0 for any u ∈ L ∞ (Z 2 ).

Lemma A.2. A linear operator obeys the maximum principle iff it is non-negative and preserves constants.

Proof. Denote by B this operator, and by u an arbitrary element of L ∞ (Z 2 ). Sups and Infs are taken over the set Z 2 .

Proof of implication ⇒. If u is constant, then from inf u ≤ Bu ≤ sup u we obtain Bu = u. If u is non-negative, then from Bu ≥ inf u ≥ 0 we obtain that Bu is non-negative.

Proof of implication ⇐. For any u ∈ L ∞ (Z 2 ), on obtains 0 ≤ B(u -inf u) = Bu -inf u, and likewise 0 ≤ B(sup u -u) = sup u -Bu, thus inf u ≤ Bu ≤ sup u as announced.

The next corollary is a states a simple necessary and sufficient condition for the maximum principle for linear and translation invariant operators.

  Definition 1.7 (Hopf-Lax operator). Let e 1 , • • • , e r ∈ Z obey the properties Ordering and Orientation of Definition 1.6. Then for any u : Z 2 → R and any x ∈ Z 2 we define Λu(x) := min 0≤i<r min t∈[0,1]

0Figure 2 :

 2 Figure 2: Voronoi cell Vor(M ), and unit ball {x ∈ R 2 ; x M ≤ 1}, for to two matrices M . Left: M = Id (non generic case). Right: matrix M with condition number κ = 1.4 and orientation θ = 0.3 (generic case), i.e. M := M κ (θ) from Definition 1.18. Voronoi region Vor(M ) shown in gray. Voronoi vectors shown as thick points, strict ones in black, others in gray. Lines (dashed for non-strict Voronoi vectors) mark the equality 2 g, M e = e 2 M ; Voronoi facets Vor(M ; e) are their intersections with the grayed Voronoi region Vor(M ).

Figure 3 :

 3 Figure 3: Unit circle {x ∈ R 2 ; x M = 1}, and collection of strict M -Voronoi vectors, for matrices M = M κ (θ) of varying condition number κ and orientation θ, see Definition 1.18. Left: κ = 5, θ ∈ [0, π/2]. Right: κ ∈ [1, 16], θ = π/3).

Figure 4 :

 4 Figure 4: Plots of R κ (θ) (left) and S κ (θ) (right), for κ ∈ {10, 50}, as a function of θ ∈ [0, π/2], see Definition 1.18. Logarithmic scale.

Figure 5 :

 5 Figure 5: Left: L p ([0, π[) norms of R κ and S κ , as a function of κ ∈ [1, 100] for different exponents p ∈ [1, ∞]. The behavior for p < 2 (dashed) and p > 2 is significantly different, see Theorem 1.19. Right: Tail distribution of R κ and S κ . Log-log scale.

  e∈V γ D Λ(e, e, e, e) + O(h 4 ), (10) with the four-linear tensor Λ = 1 2 ∇ 2 u(x) ⊗ ∇ 2 u(x) + ∇u(x) ⊗ ∇ 3 u(x). The second order O(h 2 ) contribution is thus bounded by the three term product, e∈Z γ D (e) Λ e 4 ≤ Λ e∈Z γ D (e) e 2 max e∈supp(γ D ) e 2 = Λ Tr(D)r(D) 2 . involving a contribution Λ of the PDE solution regularity, a contribution Tr(D) of the PDE data, and a contribution r(D) 2 from the scheme, where r(D) := max{ e ; e ∈ supp(γ D )} is the radius of the stencil support. Therefore, the difference between the elliptic energy E D and its discretization E h is bounded by

Figure 6 :

 6 Figure 6: Unit ball {e ∈ R 2 ; e M ≤ 1} and collection of strict M -Voronoi vectors; same graphic after a linear change of coordinates by M 1 2 . Left: diagonal matrix, Right: generic matrix. The angles of trigonometrically consecutive strict M -Voronoi vectors f, g are acute after the change of coordinates, as expected since f, M g ≥ 0 by Corollary 2.3.

  Corollary 2.5. Let (e 1 , • • • , e r ) be an M -eikonal stencil. Then there exists another M -eikonal stencil (e 1 , • • • , e r ) such that [e 1 , • • • , e r ] ⊆ [e 1 , • • • , e r ] and det(e i , e i+1 ) = 1 for all 0 ≤ i < r . Proof. Define (e 1 , • • • , e r ) by ordering trigonometrically the elements of [e 1 , • • • , e r ] ∩ Z. By construction, the points e 1 , • • • , e r are among e 1 , • • • , e r , hence for any 0 ≤ i < r there exists 0 ≤ j < r such that the vectors e j , e i , e i+1 , e j+1 are in trigonometric order. This implies e i , M e i+1 > 0 (Acuteness) and det(e i , e i+1 ) > 0 (Orientation) as required. By construction, the triangle [0, e i , e i+1 ] contains no point of Z aside from its vertices, hence | det(e i , e i+1 )| ≤ 1 by Lemma 2.4, thus det(e i , e i+1 ) = 1 as announced since this determinant is a positive integer. Proposition 2.6. Let M ∈ S + 2 , and let (e 1 , • • • , e r ) be an M -eikonal stencil. Then all strict M -Voronoi vectors belong to [e 1 , • • • , e r ].

Lemma 3. 2 . 2 M ≤ e 2 + e 1 2 M , and e 2 2 M ≤ e 1 -e 2 2 M

 22222 Let M ∈ S + 2 and let e 1 := e 1 (M ) and e 2 := e 2 (M ). Then 2| e 1 , M e 2 | ≤ e 1 2 .Proof. By construction e 2 . The result follows.Our next lemma describes the collection of strict M -Voronoi vectors, in terms of Minkowski's shortest vectors (12).Proposition 3.3. Let M ∈ S +2 and e i := e i (M ), i ∈ {1, 2}. The strict M -Voronoi vectors are:(i) if e 1 , M e 2 = 0 then {±e 1 , ±e 2 }.

Figure 7 :

 7 Figure 7: Minkowski minima λ κ 1 , λ κ 2 , and radii of Minkowski's vectors r κ 1 , r κ 2 . Logarithmic scale.

  (left): for each κ ≥ 1, e ∈ Z, α > 0 Φ κ e (α) := θ ∈ [0, π[; sin (e, e θ ) Lemma 3.7. Let κ ≥ 1, let e ∈ Z, let θ ∈ [0, π[ and let c > 0. We have the implications: e Mκ(θ) ≤ c ⇒ θ ∈ Φ κ e (c) and e ≤ c √ κ ⇒ e Mκ(θ) ≤ c √ 2.

Figure 8 :

 8 Figure 8: Left: representation of the sets Φ κ e (1/2), with e ∈ Z and e ≤ √ κ, defined in (19). Right: sets Θ κ e defined in (23) below. In this representation, we show in black the convex hull of each vector e and of the corresponding angular sector on the circle of radius κ.

( 21 )

 21 Denote by |E| the Lebesgue measure of a set E ⊆ R. Using (21) we obtain for any e ∈ Z, κ ≥ 1 and α ≤ 1, denoting σ := e (α)| = |[-arcsin σ, arcsin σ]| = 2 arcsin σ ≤ πα e √ κ .

Figure 9 :

 9 Figure 9: Cumulative distributions of λ κ 1 , r κ 1 , and tail distributions of λ κ 2 , r κ 2 . Log-log scales.

Lemma 3. 12 .

 12 Let κ ≥ 1, e ∈ Z, and C = 2/ √ π. Then Θ κ e ⊆ Φ κ e (C), and Θ κ e = ∅ ⇒ e ≤ Cκ 1 2 .Proof. Using successively (i) the definition of Θ κ e , (ii) the uniform bound on λ κ 1 (θ) obtained in Lemma 3.4, for which C = c 1 , and (iii) Lemma 3.7 (left implication) we obtain the implicationsθ ∈ Θ κ e ⇒ e Mκ(θ) = λ κ 1 (θ) ⇒ e Mκ(θ) ≤ C ⇒ θ ∈ Φ κ e (C) and e ≤ C √ κ .This establishes the announced inclusion and bound on e . Lemma 3.13. Let κ ≥ 1, e ∈ Z, and c < 1/(π √ 2). If e ≤ c √ κ then Θ κ e ⊇ Φ κ e (c). Proof. Let θ ∈ Φ κ e (c), and M := M κ (θ). By Lemma 3.7 (right implication) one has e M ≤ c √ 2. Let e ∈ Z \ {e}, and let c 0 := 1/π be the constant from Lemma 3.9. If e ≥ c 0 √ κ, then e M ≥ κ -1 2 e ≥ c 0 > c √ 2 ≥ e M . On the other hand, if e ≤ c 0 √ κ, then Φ κ e (1)∩Φ κ e (1) = ∅ by Lemma 3.9. Thus θ / ∈ Φ κ e (1) and therefore e M ≥ 1 > c √ 2 ≥ e M by Lemma 3.7 (left implication). We have shown that e M < e M for any e ∈ Z \ {e}. Thus e = e 1 (M ) = e κ 1 (θ), in other words θ ∈ Θ κ e as announced.

  Proposition 3.14 (Cumulative distribution of r κ 1 ). There exists C, c > 0 such that for any κ ≥ 1 and any λ ∈[κ -1 2 , c], one has cλ ≤ P (r κ 1 ≤ λ √ κ) ≤ Cλ.Proof. Note that for any fixed κ ≥ 1, the sets Θ κ e , e ∈ Z are pairwise disjoint. Denoting by C 0 and c 0 the constants of Lemmas 3.12 and 3.13 respectively, one thus has for λ ≤ c 0 for the first line, and denoting by c 1 , C 1 the constants of Proposition 3.6 for the second line, we conclude that when λ

Proof.

  Let κ, λ, e be as in the lemma statement. Let θ ∈ [0, π[, let M := M κ (θ), and let ϕ = (e θ , e). We define, and estimate using (17), the scalarµ := | e, M e ⊥ | e 2 M e 2 = e -2 (κ -1 κ ) cos ϕ sin ϕ κ sin 2 ϕ + 1 κ cos 2 ϕ = e -2 (κ 2 -1) κ 2 tan ϕ + cotan ϕ. tan ϕ + cotan ϕ.

  2 tan ϕ, cotan ϕ} ⇔ tan ϕ ∈ where we introduced the ratio ρ := √ κ 4λ e = 2c √ κ λ e ≥ 2. We assume that condition (27) holds and observe that since λ ≥ 2 sin ϕ ≤ tan ϕ ≤ 2c λ thus θ ∈ Φ κ e (c). By Lemma 3.13, and the assumption e ≤ c √ κ, we have Φ κ e (c) ⊆ Θ κ e . Thus θ ∈ Θ κ e , which implies e = e κ 1 (θ), and µ = µ κ (θ). Denoting by I the interval (27), the measure (24) is thus at least

  Corollary 3.19 (Tail distribution of r κ 2 . Lower bound). There exists 0 < c ≤ C < ∞ such that for any κ ≥ 1, and any λ ∈ [C, cκ 1 2 ], one has P (r κ 2 ≥ λ √ κ) ≥ cλ -2 .

  and since r κ 1 (θ) = e r κ 2 (θ) ≥ e (µ κ (θ) -1/2) by c 1 and Σ the constant and function from Proposition 3.6, we conclude assuming λ ≤ c 0 √ κ for the last inequalityP 0 λ -1 √ κ) ≥ c 2 0 c 1 λ -2 .

  Proof of proposition 1.2, on property Consistency. Consider a finitely supported γ : Z → R obeying property Symmetry of Definition 1.1, and define D := e∈Z γ(e)e ⊗ e. Consider also the linear and quadratic functions defined for all x ∈ R 2 byu 1 (x) := C + L, x , u 2 (x) := C + L, x + 1 2 x, Qx , where C ∈ R, L ∈ R 2 and Q ∈ S 2 .Then by elementary linear algebra we obtain e∈Z γ(e)(u 1 (x + e) -u 1 (x)) 2 = e∈Z γ(e) L, e 2 = e∈Z γ(e) L, e ⊗ e L = 2 L, DL . Similarly, using the symmetry of the weights, 2 e∈Z γ(e)(u 2 (x + e) -u 2 (x)) = e∈Z γ(e) (u(x + e) -2u(x) + u(x -e)) = e∈Z γ(e) e, Qe = e∈Z γ(e) Tr(e ⊗ e Q) = 2 Tr(DQ).

Note that the Perona-Malik model is often incorrectly referred to as "anisotropic diffusion", despite the fact that its diffusion tensors are proportional to the identity matrix. See[START_REF] Weickert | Theoretical foundations of anisotropic diffusion in image processing[END_REF] for a discussion on terminology.

(θ) of Minkowski's vector e κ 2 (θ), the second shortest vector with respect to the anisotropic norm • Mκ(θ) . The upper bound for the tail distribution of r κ 2 directly follows from the one obtained in Corollary

3.11 for λ κ 2 , as shown in the next lemma.

This work was partly supported by the ANR grant NS-LBR. ANR-13-JS01-0003-01

Lemma 3.15 (Tail distribution of r κ 2 . Upper bound). There exists C, c > 0 such that for any κ ≥ 1 and any λ ∈ [C, cκ Note that the estimate of the tail distribution P (r κ 2 ≥ µ) obtained in Lemma 3.15 is limited to µ ∈ [C √ κ, cκ], whereas by Lemma 3.4 the radius r κ 2 may take values in the larger interval [1, κ]. This estimate is nevertheless sufficient for our purposes; the regime in which it applies is shown in the grayed region in Figure 9.

In order to estimate r κ 2 from below, we tie it to r κ 1 , a strategy similar to the one used for the Minkowski minima λ κ 1 and λ κ 2 in §3.2. For any α ∈ R we abusively denote α+Z

Proof. Up to a linear change of coordinates, we may assume that M = Id. On the whole real line R, the minimizer of λ(l) := f -le 2 = f 2 -2l e, f + l 2 e 2 is l * := f, e / e 2 . Since λ(l * + δ) = λ(l * ) + δ 2 e 2 for any δ ∈ R, we find that k is the (an) element of (α + Z) closest to l * . Thus |k -l * | ≤ 1 2 , which concludes the proof.

Corollary 3.17. Let κ ≥ 1 and θ ∈ R. Let us denote e 1 := e κ 1 (θ), M := M κ (θ), and

Proof. Let f := e ⊥ 1 / e 1 2 , let e 2 := e κ 2 (θ), and let us note that | det(e 1 , e 2 )| = 1 = det(e 1 , f ). Up to replacing f with -f we may thus assume that det(e 1 , e 2 -f ) = 0. Therefore there exists α ∈ R such that e 2 = f -αe 1 . For any k ∈ R one has (f -ke 1 ∈ Z 2 iff k ∈ α + Z). By minimality (12) of e 2 M among elements of Z 2 non-collinear with e 1 , we thus have f -ke 1 M ≥ f -αe 1 M for any k ∈ α + Z. By Lemma 3.16

Finally, r κ 2 (θ) = e 2 = f -αe 1 ≥ |α| e 1 , since e 1 , f = 0, which concludes the proof.

Note that µ κ (θ) vanishes precisely when e κ 1 (θ) is an eigenvector of M κ (θ). We study in Lemma 3.18 the angular sectors where in contrast µ κ (θ) is large. This implies in Corollary 3.19 a lower bound on the tail distribution of r κ 2 , which is the counterpart of Lemma 3.15.

and let e ∈ Z be such that λ e ≤ c √ κ, where we denoted C = 2 and c = 1/8. Then

Moments of the Voronoi radii

We conclude in this section the proof of Theorem 1.19, we estimate the moments of the euclidean radius R κ (θ) and intrinsic radius S κ (θ) of the set of strict M κ (θ)-Voronoi vectors. As usual θ is regarded as a random variable uniformly distributed over θ ∈ [0, π[, and the the parameter κ ≥ 1 is fixed. Our first step is to estimate the tail distributions of R κ and S κ .

Proposition 3.20. There exists C, c > 0 such that for all κ ≥ 1, and all λ ∈ [C, cκ

1 2 ], one has

Proof. For readability, we omit the argument θ ∈ [0, π[ of the functions involved. The estimate for S κ immediately follows from the equivalence S κ /2 ≤ λ κ 2 ≤ S κ , see Corollary 3.5, and from the estimate of the tail distribution of λ κ 2 obtained in Corollary 3.11. Regarding R κ we have

] by Lemma 3.4. Hence, assuming C > 2c 1 , we obtain for λ ≥ C the implications

The announced estimate for P (R κ ≥ λκ 1 2 ) then follows from the estimates of the tail distribution of r κ 2 obtained in Lemma 3.15 and Corollary 3.19.

Proposition 3.20 implies in particular that P (S κ ≥ cκ

In view of the uniform bounds S κ ∈ [c 2 , 2κ

1 2 ] obtained in Corollary 3.5, and of the estimates of Proposition 3.20 which hold over a range [C, cκ 1 2 ], we cut the integration range ]0, ∞[ of (29, right) into four sub-intervals. Assuming without loss of generality that c 2 ≤ C, c ≤ 2, and κ ≥ (C/c) 2

Interval I 4 does contribute to the integral (29). Since I 1 is bounded independently of κ, its contribution to (29) also is. Since the tail distribution P (S κ ≥ λ) is decreasing in λ, the contribution to (29) of

Finally, and most importantly, the contribution of

for sufficiently large κ, with equivalence constants (9) depending only on p. This implies the estimate on S κ L p ([0,π[) announced in Theorem 1.19. The case of R κ L p ([0,π[) is analogous.

Corollary A.3 (Stability criterion).

Consider η ∈ L 1 (Z 2 ), and the operator defined by A η u(x) := e∈Z 2 η(e)u(x + e). Then A η obeys the maximum principle iff η ≥ 0 and e∈Z 2 η(e) = 1.

Proof. The non-negativity of A η is clearly equivalent to η ≥ 0, and the fact that it preserves constants to e∈Z 2 η(e) = 1.

Finally we establish Proposition 1.4. For that purpose we separate the claims on the explicit diffusion time step Id +δL γ and the implicit one (Id -δL γ ) -1 .

Proof of Proposition 1.4, on explicit diffusion. Let γ : Z → R be finitely supported, and let δ > 0. Observing that

and applying Lemma A.3 we obtain the equivalence of the two points of Proposition 1.4, as well as the optimality of the CFL condition 0 ≤ δ e∈Z γ(e) ≤ 1.

Proof of Proposition 1.4, on implicit diffusion. Let γ : Z → R be finitely supported and nonnegative, and let δ > 0. Consider the operator defined for u :

where

.

As a result, the operator associated to the implicit time step

) and non-negative. In addition the operator (Id -δL γ ) -1 preserves constants, since its inverse does, hence it obeys the maximum principle by Lemma A.2.

A.2 Eikonal stencils

We provide for completeness the proof Proposition 1.9, stating that the causality property of the Hopf-Lax update operator is equivalent to a geometrical acuteness property of the stencils. We mostly follow the steps of the appendix of [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: theory and algorithms[END_REF].

Lemma A.4. Let N ∈ S + 2 . Then the following are equivalent:

(i) N has non-negative entries.

(ii) N ξ has positive entries whenever ξ ∈ R 2 has positive entries.

Proof. Proof that (i) ⇒ (ii). Since N has non-negative coefficients, the product N ξ also has.

Observing that N has positive diagonal entries, since it is positive definite, we obtain that N ξ has positive entries. Proof that (ii) ⇒ (i). Let

where the coefficient 1 is in the i-th position, and letting ε → 0, we obtain that the i-th column of N is non-negative. The result follows.

Lemma A.5. Let N ∈ S + 2 , and for all u ∈ R 2 let λ N (u) := min

where Ξ := {ξ ∈ R 2 + ; ξ, 1 = 1} and 1 := (1, 1) ∈ R 2 . Then the following are equivalent: (i) the matrix N has non-negative entries.

(ii) for any u ∈ R 2 and any i ∈ {1, 2}, such that the minimizer in (30) obeys ξ i > 0, one has λ N (u) > u i .

Proof. Note that the quantity λ N (u) is defined as the minimum (30) of a strictly convex functional on a compact interval, hence the minimizer exists and is unique. Let u ∈ R 2 , and let us assume that the minimizer ξ in (30) has positive entries. By the Kuhn-Tucker optimality conditions, there exists λ ∈ R such that

Taking the scalar product with ξ we obtain

Thus combining the last two equations

Proof that (i) ⇒ (ii). The announced result is clear if the minimizer ξ in (30) equals (1, 0) or (0, 1), hence as above we can assume that ξ has positive entries. By assumption N has nonnegative entries, hence N ξ has positive entries by Lemma A.4, thus also λ N (u)1 -u by (32), which concludes the proof of this implication.

Proof that not (i) ⇒ not (ii). Assume that N has a negative entry. By Lemma A.4 there exists ξ with positive entries such that N ξ has at least one non-positive entry. Up to multiplication by a positive constant, we may assume that ξ ∈ Ξ. Let u := -N ξ/ ξ N . Since (30) is the optimization of a strictly convex functional over a convex set, its minimizer is characterized by the first order Kuhn-Tucker condition (31). Hence λ = λ N (u) = 0 and as announced we obtain not (ii).

Proof of Proposition 1.9, Acuteness is equivalent to Causality. By construction, one has

where N i := e i , M e i e i , M e i+1 e i , M e i+1 e i+1 , M e i+1 , and y i := x + e i . The matrices N i , 0 ≤ i < r, are positive definite by construction, and the Acuteness geometric property is equivalent to the non-negativity of their coefficients. The first implication Acuteness ⇒ Causality, is thus a direct consequence of Lemma A.5. Proof that not Acuteness ⇒ not Causality. Without loss of generality, we may assume that e 0 , M e 1 < 0, hence that N 0 has one negative entry. Choose u(y 0 ) and u(y 1 ) so as to contradict point (ii) of Lemma A.5, in other words such that defining λ := λ N 0 (u(y 0 ), u(y 1 )) one has λ < max{u(y 0 ), u(y 1 )} and the associated minimizer ξ 0 = (1 -t 0 , t 0 ) for (30) has positive entries. Note also that λ < min{ e 0 M + u(y 0 ), e 0 M + u(y 0 )}, since these values correspond to t = 0 and t = 1.

For any 0 ≤ i < r, one has λ N i (u(y i ), u(y i+1 )) ≥ min{u(y i ), u(y i+1 )}, and λ N i (u(y i ), u(y i+1 )) → e i M + u(y i ) as u(y i+1 ) → ∞. Hence choosing u(y i ) sufficiently large, for all i / ∈ {0, 1}, we obtain Λu(x) = λ, this minimum being attained solely for i = 0 and t = t 0 . This contradicts the causality property, and concludes the proof.