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

We consider discretizations of anisotropic diﬀusion and of the anisotropic eikonal equation, on two dimensional cartesian grids, which preserve their structural properties: the maximum principle for diﬀusion, and causality for the eikonal equation. These two PDEs embed geometric information, in the form of a ﬁeld of diﬀusion 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 ﬁxed 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 [Cia70], 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 [Tsi95,SV01], 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, AM11,Mir14] 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 [Dro14], or stochastic control [BOZ04], but is also fundamental in image processing [Wei98]. The Anisotropic eikonal equation is naturally required for trajectory planning, but is also relevant for seismic imaging [SV03], and extensively used in medical image segmentation [BC10]. 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 [SV03,Koc95]. The objective of this paper is to provide a complete

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 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 [Cia70] on the discrete maximum principle, [Wei98] on anisotropic diffusion in image processing, and [SV03] 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.
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.
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 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 Proposition 1.4. Let γ : Z → R be finitely supported. Then the following are equivalent: • The weights γ are non-negative.
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 [Cia70] 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 This PDE is known to admit a unique viscosity solution [CL83], 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 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 .
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.
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 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 [FF02,BR06] for the eikonal equation take the form of a fixed point problem: find u : 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) 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 [BR06], 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 [Tsi95], which is a variant of Dijkstra's algorithm on graphs, see Figure 1. This method 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 [BR06,?], 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 [SV03,Vla08,Mir13]. Proposition 1.9 (Causality is equivalent to Acuteness [SV03,Vla08]). 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  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.
• 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 [FM13] for diffusion and [Mir14] for eikonal PDEs. We refer to these works for convergence results, numerical experiments and comparisons with other approaches. The stencil constructions of [FM13,Mir14] 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 [FM13], and the M -eikonal stencil of Theorem 1.13 is the one presented in [Mir14].
Our D-diffusion stencils are also equivalent to the (two dimensional only) construction of [BOZ04]. Non-equivalent constructions can be found in [AM11] for M -eikonal stencils, and [Wei98] for D-diffusion stencils. They are, inevitably, less local and accurate, as exposed by the comparisons in [FM13,Mir14].

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.
We denote by V κ (θ) the collection of all strict M κ (θ)-Voronoi vectors, and define 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 anisotropic geometry embedded in the PDE. Several plots of R κ (θ) and S κ (θ), as a function of θ ∈ [0, π/2[ and for a fixed κ, are shown in Figure 4. The dependence of R κ (θ) and S κ (θ) 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 [Mir14] with an heuristic analysis of the accuracy of several discretizations of the eikonal PDE.
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: Theorem 1.19. For any p ∈ [1, ∞] one has uniformly in κ (L p norms are on the interval [0, π[) 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 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 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, 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 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 [Wei96] 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 [FM13]. Note that our results are irrelevant for non-linear isotropic diffusion, such as the Perona-Malik model 1 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 [Pet03] 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 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: This strategy is applied in [BOZ04] to the HJB equation of stochastic control. See also [BCM15] 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 [CCM14].
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.

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 which shows that Vor(M ; λe) = ∅, hence that λe is not a Voronoi vector. Proof. Consider x ∈ Vor(M ; f ) ∩ Vor(M ; g). Then 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. 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.

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 [Mir14], 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 [Sel74]. Obtuse superbases are an important tool in lattice classification [CS92]. They are also at the foundation of the M -eikonal stencil of [Mir14], described in Proposition 2.8 (ii), and of the D-diffusion stencil of [FM13], defined in Corollary 2.11. The proof of the following proposition is left as an easy exercise.
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.
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.

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.
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 Proof. For any i, j ∈ {0, 1, 2}, one has denoting by δ ij the Kronecker symbol Therefore denoting byD the r.h.s. of (11) we obtain for any 0 ≤ j ≤ 2 generates the three dimensional space S 2 , hence it is a basis.
The next corollary describes the D-diffusion stencil used in [FM13], and shows that it is supported on strict M -Voronoi vectors, with M := D −1 , as announced in Theorem 1.12.
Proof that γ is supported on the set of strict M -Voronoi vectors. For any f, g ∈ Z 2 one has 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 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 [FM13,Mir14] are distinct, and in particular they do not have the same number of vertices. We denote S 1 := {z ∈ R 2 ; z = 1}.
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. 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 [Ngu04], 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 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 [Ngu04], 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 [Ngu04]. 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 [Ngu04] for extensions of this result to higher dimension.
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.
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 ϕ < ϕ.
Denote by |E| the Lebesgue measure of a set E ⊆ R. Using (21) we obtain for any e ∈ Z, κ ≥ 1 and α ≤ 1, denoting σ : Lemma 3.9. Let f, g be distinct elements of Z such that such that max{ f , g } ≤ c √ κ, where c = 1/π and κ ≥ 1.
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 Assume in addition that λ √ κ ≥ 1. Then denoting by c 1 , C 1 the constants of Proposition 3.6 This concludes the proof, with c = min{c 0 √ 2, c 1 /π} and C = max{ √ 2, C 1 }.
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 This establishes the announced inclusion and bound on e .
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.
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 |Φ κ e (C 0 )| Using (22) 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 λ

Euclidean norm of Minkowski's second shortest vector
We estimate the tail distribution of the euclidean norm r κ 2 (θ) 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.

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 implies in particular that P (S κ ≥ cκ 1 2 ) > 0 for all κ ≥ (C/c) 2 , but we also know that S κ ≤ 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 I 3 = [cκ Finally, and most importantly, the contribution of I 2 is 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. 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.
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 by u 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 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.
The next corollary is a states a simple necessary and sufficient condition for the maximum principle for linear and translation invariant operators.
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 λ1 = N ξ/ ξ N + u.
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 Λu(x) = min 0≤i<r λ N i (u(y i ), u(y i+1 )), 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.