{'fraction_non_alphanumeric': 0.04258392557924539, 'fraction_numerical': 0.013173860287744843, 'mean_word_length': 4.037107101280559, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Convergence Analysis of a Generalized Full Approximation Storage Scheme for Convex Optimization Problems

Full Approximation Scheme (FAS) is a widely used multigrid method for nonlinear problems. In this paper, a new framework to analyze FAS for convex optimization problems is developed. FAS can be recast as an inexact version of nonlinear multigrid methods based on space decomposition and subspace correction. The local problem in each subspace can be simplified to be linear and one gradient decent iteration is enough to ensure a linear convergence of FAS.


INTRODUCTION
Since most real-world applications are inherently nonlinear in nature, nonlinear problems are of interest to mathematicians, physicists, biologists, engineers and many other scientists. Efficiently solving nonlinear problems is very important to many practical problems and is at the core of many numerical simulations. In this paper, we consider solving nonlinear equations arising from the minimization of a convex function in the abstract Hilbert space setting.
A traditional approach for solving nonlinear equation is the well-known Newton-Raphson method. Basically, Newton's method iteratively finds the approximate solution by linearizing the problem near the current iterate. A linear system (the Jacobian system) needs to be solved at each Newton's iteration, and linear multigrid (MG) methods are sometimes used as a solver. Practically, each linear problem can be approximately inverted by applying a few multigrid iterations. But, in this case, the quadratic rate of convergence may be sacrificed.
One alternative to Newton's method for solving nonlinear PDE is the nonlinear multigrid method, better known as the full approximation storage (FAS) scheme. This methods, developed by Brandt [1] in the 80's, often converges linearly and with optimal complexity. Recall that the success of multigrid methods relies on two ingredients: 1) high frequency components of the error will be damped by smoothers; and 2) low frequency components of the error can be approximated well on a coarse grid. The smoother used in FAS is usually the nonlinear Gauss-Seidel smoother, which solves many small-sized (typically 1-D) nonlinear problems on small patches of the mesh. For the coarse grid problem, the FAS method uses the full approximation rather than the standard defect, which makes it essentially different from linear MG methods. Due to its high efficiency, FAS method has been applied to many nonlinear PDE problems, such as in [9,11,18,22,10,13].
Although FAS is quite successful in practice, its theoretical analysis is limited. In [8], Hackbusch considered nonlinear MG methods for general nonlinear problems. By imposing conditions on the nonlinear operators and their derivatives, together with standard smoothing and approximation properties, he was able to show that the FAS converges in a sufficiently small neighborhood of the solution on a fine enough mesh. Moreover, the number of smoothing steps needs to be sufficiently large, and at least the W-cycle should be used. Later in [16,17], Reusken considered FAS for a class of second order elliptic boundary value problems with mild nonlinearity. Within this nice class of nonlinear problems, he was able to show the convergence of FAS under weaker assumptions on the nonlinear operators. We want to mention that the proofs in their work are based on the linearization of the FAS iterations, and the rate of convergence is in some sense local. For example, in [17], Reusken showed that the V-cycle FAS converges locally in a ball with radius shrinking from coarse to fine levels.
In this paper we consider a special class of nonlinear equations that can be viewed as Euler equations of certain convex objective functions. The convergence of MG methods for convex optimization problems has been studied in [20,21] under the framework of subspace correction methods [23]. In [21], Tai and Xu considered some unconstrained convex optimization problems and developed global and uniform convergence estimates for a class of subspace correction iterative methods. Their approach is based on an abstract space decomposition which is assumed to satisfy the so-called stable decomposition property and strengthened Cauchy Schwarz inequality. We want to point out that in each subspace, the original object function is used, which is only defined on the finest level, and the local problem should be solved exactly which is more expensive than what is required in the FAS scheme.
We shall borrow the theoretical framework established in [21] to analyze the FAS method. Different from the subspace correction method considered in [21] in which an exact subspace solver is used, we recast FAS as a subspace correction method with an inexact subspace solver, which reduces the computational cost, and we establish its global and uniform convergence in the framework of subspace corrections. We first show that, with a one dimensional line search approach, the FAS converges globally and uniformly under the standard assumption of the space decomposition. In addition, we borrow some techniques from the optimization literature [15] in order to properly handle the inexactness of the local solver used in FAS. Basically, we introduce a fixed step size to guarantee that the objective function is decreasing globally. For the analysis of original FAS widely used in practice, we impose an approximation property of the subspace problems and show that FAS converges globally and uniformly. We emphasize that our work represents not only a theoretical advance for the convergence analysis of the FAS scheme, but also is algorithmically simpler than the original FAS. We show that, both theoretically and numerically, each local nonlinear problem can be approximated by a linear problem, and, consequently, the computational cost is reduced significantly.
The paper is organized as follows. In Section 2, we present the optimization problem and nonlinear problem under consideration and assumptions on the problems as well as the space decomposition. Successive subspace optimization methods in a general space decomposition setting is recalled in Section 3. And it convergence analysis based on slightly weaker assumptions comparing with [21] is presented in the same section. The main global and uniform convergence analysis for FAS with line search and without line search are derived in Section 4 and 5. In Section 6, some applications are considered.

