    # Auxiliary Space Preconditioners for C^0 Finite Element Approximation of Hamilton–Jacobi–Bellman Equations with Cordes Coefficients

In the past decade, there are many works on the finite element methods for the fully nonlinear Hamilton–Jacobi–Bellman (HJB) equations with Cordes condition. The linearised systems have large condition numbers, which depend not only on the mesh size, but also on the parameters in the Cordes condition. This paper is concerned with the design and analysis of auxiliary space preconditioners for the linearised systems of C^0 finite element discretization of HJB equations [Calcolo, 58, 2021]. Based on the stable decomposition on the auxiliary spaces, we propose both the additive and multiplicative preconditoners which converge uniformly in the sense that the resulting condition number is independent of both the number of degrees of freedom and the parameter λ in Cordes condition. Numerical experiments are carried out to illustrate the efficiency of the proposed preconditioners.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Let be a bounded, open, convex polytopal domain in , where represent the dimension. In this paper, we are interested in the Hamilton–Jacobi–Bellman (HJB) equations of the following type:

 supα∈Λ(Lαu−fα)=0in Ω,u=0on ∂Ω, (1)

where is a compact metric space, and

 Lαv:=Aα:D2v+bα⋅∇v−cαv.

Here, and denote the Hessian and gradient of real-valued function , respectively. The coefficient is assumed to be uniformly elliptic, i.e., there exist constants such that

 ν––|ξ|2≤ξtAα(x)ξ≤¯¯¯ν|ξ|2∀ξ∈Rd, a.e.% in Ω, ∀α∈Λ. (2)

Further, and .

The HJB equations arise in many applications including stochastic optimal control, game theory, and mathematical finance

. In [20, 31], the HJB equations are shown to admit strong solutions under the following Cordes condition.

###### Definition 1 (Cordes condition for (1))

The coefficients satisfy that there exist and such that

 |Aα|2+|bα|2/2λ+(cα/λ)2(trAα+cα/λ)2≤1d+εa.e. in Ω, ∀α∈Λ. (3)

In the past decade, several studies have been taken on the finite element approximation of strong solutions of the HJB equations with Cordes coefficients (3). The first discontinuous Galerkin (DG) method was proposed in , which has been extended to the parabolic HJB equations in . The -interior penalty DG methods were developed in . A mixed method based on the stable finite element Stokes spaces was proposed in . Recently, the (non-Lagrange) finite element method with no stabilization parameter was proposed in , where the element is required to be -continuous at -dimensional subsimplex, e.g., -Hermite family in 2D and -Argyris family [21, 9] in 3D. The above discretizations can be naturally applied to the linear elliptic equations in non-divergence form [30, 17, 22, 12, 34]. Other related topics include the unified analysis of DGFEM and -IPDG , and the adaptivity of -IPDG [7, 19].

For all these discretizations, the discrete well-posedness is analysed under the broken -norm with possible jump terms across the boundary. This, after linearization, leads to the ill-conditioned systems with condition number on quasi-uniform meshes, where represents the mesh size. Due to the similar performance to the discrete system for fourth-order problems, it is conceivable that the linearised system from HJB equations can be effectively solved by the solvers for fourth-order problems, e.g., geometric multigrid [25, 5, 3, 33, 8] or domain decomposition [38, 6]. In , the nonoverlapping domain decomposition preconditioner was studied for the DGFEM discretization of HJB equations.

Traditional geometric multigrid methods depend crucially on the multilevel structures of underlying grids. On unstructured grids, the more user-friendly option is the algebraic multigrid method (AMG) that have been extensively studied for the second-order equations. In , the first biharmonic equation was converted to a Poisson system based on the boundary operator proposed in . Under the framework of auxiliary space preconditioning , Zhang and Xu  proposed a class of optimal solvers based on the auxiliary discretization of mixed form for the fourth-order problems. As a generalization of [25, 24, 26], it works for a variety of conforming and nonconforming finite element discretizations on both convex and nonconvex domains with unstructured triangulation.

