The Infinity Laplacian eigenvalue problem: reformulation and a numerical scheme

In this work, we present an alternative formulation of the higher eigenvalue problem associated to the infinity Laplacian, which opens the door for numerical approximation of eigenfunctions. A rigorous analysis is performed to show the equivalence of the new formulation to the traditional one. Subsequently, we present consistent monotone schemes to approximate infinity ground states and higher eigenfunctions on grids. We prove that our method converges (up to a subsequence) to a viscosity solution of the eigenvalue problem, and perform numerical experiments which investigate theoretical conjectures and compute eigenfunctions on a variety of different domains.


Introduction
The infinity Laplacian equation was introduced by G. Aronsson in [56] and has been extensively studied in the following years.It can be categorized as a nonlinear degenerate elliptic partial differential equation (PDE), which has interesting connections to Lipschitz extensions [43,44] and also to probabilistic games [33].For important contributions to the analysis of the infinity Laplacian we refer to [42,44,56].Regarding the uniqueness of Lipschitz extensions and the theory of absolute minimizers we refer the interested reader to the work of Jensen [52], and further works in [42,43,49].A numerical approximation of the infinity Laplacian equation is investigated by Oberman in [41], where he introduced a convergent finite difference scheme.In the more general context of finite weighted graphs discretizations of the infinity Laplacian operator have been studied and applied for data processing tasks by Elmoataz et al. in [17,22,25].Continuum limits of infinity Laplacian type equations, i.e., convergence of such discretizations as the graph approximates a continuum domain, have been proven using different techniques such as viscosity solutions [17], Gamma-convergence [4], and comparison principles [2,7].
In contrast, numerical methods for the eigenvalue problem of the infinity Laplacian have received much less attention so far.As we will recap in Section 2 below, the infinity Laplacian and its eigenfunctions arise as limits of p-Laplacian operators and their eigenfunctions.In the last decade nonlinear eigenvalue problems of the p-Laplacian operator have gained increasing attention [38,55] and were used for signal processing applications [15,16,21].Horak discussed numerical approximations for the two smallest eigenvalues of the p-Laplacian operator for different values of 1 < p ≤ 10 in [32], see also the work [23] of the first author of this paper.For large values of p, however, it turns out to be difficult to compute eigenvalues and corresponding eigenfunctions of the p-Laplacian due to stiffness of the discretized systems.On the other hand, the eigenvalue problem for the infinity Laplacian operator has been analytically studied for example by Juutinen, Lindqvist, Kawohl, Manfredi, and Saksman in [36,40,45,46,47].In [24] the authors investigate a inverse iteration to solve the p-Laplacian eigenvalue problem and, as a byproduct, prove convergence to infinity Laplacian eigenfunctions as p is sent to ∞ along the iteration.To the best of our knowledge a direct numerical approximation of eigenfunctions of the infinity Laplacian operator has not been investigated so far.
Let us mention that there is a large body of literature which deals methods to solve abstract nonlinear eigenvalue problems related to the minimization of a Rayleigh quotient.Some methods that are worth mentioning in this context are gradient flows of the Rayleigh quotient like [13,19], or iterative methods like for instance nonlinear power methods [11].Comprehensive reviews on the topic which also relate continuous flows and discrete iterations are [6,12].However, as we will see all of these methods are not applicable for finding eigenfunctions of the infinity Laplace operator since such functions are special minimizers of a Rayleigh quotient with infinitely many other minimizers.Furthermore, these flows and iterations cannot deal with higher eigenfunctions than ground states.
Let us now introduce the infinity Laplacian eigenvalue problem that we study in this paper.For this, let Ω be an open, bounded domain in R d .We consider the following Dirichlet eigenvalue problem of the infinity Laplacian operator as studied in [40].One looks for a function u ∈ W 1,∞ 0 (Ω) which is a viscosity solution of where u > 0, where u < 0. (1.1) Here, Λ > 0 denotes a corresponding eigenvalue of the infinity Laplacian ∆ ∞ .Positive solutions are referred to as ground states and fulfill the simpler equation Both (1.1) and (1.2) turn out to be rather challenging from a numerical point of view.One difficulty arises from the case distinction in (1.1).These different cases are based on the unknown sign of the solution itself, and thus one is not able to implement a numerical approximation scheme directly.Furthermore, both equations are non-smooth and therefore require a careful monotone discretization.Finally, another difficulty, which is a problem for ground states and higher eigenfunctions, is that solutions are non-unique.One source of non-uniqueness, namely the scaling invariance of (1.1) and (1.2), can be eliminated by searching for normalized solutions, satisfying, e.g., ∥u∥ ∞ = 1.However, the second source of non-uniqueness is intrinsic for the infinity Laplacian eigenvalue problem [28] unless on restricts oneself to very special stadium-like domains [36].
The main contributions of this work are the following: First, we give a reformulation of (1.1) as one equation which avoids the distinction of cases.Second, we define consistent monotone schemes on grids for approximating solutions of (1.1) and (1.2), and prove subsequential convergence to a viscosity solution.Third, we perform extensive numerical experiments and compute various solutions of the infinity Laplacian eigenvalue problem on different domains.The proposed numerical scheme allows to investigate open theoretical conjectures, e.g., we numerically confirm that the infinity ground state and the so-called infinity harmonic potential disagree on a square domain, proved recently in [1,5].
The structure of this paper is as follows: Section 2 recalls the mathematical background of the p-and infinity Laplacian operators and their respective eigenvalues and eigenfunctions.In Section 3 we propose the alternative formulation of the infinity Laplacian eigenfunction problem (1.1) and prove their equivalence.Based on the reformulation, we define consistent monotone schemes for approximating eigenfunctions in Section 4. In Section 5 we show numerical results using the proposed approximations.