PROBLEM AND ASSUMPTIONS
Given an energy E(v) defined on a Hilbert space V equipped with inner product (·, ·) V and norm · V , we consider the minimization problem:

2.1.
Assumptions on the Energy. We assume that the energy functional E( · ) : V → R is the continuous linear functional equal to the first Fréchet derivative at v. We further impose the following assumptions on its derivative: (E1) (Strong convexity:) There is a constant µ > 0 such that where · , · is the dual pairing between V and V.
As E is monotone, its norm E (v) V will typically be unbounded as v V → ∞. However, restricted to a bounded set, we assume that the change in E can be controlled.
(E2) (Lipschitz continuity of the first order derivative:) For fixed u 0 ∈ V, there exists a constant L such that, for all v, w ∈ B : If E satisfies (E1) -(E2), it follows [6, p.35] that there is a unique element u ∈ V with the property that and this minimizer further satisfies Euler's equation The strong convexity and the Lipschitz continuity imply the following equivalence of the norm and functional.
Furthermore the lower bound holds for all v, w ∈ V. Now, we consider the relation between the energy and the norm centered at the minimizer. The following equivalence can be easily proved using Taylor's theorem with integral remainder; see, e.g. [15].
Furthermore the lower bound holds for all v, w ∈ V. In addition, suppose u ∈ B is the minimizer of E, then for all w ∈ B, Again the lower bound holds for all w ∈ V.
Based on Assumption (E1), the upper bound can be replaced by a norm of the gradient. Since the proof is less standard, we include it here. Lemma 2.3. Suppose that E satisfies Assumption (E1) and u ∈ V is the minimizer of E; then for all v ∈ V, Proof. Fix the point v ∈ V. Now, for any w ∈ V, using the lower bound of (5), we have .
Then (7) is obtained by letting w = u in the above inequality.
We shall often use the following simple variant of Lemma 2.2.
Then J : W → R is differentiable, stongly convex, and there exists a unique element η ∈ W ∩ B such that η is a minimizer of J and Furthermore, for all w ∈ B ∩ W, The ratio L/µ sometimes is called the condition number of the nonlinear function E; see [15, page 63]. The rate of convergence of iterative methods for solving (1) usually depends on the condition number. Here we assume L/µ is uniformly bounded, i.e., it is well conditioned, as long as we remain in B. Then the Riesz map R : V → V, defined as (Rf, v) V := f, v , for any f ∈ V and v ∈ V, can be used as a preconditioner and the corresponding preconditioned gradient descent method will converge [7].
Implementing preconditioned gradient decent methods in V requires the computation of the Riesz map R which is equivalent to inverting a symmetric positive definite (SPD) operator with size dim V. In the following, we shall provide optimization methods that only require computing inverses with much smaller sizes. Of course we can also use multilevel methods to compute R and use steepest descent, nonlinear conjugate gradient, or Newton method as the outer iteration.