The propose of this work is to study the auxiliary space preconditioner to the finite element discretization of HJB equations. More specifically, the numerical scheme for fully nonlinear HJB equations leads to a discrete nonlinear problem that can be solved iteratively by a semi-smooth Newton method [31, 22, 34]. The linear system obtained from the semi-smooth Newton linearization are generally non-symmetric but coercive. To handle the non-symmetry, the existing GMRES theory  will lead to a guaranteed minimum convergence rate with a symmetric FOV-equivalent preconditioner that satisfies (20). The construction of under the auxiliary preconditioning framework follows two steps:

1. Construct appropriate auxiliary spaces and corresponding transfer operators mapping functions from original space to the auxiliary spaces;

2. Devise solvers on auxiliary spaces so that the bounds in (20) are uniform with respect to both and the parameter in the Cordes condition.

Based on the stable decomposition for auxiliary spaces, both additive and multiplicative preconditoners are shown to be efficient and -uniform for the linearised system. Further, the precondtioners only involve the Poisson-like solver which can be efficiently solved by AMG with nearly optimal complexity.

In general cases, the auxiliary space preoconditioner is additive [35, 37, 14], which usually leads to a stable but relatively large condition number in practical applications. The first contribution of this work is the construction and analysis of a multiplicative preconditioner based on the specific structure of auxiliary spaces. Having a coarse subspace, the symmetrized two-level multiplicative precondition was shown to be positive definite provided that the smoother on the fine level has contraction property 

. The condition number estimate of multiplicative precondition at the matrix level can be found in

. In this work, we show that the contracted smoother together with the stable decomposition for auxiliary spaces leads to a robust multiplicative precondititoner, which is also numerical verified with better performance than the additive version.

The parameter in the Cordes condition balances the diffusion and the constant term. We emphasis that this parameter is not involved in the monotonicity constant (2.1), which makes it possible to consider the preconditioner with uniformity on . In this work, we carefully define the norm on the auxiliary space so that the induced preconditioner is uniform with respect to . Although the preconditioner is designed for the finite element approximation, a similar idea can be applied to other discretizations.

The rest of the paper is organized as follows. In Section 2, we establish the notation and state some preliminaries results. In Section 3, we apply the FOV-equivalence preconditioner for the linear system, which can be used to solve non-symmetric systems appearing in applications to the HJB equations. In Section 4, we construct both the additive and multiplicative auxiliary space preconditioners. We also show that the condition numbers of the preconditioned systems are uniformly bounded with the stable decomposition assumption, which is verified in Section 5. Several numerical experiments are presented in Section 6 to illustrate the theoretical results.

For convenience, we use to denote a generic positive constant which may depend on , share regularity of mesh and polynomial degree, but is independent of the mesh size . The notation means . means and .

## 2 Preliminaries

In this section, we first review the strong solutions to the HJB equations (1) under the Cordes condition (3). Then we give a brief statement about the finite element scheme in .

Given an integer , let and be the usual Sobolev spaces, and denote the Sobolev norm and semi-norm. We also denote . For any Hilbert space , we denote for the dual space of , and for the corresponding dual pair. We also denote

as the Euclidian norm for vectors and the Frobenius norm for matrices.

### 2.1 H2 strong solutions to the HJB equations

We now invoke the theory of strong solutions of the HJB equations. In view of the Cordes condition (3), for each , define

 γα:=trAα+cα/λ|Aα|2+|bα|2/2λ+(cα/λ)2. (4)

And for as in (3), define a linear operator by

 Lλu:=Δu−λuu∈H2(Ω). (5)

Next, we define the operator by

 Fγ[u]:=supα∈Λ{γαLαu−γαfα}. (6)

Note that the continuity of data implies . As a consequence, it is readily seen that the HJB equation (1) is equivalent to the problem in , and on . The Cordes condition leads to the following lemma; See [31, Lemma 1] for a proof.

###### Lemma 1 (property of Cordes condition)