Mathematical background
To make this paper more self-contained, we begin by recalling the concept of viscosity solutions in Section 2.1.This is the suitable solution concept for both the infinity Laplacian equation and the eigenvalue problem (1.1).Furthermore, we recap properties of the infinity Laplacian equation, which is a substantial part of the eigenvalue problem, in Section 2.2.Finally, we summarize the analytic relationship of the p-and infinity Laplacian operators and discuss properties of their respective eigenvalues and eigenfunction in Section 2.3.

Viscosity solutions
We focus on PDEs of the following general form for a real-valued function u : Ω → R, F : Ω × R × R n × S n → R, and S n is the space of real, symmetric n × n-matrices.We further assume that F is degenerate elliptic, meaning for all u ∈ R and p ∈ R n .By N ≤ M in (2.2) we denote that the matrix M − N is positive semidefinite.Any equation fulfilling these properties is called degenerate elliptic.For a comprehensive overview on the theory of viscosity solutions we refer the interested reader to the seminal paper of Crandall, Ishii, and Lions in [53].
Definition 2.1.Any upper (respectively lower) semi-continuous function u : Ω → R is called a viscosity subsolution (respectively supersolution) of (2.1) if u(x) ≤ 0 (respectively u(x) ≥ 0) for all x ∈ ∂Ω and for all ϕ ∈ C 2 (Ω) and all x ∈ Ω such that u − ϕ has a local maximum (respectively minimum) at x, we have A continuous function u : Ω → R is said to be a viscosity solution if it is both a viscosity sub-and supersolution of (2.1).
Example 2.2.In the following, we consider the Eikonal equation on the interval Ω = (−1, 1) It is clear that there is no classical C 1 (Ω) solution to this problem.However, one can verify that there exists a unique solution in the viscosity sense given by u(x) = 1 − |x| for x ∈ [−1, 1].Any C 1 -function ϕ touching u from above in x = 0 has a slope |ϕ ′ (0)| ≤ 1 and obviously there is exists no such function touching u from below in x = 0.

The infinity Laplacian equation
The infinity Laplacian operator is defined as follows: Sometimes the operator in (2.3) is normalized by 1 |∇u| 2 , e.g., cf.[34], but this normalization does not change the equations we are considering here.A function u is said to be infinity harmonic if it solves the homogeneous infinity Laplacian equation in the viscosity sense.This equation can be derived as the limit of a sequence of p-Laplacian equations ∆ p u = div(|∇u| p−2 ∇u) = 0 under certain boundary conditions for p → ∞.
The infinity Laplacian equation is related to the absolute minimal Lipschitz extension (AMLE) problem [42,43,52,56].In this setting one searches for a continuous real-valued function which has the smallest possible Lipschitz constant in every open set whose closure is compactly contained in Ω.This interpretation has some advantages as it directly leads to numerical approximation schemes for solutions of the infinity Laplacian equation (2.4) for all open sets Ω ′ ⊂ Ω and all v such that u − v ∈ W 1,∞ 0 (Ω ′ ).The relationship between an AMLE and the infinity Laplacian is stated in Theorem 2.3 (Corollary 3.14 in [52]).A function u ∈ Lip(Ω) is an AMLE of a Lipschitz function g : ∂Ω → R if and only if u is a viscosity solution of the infinity Laplacian equation with u = g on ∂Ω.
Clearly, one can exploit this relationship to numerically solve the infinity Laplacian equation, i.e., to construct absolutely minimal Lipschitz extensions in a discrete setting [41].
Let us finally remark that infinity harmonic functions might not be C 2 differentiable in general.A well-known example from [56] is given by u(x, y) = |x| which is a C 1,1/3 infinity harmonic function.To the best of our knowledge it is still an open problem whether all viscosity solutions of the infinity Laplacian equation are in C 1 .Evans and coauthors proved C 1,α regularity of viscosity solutions for the case d = 2 in [35] and differentiability in general dimension in [31].