Assumptions on the Space Decomposition. Suppose that
is a space decomposition of V. We shall use the following assumptions on the space decomposition.
where C A > 0 is independent of v. (SS2) (Strengthened Cauchy Schwarz inequality:) There is a constant C S > 0, such that, for any

SUCCESSIVE SUBSPACE OPTIMIZATION METHODS
For k ≥ 0 and a given approximate solution u k ∈ V, one step of the Successive Subspace Optimization (SSO) method [19] is as follows: We aim to prove a linear reduction of the energy difference for one iteration of the SSO algorithm: where u is the minimizer of E and u k+1 = SSO(u k ), with a contraction factor ρ ∈ (0, 1). Ideally ρ is independent of the size of the problem. The algorithm and convergence theory has been developed in [19,21] for a nonlinear and convex energy in Banach spaces. For completeness, we include a simplified version for Hilbert space here. Let d k := E(u k ) − E(u) and δ k := E(u k ) − E(u k+1 ). The quantity d k is the difference between the current energy and the minimum energy, and δ k is the Algorithm 1: Successive Subspace Optimization Method. energy decrease associated to the k + 1 st iteration. They are connected by the trivial identity δ k = d k − d k+1 . We shall present our proof based on the following two inequalities.
Lower Bound. There exists a positive constant C L such that for k = 0, 1, 2, . . .
Upper Bound. There exists a positive constant C U such that for k = 0, 1, 2, . . .
Theorem 3.1. Assume that the discrete lower bound and upper bound hold with constants C L and C U , respectively. We then have Proof. The proof is straightforward by assumptions and rearrangement of the following inequality Verifying the lower bound is relatively easy since E is convex. Solving the convex optimization problem in each subspace will definitely decrease the energy.
Theorem 3.2. Let u k be the k-th iteration and u k+1 = SSO(u k ). If E is strongly convex in the sense of satisfying (E1), then Proof. Recalling Lemma 2.4, we observe that J i (defined in Algorithm 1) is strictly convex over V i and is Fréchet differentiable, as it inherits the structure of E. It follows that But the object on the left-hand-side is simply a directional derivative of the full energy, and it is easy to see that and consequently Remark 3.3. In view of (9), the convexity of E can be relaxed to the local convexity of J i in each subspace V i . Namely we may have a non-convex energy E which restricted to each subspace is convex and the lower bound still holds. For example, the energy used in the Optimal Delaunay Triangulation (ODT) [5] is nonconvex globally. But restricted to one vertex, it is convex, and the corresponding 1-D optimization problem has a closed form solution, which is known as ODT mesh smoothing [3]. Theorem 3.2 gurantees the energy decreasing property of ODT mesh smoothing.
The upper bound is more complicated and relies on the assumptions about the decomposition of spaces. The result is given in the following theorem.
Theorem 3.4. Let u k+1 be the k + 1-th iteration in SSO algorithm. Suppose that the space decomposition satisfies Assumptions (SS1) and (SS2) and the energy E satisfies Assumption (E1), then we have Proof. Using Lemma 2.3, with the choice v = u k+1 in (7), we have For any w ∈ V, we choose a stable decomposition w = N i=1 w i , then Here we use the fact that we solve the minimization problem on each subspace exactly and the energy decreases, therefore, v j ∈ B for all j and which finishes the proof.
Based on lower bound given in Theorem 3.2 and upper bound given in Theorem 3.4, we can conclude the convergence of SSO. Comparing with the results in [21], we use slightly weaker assumptions and the constant C U seems to be slightly better.
Corollary 3.5. Let u k be the k-th iteration and u k+1 = SSO(u k ). Suppose that the space decomposition satisfies Assumptions (SS1) and (SS2) and the energy E satisfies Assumption (E1), then we have