Under the Cordes condition (3), for any open set and , , the following inequality holds a.e. in :

 |Fγ[w]−Fγ[v]−Lλz|≤√1−ε√|D2z|2+2λ|∇z|2+λ2z2. (7)

Another key ingredient for the well-posedness of (1) is the Miranda-Talenti estimate stated as follows.

###### Lemma 2 (Miranda-Talenti estimate, [15, 20])

Suppose is a bounded convex domain in . Then, for any ,

 |v|H2(Ω)≤∥Δv∥L2(Ω)≤C|v|H2(Ω), (8)

where the constant depends only on the dimension.

Let the operator be

 ⟨M[w],v⟩:=∫ΩFγ[w]Lλvdx. (9)

By using Miranda-Talenti estimate (8) and Cordes condition (3), one can show the strong monotonicity of ,

 ⟨M[v],v⟩≥(1−√1−ε)∥v∥2λ∀v∈V,

where . Together with the Lipschitz continuity of , the compactness of and the Browder-Minty Theorem [27, Theorem 10.49], one can show the existence and uniqueness of the following problem: Find such that

 ⟨M[u],v⟩=0∀v∈V. (10)

We refer to [31, Theorem 3] for a detailed proof.

### 2.2 C0 finite element approximations of the HJB equations

Let be a conforming shape regular simplicial triangulation of polytope and be the set of all faces of . and . Let be the set of all the nodes of . Here , where is the diameter of . We also denote as the diameter of . For and , we use , respectively , to denote the -inner product over , respectively .

Following , we adopt the -Hermite finite elements () in 2D and -Argyris finite elements in 3D to solve the HJB equations (1). Define the finite element spaces as

1. For , with (cf. Fig. 1),

 Vh:={v∈H10(Ω):v|T∈Pk(T),∀T∈Th, v is C1 at all vertices},
2. For , with (cf. Fig. 2),

 Vh:={v∈H10(Ω):v|T∈Pk(T),∀T∈Th, v is C1 on all edges, v is C2 at all vertices},

where denotes set of the polynomials of degree on .

For each , we define the tangential Laplace operator as follows, where . Let be a orthogonal coordinate system on . Then, for define

 ΔTw=d−1∑i=1∂2∂t2iw.

Next, we define the jump of a vector function on an interior face as follows:

 ⟦v⟧∣∣F:=v+⋅n+|F+v−⋅n−|F,

where and is the unit outward normal vector of , respectively. For scaler function we define

 ⟦w⟧:=w|T+−w|T−.

The following lemma is critical in the design and analysis of finite element approximation of HJB equations (1).

###### Lemma 3 (discrete Miranda-Talenti identity, )

Let be a convex polytopal domain and be a conforming triangulation. For each , it holds that

 ∑T∈Th∥Δvh∥2L2(T)=∑T∈Th∥D2vh∥2L2(T)+2∑F∈Fih⟨⟦∇vh⟧,ΔTvh⟩F.

In light of (9), we define the operator by

 ⟨Mh[w],vh⟩:= ∑T∈Th(Fγ[w],Lλvh)T −(2−√1−ε)∑F∈Fih⟨⟦∇w⟧,ΔTvh−λvh⟩F.

The following finite element scheme is proposed to approximate the solutions to the HJB equations (1): Find : such that

 ⟨Mhuh,vh⟩=0∀vh∈Vh. (11)

We refer to  for the well-posedness and approximation property of discrete systems (11).

### 2.3 Semi-smooth Newton method

It is shown in  that the discretized nonlinear system (11) can be solved by a semi-smooth Newton method, which leads to a sequence of discretized linear systems. We summarized the main ideas on semi-smooth Newton here and refer  for more detials.

Following the discuss in , we define the admissible maximizers set for any ,

 Λ[v]:={\parbox57.0pt$α(⋅):Ω→Λ$\lx@parboxnewlinemeasurable∣∣∣ \parbox200.0pt$α(x)∈argmaxα∈Λ(Aα:D2hv+bα⋅∇v−cαv−fα)$\lx@parboxnewlineforalmostevery$x∈Ω$},