Eigenfunctions of the p-Laplacian and their limit
The eigenvalues and eigenfunctions of the infinity Laplacian can be constructed as the limit of respective eigenvalues and eigenfunctions of the p-Laplacian operator, see [45,46,47] for details.For this reason we shortly recall the definition of eigenvalues of the p-Laplacian in the following and also refer the interested reader to [38,55].
For 1 ≤ p < ∞ the first (smallest) eigenvalue of the p-Laplacian operator has a variational form and is given by the Rayleigh quotient for which the minimization is performed over all non-zero functions in the Sobolev space W 1,p 0 (Ω).Any minimizer of (2.5) has to satisfy the following Euler-Lagrange equation which has to be interpreted in the usual weak sense.It is well known that solutions of this equation, referred to as eigenfunctions of the p-Laplacian, are unique modulo global scaling.A particularly elegant proof of this result can be found in [37].Higher eigenfunctions are solutions of (2.6) with λ 1 (p) replaced by λ > λ 1 (p).In particular, according to [40] the second p-eigenvalue λ 2 (p) can be defined as λ 2 (p) = min {λ ∈ R : λ is a p-eigenvalue and λ > λ 1 }.
Analogously to (2.5), the first eigenvalue of the infinity Laplacian, denoted by Λ 1 , is given by where ∥φ∥ ∞ = ess sup x∈Ω |φ(x)|.It is easy to see (e.g., cf.[46]) that the Euclidean distance function d(x) = dist(x, ∂Ω) solves the minimization problem (2.7).However, solutions to the minimization problem (2.7) in W 1,∞ 0 (Ω) are in general not unique, see [8].In fact, unless the domain is a very special stadium-like domain infinitely many minimizers can be constructed.This non-uniqueness requires a more careful definition of the infinity Laplacian eigenvalue problem and correspondingly numerical methods that are based on Rayleigh quotient minimization (as discussed in the introduction) are not applicable in this context.
A first eigenfunction of the infinity Laplacian operator can be obtained through the limit of the p-Laplacian equation (2.6) for p → ∞.The limit of these equations as p → ∞ is found to be All this was proved in [46] and we subsume their results in where Moreover, any positive solution u of (2.8) realizes the minimum in (2.7).Such a function u can be constructed as a cluster point for p → ∞ of a properly normalized sequence of first eigenfunctions of the p-Laplacian operator.Furthermore, where λ 1 (p) denotes the first eigenvalue of the p-Laplacian operator given by (2.5).
An important distinction has to be made between those solutions of (2.8) which are a limit of p-Laplacian ground states and those which are not.Definition 2.5 (Variational ground states).A solution of (2.8) which is a cluster point of normalized solutions of (2.6) is called variational ground state.All other solutions are called non-variational.
Remark 1 ((Non-)uniqueness).As most eigenvalue problems, also (2.8) is invariant under scalar multiplication, meaning that if u ∈ W 1,∞ 0 (Ω)) is a solution then so is cu for any c ∈ R.However, not even normalized solutions to (2.8) are unique: In [28] Hynd, Smart, and Yu have shown the nonuniqueness of infinity ground states for a dumbbell domain.However, the ground state which was constructed there is non-variational.Yu in [36] proved that on stadium-like domains (as for instance the ball) ground states are unique up to scaling and coincide with the distance function of the domain.Whether uniqueness holds for general convex domains or variational ground states, remains an open problem.
Theorem 2.4 states that the first eigenvalue can be interpreted geometrically, i.e., Λ 1 is the reciprocal of the radius of the largest ball that fits inside the domain Ω.In general, Λ 1 cannot be detected in regions where the solution is smooth, i.e., the term |∇u(x 0 )| − Λ 1 u(x 0 ) in (2.8) is not active.According to [36] It is known that also the second eigenvalue has a geometric characterization.According to [40] the second eigenvalue of the infinity Laplacian is given by (2.9) where r 2 = sup{r > 0 : there are disjoint balls B 1 , B 2 ⊂ Ω with radius r}.Furthermore, one has Theorem 2.6 (Theorem 4.1 in [40]).Let λ 2 (p) be the second p-eigenvalue in Ω.Then it holds that and Λ 2 ∈> 0 is the second eigenvalue of the infinity Laplacian.
According to [40] higher eigenfunctions of the infinity Laplacian operator can be obtained as a viscosity solution and Λ denotes the corresponding eigenvalue.Since the sign of the solution is unknown a-priori, this is a free boundary problem and hence hard to solve numerically.This is our motivation for reformulating this eigenvalue problem in Section 3. Remark 2. The equation of the first eigenfunction in (2.8) can also be expressed through (2.10) since the first eigenfunction does not change sign.