FAS WITH LINE SEARCH
In the SSO method, cf. Algorithm 1, we need to solve the optimization problem min w∈V i E(v i−1 + w) in each subspace, which requires evaluation of the global energy E and its derivative E in the space V. Although the size of the optimization problem is reduced to dim V i , such evaluations are still in the original space of size dim V, which is still expensive.
More precisely, denote by I i : V i → V the natural inclusion and R i = I i : V → V i the natural restriction of functionals. We need to evaluate the gradient and its inverse, if Newton's method is used, several times. This is practical only if the natural inclusion I i is efficient to realize, e.g., a one-dimensional subspace generated by one basis function of V and the resulting method is the so-called non-linear Gauss-Sedel iteration.
For nonlinear elliptic partial differential equations discretized by finite difference or finite element methods, a stable decomposition V = J k=1 V k is usually associated to multi-level grids and dim V k = γ k dim V, γ < 1. For each subspace V k , the prolongation should be applied all the way to the whole space V. Even with this stable and multilevel decomposition, the computational cost will be N log N , In contrast, optimal multilevel methods only requires the prolongation and restriction among consecutive levels which will make the evaluation of gradient and Hessian easier and results in optimal complexity, i.e., O(N ) algorithms. The most commonly used nonlinear multigrid is called the full approximation scheme (FAS) [1,2] because the problem in the coarse grid is solved for the full approximation rather than the correction.
We shall view FAS as an SSO iteration with an inexact local solver. We will utilize a locally defined energy E i in V i and solve a perturbed optimization problem. In addition to prolongation and restriction operators, we also need a projection op- In Step 4, the orthogonalization step, we perform a line search to find the optimal step size which still requires the evaluation of some of the "fine level" functions The computational cost is reduced comparing with evaluation at v i−1 + w for multiple w ∈ V i . Algorithm 2 is an intermediate step towards the convergences proof of original FAS. In the next section, we shall further relax this constraint by considering a simple choice of step size.
Again verifying the lower bound is relatively easy due to the line search and the convexity of E.
Solve the subspace problem residual problem: Compute the search direction: Orthogonalize the subspace correction via line search: Apply the subspace correction: Proof. We apply a similar technique as in the proof of Theorem 3.2. Due to the line search, we still have an orthogonality property that can be utilized, namely, Then, applying Lemma 2.4, with the subspace W = span{s i }, and noting that The upper bound is more complicated. We always assume E i : V i → R is Fréchet differentiable for all points v ∈ V i . In addition, we introduce the following assumptions on the local energy which is just the local version of (E1)-(E2).
(E3) (Strong convexity of the first order derivative:) There exists a constant µ i such that for all v, w ∈ V i As Theorem 4.1 implies the energy is always decreasing and iterates will remain in the ball B, but the search region, e.g., ξ i + s i may not be in B. To be able to use Lipschitz continuity, we introduce a larger domain We then introduce a Lipschitz continuity of E i with respect to this larger domain: (E4) (Lipschitz continuity of the first order derivative:) There exists a constant L i such that Later, we will show that ξ i + s i ∈ B i so that we can still take advantage of the Lipschitz continuity of E i (·) in our analysis. Notice that the Lipschitz continuity of E i is imposed for the ballB which is related to B used in (E2) and there is no relation between E and E i is explicitly imposed so far. Indeed E and E i are just related through the upper and lower bound of the first derivatives and norms. In general, based on the assumptions (E3) and (E4), we have the following lemma which gives the similar results as Lemma 2.1 and 2.2.
In order to better understanding the choice of the step size, we introduce a scalar function Inherit from E, f is also strong convex. In Step 2 of the local problem, we are not minimizing an approximated energy E i , i.e., not solving E i (Q i v i−1 +s i ) = 0. Instead a so-called τ -perturbation is added to the right hand side. We now give our motivation for adding the τperturbation. Lemma 4.3. Let s i be computed as in Algorithm 2 and suppose that E i satisfies Assumption (E3). Then s i is a decent direction in the sense that Proof. The local problem (11) can be rewritten as follows: find η i ∈ V i s.t.
Here recall that ξ i = Q i v i−1 ∈ V i and η i = ξ i + s i . Choosing w = s i and using the strong convexity of E i , we obtained the inequality In terms of the scalar function f (α), Lemma 4.3 implies f (0) < 0. As f is continuous, we conclude that the optimal point is positive, α * i > 0, and for all α ∈ (0, α * ], f (α) < f (0). Inherited from E, f is strongly convex, and f is Lipschitz continuous with a Lipschitz constant that is scaled by s i 2 . More precisely, we have the following inequalities on f .
Furthermore, f is Lipschitz in the following sense: for all 0 < α, β ≤ α L,i , The proof is based on the following identity Then, by Assumption (E1), To use the Lipschitz inequality involving E , we need to ensure that the points of evaluation are inside the set B, which imposes an upper bound on α and β. As . This implies v i−1 + αs i ∈ B, since the energy is decreased for such values of α. We now estimate α o L,i . As f (α * i ) = 0 and f is Lipschitz in (0, α o L,i ), we have the bound On the other hand, as in Theorem 4.1, The desired bound on the interval (0, α o L,i ), it is also Lipschitz with the same constant on the smaller interval (0, α L,i ) ⊆ (0, α o L,i ). The proof is complete. To use the Lipschitz continuity of E i , we require ξ i + s i ∈ B i = Q iB which will be proved by a lower bound of the optimal step size.
The lower bound is obtained by the strong convexity of E i and Lipschitz continuity of E : . Note that v i−1 + ε i ∈ B by Lemma 4.4 so that we can use Lipschitz continuity of E .
Next we show the norm of s i is bounded and thus ξ i + s i ∈ B i .
To start, we know that v i−1 ∈ B; so by the definition ofB cf. (16), it suffices to prove that s i 2 V ≤ χ. By Theorem 4.1 and Lemma 4.5, we have and the result is proven.
Remark 4.7. If necessary, we can improve our bound above to some degree. In particular, we know that v i−1 + α L,i s i ∈ B. So we only need to estimate the size of Thus, if 1 > α L,i , we only need to estimate the V-norm of (1 − α L,i )s i : Next, we present a refined upper bound on the step size α * i . Lemma 4.8. Assume the energy E satisfies the strong convexity Assumption (E1) and the local energy E i satisfies the Lipschitz assumptions (E4), then we have the upper bound Proof. The upper bound can be proved using the strong convexity of E and the Lipschitz continuity of E i , with the knowledge that ξ i , ξ i + s i ∈ B i : Now we are ready to prove the upper bound.
Proof. Note, for any w ∈ V, we choose a stable decomposition w = N i=1 w i , then I 1 and I 2 defined in respective order. I 1 can be estimated in exactly the same way as in Theorem 3.4. Therefore, In the last step, we have used the relation s i = α −1 i ε i and the lower bound of α i . Then use inequality (7) in Lemma 2.3 with v = u k+1 and following the same procedure as in Theorem 3.4, we complete the proof.
, where R i is the Riesz map V i → V i and its realization is the inverse of an SPD matrix of size dim V i . Solving a linear local problem will save the computational cost of FAS dramatically; see Section 6. In this setting, FAS coincises with the coordinate descent methods analyzed in [14].
Corollary 4.11. In addition to the hypotheses of the last theorem, let us assume that E i is quadratic, chosen as in (18). Then, Using Theorem 4.1, 4.9, and 3.1, we obtain the following linear convergence result.
Furthermore if E i is quadratic, chosen as in (18), then