where denotes the broken Hessian of . As shown in [31, Lemma 9 & Theorem 10], the set is not empty for any .

The semi-smooth Newton method is now stated as follows. Start by choosing an initial iterate . Then, for each nonnegative integer , given the previous iterate , choose an . Next the function is defined by ; the functions , , and are defined in a similar way. Then find the solution of the linearised system

 bjλ,h(uj+1h,vh)=∑T∈Th(γαjfαj,Δvh)T∀vh∈Vh, (12)

where the bilinear form is defined by

 bjλ,h(wh,vh):= ∑T∈Th(γαjLαjwh,Lλvh)T −(2−√1−ε)∑F∈Fih⟨⟦∇wh⟧,ΔTvh−λvh⟩F.

Following , we define inner product on as

 (wh,vh)λ,h:=∑T∈Th(D2wh,D2vh)T+2λ(∇wh,∇vh)Ω+λ2(wh,vh)Ω,

and the norm . It is also shown in  that the bilinear forms are uniformly coercive and bounded on with norm , with constants independent of iterates. Since the preconditioners in this work take advantage on the coercivity and boundedness of , we summarize the relevant results in the following lemma (see [34, Lemmas 4.1 & 4.2] for a detailed proof).

###### Lemma 4 (coercivity and boundedness of bilinear form)

For every , we have equationparentequation

 bjλ,h(vh,wh) ≤C∥vh∥λ,h∥wh∥λ,h, bjλ,h(vh,vh) ≥(1−√1−ε)∥vh∥2λ,h.

Here, the constant depends only on , shape regularity of the grid and polynomial degree .

## 3 FOV-equivalent preconditioners for GMRES methods

The preconditioned GMRES (PGMRES) methods are among the most effective iterative methods for non-symmetric linear systems arising from discretizations of PDEs. Our study will start by discussing PGMRES methods in an operator form. Let be a linear operator which may be non-symmetric or indefinite, defined on a finite dimensional space , and be a given functional in its dual space . The linear equation considered here is of the following form

 Gx=g. (14)

Let be an inner product on , and be the preconditioner. The PGMRES method for solving (14) is stated as follows: Begin with an initial gauss and denote the initial residual, the -th steps of PGMRES method seeks such that

 xk=argmin~xk∈Kk(PG,Pr0)+x0∥PG(x−~xk)∥M,

where is the Krylov subspace of dimension generated by and .

In the semi-smooth Newton steps, the discrete linear equations (12) have a common form: Find such that

 bλ,h(uh,vh)=fh(vh)∀vh∈Vh, (15)

where we shall omit to denote independence of the bilinear form and of the right-hand side on the iteration number of the semi-smooth Newton method. Define the operator by

 ⟨Bλ,huh,vh⟩:=bλ,h(uh,vh)∀uh,vh∈Vh, (16)

then the discrete system (15) can be written in an operator form, namely

 Bλ,huh=fh. (17)

Moreover, a general operator is used to denote the preconditioner. Given an inner product , we can estimate the convergence rate of the PGMRES method. It is proved in [10, 28] that if is the -iteration of PGMRES method and is the exact solution of (17), then

 ∥Pλ,hBλ,h(uh−umh)∥Mλ,h∥Pλ,hBλ,h(uh−u0h)∥Mλ,h≤(1−γ2Γ2)m/2,

where

 γ≤(vh,Pλ,hBλ,hvh)Mλ,h(vh,vh)Mλ,h,∥Pλ,hBλ,hvh∥Mλ,h∥vh∥Mλ,h≤Γ∀vh∈Vh. (18)

Therefore, we conclude that as long as we find an operator and a proper inner product such that condition (18) is satisfied with constants and independent of the discretization parameter and the Cordes condition parameter , then is a uniform preconditioner for GMRES method. Such preconditioners are usually referred to as FOV-equivalent preconditioners. In what follows, we always take

 Pλ,h to be an SPD operator,andMλ,h=P−1λ,h.