Reformulation of the infinity Laplacian eigenvalue problem
In the following, we present an equivalent formulation of the higher infinity Laplacian eigenvalue problem, which allows us to avoid the distinction of cases in (2.10).To this end, we introduce the function H Λ : R × R × S n → R, defined as and consider the associated problem of finding a viscosity solution to the equation The following is our main theorem and states that the formulations through F Λ and H Λ are equivalent.
Theorem 3.1 (Equivalent formulation of the eigenvalue problem).
Proof.Assume F Λ (u, ∇u, D 2 u) = 0 in the viscosity sense.We need to make a case distinction on the sign of the solution u to prove that u is a subsolution of H λ (u, ∇u, D 2 u) ≤ 0. Showing that it is also a supersolution works in an analogous way.
Case 1.1 Let φ be a C 2 function touching u from above in x such that u(x) > 0. Then we have min and hence and hence Combining these two we infer Hence, we have shown that u is a subsolution of H Λ (u, ∇u, D 2 u) ≤ 0.
Case 1.2 Let now φ be a C 2 function touching u from above in x such that u(x) < 0. Then by assumption we have max( Case 1.3 Let φ be a C 2 function touching u from above in x such that u(x) = 0 meaning φ(x) = 0. Then we have −∆ ∞ φ(x) ≤ 0 which implies the subsolution property Now we prove the converse statement and assume that H Λ (u, ∇u, D 2 u) = 0 in the viscosity sense.Again, we consider the different possible signs of u and shall prove that u is a subsolution of F Λ (u, ∇u, D 2 u) ≤ 0. Analogously, one shows that u is a supersolution, as well.
Case 2.1 Let φ be a C 2 function touching u from above in x such that u(x) > 0. Then it holds that H Λ (φ(x), ∇φ(x), D 2 φ(x)) ≤ 0 and we must show that On the other hand, if , we obtain the subsolution property In this case we get and hence from above we see that This implies that |∇φ(x 0 )| = 0 and hence also −∆ ∞ φ(x 0 ) = 0. Thus, in both cases −∆ ∞ φ(x 0 ) ≤ 0 such that u is a viscosity subsolution.
For completeness we also prove that the equation H Λ (u, ∇u, D 2 u) = 0 is degenerate elliptic.Proof.To prove degenerate ellipticity it suffices to argue that for a, b ∈ R with a ≥ b, the function t → f (t) := min(a, −t) + max(b, −t) + t is non-increasing.To see this, a simple case distinction shows that one can express f as which is obviously a non-increasing function.

Numerical method
In this section we propose methods to approximate eigenfunctions of the infinity Laplacian.First, we recall the concept of monotone schemes in Section 4.1 as these are needed to construct numerical schemes which approximate eigenfunctions of the infinity Laplacian on general unstructured grids.Then, we sketch the approximation of the distance function and the first infinity Laplacian eigenvalue in Section 4.2.Finally, our main contribution in this section is that we define consistent monotone schemes to approximate ground states and higher eigenfunctions of the infinity Laplacian in Section 4.3 and Section 4.4, respectively.

Monotone schemes and convergence without comparison principle
In order to numerically compute approximate viscosity solutions to the abstract degenerate elliptic equation (2.1), which in particular allows us to solve the infinity eigenvalue problems, we make use of monotone schemes and follow the description by Oberman in [39].We first define an unstructured grid on the domain Ω as a graph consisting of a set of vertices V = V inn ∪ V bdry where V inn := {x i ∈ Ω : i = 1, . . ., M } for M ∈ N are inner vertices and V bdry := {x i ∈ ∂Ω : i = M +1, . . ., M +N } for N ∈ N are boundary vertices.The number of total vertices is abbreviated by K := M + N ∈ N. To each point x i ∈ V we associate a list of global neighbors indices given by N i = {i 1 , . . ., i ki } ⊂ {1, . . ., K} for some k i ∈ N. We assume the symmetry condition j ∈ N i if and only if i ∈ N j .
A grid function F : V inn → R is a real-valued function defined on V inn which is based on values u i = u(x i ) of a function u : Ω → R and is given by: where the functions Fi on the right are possibly different for every grid point x i ∈ V inn .Then, a discrete solution of (2.1) on the unstructured grid introduced above is a grid function û : V → R which satisfies Any such function can be trivially turned into a function defined on Ω using a closest point projection which allows us to define u K := û • π K : Ω → R. Note that the approximation quality of the grid depends on the choice of the neighborhood and, in particular, it comes with the intrinsic errors ) ) where S n denotes the unit sphere in R n .The first error term describes the grid resolution in the interior of Ω, the second one the resolution of the boundary, and the third error is the directional resolution of inner vertices.
Remark 3.For convergence of our numerical scheme we shall require that dx M , dx N , dθ M → 0.
We emphasize that for standard finite difference schemes (for instance, five point stencils in two dimensions) do not satisfy that dθ M → 0. For discretizing the linear Laplace equation, which just depends on second derivatives in the coordinate directions, this is not necessary.However, for the infinity of even p-Laplacians, which depend on second derivatives in the direction of the gradient, it is required that the direction resolution goes to zero, cf.[2,4,7,9,17,29,39,41].In order for a grid to satisfy dθ M → 0 the number of points in the computational stencil N i has to tend to infinity.For quantitative rates how fast this growth has to be for the infinity Laplace equation we refer to [2,7].• consistent with respect to F in (2.1) if for every x ∈ Ω and ϕ ∈ C 2 (Ω) it holds The following proposition gives a convergence criterium for degenerate elliptic schemes.It is a straightforward generalization of the classical Barles-Souganidis theorem [54] to the case where the equation does not admit a comparison principle or unique solutions, see also [20].The proposition shows that any continuous cluster point of a sequence of solutions of (4.1) solves (2.1).Proposition 4.2.Assume that the grid function F is degenerate elliptic and consistent according to Definition 4.1, and let ûK : and a subsequence of u K converges uniformly to a continuous function u : Ω → R then u is a viscosity solution of (2.1).
Proof.We only show the subsolution property since showing the supersolution property works analogously.
Let x ∈ Ω and ϕ ∈ C 2 (Ω) such that u − ϕ has a strict maximum at x.Then, since dx M → 0, there exist sequences M k , K k → ∞ as and, after suitable numbering, a sequence of grid points x k → x such that u K k − ϕ has a maximum at x k in the grid.Hence, for all neighbors x i k of x k it holds Using dθ M → 0 and degenerate ellipticity of For x ∈ ∂Ω, since dx N → 0, there exists a sequence of grid points (x k ) k∈N ⊂ ∂Ω converging to x and hence we can use the uniform convergence of u K to u and continuity of u to infer this is equivalent to u(x) = 0 and we can conclude the proof.
Because of the lack of uniqueness of the infinity Laplacian eigenvalue problems discussed in Section 2.3 we can only hope for a numerical scheme which possesses convergent subsequences.Still, thanks to Proposition 4.2 any continuous subsequential limit is a solution to our infinity Laplacian eigenvalue problems.Note, furthermore, that in the absence of a comparison principle stability of the numerical scheme F [u] = 0, as assumed in [54], does not suffice to prove convergence.Instead we will need a stronger form of stability, namely uniform Lipschitz continuity of discrete solutions, which allows us to obtain the existence continuous subsequential limits and lets us apply Proposition 4.2.This approach is also mentioned in [20].
In the following proposition, which can be interpreted as discrete-to-continuum version of the Arzela-Ascoli theorem, we state that a uniformly bounded sequence of grid functions with uniformly bounded Lipschitz constants gives rise to a precompact sequence of continuum functions u K := ûK • π K .The proof is contained in [4] and we only sketch the necessary ingredients here.Proposition 4.3.Let Ω be a Lipschitz domain and assume that the sequence of grid functions ûK : Assume that dx M , dx N , dθ M → 0 as K → ∞.Then the sequence of functions u K := ûK • π K : Ω → R for K ∈ N admits a subsequence converging to a Lipschitz continuous function u : Ω → R.
Proof.The statement of the proposition was proved in large generality in [4].However, the authors used slightly different assumptions which are necessary for other results there but for the compactness statement alone they are not.Following the proofs of [4,  verbatim and utilizing that since Ω is a Lipschitz domain there exists a constant C > 0 such that d g (x, y) ≤ C|x − y| for all x, y ∈ Ω, where d g (•, •) denotes the geodesic distance in Ω, we obtain the existence of a convergent subsequence.The Lipschitz continuity of the limit follows from [4, Lemma 5].
Hence, for showing convergence of our numerical scheme it will suffice to show that it is degenerate elliptic, consistent, and allows for solutions with uniformly bounded norms and Lipschitz constants.Then Propositions 4.2 and 4.3 imply subsequential convergence to a viscosity solution.

Approximation of the distance function and infinity eigenvalues
As we have already seen in Theorem 2.4 the first eigenvalue Λ 1 is directly linked to the geometry of the domain, i.e., For simple domains Ω ⊂ R 2 , such as a circle, square, or triangle, the so-called in-radius r 1 and hence also the first eigenvalue Λ 1 can be easily calculated by geometric reasoning.In general, for a more complicated domain Ω we have to compute the distance function d(x) = dist(x, ∂Ω) which is the unique solution of the following Eikonal equation From the solution of (4.6) we thus obtain the in-radius r 1 together with the set of points in Ω where this maximal distance to the boundary is attained.The solution of the Eikonal equation on a discrete grid can be approximated with different methods, the best-known of which is the fast marching method [48].Originally formulated on structured grids, it was generalized to weighted graphs in [27].Alternatively, it was shown in [26], that the solution of (4.6) coincides with the solution to the optimization problem max and it was characterized as nonlinear eigenfunction of a subdifferential operator in [14].There, the same was shown for a graph analogue of (4.7).Therefore, one can also use the gradient flow based methods [10,13,19] to solve discrete versions of (4.7) or (4.6), respectively.Hence, by employing any of these methods, one obtains a discrete distance function and an associated first eigenvalue Λ1 (cf.(4.5)) subordinate to the discrete grid V , defined in Section 4.1.Uniform convergence of distance functions on general graphs to continuum distance functions was proved in [2,3,4].A fortiori, this also implies the convergence of the first eigenvalue Λ 1 in (4.5).
One should remark that the second eigenvalue Λ 2 cannot be approximated as easily.Remember that it has a geometric characterization as reciprocal of the maximal radius of two equal nonintersecting balls which fit into the domain (cf.(2.9)).For many symmetric domains (e.g.circle, square, isosceles right triangle, L-shape, etc.) the solution of this sphere packing problem can be derived using elementary geometric reasoning.However, we could not find a circle / sphere packing algorithm in the literature which works for general domains.
Furthermore, higher infinity-eigenvalues have not yet been characterized.Only in some special cases one knows that they are given by the reciprocal of the maximal radius of k equal non-overlapping spheres which fit into the domain [40].In these cases one can use known solutions of the general sphere packing problem to obtain the eigenvalue.

Approximation of the first eigenfunction
Let us consider the first eigenfunction problem: As in Section 4.1 we subdivide the set of vertices in the discrete grid into V = V inn ∪ V bdry where V inn denotes the inner nodes of the grid and V bdry denotes the boundary vertices.We approximate (4.8) by the discrete scheme (4.1)where and F + 1 and F2 are degenerate elliptic and consistent grid functions, which implies the same for F .Note that the superscript in F + 1 serves to distinguish this scheme from a similar one for higher eigenfunctions, introduced in the next section.
We first discuss the grid function F + 1 which is the novelty of our approach.Taking (4.8) into account we have to approximate the term |∇u| − Λ 1 u.In the following, we fix a vertex x i and denote the distances to its neighbors by where the index j max is chosen such that The number Λ1 is the reciprocal of the maximal value of the distance function on the grid.By definition of j max the quantity is non-decreasing in the differences (u i − u j )/d ij and hence degenerate elliptic.This is in strong contrast to the naïve approximation |∇u(x i )| ≈ max j∈Ni |u i − u j |/d ij which is not monotone.Furthermore, using a Taylor expansion for any function u Choosing j = j max we observe which together with convergence of the eigenvalue Λ1 (see Section 4.2) implies the consistency of (4.12).Note that thanks to the homogeneous nature of (4.1) we can readily work with F + 1 , which enjoys the numerically convenient property of a Lipschitz constant which is close to one.
Next we recap the approximation of the infinity Laplacian due to Oberman in [41].One defines a discrete Lipschitz constant L(u i ) of u in x i for i ∈ {1, . . ., M } as In [41,Theorem 5] it was proved that the minimizer of this discrete Lipschitz constant with respect to u i is given by where the indices r, s ∈ N i are chosen such that Furthermore, u * i is non-decreasing as a function of {u j : j ∈ N i } and under some mild technical conditions on the grid it holds where dx M and dθ M denote the errors from Equation (4.3).Hence, we can define the grid function F2 in (4.9), evaluated in a grid point x i as Note that the following quantity is consistent: Furthermore, it holds that Hence, there is a multiple of F2 which is consistent and another one which is degenerate elliptic.However, since (4.1) is an homogeneous equation we can simply work with F2 which has the convenient property of having unit Lipschitz constant (cf.[41]).
We are now able to state our main result in this section, addressing subsequential convergence of solutions to (4.1) to solutions of (4.8).For this we heavily rely on Propositions 4.2 and 4.3.To be able to apply them, we have to normalize the discrete appropriately which is indicated by the homogeneous nature of the considered equations, anyway.Theorem 4.4.Let Ω be a Lipschitz domain and assume that for every i ∈ {1, . . ., M } and every r ∈ N i there exists s ∈ N i such that x i − x r and x s − x i are parallel vectors.Let the grid functions ûK : V → R solve (4.1) with F given by (4.9), normalized such that If dx M , dx N , dθ M → 0 as K → ∞, then a subsequence of the sequence u K := ûK • π K converges uniformly to a continuous function u : Ω → R which is a viscosity solution of (4.8).
Remark 4. Let us briefly comment on the assumption on the grid.We just pose it in order to cite the consistency result [41,Theorem 6] which states that (4.14) is valid.The assumption requires some symmetry of the grid and the stencil N i which is satisfied for instance, by a uniform rectangular or hexagonal grid with a distance-based stencil.We would like to remark that in [41] the author also poses the assumption that all points in the stencil N i have distance of the same order to the central point x i which essentially requires the stencil to be a ring of points.However, carefully redoing the proof of [41,Theorem 6] one observes that actually the condition is not needed and furthermore that the factor of 2 in (4.14) is missing in the original reference.
Proof.Using [41, Theorem 6] we see that (4.14) holds true .Since we have already established degenerate ellipticity and consistency of F + 1 and F2 (and hence also of the combined scheme F ), it only remains to prove that which will allow us to conclude the proof by applying Propositions 4.2 and 4.3.
For this we simply note for any x j ∈ V bdry we have the trivial discrete Lipschitz estimate where d denotes the graph distance function from x i to V bdry defined through the dynamic programming principle Consequently, taking the maximum over i and using that Λ1 is the maximum value of the distance function, it holds Using the convergence of the eigenvalue Λ1 to Λ 1 (see Section 4.2), we can take the supremum over This concludes the proof.
Remark 5 (Existence of discrete solutions).Note that we do not prove existence of roots of the grid function (4.9) since this requires some lengthy discrete theory which is beyond the scope of this paper.The pipeline for proving existence should follow the approach in the continuum: Starting with discrete p-Laplacian problems for which existence can be proved on a variational level, one proves compactness as p → ∞ and passes to the limit in the optimality conditions of the discrete p-Laplacian problems to arrive at our discretization F [û] = 0.
Remark 6 (Local monotonicity).One might ask whether the grid function F + 1 , and hence also the combined function F , is monotone in the nodal values u i .In the theory of monotone schemes this would ensure that the scheme F [u] = 0 possesses a unique solution, which cannot be expected.However, F + 1 can be rewritten as is non-negative if the stencil size dx M , and hence in particular d ijmax , is sufficiently small.Since the term u jmax does not change for sufficiently small changes in u i , function F + 1 is at least locally monotone.Furthermore, from this representation it can be seen that F + 1 is locally 1-Lipschitz.
Remark 7 (Concave approximations of the ground state problem).As already discussed in Section 1 the solution of the ground state problem (4.8) are in general not unique.To alleviate this problem one can slightly modify the original problem and study a family of concave approximations, parameterized by a parameter α ∈ (0, 1), of the following form: This equation has been introduced and analyzed in [18], admits a comparison principle and hence has a unique solution.Furthermore, choosing α in (4.16) such that α ↗ 1 it was shown that the solutions of the corresponding problems converge towards the pointwise maximal solution of (4.8).It is trivial to see that changing F + 1 [u], which was defined in (4.10), to yields a convergent discretization in the sense that Theorem 4.4 is valid.Even more, since (4.16) has a unique solution, the convergence is not only subsequential.

Approximation of higher eigenfunctions
Similarly as before, we would like to approximate our reformulation for higher eigenfunctions as monotone scheme.Analogously to the previous section we approximate this equation with F [u] = 0 where F [u](x i ) for x i ∈ V inn is now given by and F + 1 and F2 are as in (4.10) and (4.15), respectively.The function F − 1 is given by Here the index j min is chosen such that and hence (4.19) approximates a multiple of −|∇u| − Λu.As for the first eigenfunction, one can easily see the consistency of the grid function.Furthermore, monotonicity in the differences u i − u j is proved as in Proposition 3.2, using that the non-monotone term − F2 [u](x i ) in (4.18) always vanishes and the first two terms are obviously monotone.

Numerical solution of the schemes
Now we describe how we solve F [u] = 0 where F given by (4.9) or (4.18).Due to the non-smoothness of the grid functions F , Newton-type methods are not applicable to compute a root of F .Also quasi-Newton methods require some degree of (directional) differentiability in the root x * such that F [u * ] = 0 in order to converge (cf.e.g.[50,51]).Due to the strong non-smoothness of F given by (4.9) or (4.18), this is too much of an assumption.Hence, it seems natural to study the simple fixed-point iteration where is referred to as Euler map.Obviously, roots of F correspond to fixed-points of Ê.The terminology "Euler map" stems from the obvious fact that (4.21) can be seen as explicit Euler discretization of the ODE u(t) = − F [u(t)] with time step size ρ > 0. Is is well-known (cf.[41], for instance) that if ρ > 0 is smaller than the reciprocal Lipschitz constant of F and F is monotone in the sense that u ≥ v implies F [u] ≥ F [v] in the partial order in R M , then the Euler map E is a contraction.Since this would in particular imply a unique fixed point of E and hence a unique root of F , we cannot expect this in our case.However, due to the "local monotonicity" of F (cf. Remark 6) one can expect that in the proximity of a root the map F is monotone and hence Ê is a contraction there.In practice, the fixed point iteration (4.21) converges very reliably on our numerical experiments.For designing a stopping criterion we utilize both the relative changes of the iterates and the accuracy of the root.The detailed algorithm to find a root of F , and hence an infinity Laplacian eigenfunction, is given in Algorithm 1.

Experimental results
In the following, we present numerical results which use the schemes and algorithms from Section 4.Many of the experiments deal with open questions and conjectures regarding infinity eigenfunctions and, thereby, we hope to shed some light on the theory.
The computations take place on a regular grid which discretizes the unit square [−1, 1] 2 .In order to compute on more general domains we simply restrict the computations on those grid nodes which belong to the domain of interest (see, for instance, Section 5.1.3below).
Algorithm 1 Root finding of F [u] = 0 where F is given by (4.9) or (4.18) k ← k + 1 8: end while 9: return u In all experiments apart from the very first one we choose the number of local neighbors k i of a generic node x i ∈ V inn -which appears in (4.11) and (4.13), for instance-as k i = 120.In our regular grid this corresponds to a quadratic stencil of 11 × 11 around the node of interest.If parts of the stencil leave the computational domain-which happens close to the boundary, for instance-we simply reduce the number of neighbors of the corresponding node.We discretize the unit square including its boundary with 97 × 97 nodes.
If not stated differently, the inputs in Algorithm 1 were chosen as follows.When computing infinity ground states (cf.Section 5.1 below) the initial guess u 0 ∈ R M is chosen as discrete distance function of the domain1 .The constant ρ > 0-which should be chosen smaller than the reciprocal Lipschitz constant of F -is chosen as ρ = 0.9.Note that the Lipschitz constants of F given by (4.9) or (4.18) are close to one.We allow for a maximum of K = 5000 iterations and choose the tolerance TOL = 10 −7 .Let us remark that in almost all our experiments the algorithm required only a few hundred iterations in order to reach the tolerance.Furthermore, the tolerance should scale with the square of the characteristic grid size which can be seen from (4.14).
Our implementation uses MathWorks MATLAB ® R2018b and a typical test-case requires a few minutes of computation on a standard laptop computer.Code is available on GitHub. 2

Infinity ground states
In this section we perform numerical experiments for infinity ground states, by computing a root of (4.9).The eigenvalue occurring in (4.10) is chosen as maximum of the distance function on the grid.

Influence of the number of local neighborhood size
First, we would like to investigate the influence of the number of local neighbors k i on the computed ground state.Remember that due to Theorem 4.4 one can expect more accurate results as the number increases.In Figure 1 we show the level lines of the ground state on the unit square [−1, 1] 2 , computed using neighborhoods of size 3 × 3, 5 × 5, 7 × 7, and 11 × 11.Looking at the level lines, one can observe that the smoothness of the ground state increases as the neighborhood size grows.This can be explained by a more accurate approximation of ∇u and its norm.Further experiments show the same behavior for the infinity harmonic function on the punctured square [−1, 1] 2 \ {0} (see also [41] for similar observations).In the following experiments we will use the 11 × 11 stencil in order to produce accurate results.

Infinity ground state and infinity harmonic on the square
In this experiment we numerically investigate the long-standing conjecture that the infinity harmonic function on the punctured square, i.e., the square without its center point, is a ground state (cf.e.g [47]).Note that the analogue of this statement is known to be true on stadium-like domains [36] (like for instance the ball) but is false in general [30].In fact, only recently the conjecture was proved false by constructing the explicit unique solution of the infinity Laplace equation on the punctured square which does not solve the ground state problem [1,5].
In Figure 2 we show the infinity ground state on the square, computed with our method, the infinity harmonic function on the punctured square, and their pointwise difference.Note that we compute the infinity harmonic function by simply solving the scheme F2 [u] = 0 together with appropriate boundary conditions, where F2 is given by (4.15).We compute both solutions with high accuracy such that they solve their respective equations up to an L ∞ -error of at most 10 −9 .Although not visible from the solutions themselves, their difference, plotted on the right, exhibits an L ∞ -norm of order 10 −3 which confirms that ground state and infinity harmonic do not coincide. 3Furthermore, as predicted by theory [47], the ground state is pointwise larger or equal than the infinity harmonic function which is reflected by the difference being non-negative.

Infinity ground states on different domains
With this test-case we demonstrate the aptness of our algorithm to compute infinity ground state also on more complicated and in particular non-convex domains.Figure 3 shows the computed ground states on six different convex (top row) and non-convex (bottom row) domains.The shapes of the ground states are similar to p-Laplacian ground states for large values of p, see for instance [23,32].

Regularity of ground states
Next we study the regularity of ground states which is still an open problem from the theoretical perspective.Some results on singular sets of ground states were proven in [36], among which are the statements that ground states are non-differentiable in the maximal set of the distance function and that singular points of the gradient are not isolated.Furthermore, in two dimensions ground states are C 1 away from their maximal set if and only if they are infinity harmonic there.However, general statements on the regularity of ground states outside the maximal set are still pending.
Here we recap the ground state computed for the dumbbell shape (see also bottom center in Figure 3).Figure 4 shows the level lines of the ground state and exhibits a non-smoothness along the line segment which connects to the two maxima.Here the level lines show kinks and even touch in the center of the domain.This suggests that ground states are non-differentiable, in general.3).The gradient looks singular between the two maxima.

