{'fraction_non_alphanumeric': 0.07370477392047536, 'fraction_numerical': 0.02376785641610028, 'mean_word_length': 3.624882364012799, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
THE HESSIAN RIEMANNIAN FLOW AND NEWTON’S METHOD FOR EFFECTIVE HAMILTONIANS AND MATHER MEASURES

Effective Hamiltonians arise in multiple problems, including homogenization of Hamilton-Jacobi equations, nonlinear control systems, Hamiltonian dynamics, and Aubry-Mather theory. In Aubry-Mather theory, related objects, Mather measures, are also of great importance. Here, we combine ideas from mean-field games with the Hessian Riemannian flow to compute effective Hamiltonians and Mather measures simultaneously. We prove the existence and convergence of the Hessian Riemannian flow. In addition, we explore the relation between the Hessian Riemannian flow and a variant of Newton’s method that greatly improves the performance of our numerical methods. In our numerical experiments, we see that our algorithms preserve the non-negativity of Mather measures and are more stable than previous methods in problems that are close to singular.


Introduction
Let T d be the unit d-dimensional torus. Given P ∈ R d and a Hamiltonian H : T d × R d → R, the effective Hamiltonian, H(P ), is the unique real number for which there exists a periodic viscosity solution, u : T d → R, of the Hamilton-Jacobi equation (1.1) If H and u solve (1.1), the infimum of (1.2) is −H(P ), µ is supported on the graph (x, −D p H(x, P + D x u(x))) and the x-projection of µ, m(x), is a weak solution of − div(mD p H) = 0. (1.3) Here, we develop an algorithm that computes u, H, and the projected Mather measure, m, simultaneously. For that, we draw inspiration from numerical methods for mean-field games (MFGs).
MFGs study the behavior of rational non-cooperating agents or players in a large population [6,18]. A typical MFG model comprises a system of a Hamilton-Jacobi (HJ) equation and a transport or Fokker-Plank (FP) equation. In MFGs, the FP equation is the adjoint operator of the HJ equation. The following is a first-order stationary MFG: H(x, P + D x u) = H + g(m), − div(mD p H) = 0, (1.4) where u(x) determines the cost for an agent at x ∈ R d , m is a probability density that gives the agents' distributions, and g determines the interaction between agents. When g = 0, (1.4) is the cell problem, H the effective Hamiltonian, and m the projected Mather measure.
To approximate the cell problem, we consider a MFG that also arises in the study of entropy penalized Mather measures [8]: H(x, P + D x u k ) = H k (P ) + 1 k ln m k − div(m k D p H(x, P + D x u k )) = 0, (1.5) where k > 0 is an integer and Under Assumptions 1 and 2 in Section 3, as k → ∞, H k (P ) converges to H, u k converges to a viscosity subsolution of (1.1), and m k converges to a projected Mather measure [8]. To solve (1.5), we construct the Hessian Riemannian flow, that preserves the non-negativity of m. More precisely, we consider the system of PDEs: (1. 8) In Section 3, we establish the following convergence theorem: Theorem 1.1. Suppose that Assumptions 1 and 2 (see Section 3) hold, and that (1.7) has a solution, (m, u) ∈ C [0, ∞); W 1,2 + (T d ) × W 1,∞ (T d ) . Assume further that (m * , u * ) is the periodic smooth solution of (1.5). Then, there exists a sequence {t i } such that u(t i ) → u * in W 1,2 T d as i → +∞. Moreover, we have u(t) → u * in L 2 T d and m(t) → m * in L 1 (T d ), as t → ∞.
Next, in Section 4, we discretize (1.7) in space and obtain a system of ODEs: Ṁ where (M, U ) ∈ R N × R N , N is the number of grid points, and F is defined in (4.7). There, we prove the following theorem: In Section 5, we explore the connection with a variant of Newton's method that is equivalent to the Crank-Nicolson scheme for (1.9). Numerical results and performance comparisons follow in Section 6. In our numerical experiments, we get an accurate approximation for u, H and m. In particular, our methods are stable for problems that are nearly singular. Finally, conclusions and a brief discussion of future work are presented in Section 7.
Notation. We use |·| to represent the l 2 -norm of a matrix or a vector, and · to represent the L 2 -norm of a function. Denote by W 1,p * (T d ) and W 1,p + (T d ), respectively, the spaces of nonnegative and strictly positive functions in Besides, we also denote the inner product of two vectors in a Euclidean space by ·, · . We identify the d-dimensional torus, T d , with [0, 1] d . Finally, R N + is the subset of R N of vectors with positive components.

Previous work
Multiple authors studied and proposed numerical methods for the computation of effective Hamiltonians. Here, we give a brief overview of the various approaches in the literature.
Two approaches described in [25] use the asymptotic behavior of Hamilton-Jacobi equations to compute H. The first approach, called small-δ method, introduces a parameter, δ > 0, and considers the stationary equation According to [19], −δu δ converges uniformly to H(P ) on R d as δ → 0. Thus, we can choose a small δ and solve (2.1) numerically to get an approximation for H(P ). The second method, called in [25] the large-T method, uses a large-time approximation where v is a continuous, periodic function. Under suitable assumptions, (2.2) has a unique viscosity solution on T d × [0, T ], see [26], and [25] established that −u(x, t)/t → H(P ) for a general, not necessarily convex, Hamiltonian, H. Alternatively, the effective Hamiltonian can be computed using a representation formula that arises as a dual problem of an infinite-dimensional linear programming problem [12,14]. This is the idea used in [16], where H(P ) is computed through the formula, by discretizing the spatial variable and solving the minimax problem.
The preceding approaches are slow from the computational point of view. Thus, significant efforts were spent on developing fast algorithms. These include solving a homogenization problem directly [20,24] and employing a Newton-type method [5] to (1.1).
In [20,24], given a function, f , the authors of [20] considered the oscillatory equation Then, the value of f at point x 0 , which is close enough to the minimum of u − P · x, yields an approximation of H(P ) [20]. However, f (x) has a formula that involves the minimum of H, which may be hard to compute for general Hamiltonians.
The generalized Newton method in [5] uses a novel approach to compute the effective Hamiltonian. There, (1.1) is discretized directly into a nonlinear system, F (X) = 0, where X encodes a discretized version of u and H. Then, the resulting system is solved by the Newton method.
The focus of the preceding methods is the computation of the effective Hamiltonian and the viscosity solution. Mather measures do not play a role. In contrast, the variational method [8] approximates the projected Mather measure and the effective Hamiltonian by where k ∈ N and u k is the minimizer of We observe that (1.5) is the Euler-Lagrange equation corresponding to the functional in (2.4). If H(x, p) satisfy Assumptions 1 and 2, the results in [8] imply that m k → m and H k (P ) → H (P ) as k → ∞. Inspired by this, [13] proposed a numerical method, which solves the Euler-Lagrange equation of (2.4) by finite-difference methods and gets H k (P ) using (2.3). Numerical experiments in [13] show that this approximation is more efficient than the algorithm in [16] but with less accuracy. Besides, as pointed out in [13], this scheme is unstable when k is too large for a fixed mesh. In contrast, our methods seem to be stabler as illustrated in Section 6.

MFGs and Effective Hamiltonians
To solve the cell problem and compute the projected Mather measure, we combine (1.1) and (1.3) into the system where m 0 is a probability measure. Taking into account that .
We notice that if u is the viscosity solution of (3.1), so is u + C, where C is an arbitrary constant. So, for the uniqueness of u, we require T d u = 0. Hence, our goal is to solve The previous equation may not have a solution, (m, u), in W 1,2 * T d × W 1,∞ T d . For example, m may be singular. We tackle this matter by introducing various approximation procedures. First, we attempt to use a monotone flow as in [1] to approximate the solution of (3.3). However, we observe that this flow may not preserve the non-negativity of m. This leads us to introduce the Hessian Riemannian flow. Under the assumption of the existence of a solution to (3.1) with m > 0, we prove the convergence for u. Unfortunately, the convergence for m may not hold due to the non-uniqueness of solutions of (3.1) and the possibility of m vanishing. Hence, we add an entropy penalization term to the Hessian Riemannian flow that gives both the positivity and the convergence for m.
3.1. The monotone flow. A way to compute the solution of (3.3) is the monotone flow method introduced in [1]. First, we recall that the operator F defined in (3.2 The monotonicity of F suggests the monotone flow, provided t 0, m 0 and m 0. Suppose that (m * , u * ) is the solution of (3.3). Then, (m * , u * ) also solves (3.5), since d dt m * = d dt u * = 0. Thus, if we suppose further that m ∈ C [0, ∞); according to (3.6). In this case, (3.5) defines a contraction in the region where m is nonnegative. However, there are several issues about the monotone flow. First, we do not know if it is globally defined. Besides, the projected Mather measure may be singular. Finally, the convergence is not guaranteed either. In Example 3.1 below, we show that the monotone flow may not preserve the non-negativity of m. Hence, (3.5) may not give a global contraction.
Let (m 0 , 0) to be the initial point and 1 0 m 0 dx = 1. It is easy to check that (m, u) = (m 0 + sin (2πx) t, 0) is the solution for (3.7). However, m(t) becomes negative in some regions as t → +∞.
Another reason why the convergence may fail is that the solution of (3.1) may not be unique, as the next example illustrates.
It is easy to see that H = 1 2 , u = 0 and m = f (x 2 ), where f is any function that depends only on the second component of x solving (3.8). Thus, m is not unique.
To guarantee the non-negativity of m in the monotone flow, we use the Hessian Riemannian gradient flow introduced in [ where E is the closure of an open, nonempty, convex set E ⊂ R n , A ∈ R m×n with m n, b ∈ R n , and f ∈ C 1 (R n ). To solve this problem, the authors introduced a Riemannian metric, g, derived from the Hessian matrix, ∇ 2 h, of a Legendre-type convex function [2], h, on E. Then, they used the steepest descent flow to generate trajectories in the relative interior of the feasible set, F := E ∩ {x|Ax = b}. In the steepest descent method, the authors sought a trajectory x(t) solving where ∇ H f |F (x) is the projection w.r.t. g of the gradient of f into the admissible directions. According to [2], (3.9) is well-posed. Moreover, this steepest descent flow never leaves the admissible set and leads to a local minimum. A similar idea can be used for monotone operators and lead us to the Hessian Riemannian flow.

The Hessian Riemannian flow.
To guarantee the non-negativity of m, we introduce the Hessian Riemannian flow. More precisely, we define the convex function, h : Thus, the Hessian of h corresponds to the matrix Then, we consider the Hessian Riemannian flow, which can be rewritten as . (3.10) The mass of m is preserved by this flow, because Before proving the convergence of (3.10), we impose a convexity assumption on H.
Then, we have the following convergence results: Suppose that Assumption 1 holds and that (3.10) has a solution, (m, u) ∈ In addition, by the periodicity of u(t) and m(t).
For the convergence, we define a Lyapunov function for t > 0, We know that φ(t) 0 since the mapping z → z ln z is convex for all z 0. Next, by differentiating φ(t) and using the fact thatṁ * = 0,u * = 0, and where we apply Assumption 1 in the last inequality. Then, we have d dt Proof. Suppose that lim inf t→+∞ g(t) = 0. Then, we can find t 0 0, > 0, such that for any t > t 0 , we have g(t) > . This contradicts the fact that +∞ 0 g(t) < +∞.
Unfortunately, the convergence of m for (3.10) may not hold since solutions of F (m, u) = 0 may not be unique and m * may fail to be positive as shown in Example 3.2. This observation motivates us to introduce an entropy penalization that we discuss next.
3.4. Entropy penalization. To obtain uniqueness for the projected Mather measure, we consider the entropy penalized model given by (1.5). Combining (1.6) and the first equation of (1.5), we get Thus, H k (P ) can be rewritten as To prove the convergence and the existence of the solutions of (1.5), we introduce another assumption on H, in addition to Assumption 1.
Assumption 2. The Hamiltonian, H, satisfies the following natural growth conditions: for some constant C > 0.
According to [8], under Assumptions 1 and 2, for each k, there exists a unique smooth solution, (u k , m k ), of (1.5). Besides, we have m k m weakly as a measure on T d , and, for each 1 q < ∞, In addition, Moreover, Therefore, as k → +∞, we get that u is a subsolution for (1.1). The monotone flow for (1.5) may not preserve the mass of m. Instead, we explore its Hessian Riemannian flow, which is given in (1.7). We notice that the mass of m is constant since Next, we prove the convergence for both m and u, We also define the same Lyapunov function for t > 0.
where we use Assumption 1 in the last inequality.
Then, by Lemma 3.4, we have Since m * is strictly positive on T d , we have By the Poincaré inequality, we obtain Thus, we conclude that we get, Accordingly, it follows that u(t) → u * in L 2 T d . Besides, by rewriting (3.12), we obtain By Lemma 3.5 below, we get for any 0 < < 1, Then, using (3.13), we obtain (3.14) Because (3.14) holds for any ∈ (0, 1), we consider the limit → 0, and get Lemma 3.5. Suppose that a, ∈ R with a > 0, 0 < < 1, then for any z > 0, we have Proof. When z a, we define So, we have Thus, g achieves its minimum when z = a 1− . Since g a
We differentiate f with respect to z, and get Thus, f achieves its minimum at z = a 1+ . Evaluating f at z = a 1+ , we obtain So, f (z) 0 when z < a. Therefore, we conclude that for any z > 0, a ln a z − (a − z) − a ( + ln(1 − )) |z − a| .

A numerical scheme for the Hessian Riemannian flow
Let T d ∆ be a uniform grid on T d , and x i , i = 1, . . . , N the grid points. We discretize u and m at the grid by U = (u 1 , . . . , u N ) and M = (m 1 , . . . , m N ). In addition, we impose periodicity of u and m using a straightforward convention; for d = 1, we set u 0 = u N and m 0 = m N . Our difference scheme for H is An example of G i is given in Section 6. Let L U : R N → R N be the linearized operator of G at U ∈ R N and L * U its adjoint operator. We define F : Then, the space-discretized version of (1.5) is where H : R → R is a numerical approximation of the effective Hamiltonian. Before solving (4.3), we lay out the assumptions on G. To ensure the monotonicity of F , we require each component of G(U ) to be convex.
To prove the local existence of the Hessian Riemannian flow, we need to assume that each partial derivative of G i (U ), 1 i N , is locally Lipschitz.
Finally, to guarantee the convergence of the Hessian Riemannian flow, we require the following property: Under Assumption 3, F is a monotone operator, as we prove below in Lemma 4.2.
Proof. Let (M, U ) , (Θ, V ) be as above. We have because of the positivity of m j and of the convexity of G j . Similarly, Combining (4.5) and (4.6), we conclude that To construct the Hessian Riemannian flow corresponding to (4.3), we define F : Accordingly, the Hessian Riemannian flow is Under Assumptions 3 and 4, F is locally Lipschitz continuous on R N + × R N . Moreover, since F depends only on (M, U ), we have local existence of the solution for (4.8); that is, given (M 0 , U 0 ) ∈ R N + × R N , the initial value problem in (4.8) has a unique solution on t ∈ (0, T ), for some 0 < T +∞.
is bounded as t → T . Let (M * , U * ) = (m * 1 , . . . , m * N , u * 1 , . . . , u * N ) be the solution of (4.3). Define By the convexity of the mapping z → z ln z, z ∈ R, we have (4.10) Differentiating φ and usingṁ * j =u * j = 0, we get using, as before, the identity Thus, using the definition of L and L * , we have Due to the convexity of G j and the monotonicity of z → ln z, φ is decreasing. In addition, due to (4.10), we have Next, we prove (4.8) is well-posed. We finish this section with the lemma used in the proof of the previous theorem. Let Proof. We prove Ω ⊂ R N + ×R N by contradiction. Suppose that Ω ⊂ R N + ×R N . We can find a sequence t i such that (M (t i ), (m 1 (t), . . . , m N (t)) , U (t) = (u 1 (t), . . . , u N (t)) . From (4.8), we know that Since (M (t i ), U (t i )) → (M * , U * ), the left-hand side of (4.13) converges to −∞. However, by Proposition 4.3, (M * , U * ) is bounded. So, the right-hand side of (4.13) is finite, which gives a contradiction.
To guarantee that   Proof. Let Because of (4.10), φ 0. According to (4.11), we have Using Lemma 3.4 and Assumption 5, we can find a sequence, Under Assumption 6, we have N j=1 u j (t i ) − u * j = 0. Then, combining (4.14), and Lemma 4.8 below, we conclude that Since φ is decreasing, we know Thus, Proof. Assume that a N +1 = a 1 . We consider the linear subspace, equipped with the standard l 2 -norm. Then, D is isomorphic to the quotient space R N /R. We notice that R N /R has another norm given by In addition, because all norms in a finite-dimensional linear space are equivalent, we conclude that (4.15) holds.
Finally, we record the proof of Theorem 1.2.
Proof (of Theorem 1.2). The global existence is given by Proposition 4.4 and the convergence follows from Proposition 4.7.

Newton's Method for Effective Hamiltonians
Here, we explore the connection between the Hessian Riemannian flow and Newton's method and construct a numerical scheme that, in our numerical tests, improves substantially the speed of the Hessian Riemannian flow. To motivate our method, we begin by discretizing (4.8) using the implicit Euler method. Let (M j , U j ) represent the result of the j−th literation, and (M 0 , U 0 ) be the initial value. The implicit Euler method computes M j+1 , U j+1 implicitly using the equation where ξ is the step size. Adding ξ 2 F M j+1 , U j+1 − ξ 2 F M j , U j to both sides of (5.1), we get The prior identity is the Crank-Nicolson scheme with a step size ξ for the following ODE: Let J F (M, U ) be the Jacobi matrix of F at (M, U ). Then, (5.2) is equivalent to Here, we fix two non-negative parameters κ and τ and consider the following generalization of the explicit Euler method for (5.3): When κ = 1 and τ = 0, we obtain Newton's method. When κ = 0 and τ > 0, we obtain the explicit Euler scheme for (1.9).

Numerical results
In this section, we discuss numerical results. Our algorithms were implemented in Mathematica 10 on MacBook Air (CPU: 1.6 GHz Intel Core i5; Memory: 4 GB 1600 MHz DDR3). For the Hessian Riemannian flow, we use the built-in routine, NDSolve, of Mathematica to solve (4.8) with default settings. For Newton's method, we use the iteration given in (5.4).  This formulation has the advantage that M = e W is automatically positive. While the Hessian Riemannian flow preserves positivity, small truncation errors sometimes give rise to negative values of M . Without the above transformation, this would cause serious numerical difficulties.
6.1. One-dimensional case. For d = 1, p ∈ R, and x ∈ T 1 , we set Thus, (3.1) becomes In this case, we compute the analytic value of the effective Hamiltonian using the explicit formula given in [5]: . , u N ) be the approximation of (m, u) at the grid points. Then, we approximate |Du(x) + P | at x i by Accordingly, we get We also use similar schemes for H in other examples. In the algorithms, we set the initial value, (M 0 , U 0 ) = (m 0 1 , . . . , m 0 N , u 0 1 , . . . , u 0 N ), where m 0 i = 1 + 0.9 cos(2πx i ) and u 0 i = 0.2 cos(2πx i ). For Newton's method, we choose τ = κ = 1. Figure 1 plots the effective Hamiltonians versus their approximated values calculated using the Hessian Riemannian flow (HRF) and Newton's method (NM). Figure 2 shows the evolution of H, m, and u for P = 0.5 and k = 100. In Figure 1, we see that our method is extremely accurate away from the flat part of the effective Hamiltonian. In the flat part, the Mather measure corresponding to the different values of P is not strictly positive (see, for example, Figures 2b and 2e at the terminal time), and the logarithmic term seems to slow the convergence speed.
To illustrate the convergence of our methods, we introduce error functions measuring the difference between the numerical result and the exact solution of (1.5). Let (m(t), u(t), H(t)) denote either the solution of the Hessian Riemannian flow or Newton's method. Besides, (m * , u * ) is the exact solution of (1.5) and H * be the corresponding effective Hamiltonian.
Inspired by Theorem 1.1, we define errors:  for the Hessian Riemannian flow and for Newton's method. We see that the errors decrease exponentially.
6.2. Higher Dimensions. Let d = 2, p = (p 1 , p 1 ) ∈ R 2 , and x = (x 1 , x 2 ) ∈ T 2 . We consider the two-dimensional Hamiltonian discussed in [16]: Here, we choose P = (1.5, 2.5) for which H=4.4099660 according to [16]. For Newton's method, we set τ = 2 and κ = 1. Fixing N = 144, Table 1 shows computed H(P ) at t = 14  for different values of k. We see that when k = 10000, we get a very accurate approximation for H. Figure 4 plots m and u at t = 14 when k = 10000.   (m i (t) − m i (t)) 2 + (u i (t) − u i (t)) 2 versus time in Figure   5a, which shows that the curve is not strictly decreasing. Thus, F fails to be monotone. In contrast, we plot φ defined in (4.9) (see Figure 5b), which shows that φ(t) is decreasing, as proven.
6.4. Speed Comparison between the Hessian Riemannian flow and Newton's method. Here, we compare the speed of the Hessian Riemannian flow in (4.8), which is solved by NDSolve, with the speed of Newton's method in (5.4). We consider the Hamiltonian, H(x, p) = |p| 2 2 − sin(2πx).
In the numerical experiment, we set P = 0.5 and k = 100. In this case, H(P ) = 1. The initial point is given by (M 0 , U 0 ) = (m 0 1 , . . . , m 0 N , u 0 1 , . . . , u 0 N ), where m 0 i = 1 + 0.9 cos(2πx i ) and u 0 i = 0.2 cos(2πx i ). For Newton's method, we choose τ = κ = 1. For each value of N , we compute H by the Hessian Riemannian flow for a large time, T , and use it as a benchmark, named H . Then, we use the Hessian Riemannian flow and Newton's method to compute H(T ) such that |H − H | < and record T and the corresponding CPU time (measured in seconds) in Table 2. Here, we choose T = 50 and = 0.001. We see that Newton's method is substantially faster than the Hessian Riemannian flow as implemented using the built-in routine, NDSolve, in Mathematica. Besides, as N increases, the execution time of Newton's method grows slower than that of the Hessian Riemannian flow.
6.5. Stability of the Hessian Riemannian flow and Newton's method. Though not stated explicitly in [13], the algorithm described there computes both H and the projected Mather measure. However, as stated in that paper, that scheme becomes unstable if k is too large compared to the mesh size, N . Here, we show that the Hessian Riemannian flow and Newton's method overcome this issue.  To illustrate the stability of our methods, we consider H(x, p) = |p| 2 2 − sin(2πx).
For the implementation, we choose P = 0.5, k = 100000, and N = 20. We use the initial value, (M 0 , U 0 ) = (m 0 1 , . . . , m 0 N , u 0 1 , . . . , u 0 N ), where m 0 i = 1 + 0.9 cos(2πx i ) and u 0 i = 0.2 cos(2πx i ). Figure 6 shows the evolution of H(P ), m and u by the Hessian Riemannian flow and Newton's method, which illustrates that the Hessian Riemannian flow and Newton's method are stable for nearly singular equations, corresponding to a large value of k.

Conclusion
In this paper, we suggested two methods to calculate simultaneously the effective Hamiltonian and the Mather measure: the Hessian Riemannian flow and Newton's method. We proved the existence and convergence of the Hessian Riemannian flow. We showed that this method guarantees the non-negativity of m. Besides, we pointed out the relation between the implicit discretization of the Hessian Riemannian flow and Newton's method. In our numerical experiments, Newton's method is faster than the Hessian Riemannian flow. Both methods preserve the positivity of the Mather measure. Moreover, the Hessian Riemannian flow and Newton's method seem to be stable for large k, a case where the variational method in [13] faces difficulties.