Next, we give a general principle for constructing . Define an SPD operator by

 ⟨Aλ,hwh,vh⟩:=(wh,vh)λ,h∀wh,vh∈Vh. (19)

Recalling Lemma 4 (coercivity and boundedness of bilinear form), is coercive and bounded on with the inner product . It is therefore that an efficient preconditioner for can also be used as an FOV-preconditioner for the GMRES algorithm applied to , which is shown in the following lemma.

###### Lemma 5 (FOV-equivalent preconditioner)

Let and be the operators defined in (19) and (16), respectively. If an SPD operator satisfies that

 α⟨P−1λ,hvh,vh⟩≤⟨Aλ,hvh,vh⟩≤β⟨P−1λ,hvh,vh⟩∀vh∈Vh, (20)

with constants independent of both and , then is a uniform FOV-equivalent preconditioner of .

###### Proof

From (16), (19) and Lemma 4 (coercivity and boundedness of bilinear form), we see that for any equationparentequation

 (1−√1−ε)⟨Aλ,huh,uh⟩ ≤⟨Bλ,huh,uh⟩ (21a) ⟨Bλ,huh,vh⟩ ≤C⟨Aλ,huh,uh⟩1/2⟨Aλ,hvh,vh⟩1/2, (21b)

where is independent of both and .

Recalling , then for any , we have

 ∥Pλ,hBλ,huh∥P−1λ,h =supvh∈Vh,vh≠0(Pλ,hBλ,huh,vh)P−1λ,h(vh,vh)1/2P−1λ,h ≤β1/2supvh∈Vh,vh≠0⟨Bλ,huh,vh⟩⟨Aλ,hvh,vh⟩1/2    (by (???)) ≤Cβ1/2⟨Aλ,huh,uh⟩1/2                (by (???)) ≤Cβ∥uh∥P−1λ,h,                           (by (???))

which yields the second inequality of (18) with . On the other side, we have

 (uh,Pλ,hBλ,huh)P−1λ,h =⟨Bλ,huh,uh⟩ ≥(1−√1−ε)⟨Aλ,huh,uh⟩        (by (???)) ≥(1−√1−ε)α(uh,uh)P−1λ,h,      (by (???))

which yields the first inequality of (18) with .

## 4 Fast auxiliary space preconditioners

In this section, we construct both additive and multiplicative auxiliary space preconditioners for SPD operator . From Lemma 5 (FOV-equivalent preconditioner), those preconditioners can be applied to the discrete linearised systems (12) arising from each semi-smooth Newton step of solving the HJB equations.

### 4.1 Space decomposition

For the purpose of constructing auxiliary space preconditioners, we give the following space decomposition of as

 Vh=Vh+Π0V0,

where the auxiliary space denotes the continuous piecewise linear element space on with homogeneous Dirichlet boundary condition, and is a linear injective map which will be defined later. As a result, the induced operator is also SPD, and hence we define on . We also introduce a projection by

 ⟨A0P0vh,wh⟩:=(vh,Π0wh)λ,h∀vh∈Vh,wh∈V0.

A direct calculation shows the following identity

 Π′0Aλ,h=A0P0. (22)

#### Smoother and norm on V0

Define the discrete Laplacian operator by

 (−Δhwh,vh)Ω:=(∇wh,∇vh)Ω∀wh,vh∈V0. (23)

Then, the smoother on , denoted by , is defined by

 ⟨R−10uh,vh⟩=(λuh−Δhuh,λvh−Δhvh)Ω∀uh,vh∈V0. (24)

Note that for any given , can be obtained by solving the following two discrete Poisson-like equations equationparentequation

 ⟨fh,vh⟩ =λ(zh,vh)Ω+(∇zh,∇vh)Ω∀vh∈V0, (25a) (zh,wh)Ω =λ(uh,wh)Ω+(∇uh,∇wh)Ω∀wh∈V0. (25b)