Discrete non-uniqueness on the rectangle
In this experiment we address the question of uniqueness of infinity ground states, computed with our method.As mentioned in Remark 1, it is known that ground states are in general not unique, see, e.g., [28] in which a non-convex dumbbell domain was constructed where uniqueness fails.However, for convex domains uniqueness is neither proved nor disproved apart from the case of stadium-like domains [36].
In line with these observations for the continuous problem, non-uniqueness can also be observed when numerically computing ground states with our discrete scheme.Depending on the initialization of the scheme one may get two substantially different solutions of the discrete problem.
In Figure 5 we show surface and contour plots of two different discrete ground states on the rectangle [−1, 1] × [−0.5, 0.5], computed with our method.Both results fulfill F [u] = 0 with very high accuracy.In this experiment we fix the value u(0, 0) = 0.5 which is no loss of generality due to the homogeneity of the eigenvalue problem (4.9).The first result was computed by initializing Algorithm 1 with the distance function, whereas the second one was initialized with zero.The two results differ significantly: the first one attains its maximum on the so-called high ridge of the rectangle, given by [−0.5, 0.5] × {0}.In contrast, the second result attains its maximum only in the point (0, 0) and its level lines are not parallel to the long sides of the rectangle as it is the case for the first ground state.However, for finer directional and spatial resolutions of the grid the two computed ground states do not differ as severely anymore, as can be seen in the two bottom rows of Figure 5.This is, of course, in Figure 5: Discrete non-uniqueness of ground states on the rectangle for different grid resolutions.Surface and contour plots of computed results, initialized with distance function (left) and zero (right).The pointwise difference at a low resolution is substantial (top rows).For high resolutions they are less different (bottom rows).line with our convergence result Theorem 4.4 and the fact that the continuum ground state attains its maximum on the entire high ridge [36].Note that in particular the high directional resolution, which is guaranteed through large computational stencils, is the limiting factor in computing solutions on even finer grids.
Alternatively, instead of computing high resolution ground states, which is numerically cumbersome, we suggest to solve the concave approximation (4.16) of the original discrete problem in order Figure 6: The corresponding results to Figure 5 for the concave appoximation model (4.16).
to gain a unique solution.Naturally, we choose α ∈ (0, 1) close to 1, as discussed in Remark 7, to numerically approximate the original discrete problem well.In Figure 6 we show that, even for small stencils and a low resolution, this slight modification effectively selects a ground state which is independent of the initialization of the scheme and has the desired properties of the high ridge as discussed above.