FAS WITHOUT LINE SEARCH
In this section, we consider the FAS algorithm without line search. Basically, it is similar with Algorithm 2 without the line search step but with a fixed step size α i as shown below. Notice that, now, there is no need to repeatedly evaluate E and its derivatives in the subspace but only compute R i E (v i−1 ) once for the local problem. We shall also consider the original FAS which corresponds to α i = 1 and prove its convergence based on an approximation property.
Compute the subspace residual functional: Solve the subspace problem residual problem: Find η i ∈ V i , such that Compute the search direction and the step size: Apply the subspace correction:

Algorithm 3: FAS algorithm without line search
Recall the scalar function f (α) : The optimal step size α * i = arg min α f (α). Our choice of α i = arg min α f Q (α). Lemma 5.1. Assume the energy E satisfies the Lipschitz continuity assumption (E2) and the local energy E i satisfies the strong convexity assumptions (E3), then Proof. The lower bound is obtained by the definition of α i and Lemma 4.3. To prove the upper bound, we notice that f (α * i ) = 0 and thus Theorem 5.2. Let u k be the k-th iteration and u k+1 = FAS(u k ). Suppose that E satisfies Assumption (E1)-(E2) and the local energy E i is strong convexity with assumptions (E3), we have In the last step, we have used the definition of α i and this completes the proof.
Since α i has the same lower bound of α * i , we can derive the upper bound exactly the same as the proof of Theorem 4.9 by replacing ξ i = α * i s i by α i s i . Thus, we only state the theory below and the proof is omitted.
. We summarize the linear convergence result below.
Corollary 5.4. Let u k be the k-th iteration and u k+1 = FAS(u k ). Suppose that the space decomposition satisfies Assumptions (SS1) and (SS2), the energy E satisfies Assumption (E1)-(E2), and the energy E i satisfies Assumption (E3)-(E4), then we have The Lipschitz constant L is used in the step size α i which can be replaced by a local Lipschitz constant for the scalar function f (α) for α ∈ (0, α * i ) and popular line search algorithms can be used.
Notice that the original FAS does not contain the computation of step size α i using the current residual and search direction. It is simply α i = 1. Our choice of step size is motivated by the choice of step size in the gradient decent method [15]. We shall prove α i = 1 is also allowed provide the following approximation property is satisfied.
(AP) Both E and E i are twice Fréchet differentiable. Furthermore there exists a constant < µ/2 so that for all w ∈B and all For quadratic energy, R i E I i is the matrix formed by the triple product and E i is the matrix obtained using the bilinear form restricted to the subspace V i . They should be close in certain norm.
Theorem 5.5. Let u k be the k-th iteration and u k+1 = FAS(u k ) with local step size α i = 1. Suppose that E satisfies Assumption (E1) and Approximation property (AP) holds with < µ/2, we have If in addition (E2) holds, then the lower bound holds with C L = L/2 − .
We shall use the auxiliary function f Q (α) to get an improved lower bound. First we can rewrite the quadratic function as where α i is given in (22). Substitute α = 1 and expand the term α i s i Due to (24), we conclude v i ∈ B and thus by Assumption (E2), , which leads to the estimate Hence by (26), we get Note that which completes the proof.
The upper bound for α i = 1 is easy as no need to have a upper bound of the step size.
Theorem 5.6. Let u k be the k-th iteration and u k+1 = FAS(u k ) with local step size α i = 1. Suppose the space decomposition satisfies (SS1) and (SS2), the energy E satisfies (E1) -(E2), and assumption (AP) holds. Then we have the upper bound The first term is bounded as before. Therefore, For the second term, we insert which is zero in V i and use assumption (AP) to get Remark 5.7. When (AP) holds, we can chose the local quadratic energy Then the local problem becomes one Newton's iteration in subspace V i When (AP) holds, one can easily show (E3) and (E4) holds with constant µ − and L + by a simple perturbation argument. In this setting, the block Newton's method proposed in [12] can be interpreted as a FAS with space decomposition. We will investigate the randomized version somewhere else.
Corollary 5.8. Let u k be the k-th iteration and u k+1 = FAS(u k ). Suppose that the space decomposition satisfies Assumptions (SS1) and (SS2), the energy E satisfies Assumption (E1)-(E2), and the energy E i satisfies Assumption (AP) with < µ/2, then we have .