It can be shown that the above two equations can be solved within operations, where denotes the number of degrees of freedom. We will give a detailed explanation in Remark 4 (computational complexity). The smoother induces a norm on , i.e., . By using (23) and (24), it is straightforward to show that

 ∥vh∥2R−10= ∥Δhvh∥2L2(Ω)+2λ∥∇vh∥2L2(Ω)+λ2∥vh∥2L2(Ω).

The relationship between and is shown in the following lemma, whose proof is postponed to Section 5.

###### Lemma 6 (spectral equivalence of R0)

Let be the operator defined in (24) and . Then,

 ∥v0∥R−10≃∥v0∥A0∀v0∈V0,

with hidden constants independent of both and .

#### Smoother on Vh

Let denote the Gauss-Seidel smoother for and be the symmetric Gauss-Seidel smoother, i.e.,

 I−¯RhAλ,h=(I−R′hAλ,h)(I−RhAλ,h).

We also define as a norm on . Note that is an SPD operator, thus has the following contraction property

 ∥I−RhAλ,h∥λ,h<1. (26)

#### Transfer operator

We now give the definition of . To this end, we first give some notation on the degrees of freedom of finite element space . We denote the degrees of freedom as

 Nα(φ)=\fintDα∇kα(φ)(t1,…,tkα)

where is the domain of the integral with respect to the degree of freedom, denotes the integral average on . In general, is a subsimplex of the triangulation. When is a point, the average of the integral is reduced to the evaluation on the point. are identical or different unit vectors to denote the direction of the derivative where . When only one direction is involved for the derivative, and the direction is denoted by . Let be the nodal basis function corresponding to . Define , and . We are now ready to give the definition of as follows: For any

 Nα(Π0ph) =1#ωα∑T⊂ωα\fintDα∂n(ph|T)(n⋅tα) (27) when Dα⊂∂Ω and kα=1, Nα(Π0ph) =1#ωα∑T⊂ωαNα(ph|T), else,

where is the unit outer normal vector of . We note that the degrees of freedom corresponding to the second-order derivative vanish since is piecewise linear.

The general theory of auxiliary space preconditioning simplifies the analysis of preconditioners to the verification of the following two key assumptions.

[stable decomposition] There exists a uniform constant independent of both and , such that for any , there exist and satisfy equationparentequation

 v =vh+Π0v0, (28a) ∥vh∥2¯R−1h+∥v0∥2R−10 ≤c20∥v∥2λ,h. (28b)

[boundedness] There exist uniform constants and independent of both and such that equationparentequation

 ∥v0∥A0 ≤c1∥v0∥R−10∀v0∈V0, (29a) ∥vh∥λ,h ≤c2∥vh∥¯R−1h∀vh∈Vh. (29b)

Firstly, we introduce the additive preconditioner as

 Pa:=¯Rh+Π0R0Π′0. (30)

The following theorem plays a fundamental role in the theory of auxiliary space preconditioning .

###### Theorem 1 (spectral equivalence of additive preconditioner)

Let be the preconditioner defined in (30). If Assumption 4.1 (stable decomposition) and Assumption 4.1 (boundedness) hold, then we have

 c−20(vh,vh)λ,h≤(PaAλ,hvh,vh)λ,h≤(c21+c22)(vh,vh)λ,h∀vh∈Vh,

That is, is a uniform spectral equivalence preconditioner of .

In light of the above theorem and Lemma 5 (FOV-equivalent preconditioner), one can see that is a uniform FOV-equivalent preconditioner of as long as Assumption 4.1 (stable decomposition) and Assumption 4.1 (boundedness) are verified. We postpone those verifications to Section 5.

###### Remark 1 (additive preconditioner with Jacobi smoother)

Let be the Jacobi smoother of . From the norm equivalence between and [36, Lemma 4.6], the additive preconditioner

 ~Pa:=D−1h+Π