Higher eigenfunctions
Until now we numerically computed first eigenfunctions of the infinity Laplacian.In the following we concentrate on the computation of higher eigenfunctions, which comes with two major challenges.
First, the computation of higher eigenvalues is difficult, as already explained in Section 4.2.While the second eigenvalue has a variational characterization in terms of a sphere packing problem (2.9), such a property is unknown for higher eigenvalues so far.However, for certain symmetric domains like the square one can compute higher eigenvalues explicitly.In this section we concentrate on symmetric domains for our experiments, such that we are able to deduce the unknown eigenvalue from its geometry.
This leaves the problem of finding a suitable initialization for Algorithm 1.Here, we investigate three different strategies: The first one is to use the symmetry of the domain to fix the maximal and minimal value of the eigenfunction in two designated points and initialize the rest with zero.The second one consists in initializing with higher eigenfunctions of the linear Laplacian.The third strategy is adapted to second eigenfunctions and consists in initializing randomly and slightly modifying Algorithm 1 by adding the normalization step after the fixed-point iteration.This assures that the maximum of the positive and negative parts of u are equal, which is a necessary condition for second eigenfunctions [40].
Figure 7 shows second infinity eigenfunctions on three different domains, for which the second eigenvalue can be computed.The algorithm was initialized with zero and the peak values of the eigenfunctions were fixed.Again the results are similar to second eigenfunctions of the p-Laplacian for large p (e.g.[32]).In Figure 8 we show three different eigenfunctions on the square where we initialized with the first three eigenfunctions of the standard Laplacian on the square which can be computed with standard linear algebra tools 4 .Finally, we also show a result which was computed with the normalizations steps (5.1) of the positive and negative parts of the solution.The initialization was chosen as random noise and also this method converges to a second eigenfunction nicely, as visualized in Figure 9, which shows the solution after 0, 300, and 655 iterations of Algorithm 1 with the normalizations (5.1).