APPLICATION AND NUMERICAL EXPERIMENTS
In this section we shall apply our theory to a model nonlinear problem with polynomial nonlinearity and provide numerical examples to illustrate the efficiency of a variant of FAS with a local quadratic energy.
We then verify that E satisfies our assumptions. The space V = H 1 0 (Ω) with norm ∇v for v ∈ V. Obviously, E is strictly convex with constant ε 2 . By Hölder inequality and Sobolev embedding theorem H 1 0 (Ω) → L p (Ω), assuming 2 ≤ p < ∞, for d = 2, and 2 ≤ p ≤ 6, for d = 3, we have the bound Without loss of generality, which implies By the mean value theorem, there exists a z = tv If w, v ∈ B, then by Minkowski inequality and (30) z L p ≤ C 1 , and thus Namely (E2) holds with L := (p − 1)C(Ω, p)C p−2 1 + ε 2 . Now, suppose that Ω ⊂ R 2 is a polygonal domain and T H is a conforming triangulation of Ω. Let T h be the triangulation obtained by quadri-secting T H . Specifically, if K i ∈ T h is one of the four daughter triangles (i = 1, · · · , 4) obtained by quadri-secting K ∈ T H -that is by connecting the midpoints of K -then h K i = H K /2, i = 1, · · · , 4. Define With a similar definition for S H . Then, S H ⊂ S h , and the containment is proper. We shall consider the minimization of energy E restricted to S h which is a subspace of and thus now V = S h with norm |v| 1 = ∇v . Notice that (E1) and (E2) still holds as S h ⊂ H 1 0 (Ω). Next we give a two-level space decomposition of V as follows. Let N = {x i } N i=1 ⊂ R 2 be the set of interior nodes of T h and define the Lagrange nodal basis B h is a bona fide basis for S h , and we may use the following decomposition The fact that this forms a stable decomposition is well known, i.e., Assumption (SS1) holds and will be briefly proved below.
Let Q H be the L 2 -projection to S H . For any v ∈ S h letṽ = (I − Q H )v and v = N i=1ṽ i be the nodal decomposition. By the approximation property of Q H on quasi-uniform grids and the stability of nodal decomposition in the L 2 -norm, we have By the H 1 -stability of Q H on quasi-uniform grids, we also have |Q H v| 1 |v| 1 . In conclusion, Assumption (SS1) holds.
We now show that assumption (SS2) holds. Suppose that for some z i,j ∈ S h between w i,j ∈ B and w i,j + u i ∈ B which satisfies the bound (30). The functions u i , 1 ≤ i, j ≤ N are local. In particular, the support of v i , denoted S i is exactly equal to the union of those triangles that have the node x i as a vertex. Let In a similar fashion, using the fact that ∇z i,j ≤ C 1 and the embedding theorem, where C 1 is independent of N , we can prove that N i,j=0 with a constant independent of N . Therefore, there is a C S such that where C S is independent of N , and assumption (SS2) is confirmed.