Conclusion
In the first part of this work we have presented a reformulation of the eigenvalue problem for the infinity Laplacian, such that higher eigenfunctions can be computed, avoiding the distinction of cases on their sign.In the second part, we proposed consistent and monotone schemes for approximating ground states and higher eigenfunctions, building on our reformulation from the first part.The discrete problem is solved by a fixed-point iteration.Our numerical results show the aptness of the proposed numerical scheme to approximate eigenfunctions even on complicated domains.This appears to be the first numerical method in the literature to compute infinity Laplacian eigenfunctions.
There are four open problems related to our work, which will be subject to future research.First, our aim is to find an algorithm to compute the second eigenvalue on general domains.While on symmetric domains one can easily compute it using the distance function on the expected nodal domains, a theoretically sound approach for more complicated domains has to be investigated.We believe that the several variational characterizations of the second eigenvalue in [40] can be a promising route for that.
Second, the fixed-point iteration, which we currently use to compute roots of the grid functions, can possibly be replaced by a more involved non-smooth Newton-type method.However, currently we are not aware of an alternative method which can handle the strong non-smoothness of the problem.
Third, to replace subsequential by full convergence as the grid gets finer, a promising avenue might be to couple the concave approximation parameter α in (4.16) with the grid parameters in order to ensure convergence to the unique pointwise maximal ground state.For higher eigenfunctions, however, no such selection principle is known yet, which is another interesting avenue of research.
Finally, proving convergence rates for our numerical scheme is an especially challenging endeavour due to the strong non-smoothness of infinity ground states.In parts of the domain where a solution of min(|∇u| − Λu, −∆ ∞ u) = 0 is C 2 , it is know [47] that it in fact solves ∆ ∞ u = 0.For this infinity Laplacian equation convergence rates in the supremum norm were established recently [2,7].The techniques used there strongly rely on comparison principles which are not available for the infinity eigenfunction problem.Furthermore, if rates can be proved, they will be extremely slow which can be seen from the example in Section 5.1.5.
or non-financial interest in the subject matter or materials discussed in this manuscript.The authors have no financial or proprietary interests in any material discussed in this article.

Theorem 2 . 4 .
Let Ω be an open, bounded domain.Then there exists a positive viscosity solution

Definition 4 . 1 (
Properties of grid functions).The grid function F is called • degenerate elliptic if for i = 1, . . ., M the functions Fi are non-decreasing in the variables 2, . . ., k i .

Figure 1 :
Figure 1: From left to right: level lines of infinity ground state on the square for stencils of size 3 × 3, 5 × 5, 7 × 7, and 11 × 11.Smoothness of the level lines increases with larger neighborhoods.

Figure 3 :
Figure 3: Infinity ground states on different domains.All results were initialized with the distance function of the domain.

Figure 4 :
Figure 4: Level lines of the dumbbell ground state (cf.bottom center in Figure3).The gradient looks singular between the two maxima.

Figure 8 :
Figure 8: Infinity eigenfunction on the square, initialized with first three Laplacian eigenfunctions.