Numerical Examples.
In this subsection, we present some numerical results for the nonlinear problems described in the previous two subsections to illustrate our theoretical results. For both problems, we will use linear finite element method to discretize and use FAS to solve the discretized nonlinear equations. Our algorithms are implemented in Matlab based on the software package iFEM [4]. The numerical experiments are conducted on a System76 Galago with an Intel Core i7-8550U CPU and 32GB RAM.
We mainly focus on three different implementations of FAS (Algorithm 3), based on different choices of space decomposition and local energy. Geometric multigrid setting is considered here, i.e., we have a set of uniformly refined meshes and nested linear finite element spaces V 1 ⊂ V 2 ⊂ · · · ⊂ V J , where V = span{φ 1 , φ 2 , · · · , φ N } with φ i being the nodal linear finite element basis on level .
(1) First implementation is the original FAS. We consider standard multilevel nodal based space decomposition V = J =1 N i=1 φ i and the local energy E i is defined as the restriction of E on the subspace span{φ i }. Newton's method is used to solve the local nonlinear problem and we set the tolerance to be 10 −10 and at most 100 iterations are allowed (in general, less than 5 iterations are needed for solving the local problems in all of our numerical tests). We use small tolerance to make sure each local problem is solved exactly in order to be consistent with our theoretical analysis.
(2) Second implementation is a simplified version of FAS based on Remark 4.10 and we refer it as "FASq1". We again consider the multilevel nodal based (18) is used, which leads to solve a linear problem locally.
(3) Third implementation is a further simplified version and we refer it as "FASq2". In this case, we use space decomposition V = J =1 V and consider quadratic energy (18). As mentioned in Remark 4.10, this involves the Riesz map which can be computed by inverting an SPD matrix defined on V . For our example, this is equivalent to solving a discrete Laplacian matrix on each level, which is still expensive. Therefore, we solve the discrete Laplacian matrix approximately by just applying on step of symmetric Gausse Seidel method. In all of our numerical experiments, we use Newton's method to solve the nonlinear problem on the coarest level. We use 10 −10 as the tolerance and maximal number of iterations is 100, which means that the coarse problem is solved exactly. Moreover, we use α i = 1 in the tests to make sure our implementation is simple and practical. The overall stopping criterion of FAS is 10 −10 .   Table 1, 2, and 3, we report the numerical results of FAS, FASq1, and FASq2, respectively. Here, we fix the finest mesh size h = 1/64 but change p and ε to adjust the nonlinearity. In this case, bigger p and/or smaller ε lead to stronger nonlinearity.
Number of iterations and convergence rates (in the parenthesis) are listed in Table 1, 2, and 3. "-" means that it takes more than 500 iterations and we simply stop the program after 500 iterations. As we can see, FAS is the most robust one and converges for all the choices of our parameters. The number of iterations are quite stable, ranging from 10 − 21 iterations, and the convergence rate is about 0.2. For FASq1 and FASq2, both implementations perform well when p is relatively small and/or ε is relative big. We can clearly see that the number of iterations grows when p gets larger or ε gets smaller. Both implementations fail to converge within 500 iterations when nonlinearity is strong, while FASq1 seems to be slightly more robust than FASq2 since it converges for slight larger set of parameters. This is consistent with our theoretical results because when the problem becomes more nonlinear, quadratic energy is not a good approximation of the original energy E any more. However, the advantage of using quadratic energy on local subspaces is that we only need to solve linear problems locally, which could save computational cost considerably.
Next, we compare the CPU time of FAS and FASq2. The reason we choose FASq2 to compare is that FASq2 only involves symmetric Gauss-Seidel smoother on each level, which basically has the same cost as multigrid method for soling linear problems. This could dramatically improve the computational complexity for solving our model problem (28) and the results are shown in Table 4. In Table 4, we fix ε = 1 and p = 6 and change h. As we can see, for these choices of p and ε, the quadratic energy provides a good approximation of the global energy restricted to the subspace, therefore, the number of iterations of FASq2 is similar with the number of iterations of FAS and remain robust with respect to the mesh size h. The CPU time of FAS grows faster than linear, which is due to the inefficiency of large for loops in MATLAB.
In contrast, FASq2 is significantly faster than FAS and scales linearly. This demonstrate that, when nonlinearity is mild, we can use a simple quadratic energy and save computational cost.
On the other hand, we want to point out that FAS is more robust than FASq2 as shown before. We have also tested the quadratic energy defined by the Hessian at the previous iteration c.f., (27), which is more or less equivalent to using one approximated Newton's iteration, and the results are similar. Therefore, in practice, we should consider the trade-off between robustness and efficiency in order to decide which kind of local energy should be used on each subspace.