ON CONVERGENCE OF THE ADDITIVE ... - Semantic Scholar
SIAM J. NUMER. ANAL. Vol. 43, No. 5, pp. 1850–1871
c 2005 Society for Industrial and Applied Mathematics
ON CONVERGENCE OF THE ADDITIVE SCHWARZ PRECONDITIONED INEXACT NEWTON METHOD∗ HENG-BIN AN† Abstract. The additive Schwarz preconditioned inexact Newton (ASPIN) method was recently introduced [X.-C. Cai and D. E. Keyes, SIAM J. Sci. Comput., 24 (2002), pp. 183–200] to solve the systems of nonlinear equations with nonbalanced nonlinearities. Although the ASPIN method has successfully been used to solve some difficult nonlinear equations, its convergence property has not been studied since it was proposed. In this paper, the convergence property of the ASPIN method is studied, and the obtained result shows that this method is locally convergent. Furthermore, the convergence rate for the ASPIN method is discussed and the obtained result is similar to that of the inexact Newton method. Key words. nonlinear systems, inexact Newton methods, nonlinear preconditioning, additive Schwarz, local convergence, convergence rate AMS subject classifications. 65H10, 65N12, 65N22 DOI. 10.1137/040611653
1. Introduction. Consider the nonlinear system of equations (1.1)
F (u) = 0,
where F : Rn → Rn is a continuously differentiable function. For convenience of discussion, let F = (F1 , F2 , . . . , Fn )T , u = (u1 , u2 , . . . , un )T , and J(u) = F (u). A numerical solution for (1.1) is often required in many scientific and engineering computing areas such as the discretization of nonlinear partial differential equations; see [7, 16]. The inexact Newton method [8] is one of the most important and effective tools for solving such systems, in particular, when the problem is large and sparse. In applications, some global strategies, such as linesearch or trust region techniques, are often needed because the inexact Newton method is locally convergent [1, 2, 3, 4, 12]. In particular, if the linesearch backtracking technique is augmented in the inexact Newton method, then the inexact Newton with backtracking (INB) method is obtained [12, 13, 19]. This method is more robust and it can be briefly described here. Suppose u(0) is a given initial guess and let u(k) be the current approximate solution; the next approximate solution u(k+1) can be obtained through the following steps. Algorithm 1.1 (INB [12]). 1. Inexactly solve the system (1.2) J u(k) p = −F u(k) , and obtain an inexact Newton direction p(k) such that (1.3) F u(k) + J u(k) p(k) ≤ ηk F u(k) . 2. Compute the new approximate solution u(k+1) = u(k) + λk p(k) . ∗ Received by the editors June 24, 2004; accepted for publication (in revised form) February 11, 2005; published electronically November 22, 2005. http://www.siam.org/journals/sinum/43-5/61165.html † State Key Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, P.O. Box 8009, Beijing 100088, People’s Republic of China ([email protected]).
1850
ON CONVERGENCE OF THE ASPIN METHOD
1851
Here ηk ∈ [0, 1) is the forcing term that controls how accurately system (1.2) should be solved, and p(k) is the inexact Newton direction of F at u(k) . Step 2 in Algorithm 1.1 is a linesearch procedure that is used to find a satisfied step factor λk ∈ (0, 1] and then form the next approximate solution. Usually, we use linear iterative methods, such as the classical splitting method or the modern Krylov subspace method, to inexactly solve system (1.2). Thus, the inexact Newton method is an inner-outer iterative method. In particular, when the Krylov subspace method is used in an inner iteration, we obtain the Newton–Krylov subspace method, which has been used successfully in many areas [1, 2, 3, 4, 16]. Although the inexact Newton method works very well for most nonlinear equations, this may often fail when it is used to solve some difficult problems. Many numerical experiments show that most failed cases in the inexact Newton method result from stagnation, particularly when it is used to solve some problems with nonbalanced nonlinearities [1, 6]. Usually, the stagnation phenomenon is caused by the lack of a good initial guess and/or problematic regions such as boundary layers, singularities in the domain, and/or multiphysics domain, etc. See [17]. Considering this, Cai and Keyes [6] recently proposed a nonlinearly preconditioned inexact Newton algorithm: first convert system (1.1) into another nonlinear system F(u) = 0 such that the two systems have the same solution u∗ ∈ Rn ; then use Algorithm 1.1 to solve F(u) = 0. F and F may have completely different forms, but they must have the same solution. Usually, F has more uniform nonlinearities, so it is relatively easy to solve. In [6], an especially preconditioned case, where F is obtained by the single-level nonlinear additive Schwarz method, is discussed in detail. The corresponding method is the additive Schwarz preconditioned inexact Newton (ASPIN) method. Numerical results in [6] show that the ASPIN method can solve some difficult problems where the traditional inexact Newton method fails. Although the ASPIN method has better numerical results than the traditional inexact Newton method, it is unfortunate that until now the convergence property for the ASPIN method has not been given much importance except for some preliminary convergence analysis in the context of semilinear PDEs in paper [17]. In this paper, we show that the ASPIN method is locally convergent; thus we give theoretical support for the ASPIN method. Moreover, we will discuss the convergence rate of the ASPIN method. The rest of the paper is organized as follows. In section 2, we briefly discuss the ASPIN method, and some of its properties are listed. In section 3, we show that the ASPIN method is locally convergent, and its convergence rate is discussed in section 4. Finally, in section 5, some brief conclusions are given. 2. The ASPIN method. Assume that F (u∗ ) = 0 and J(u∗ ) is invertible. To find the solution u∗ of system (1.1), the ASPIN method solves another nonlinear system F(u) = 0, which is obtained from (1.1) through the additive Schwarz preconditioning technique. Specifically, the ASPIN method can be described as follows. Let S = {1, 2, . . . , n} be an index set, i.e., one integer corresponds to each ui and Fi . Assume that S has a partition {S1 , S2 , . . . , SN } such that N i=1
Si = S
and Si ⊂ S.
1852
HENG-BIN AN
Here the subsets may overlap. Let ni = |Si | be the dimension of Si ; then N
ni ≥ n.
i=1
Assume that Si = {i1 , i2 , . . . , ini }, where i1 < i2 < · · · < ini . For i = 1, 2, . . . , N , define matrices Ei ∈ Rni ×n by 1, l = ik , (Ei )k,l = 0, l = ik . Let Pi = EiT Ei and Vi = Pi Rn ,
FSi = Pi F.
It is easy to see that Pi is the orthogonal projection from Rn onto Vi . For each u ∈ Rn , we define Ti (u) ∈ Vi such that (2.1)
FSi (u − Ti (u)) = 0,
i = 1, 2, . . . , N,
and let F(u) =
N
Ti (u),
i=1
which is referred to as the additive Schwarz preconditioned nonlinear function. The ASPIN method tries to find the solution u∗ of (1.1) by solving the nonlinear system (2.2)
F(u) = 0
with the inexact Newton method. About the solvability of (2.1), we have the following proposition. Proposition 2.1. If Ei J(u∗ )EiT is invertible for each i, then there exists a neighborhood U of u∗ and a unique continuously differentiable function Ti : Rn → Vi for each i such that (2.1) holds for each u ∈ U , and also Ti (u∗ ) = 0. Moreover, (2.3)
−1 Ti (u) = EiT Ei J(u − Ti (u))EiT Ei J(u − Ti (u)).
Proof. Theorem 1.1 in [10] shows that there exist a neighborhood U1 of u∗ and a unique continuous function Ti : Rn → Vi for each i such that (2.1) holds for each u ∈ U1 , and Ti (u∗ ) = 0. Also, we know from [10] that Ti satisfies (2.4)
Ti (v) − Ti (u) = DTi (v, u)(v − u),
v, u ∈ U1 ,
where −1 DTi (v, u) = EiT Ei DF (v − Ti (v), u − Ti (u))EiT Ei DF (v − Ti (v), u − Ti (u)),
ON CONVERGENCE OF THE ASPIN METHOD
while
1853
1
DF (v − Ti (v), u − Ti (u)) =
J([v − Ti (v)] + t[(u − Ti (u)) − (v − Ti (v))]) dt. 0
Since J(u) is continuous and Ei J(u∗ )EiT is invertible for each i, Lemma 2.3.3 in [20] shows that there exists a neighborhood U2 of u∗ such that Ei J(u)EiT is invertible for each u ∈ U2 . Now let U ⊂ U1 ∩U2 be a neighborhood of u∗ such that u−Ti (u) ∈ U2 for each u ∈ U and for each i. For u ∈ U , let A(u) = EiT [Ei J(u − Ti (u))EiT ]−1 Ei J(u − Ti (u)); then, according to (2.4) and Lemma 2.3.3 in [20], it is easy to verify that Ti (u + h) − Ti (u) − A(u)h = 0. h h→0 lim
Thus, Ti (u) is Fr´echet-differentiable, and Ti (u) = A(u) = EiT [Ei J(u − Ti (u))EiT ]−1 Ei J(u − Ti (u)). In addition, it is easy to see from Lemma 2.3.3 in [20] that Ti (u) is continuous. It should be pointed out that formula (2.3) has been given in [6], but it has been obtained in a different way. In addition, it should be noted that the condition of Proposition 2.1 is satisfied for any partition of S if J(u∗ ) is positive definite. Since the inexact Newton method concerns the Jacobian of the system, an analysis of the basic property of the Jacobian F (u) is necessary. Because Ti (u) is continuous and Ti (u∗ ) = 0, we know that when u is sufficiently close to u∗ , Ti (u) will be sufficiently close to 0, and as a result, u−Ti (u) will be close to u. Since J(u) is continuous, we may replace J(u − Ti (u)) by J(u); therefore, −1 Ei J(u − Ti (u)) Ti (u) = EiT Ei J(u − Ti (u))EiT −1 ≈ EiT Ei J(u)EiT Ei J(u) ≡ Ri (u). Let J (u) = F (u); then J (u) =
N
Ti (u)
i=1
≈
N
Ri (u) ≡ B(u).
i=1
In implementation of the ASPIN method, the Jacobian J (u) is replaced by B(u), since the latter is easier to use. Remark 2.1. From the proof of Proposition 2.1, we know that if Ei J(u∗ )EiT is invertible for each i, then Ei J(u − Ti (u))EiT and Ei J(u)EiT are invertible in the neighborhood U of u∗ . In addition, (Ei J(u − Ti (u))EiT )−1 and (Ei J(u)EiT )−1 are continuous in U , so J (u) and B(u) are all continuous in U . Furthermore, it is easy to see that (2.5)
lim∗ J (u) = lim∗ B(u) =
u→u
u→u
N
−1 EiT Ei J(u∗ )EiT Ei J(u∗ ) = J (u∗ ).
i=1
From [6] and [10], we can obtain the following result.
1854
HENG-BIN AN
Proposition 2.2. If Ei J(u∗ )EiT is invertible for each i, then there exists a neighborhood D ⊂ U of u∗ with U determined in Proposition 2.1 such that (i) J (u) and B(u) are nonsingular in D; (ii) the nonlinear systems (1.1) and (2.2) are equivalent in the sense that they have the same solution in D. Remark 2.2. Remark 2.1 points out that J (u) and B(u) are all continuous in D ⊂ U ; thus Lemma 2.3.3 in [20] and Proposition 2.2 show that J (u)−1 and B(u)−1 are continuous in D. For convenience of discussion, we describe the ASPIN algorithm here. Let f (u) =
1 F(u)T F(u), 2
which will be used in linesearch in the inexact Newton method. Assume that u(0) is a given initial guess and u(k) is the current approximate solution; the next approximate solution u(k+1) for system (2.2) can be computed through the following steps. Algorithm 2.1 (ASPIN [6]). 0. Let ηmax ∈ (0, 1), α ∈ (0, 1), 0 < θmin < θmax < 1 be given. 1. Compute the nonlinear residual g (k) = F(u(k) ) through the following steps. (k) 1.1. Find gi = Ti (u(k) ) by solving the local subdomain nonlinear systems (k) = 0, FSi u(k) − gi
i = 1, 2, . . . , N,
(k)
with the initial point gi = 0. 1.2. Form the global residual g (k) =
N
(k)
gi .
i=1
1.3. Check the stopping conditions on g (k) . 2. Find the approximate inexact Newton direction p(k) by solving the system B(u(k) )p = −F u(k) such that
(k)
F u + B u(k) p(k) ≤ ηk F u(k) , where ηk ∈ [0, ηmax ] is the forcing term. 3. Perform linesearch along p(k) : 3.1. Let λk = 1. 3.2. While f (u(k) + λk p(k) ) > f (u(k) ) + αλk F(u(k) )T B(u(k) )p(k) , do • choose θ ∈ [θmin , θmax ], • let λk = θλk . 3.3. Let u(k+1) = u(k) + λk p(k) . In step 1.1 of Algorithm 2.1, N subdomain nonlinear systems have to be solved in order to evaluate the preconditioned function F at a given point. Step 3 of Algorithm 2.1 is the linesearch procedure to find a satisfied step. For more details about the ASPIN method, see [6]. We point out that Algorithm 2.1 can be implemented in parallel. For details about implementation, see [6, 7].
ON CONVERGENCE OF THE ASPIN METHOD
1855
3. Local convergence of the ASPIN method. We will prove in this section that the ASPIN method is locally convergent. Note that by Propositions 2.1 and 2.2, the ASPIN method is based on the local property of F (u) at the solution u∗ of system (1.1), so it seems impossible to obtain a global convergence result for this method. In this section and the following, · always denotes the Euclidean norm for both vectors and matrices, and N (u, ρ) = {v | v − u < ρ} represents the open ball with center u and radius ρ. Since the analysis of secondary iteration would complicate the discussion without gaining more insight into the method, we assume that (A1 ) the value of F at each iterative point is evaluated exactly, i.e., (2.1) holds with u replaced by u(k) + λk p(k) for each k; moreover, from now on, we assume that (A2 ) Ei J(u∗ )EiT is invertible for each i; (A3 ) D represents the neighborhood determined in Proposition 2.2; and (A4 ) δ > 0 is a fixed small number such that N (u∗ , δ) ⊂ D; in addition, the following inequalities hold for any u ∈ N (u∗ , δ): (I1 ) B(u) ≤ 2M ; (I2 ) B(u)−1 ≤ 2M ; (I3 ) J (u) ≤ 2M ; (I4 ) J (u)−1 ≤ 2M ; max (I5 ) J (u) − B(u) ≤ 4M1−η (1+ηmax ) ; 1 ∗ ∗ u − u∗ , (I6 ) F(u) − F(u ) − J (u )(u − u∗ ) ≤ 2M where M := max{J (u∗ ), J (u∗ )−1 }. It is easy to see from Remarks 2.1 and 2.2 that inequalities (I1 )–(I5 ) may hold with δ small enough. The last inequality may hold by Lemma 3.2.10 in [20]. In addition, we assume that the parameter α in Algorithm 2.1 is small enough so that 64αM 4 ≤
(3.1)
1 − ηmax . 3 + ηmax
It should be pointed out that the above assumptions are not so strict; see the appendix, where an example is given. Now we show that the inexact Newton direction computed in Algorithm 2.1 is also a regular inexact Newton direction for F in the sense that (1.3) holds. Proposition 3.1. Assume that u ∈ N (u∗ , δ). If (3.2)
F(u) + B(u)p ≤ ηF(u),
η ∈ [0, ηmax ],
then (3.3)
F(u) + J (u)p ≤
1+η F(u). 2
Proof. By (I2 ) and (3.2), p = B(u)−1 [B(u)p + F(u) − F(u)] ≤ (1 + η)B(u)−1 F(u) ≤ 2M (1 + η)F(u).
1856
HENG-BIN AN
Thus, according to (I5 ) and (3.2), J (u)p + F(u) = [J (u) − B(u)]p + B(u)p + F(u) ≤ J (u) − B(u)p + B(u)p + F(u) 1 − ηmax · 2M (1 + η)F(u) + ηF(u) ≤ 4M (1 + ηmax ) 1−η ≤ F(u) + ηF(u) 2 1+η ≤ F(u). 2 Thus we obtain the required inequality. Remark 3.1. If u ∈ N (u∗ , δ) and (3.2) holds, then we have F(u)T B(u)p = F(u)T [B(u)p + F(u) − F(u)] ≤ −F(u)2 + F(u)B(u)p + F(u) ≤ −(1 − η)F(u)2 < 0. In the same way, (3.3) shows that F(u)T J (u)p < 0.
(3.4)
In particular, (3.4) shows that if u ∈ N (u∗ , δ) and p is computed by the ASPIN method, then p is a descent direction for the function f (u) = 12 F(u)2 at point u.1 The following lemma is needed in our analysis. Lemma 3.2 (see [4, Lemma 3.4]). Let u ∈ Rn and H : Rn → Rn be continuously differentiable in a neighborhood of u. Assume that H(u) = 0 and H (u) is nonsingular. If p ∈ Rn such that H(u) + H (u)p ≤ ηH(u),
η ∈ [0, 1),
then |∇h(u)T p| 1−η ≥ ∇h(u) > 0, p (1 + η)κ(H (u)) where κ(H (u)) is the condition number for H (u) and h(u) = 21 H(u)T H(u). Based on Lemma 3.2, we have the following result. Lemma 3.3. Assume that u ∈ N (u∗ , δ). If F(u) + B(u)p ≤ ηF(u),
η ∈ [0, ηmax ],
then it holds that |F(u)T J (u)p| ≥ 4α|F(u)T B(u)p|. Proof. According to Proposition 3.1, we have F(u) + J (u)p ≤ 1p
1+η F(u); 2
is a descent direction for f (u) if ∇f (u)T p < 0. In addition, note that ∇f (u) = J (u)T F (u).
ON CONVERGENCE OF THE ASPIN METHOD
1857
therefore, by Lemma 3.2, (I3 ), and (I4 ), 1−η F(u)T J (u)p (3 + η)κ(J (u)) 1−η F(u)p ≥ (3 + η)κ(J (u))J (u)−1 1−η ≥ F(u)p. 8M 3 (3 + η)
|F(u)T J (u)p| ≥
(3.5)
Thus, by (I1 ), (3.5), and (3.1), 4α|F(u)T B(u)p| ≤ 4αB(u)F(u)p ≤ 8αM F(u)p 8M 3 (3 + η) ≤ 8αM |F(u)T J (u)p| 1−η 3 + ηmax ≤ 64αM 4 |F(u)T J (u)p| 1 − ηmax ≤ |F(u)T J (u)p|. This concludes the proof. The following lemma shows that if u is sufficiently close to u∗ , then the direction p obtained in the ASPIN method will not be too long. δ Lemma 3.4. Assume that u ∈ N (u∗ , 2δ ) with F(u) ≤ 8M . If p ∈ Rn such that F(u) + B(u)p ≤ ηF(u),
η ∈ [0, ηmax ],
then [u, u + p] ⊂ N (u∗ , δ), where [u, u + p] represents the line segment between u and u + p. δ Proof. Since F(u) ≤ 8M , by (I2 ), we have p = B(u)−1 [B(u)p + F(u) − F(u)] ≤ B(u)−1 · [B(u)p + F(u) + F(u)] ≤ 2M (1 + η)F(u) ≤ 2M (1 + ηmax )F(u) δ(1 + ηmax ) ≤ 4 δ < . 2 Therefore, by u ∈ N (u∗ , 2δ ), (u + p) − u∗ ≤ u − u∗ + p
0 such that (3.6)
∇f (v) − ∇f (w) ≤ γv − w
∀ v, w ∈ N (u∗ , δ).
If p ∈ Rn such that F(u) + B(u)p ≤ ηF(u),
η ∈ [0, ηmax ],
then the linesearch procedure along p in Algorithm 2.1 will terminate in finite iterations and the obtained λ satisfies αθmin |F(u)T J (u)p| λ ≥ min 1, . γp2 Proof. Because u ∈ N (u∗ , 2δ ), Lemma 3.3 shows that |F(u)T J (u)p| ≥ 4α|F(u)T B(u)p|.
(3.7)
Since F(u)T J (u)p < 0, F(u)T B(u)p < 0, (3.7) shows that F(u)T J (u)p ≤ 4αF(u)T B(u)p, so (3.8)
F(u)T J (u)p − αF(u)T B(u)p ≤ 3αF(u)T B(u)p < 0.
δ Because u ∈ N (u∗ , 2δ ) and F(u) ≤ 8M , Lemma 3.4 shows that [u, u + p] ⊂ ∗ N (u , δ). Thus, by the mean value theorem, there exists ξ ∈ [u, u + p] such that
f (u + λp) = f (u) + λ∇f (ξ)T p. Therefore, f (u + λp) = f (u) + λ∇f (ξ)T p = f (u) + αλF(u)T B(u)p − αλF(u)T B(u)p + λ∇f (ξ)T p = f (u) + αλF(u)T B(u)p + λ{[∇f (ξ)T p − ∇f (u)T p] +[∇f (u)T p − αF(u)T B(u)p]} = f (u) + αλF(u)T B(u)p + λ{λζ + [F(u)T J (u)p − αF(u)T B(u)p]}, where ζ=
∇f (ξ)T p − ∇f (u)T p . λ
By (3.6), we have |ζ| ≤ γp2 , so
f (u + λp) ≤ f (u) + αλF(u)T B(u)p + λ λγp2 + [F(u)T J (u)p − αF(u)T B(u)p] . Thus, if λγp2 + [F(u)T J (u)p − αF(u)T B(u)p] ≤ 0,
ON CONVERGENCE OF THE ASPIN METHOD
1859
then λ is acceptable. Since λ is reduced by a factor θ ≤ θmax < 1 at each iteration of the while-loop, it follows from (3.8) that the while-loop will terminate in finite steps. Let λ be the ultimate step factor. If λ = 1, then the needed conclusion trivially holds. Now suppose that the linesearch procedure is implemented at least once, and let λ− be the penultimate value; then the above argument shows that λ− γp2 + [F(u)T J (u)p − αF(u)T B(u)p] > 0. Consequently, λ− >
|F(u)T J (u)p| − α|F(u)T B(u)p| . γp2
Therefore, it follows from (3.7) that λ ≥ θmin λ− |F(u)T J (u)p| − α|F(u)T B(u)p| > θmin γp2 T 4α|F(u) B(u)p| − α|F(u)T B(u)p| ≥ θmin γp2 T αθmin |F(u) B(u)p| ≥ . γp2 Thus, we have obtained the required conclusion. Lemma 3.6. If u ∈ N (u∗ , 2δ ) and F(u)
F u(k+1) . Thus, by induction, the ASPIN method can generate a sequence {u(k) } ⊂ N (u∗ , 2δ ). Next we show that f u(k) → 0. By (3.10),
(3.11)
T T F u(k) B u(k) p(k) = F u(k) [B u(k) p(k) + F u(k) − F u(k) ] ≤ −2(1 − ηmax )f u(k) ,
and it follows that T f (u(k) + λk p(k) ) ≤ f u(k) + αλk F u(k) B u(k) p(k) ≤ [1 − 2αλk (1 − ηmax )]f u(k) . Therefore, by Theorem 3.5, we have f (u
(k)
+ λk p
(3.12)
(k)
T αθmin |F u(k) B u(k) p(k) | (k) f u ) ≤ 1 − 2α(1 − ηmax ) γp(k) 2 ≡ (1 − ctk )f u(k) ,
where we set c=
2α2 θmin (1 − ηmax ) γ
and T |F u(k) B u(k) p(k) | . tk = p(k) 2 Because {f u(k) } is nonnegative and strictly decreased, limk→∞ f u(k) exists. (k) If limk→∞ f u > 0, then (3.12) shows that f u(k+1) ≤ lim (1 − ctk ) ≤ 1, lim k→∞ f u(k) k→∞ so lim (1 − ctk ) = 1
k→∞
or, equivalently, lim tk = 0.
k→∞
1862
HENG-BIN AN
Since (3.11) shows that T |F u(k) B u(k) p(k) | 2(1 − ηmax )f u(k) ≤ = tk , p(k) 2 p(k) 2 we have
f u(k) lim = 0. k→∞ p(k) 2
Because limk→∞ f u(k) > 0, it follows that (3.13)
p(k) → ∞ (k → ∞).
But on the other hand, by (I2 ) and (3.10), we have −1 (k) (k) p(k) = B u(k) [B u p + F u(k) − F u(k) ] −1 ≤ B u(k) · [B u(k) p(k) + F u(k) + F u(k) ] ≤ 2M (1 + ηk )F u(k) ≤ 2M (1 + ηmax )F(u(0) ) δ(1 + ηmax ) , ≤ 4 which contradicts (3.13). Thus, we must have f (u(k) ) → 0. Since {u(k) } ⊂ N (u∗ , δ), by (I6 ), u(k) − u∗ = 2u(k) − u∗ − u(k) − u∗
≤ 2J (u∗ )−1 J (u∗ )(u(k) −u∗ )−2M F u(k) −F(u∗ )−J (u∗ )(u(k) −u∗ ) ≤ 2M [J (u∗ )(u(k) − u∗ ) − F u(k) − F(u∗ ) − J (u∗ )(u(k) − u∗ )] ≤ 2M F u(k) . Because F(u(k) ) = 2f (u(k) ) → 0, we have u(k) → u∗ . Now, we complete the discussion for local convergence of the ASPIN method. 4. Convergence rate of the ASPIN method. In this section, we discuss the convergence rate of the ASPIN method. We will show that the ASPIN method is quadratically convergent under suitable conditions. The following theorem shows that under suitable assumptions, λk = 1 is acceptable for all k sufficiently large. Theorem 4.1. Assume that f (u) is twice continuously differentiable in N (u∗ , δ), and there exists γ > 0 such that for any v, w ∈ N (u∗ , δ), ∇f (v) − ∇f (w) ≤ γv − w, ∇2 f (v) − ∇2 f (w) ≤ γv − w. Let {u(k) } be the sequence generated by the ASPIN method such that u(k) → u∗ . If ηk → 0, then u(k+1) = u(k) + p(k) for all sufficiently large k. Proof. Because u(k) → u∗ , without loss of generality, we may assume that (k) {u } ⊂ N (u∗ , 2δ ). Thus, it follows from Proposition 3.1 that (4.1)
1 + ηmax (k) F u(k) + J u(k) p(k) ≤ F u . 2
ON CONVERGENCE OF THE ASPIN METHOD
1863
Therefore, by (I4 ) and (4.1),
(4.2)
−1 (k) (k) p + F u(k) − F u(k) ] [J u p(k) = J u(k) −1 ≤ J u(k) · [J u(k) p(k) + F u(k) + F u(k) ] 1 + ηmax + 1 F u(k) ≤ 2M 2 = M (3 + ηmax )F u(k) .
Because F(u(k) ) → 0, the above inequality shows that p(k) → 0 (k → ∞). By (I4 ), we have T −1 −1 (k) 1 ∇f u(k) = J u(k) F u(k) ≥ J u(k) F u(k) ; ≥ F u 2M thus, Lemma 3.2 in connection with (I3 ), (I4 ), and (4.1) shows that T |∇f u(k) p(k) | 1 − 1+η2max 1 · ∇f u(k) ≥ 1+ηmax · (k) (k) p ) κ(J u 1+ 2 1 − ηmax 1 1 ≥ · · F u(k) 2 3 + ηmax 4M 2M 1 − ηmax = F u(k) , 8M 3 (3 + ηmax ) or we have (4.3)
T 8M 3 (3 + ηmax ) F u(k) p(k) ≤ · |∇f u(k) p(k) | 1 − ηmax T ≡ a|∇f u(k) p(k) |,
where a=
8M 3 (3 + ηmax ) . 1 − ηmax
From (4.2) and (4.3), we obtain p(k) 2 ≤ M (3 + ηmax )F u(k) p(k) T 8M 4 (3 + ηmax )2 ≤ · |∇f u(k) p(k) | 1 − ηmax T ≡ b|∇f u(k) p(k) |, where b=
8M 4 (3 + ηmax )2 . 1 − ηmax
1864
HENG-BIN AN
Next we show that λk = 1 is acceptable for all k sufficiently large. First note that if F = (F1 , F2 , . . . , Fn )T , then ∇2 f (u) = J (u)T J (u) +
n
Fi (u)∇2 Fi (u)
i=1
≡ J (u) J (u) + S(u). T
Since u(k) → u∗ and F(u∗ ) = 0, it follows that S(u(k) ) → 0. For each k, by the mean value theorem, there exists a ξ (k) on the line segment between u(k) and u(k) + p(k) such that 1 T T 1 f u(k) + p(k) − f u(k) − ∇f u(k) p(k) = ∇f u(k) + ∇2 f ξ (k) p(k) p(k) . 2 2 This gives (k) f u + p(k) − f u(k) − 1 ∇f u(k) T p(k) 2 1 T = ∇f u(k) + ∇2 f (ξ (k) )p(k) p(k) 2 T 1 = ∇f u(k) + ∇2 f u(k) p(k) p(k) + (p(k) )T ∇2 f (ξ (k) ) − ∇2 f u(k) p(k) 2
2 T (k)
1 ≤ J u(k) F u + J u(k) p(k) · p(k) + S u(k) + γ p(k) p(k)
2
2
1 ≤ ηk J u(k) · F u(k) · p(k) + S u(k) + γ p(k) p(k)
2
T 1 ≤ − 2aM ηk + b S u(k) + γ p(k) ∇f u(k) p(k) 2 T 1 ≡ − k ∇f u(k) p(k) ; 2 therefore, (4.4)
1 T f (u(k) + p(k) ) − f u(k) ≤ (1 − k )∇f u(k) p(k) . 2
Since ηk , S(u(k) ), and p(k) all converge to zero, it follows that k → 0. Thus, for all k sufficiently large, we have k
0 such that ∇f (v) − ∇f (w) ≤ γv − w, ∇2 f (v) − ∇2 f (w) ≤ γv − w for any v, w ∈ N (u∗ , δ). If {u(k) } is a sequence generated by the ASPIN method such that u(k) → u∗ and u(k+1) = u(k) + p(k) for all sufficiently large k, then (i) u(k) → u∗ superlinearly if and only if rk = o(F u(k) ), where rk = F u(k) + B u(k) p(k) ; (ii) u(k) → u∗ quadratically if and only if rk = O(F u(k) 2 ).
ON CONVERGENCE OF THE ASPIN METHOD
1867
Proof. Since u(k+1) = u(k) + p(k) for all sufficiently large k, without loss of generality, we assume that u(k+1) = u(k) + p(k) for all k in the following argument. Assume that u(k) → u∗ superlinearly. Since rk = F u(k) + B u(k) p(k) = F u(k) − F(u∗ ) − J (u∗ ) u(k) − u∗ − B u(k) − J (u∗ ) u(k) − u∗ + B u(k) u(k+1) − u∗ , (4.10) by Lemma 3.2.10 in [20], continuity of B(u) at u∗ and the assumption that u(k) → u∗ superlinearly, we have rk ≤ o(u(k) − u∗ ) + o(1)u(k) − u∗ + o(u(k) − u∗ ). Therefore, by Lemma 4.3, it follows that rk = o(u(k) − u∗ ) = o(F u(k) ). Conversely, assume that rk = o(F u(k) ). Since
(4.11)
u(k+1) − u∗ = (u(k) − u∗ ) + p(k) −1 = (u(k) − u∗ ) + B u(k) [rk − F u(k) ] −1 (k) B u − J (u∗ ) (u(k) − u∗ ) = B u(k) + rk − F u(k) − F(u∗ ) − J (u∗ )(u(k) − u∗ ) ,
thus, by (I2 ), the continuity of B(u) at u∗ , the assumption that rk = o(F u(k) ), and Lemma 3.2.10 in [20], we have u(k+1) − u∗ ≤ 2M [o(1)u(k) − u∗ + o(F u(k) ) + o(u(k) − u∗ )]. Therefore, by Lemma 4.3, u(k+1) − u∗ = o(u(k) − u∗ ) + o(F u(k) ) = o(u(k) − u∗ ). Since F is twice continuously differentiable in N (u∗ , δ), Lemma 4.2 shows that both J (u) and B(u) are Lipschitz continuous in a neighborhood V ⊂ N (u∗ , 2δ ) of u∗ . Thus, there exists L > 0 such that for any u ∈ V , J (u) − J (u∗ ) ≤ Lu − u∗
(4.12) and
B(u) − B(u∗ ) ≤ Lu − u∗ .
(4.13)
By (4.12) and Lemma 3.2.12 in [20], (4.14)
F(u) − F(u∗ ) − J (u∗ )(u − u∗ ) ≤
L u − u∗ 2 2
∀ u ∈ V.
In a similar way, it follows from (4.10), (4.11), (4.13), (4.14), and Lemma 4.3 that u(k) → u∗ quadratically if and only if rk = O(F u(k) 2 ). This concludes the proof.
1868
HENG-BIN AN
From Theorems 4.1 and 4.4, we can obtain the following result. Corollary 4.5. Assume that both F (u) and f (u) are twice continuously differentiable in N (u∗ , δ), and there exists γ > 0 such that for any v, w ∈ N (u∗ , δ), ∇f (v) − ∇f (w) ≤ γv − w, ∇2 f (v) − ∇2 f (w) ≤ γv − w. Let {u(k) } be a sequence generated by the ASPIN method such that u(k) → u∗ . Then (i) u(k) → u∗ superlinearly if ηk → 0; (ii) u(k) → u∗ quadratically if ηk = O(F u(k) ). Corollary 4.5 reflects how the forcing term influences the convergence rate of the ASPIN method. This result is similar to Corollary 3.5 in [8] for the inexact Newton method. In particular, by Corollary 4.5, we can determine the convergence rate of the ASPIN method by choosing proper forcing terms. 5. Conclusion. The inexact Newton method is one of the effective tools for solving large sparse systems of nonlinear equations. By using nonlinear additive Schwarz preconditioning technique, Cai and Keyes [6] introduced the ASPIN method. This method is very effective for solving some nonlinear problems with strong nonbalanced nonlinearities. However, the convergence of the ASPIN method is not discussed by them or others. In this paper, we discussed the convergence property of the ASPIN method and thus we provided a theoretical support for the ASPIN method. The convergence result is local since the design of the ASPIN only concerns the local properties of the original function. Furthermore, we discussed the convergence rate for the ASPIN method, and the result shows that the convergence rate of the ASPIN method is similar to that of the inexact Newton method. Thus, we can obtain the desired convergence rate by choosing proper forcing terms. Appendix. We give a simple example to show that our main assumptions, inequalities (I1 )–(I6 ) in section 3, are not so strict. Consider the nonlinear equations 2u1 − u2 + λeu1 − λ = 0, −u1 + 2u2 + λeu2 − λ = 0, where λ > 0. It is obvious that u∗ = (0, 0)T is a solution of the system. Let S = {1, 2},
S1 = {1},
S2 = {2}
and T1 (u) =
T11 (u) , 0
T2 (u) =
0 , T22 (u)
where T11 , T22 : R2 → R. By FSi (u − Ti (u)) = 0, we have 2(ui − Tii (u)) − u3−i + λeui −Tii (u) − λ = 0,
(A.1)
i = 1, 2.
It is easy to see that Tii (u) are continuous functions. Furthermore, we can obtain T1 (u)
=
−1 1 − 2 + λeu1 −T11 (u) , 0 0
T2 (u)
=
0
− 2 + λeu2 −T22
(u) −1
0 . 1
ON CONVERGENCE OF THE ASPIN METHOD
Thus,
1
J (u) =
−1 − 2 + λeu2 −T22 (u)
Besides, it is easy to see that B(u) =
1 −1 − (2 + λeu2 )
Since B(u∗ ) = J (u∗ ) =
1869
−1 − 2 + λeu1 −T11 (u) . 1 −1
− (2 + λeu1 ) 1
1 −(2 + λ)−1
.
−(2 + λ)−1 , 1
so M = max{J (u∗ ), J (u∗ )−1 } =
2+λ . 1+λ
Now we will choose some proper δ such that inequalities (I1 )–(I6 ) hold. Let δ ∈ (0, 12 ) and assume that u = (u1 , u2 )T ∈ N (u∗ , δ). Note that if 1 a A= , b 1 then A ≤ 1 + max{|a|, |b|}. Therefore, B(u) ≤ 1 + max (2 + λeui )−1 < 1 + (2 + λe−δ )−1 < 2M. 1≤i≤2
In addition, since B(u)−1 =
1 u 1 1 − (2 + λe )−1 (2 + λeu2 )−1
1 (2 + λeu2 )−1
(2 + λeu1 )−1 , 1
one can easily verify that B(u)−1 ≤
2 + λe−δ < 2M. 1 + λe−δ
Because u ∈ N (u∗ , δ), by (A.1), we have |ui − Tii (u)|
1 + x,
x ∈ (0, 1),
1 and ex < 1 + x, 2
x ∈ (−1, 0),
we have |ui − Tii (u)| ≤
2 |u3−i |, 4+λ
i = 1, 2.
Therefore,
T11 (u) − u1 + (2 + λ)−1 u2
F(u) − F(u∗ ) − J (u∗ )(u − u∗ ) =
T22 (u) − u2 + (2 + λ)−1 u1
|T11 (u) − u1 | + (2 + λ)−1 |u2 |
≤
|T22 (u) − u2 | + (2 + λ)−1 |u1 |
3 ≤ u < 2M u − u∗ . 2+λ Thus, inequality (I6 ) holds. Summing up the above discussion, we know that if 1 (1 + λ)2 (1 − ηmax ) , , δ ≤ min 2 3eλ(2 + λ)(1 + ηmax ) then for any u ∈ N (u∗ , δ), inequalities (I1 )–(I6 ) hold. For this example, other assumptions in the paper can also be checked to be true. Acknowledgments. The author thanks the anonymous referees and Carol S. Woodward for many valuable suggestions. REFERENCES [1] S. Bellavia and B. Morini, A globally convergent Newton-GMRES subspace method for systems of nonlinear equations, SIAM J. Sci. Comput., 23 (2001), pp. 940–960. [2] S. Bellavia, M. Macconi, and B. Morini, A hybrid Newton-GMRES method for solving nonlinear equations, in Numerical Analysis and Its Applications, Lecture Notes in Comput. Sci. 1988, L. Vulkov, J. Wasniewski, and P. Yalamov, eds., Springer-Verlag, Berlin, 2000, pp. 68–75. [3] P. N. Brown and Y. Saad, Hybrid Krylov methods for nonlinear systems of equations, SIAM J. Sci. Statist. Comput., 11 (1990), pp. 450–481. [4] P. N. Brown and Y. Saad, Convergence theory of nonlinear Newton-Krylov algorithms, SIAM J. Optim., 4 (1994), pp. 297–330. [5] X.-C. Cai, M. Dryja, and M. Sarkis, Restricted additive Schwarz preconditioners with harmonic overlap for symmetric positive definite linear systems, SIAM J. Numer. Anal., 41 (2003), pp. 1209–1231. [6] X.-C. Cai and D. E. Keyes, Nonlinearly preconditioned inexact Newton algorithms, SIAM J. Sci. Comput., 24 (2002), pp. 183–200. [7] X.-C. Cai, D. E. Keyes, and L. Marcinkowski, Nonlinear additive Schwarz preconditioners and applications in computational fluid dynamics, Internat. J. Numer. Methods Fluid Mech., 40 (2002), pp. 1463–1470.
ON CONVERGENCE OF THE ASPIN METHOD
1871
[8] R. S. Dembo, S. C. Eisenstat, and T. Steihaug, Inexact Newton methods, SIAM J. Numer. Anal., 19 (1982), pp. 400–408. [9] J. E. Dennis and R. B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice–Hall, Englewood Cliffs, NJ, 1983. [10] M. Dryja and W. Hackbusch, On the nonlinear domain decomposition method, BIT, 37 (1997), pp. 296–311. [11] M. Dryja and O. B. Widlund, Domain decomposition algorithms with small overlap, SIAM J. Sci. Comput., 15 (1994), pp. 604–620. [12] S. C. Eisenstat and H. F. Walker, Globally convergent inexact Newton methods, SIAM J. Optim., 4 (1994), pp. 393–422. [13] S. C. Eisenstat and H. F. Walker, Choosing the forcing term in an inexact Newton method, SIAM J. Sci. Comput., 17 (1996), pp. 16–32. [14] D. R. Fokkema, G. L. G. Sleijpen, and H. A. van der Vorst, Accelerated inexact Newton schemes for large systems of nonlinear equations, SIAM J. Sci. Comput., 19 (1998), pp. 657–674. [15] I. E. Kaporin and O. Axelsson, On a class of nonlinear equation solvers based on the residual norm reduction over a sequence of affine subspaces, SIAM J. Sci. Comput., 16 (1995), pp. 228–249. [16] D. A. Knoll and D. E. Keyes, Jacobian-free Newton-Krylov methods: A survey of approaches and applications, J. Comput. Phys., 193 (2004), pp. 357–397. [17] S. H. Lui, Nonlinearly preconditioned Newton’s method, in Domain Decomposition Methods in Science and Engineering, Fourteenth International Conference on Domain Decomposition Methods, Cocoyoc, Mexico, I. Herrera, D. E. Keyes, O. B. Widlund, and R. Yates, eds., National Autonomous University of Mexico (UNAM), Mexico City, Mexico, 2003, pp. 95– 105. [18] L. Luksan, Inexact trust region method for large sparse systems of nonlinear equations, J. Optim. Theory Appl., 81 (1994), pp. 569–591. [19] M. Pernice and H. F. Walker, NITSOL: A Newton iterative solver for nonlinear systems, SIAM J. Sci. Comput., 19 (1998), pp. 302–318. [20] J. M. Ortega and W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Classics Appl. Math. 30, SIAM, Philadelphia, PA, 2000. [21] Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Comput., 7 (1986), pp. 856–869.
c 2005 Society for Industrial and Applied Mathematics
ON CONVERGENCE OF THE ADDITIVE SCHWARZ PRECONDITIONED INEXACT NEWTON METHOD∗ HENG-BIN AN† Abstract. The additive Schwarz preconditioned inexact Newton (ASPIN) method was recently introduced [X.-C. Cai and D. E. Keyes, SIAM J. Sci. Comput., 24 (2002), pp. 183–200] to solve the systems of nonlinear equations with nonbalanced nonlinearities. Although the ASPIN method has successfully been used to solve some difficult nonlinear equations, its convergence property has not been studied since it was proposed. In this paper, the convergence property of the ASPIN method is studied, and the obtained result shows that this method is locally convergent. Furthermore, the convergence rate for the ASPIN method is discussed and the obtained result is similar to that of the inexact Newton method. Key words. nonlinear systems, inexact Newton methods, nonlinear preconditioning, additive Schwarz, local convergence, convergence rate AMS subject classifications. 65H10, 65N12, 65N22 DOI. 10.1137/040611653
1. Introduction. Consider the nonlinear system of equations (1.1)
F (u) = 0,
where F : Rn → Rn is a continuously differentiable function. For convenience of discussion, let F = (F1 , F2 , . . . , Fn )T , u = (u1 , u2 , . . . , un )T , and J(u) = F (u). A numerical solution for (1.1) is often required in many scientific and engineering computing areas such as the discretization of nonlinear partial differential equations; see [7, 16]. The inexact Newton method [8] is one of the most important and effective tools for solving such systems, in particular, when the problem is large and sparse. In applications, some global strategies, such as linesearch or trust region techniques, are often needed because the inexact Newton method is locally convergent [1, 2, 3, 4, 12]. In particular, if the linesearch backtracking technique is augmented in the inexact Newton method, then the inexact Newton with backtracking (INB) method is obtained [12, 13, 19]. This method is more robust and it can be briefly described here. Suppose u(0) is a given initial guess and let u(k) be the current approximate solution; the next approximate solution u(k+1) can be obtained through the following steps. Algorithm 1.1 (INB [12]). 1. Inexactly solve the system (1.2) J u(k) p = −F u(k) , and obtain an inexact Newton direction p(k) such that (1.3) F u(k) + J u(k) p(k) ≤ ηk F u(k) . 2. Compute the new approximate solution u(k+1) = u(k) + λk p(k) . ∗ Received by the editors June 24, 2004; accepted for publication (in revised form) February 11, 2005; published electronically November 22, 2005. http://www.siam.org/journals/sinum/43-5/61165.html † State Key Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, P.O. Box 8009, Beijing 100088, People’s Republic of China ([email protected]).
1850
ON CONVERGENCE OF THE ASPIN METHOD
1851
Here ηk ∈ [0, 1) is the forcing term that controls how accurately system (1.2) should be solved, and p(k) is the inexact Newton direction of F at u(k) . Step 2 in Algorithm 1.1 is a linesearch procedure that is used to find a satisfied step factor λk ∈ (0, 1] and then form the next approximate solution. Usually, we use linear iterative methods, such as the classical splitting method or the modern Krylov subspace method, to inexactly solve system (1.2). Thus, the inexact Newton method is an inner-outer iterative method. In particular, when the Krylov subspace method is used in an inner iteration, we obtain the Newton–Krylov subspace method, which has been used successfully in many areas [1, 2, 3, 4, 16]. Although the inexact Newton method works very well for most nonlinear equations, this may often fail when it is used to solve some difficult problems. Many numerical experiments show that most failed cases in the inexact Newton method result from stagnation, particularly when it is used to solve some problems with nonbalanced nonlinearities [1, 6]. Usually, the stagnation phenomenon is caused by the lack of a good initial guess and/or problematic regions such as boundary layers, singularities in the domain, and/or multiphysics domain, etc. See [17]. Considering this, Cai and Keyes [6] recently proposed a nonlinearly preconditioned inexact Newton algorithm: first convert system (1.1) into another nonlinear system F(u) = 0 such that the two systems have the same solution u∗ ∈ Rn ; then use Algorithm 1.1 to solve F(u) = 0. F and F may have completely different forms, but they must have the same solution. Usually, F has more uniform nonlinearities, so it is relatively easy to solve. In [6], an especially preconditioned case, where F is obtained by the single-level nonlinear additive Schwarz method, is discussed in detail. The corresponding method is the additive Schwarz preconditioned inexact Newton (ASPIN) method. Numerical results in [6] show that the ASPIN method can solve some difficult problems where the traditional inexact Newton method fails. Although the ASPIN method has better numerical results than the traditional inexact Newton method, it is unfortunate that until now the convergence property for the ASPIN method has not been given much importance except for some preliminary convergence analysis in the context of semilinear PDEs in paper [17]. In this paper, we show that the ASPIN method is locally convergent; thus we give theoretical support for the ASPIN method. Moreover, we will discuss the convergence rate of the ASPIN method. The rest of the paper is organized as follows. In section 2, we briefly discuss the ASPIN method, and some of its properties are listed. In section 3, we show that the ASPIN method is locally convergent, and its convergence rate is discussed in section 4. Finally, in section 5, some brief conclusions are given. 2. The ASPIN method. Assume that F (u∗ ) = 0 and J(u∗ ) is invertible. To find the solution u∗ of system (1.1), the ASPIN method solves another nonlinear system F(u) = 0, which is obtained from (1.1) through the additive Schwarz preconditioning technique. Specifically, the ASPIN method can be described as follows. Let S = {1, 2, . . . , n} be an index set, i.e., one integer corresponds to each ui and Fi . Assume that S has a partition {S1 , S2 , . . . , SN } such that N i=1
Si = S
and Si ⊂ S.
1852
HENG-BIN AN
Here the subsets may overlap. Let ni = |Si | be the dimension of Si ; then N
ni ≥ n.
i=1
Assume that Si = {i1 , i2 , . . . , ini }, where i1 < i2 < · · · < ini . For i = 1, 2, . . . , N , define matrices Ei ∈ Rni ×n by 1, l = ik , (Ei )k,l = 0, l = ik . Let Pi = EiT Ei and Vi = Pi Rn ,
FSi = Pi F.
It is easy to see that Pi is the orthogonal projection from Rn onto Vi . For each u ∈ Rn , we define Ti (u) ∈ Vi such that (2.1)
FSi (u − Ti (u)) = 0,
i = 1, 2, . . . , N,
and let F(u) =
N
Ti (u),
i=1
which is referred to as the additive Schwarz preconditioned nonlinear function. The ASPIN method tries to find the solution u∗ of (1.1) by solving the nonlinear system (2.2)
F(u) = 0
with the inexact Newton method. About the solvability of (2.1), we have the following proposition. Proposition 2.1. If Ei J(u∗ )EiT is invertible for each i, then there exists a neighborhood U of u∗ and a unique continuously differentiable function Ti : Rn → Vi for each i such that (2.1) holds for each u ∈ U , and also Ti (u∗ ) = 0. Moreover, (2.3)
−1 Ti (u) = EiT Ei J(u − Ti (u))EiT Ei J(u − Ti (u)).
Proof. Theorem 1.1 in [10] shows that there exist a neighborhood U1 of u∗ and a unique continuous function Ti : Rn → Vi for each i such that (2.1) holds for each u ∈ U1 , and Ti (u∗ ) = 0. Also, we know from [10] that Ti satisfies (2.4)
Ti (v) − Ti (u) = DTi (v, u)(v − u),
v, u ∈ U1 ,
where −1 DTi (v, u) = EiT Ei DF (v − Ti (v), u − Ti (u))EiT Ei DF (v − Ti (v), u − Ti (u)),
ON CONVERGENCE OF THE ASPIN METHOD
while
1853
1
DF (v − Ti (v), u − Ti (u)) =
J([v − Ti (v)] + t[(u − Ti (u)) − (v − Ti (v))]) dt. 0
Since J(u) is continuous and Ei J(u∗ )EiT is invertible for each i, Lemma 2.3.3 in [20] shows that there exists a neighborhood U2 of u∗ such that Ei J(u)EiT is invertible for each u ∈ U2 . Now let U ⊂ U1 ∩U2 be a neighborhood of u∗ such that u−Ti (u) ∈ U2 for each u ∈ U and for each i. For u ∈ U , let A(u) = EiT [Ei J(u − Ti (u))EiT ]−1 Ei J(u − Ti (u)); then, according to (2.4) and Lemma 2.3.3 in [20], it is easy to verify that Ti (u + h) − Ti (u) − A(u)h = 0. h h→0 lim
Thus, Ti (u) is Fr´echet-differentiable, and Ti (u) = A(u) = EiT [Ei J(u − Ti (u))EiT ]−1 Ei J(u − Ti (u)). In addition, it is easy to see from Lemma 2.3.3 in [20] that Ti (u) is continuous. It should be pointed out that formula (2.3) has been given in [6], but it has been obtained in a different way. In addition, it should be noted that the condition of Proposition 2.1 is satisfied for any partition of S if J(u∗ ) is positive definite. Since the inexact Newton method concerns the Jacobian of the system, an analysis of the basic property of the Jacobian F (u) is necessary. Because Ti (u) is continuous and Ti (u∗ ) = 0, we know that when u is sufficiently close to u∗ , Ti (u) will be sufficiently close to 0, and as a result, u−Ti (u) will be close to u. Since J(u) is continuous, we may replace J(u − Ti (u)) by J(u); therefore, −1 Ei J(u − Ti (u)) Ti (u) = EiT Ei J(u − Ti (u))EiT −1 ≈ EiT Ei J(u)EiT Ei J(u) ≡ Ri (u). Let J (u) = F (u); then J (u) =
N
Ti (u)
i=1
≈
N
Ri (u) ≡ B(u).
i=1
In implementation of the ASPIN method, the Jacobian J (u) is replaced by B(u), since the latter is easier to use. Remark 2.1. From the proof of Proposition 2.1, we know that if Ei J(u∗ )EiT is invertible for each i, then Ei J(u − Ti (u))EiT and Ei J(u)EiT are invertible in the neighborhood U of u∗ . In addition, (Ei J(u − Ti (u))EiT )−1 and (Ei J(u)EiT )−1 are continuous in U , so J (u) and B(u) are all continuous in U . Furthermore, it is easy to see that (2.5)
lim∗ J (u) = lim∗ B(u) =
u→u
u→u
N
−1 EiT Ei J(u∗ )EiT Ei J(u∗ ) = J (u∗ ).
i=1
From [6] and [10], we can obtain the following result.
1854
HENG-BIN AN
Proposition 2.2. If Ei J(u∗ )EiT is invertible for each i, then there exists a neighborhood D ⊂ U of u∗ with U determined in Proposition 2.1 such that (i) J (u) and B(u) are nonsingular in D; (ii) the nonlinear systems (1.1) and (2.2) are equivalent in the sense that they have the same solution in D. Remark 2.2. Remark 2.1 points out that J (u) and B(u) are all continuous in D ⊂ U ; thus Lemma 2.3.3 in [20] and Proposition 2.2 show that J (u)−1 and B(u)−1 are continuous in D. For convenience of discussion, we describe the ASPIN algorithm here. Let f (u) =
1 F(u)T F(u), 2
which will be used in linesearch in the inexact Newton method. Assume that u(0) is a given initial guess and u(k) is the current approximate solution; the next approximate solution u(k+1) for system (2.2) can be computed through the following steps. Algorithm 2.1 (ASPIN [6]). 0. Let ηmax ∈ (0, 1), α ∈ (0, 1), 0 < θmin < θmax < 1 be given. 1. Compute the nonlinear residual g (k) = F(u(k) ) through the following steps. (k) 1.1. Find gi = Ti (u(k) ) by solving the local subdomain nonlinear systems (k) = 0, FSi u(k) − gi
i = 1, 2, . . . , N,
(k)
with the initial point gi = 0. 1.2. Form the global residual g (k) =
N
(k)
gi .
i=1
1.3. Check the stopping conditions on g (k) . 2. Find the approximate inexact Newton direction p(k) by solving the system B(u(k) )p = −F u(k) such that
(k)
F u + B u(k) p(k) ≤ ηk F u(k) , where ηk ∈ [0, ηmax ] is the forcing term. 3. Perform linesearch along p(k) : 3.1. Let λk = 1. 3.2. While f (u(k) + λk p(k) ) > f (u(k) ) + αλk F(u(k) )T B(u(k) )p(k) , do • choose θ ∈ [θmin , θmax ], • let λk = θλk . 3.3. Let u(k+1) = u(k) + λk p(k) . In step 1.1 of Algorithm 2.1, N subdomain nonlinear systems have to be solved in order to evaluate the preconditioned function F at a given point. Step 3 of Algorithm 2.1 is the linesearch procedure to find a satisfied step. For more details about the ASPIN method, see [6]. We point out that Algorithm 2.1 can be implemented in parallel. For details about implementation, see [6, 7].
ON CONVERGENCE OF THE ASPIN METHOD
1855
3. Local convergence of the ASPIN method. We will prove in this section that the ASPIN method is locally convergent. Note that by Propositions 2.1 and 2.2, the ASPIN method is based on the local property of F (u) at the solution u∗ of system (1.1), so it seems impossible to obtain a global convergence result for this method. In this section and the following, · always denotes the Euclidean norm for both vectors and matrices, and N (u, ρ) = {v | v − u < ρ} represents the open ball with center u and radius ρ. Since the analysis of secondary iteration would complicate the discussion without gaining more insight into the method, we assume that (A1 ) the value of F at each iterative point is evaluated exactly, i.e., (2.1) holds with u replaced by u(k) + λk p(k) for each k; moreover, from now on, we assume that (A2 ) Ei J(u∗ )EiT is invertible for each i; (A3 ) D represents the neighborhood determined in Proposition 2.2; and (A4 ) δ > 0 is a fixed small number such that N (u∗ , δ) ⊂ D; in addition, the following inequalities hold for any u ∈ N (u∗ , δ): (I1 ) B(u) ≤ 2M ; (I2 ) B(u)−1 ≤ 2M ; (I3 ) J (u) ≤ 2M ; (I4 ) J (u)−1 ≤ 2M ; max (I5 ) J (u) − B(u) ≤ 4M1−η (1+ηmax ) ; 1 ∗ ∗ u − u∗ , (I6 ) F(u) − F(u ) − J (u )(u − u∗ ) ≤ 2M where M := max{J (u∗ ), J (u∗ )−1 }. It is easy to see from Remarks 2.1 and 2.2 that inequalities (I1 )–(I5 ) may hold with δ small enough. The last inequality may hold by Lemma 3.2.10 in [20]. In addition, we assume that the parameter α in Algorithm 2.1 is small enough so that 64αM 4 ≤
(3.1)
1 − ηmax . 3 + ηmax
It should be pointed out that the above assumptions are not so strict; see the appendix, where an example is given. Now we show that the inexact Newton direction computed in Algorithm 2.1 is also a regular inexact Newton direction for F in the sense that (1.3) holds. Proposition 3.1. Assume that u ∈ N (u∗ , δ). If (3.2)
F(u) + B(u)p ≤ ηF(u),
η ∈ [0, ηmax ],
then (3.3)
F(u) + J (u)p ≤
1+η F(u). 2
Proof. By (I2 ) and (3.2), p = B(u)−1 [B(u)p + F(u) − F(u)] ≤ (1 + η)B(u)−1 F(u) ≤ 2M (1 + η)F(u).
1856
HENG-BIN AN
Thus, according to (I5 ) and (3.2), J (u)p + F(u) = [J (u) − B(u)]p + B(u)p + F(u) ≤ J (u) − B(u)p + B(u)p + F(u) 1 − ηmax · 2M (1 + η)F(u) + ηF(u) ≤ 4M (1 + ηmax ) 1−η ≤ F(u) + ηF(u) 2 1+η ≤ F(u). 2 Thus we obtain the required inequality. Remark 3.1. If u ∈ N (u∗ , δ) and (3.2) holds, then we have F(u)T B(u)p = F(u)T [B(u)p + F(u) − F(u)] ≤ −F(u)2 + F(u)B(u)p + F(u) ≤ −(1 − η)F(u)2 < 0. In the same way, (3.3) shows that F(u)T J (u)p < 0.
(3.4)
In particular, (3.4) shows that if u ∈ N (u∗ , δ) and p is computed by the ASPIN method, then p is a descent direction for the function f (u) = 12 F(u)2 at point u.1 The following lemma is needed in our analysis. Lemma 3.2 (see [4, Lemma 3.4]). Let u ∈ Rn and H : Rn → Rn be continuously differentiable in a neighborhood of u. Assume that H(u) = 0 and H (u) is nonsingular. If p ∈ Rn such that H(u) + H (u)p ≤ ηH(u),
η ∈ [0, 1),
then |∇h(u)T p| 1−η ≥ ∇h(u) > 0, p (1 + η)κ(H (u)) where κ(H (u)) is the condition number for H (u) and h(u) = 21 H(u)T H(u). Based on Lemma 3.2, we have the following result. Lemma 3.3. Assume that u ∈ N (u∗ , δ). If F(u) + B(u)p ≤ ηF(u),
η ∈ [0, ηmax ],
then it holds that |F(u)T J (u)p| ≥ 4α|F(u)T B(u)p|. Proof. According to Proposition 3.1, we have F(u) + J (u)p ≤ 1p
1+η F(u); 2
is a descent direction for f (u) if ∇f (u)T p < 0. In addition, note that ∇f (u) = J (u)T F (u).
ON CONVERGENCE OF THE ASPIN METHOD
1857
therefore, by Lemma 3.2, (I3 ), and (I4 ), 1−η F(u)T J (u)p (3 + η)κ(J (u)) 1−η F(u)p ≥ (3 + η)κ(J (u))J (u)−1 1−η ≥ F(u)p. 8M 3 (3 + η)
|F(u)T J (u)p| ≥
(3.5)
Thus, by (I1 ), (3.5), and (3.1), 4α|F(u)T B(u)p| ≤ 4αB(u)F(u)p ≤ 8αM F(u)p 8M 3 (3 + η) ≤ 8αM |F(u)T J (u)p| 1−η 3 + ηmax ≤ 64αM 4 |F(u)T J (u)p| 1 − ηmax ≤ |F(u)T J (u)p|. This concludes the proof. The following lemma shows that if u is sufficiently close to u∗ , then the direction p obtained in the ASPIN method will not be too long. δ Lemma 3.4. Assume that u ∈ N (u∗ , 2δ ) with F(u) ≤ 8M . If p ∈ Rn such that F(u) + B(u)p ≤ ηF(u),
η ∈ [0, ηmax ],
then [u, u + p] ⊂ N (u∗ , δ), where [u, u + p] represents the line segment between u and u + p. δ Proof. Since F(u) ≤ 8M , by (I2 ), we have p = B(u)−1 [B(u)p + F(u) − F(u)] ≤ B(u)−1 · [B(u)p + F(u) + F(u)] ≤ 2M (1 + η)F(u) ≤ 2M (1 + ηmax )F(u) δ(1 + ηmax ) ≤ 4 δ < . 2 Therefore, by u ∈ N (u∗ , 2δ ), (u + p) − u∗ ≤ u − u∗ + p
0 such that (3.6)
∇f (v) − ∇f (w) ≤ γv − w
∀ v, w ∈ N (u∗ , δ).
If p ∈ Rn such that F(u) + B(u)p ≤ ηF(u),
η ∈ [0, ηmax ],
then the linesearch procedure along p in Algorithm 2.1 will terminate in finite iterations and the obtained λ satisfies αθmin |F(u)T J (u)p| λ ≥ min 1, . γp2 Proof. Because u ∈ N (u∗ , 2δ ), Lemma 3.3 shows that |F(u)T J (u)p| ≥ 4α|F(u)T B(u)p|.
(3.7)
Since F(u)T J (u)p < 0, F(u)T B(u)p < 0, (3.7) shows that F(u)T J (u)p ≤ 4αF(u)T B(u)p, so (3.8)
F(u)T J (u)p − αF(u)T B(u)p ≤ 3αF(u)T B(u)p < 0.
δ Because u ∈ N (u∗ , 2δ ) and F(u) ≤ 8M , Lemma 3.4 shows that [u, u + p] ⊂ ∗ N (u , δ). Thus, by the mean value theorem, there exists ξ ∈ [u, u + p] such that
f (u + λp) = f (u) + λ∇f (ξ)T p. Therefore, f (u + λp) = f (u) + λ∇f (ξ)T p = f (u) + αλF(u)T B(u)p − αλF(u)T B(u)p + λ∇f (ξ)T p = f (u) + αλF(u)T B(u)p + λ{[∇f (ξ)T p − ∇f (u)T p] +[∇f (u)T p − αF(u)T B(u)p]} = f (u) + αλF(u)T B(u)p + λ{λζ + [F(u)T J (u)p − αF(u)T B(u)p]}, where ζ=
∇f (ξ)T p − ∇f (u)T p . λ
By (3.6), we have |ζ| ≤ γp2 , so
f (u + λp) ≤ f (u) + αλF(u)T B(u)p + λ λγp2 + [F(u)T J (u)p − αF(u)T B(u)p] . Thus, if λγp2 + [F(u)T J (u)p − αF(u)T B(u)p] ≤ 0,
ON CONVERGENCE OF THE ASPIN METHOD
1859
then λ is acceptable. Since λ is reduced by a factor θ ≤ θmax < 1 at each iteration of the while-loop, it follows from (3.8) that the while-loop will terminate in finite steps. Let λ be the ultimate step factor. If λ = 1, then the needed conclusion trivially holds. Now suppose that the linesearch procedure is implemented at least once, and let λ− be the penultimate value; then the above argument shows that λ− γp2 + [F(u)T J (u)p − αF(u)T B(u)p] > 0. Consequently, λ− >
|F(u)T J (u)p| − α|F(u)T B(u)p| . γp2
Therefore, it follows from (3.7) that λ ≥ θmin λ− |F(u)T J (u)p| − α|F(u)T B(u)p| > θmin γp2 T 4α|F(u) B(u)p| − α|F(u)T B(u)p| ≥ θmin γp2 T αθmin |F(u) B(u)p| ≥ . γp2 Thus, we have obtained the required conclusion. Lemma 3.6. If u ∈ N (u∗ , 2δ ) and F(u)
F u(k+1) . Thus, by induction, the ASPIN method can generate a sequence {u(k) } ⊂ N (u∗ , 2δ ). Next we show that f u(k) → 0. By (3.10),
(3.11)
T T F u(k) B u(k) p(k) = F u(k) [B u(k) p(k) + F u(k) − F u(k) ] ≤ −2(1 − ηmax )f u(k) ,
and it follows that T f (u(k) + λk p(k) ) ≤ f u(k) + αλk F u(k) B u(k) p(k) ≤ [1 − 2αλk (1 − ηmax )]f u(k) . Therefore, by Theorem 3.5, we have f (u
(k)
+ λk p
(3.12)
(k)
T αθmin |F u(k) B u(k) p(k) | (k) f u ) ≤ 1 − 2α(1 − ηmax ) γp(k) 2 ≡ (1 − ctk )f u(k) ,
where we set c=
2α2 θmin (1 − ηmax ) γ
and T |F u(k) B u(k) p(k) | . tk = p(k) 2 Because {f u(k) } is nonnegative and strictly decreased, limk→∞ f u(k) exists. (k) If limk→∞ f u > 0, then (3.12) shows that f u(k+1) ≤ lim (1 − ctk ) ≤ 1, lim k→∞ f u(k) k→∞ so lim (1 − ctk ) = 1
k→∞
or, equivalently, lim tk = 0.
k→∞
1862
HENG-BIN AN
Since (3.11) shows that T |F u(k) B u(k) p(k) | 2(1 − ηmax )f u(k) ≤ = tk , p(k) 2 p(k) 2 we have
f u(k) lim = 0. k→∞ p(k) 2
Because limk→∞ f u(k) > 0, it follows that (3.13)
p(k) → ∞ (k → ∞).
But on the other hand, by (I2 ) and (3.10), we have −1 (k) (k) p(k) = B u(k) [B u p + F u(k) − F u(k) ] −1 ≤ B u(k) · [B u(k) p(k) + F u(k) + F u(k) ] ≤ 2M (1 + ηk )F u(k) ≤ 2M (1 + ηmax )F(u(0) ) δ(1 + ηmax ) , ≤ 4 which contradicts (3.13). Thus, we must have f (u(k) ) → 0. Since {u(k) } ⊂ N (u∗ , δ), by (I6 ), u(k) − u∗ = 2u(k) − u∗ − u(k) − u∗
≤ 2J (u∗ )−1 J (u∗ )(u(k) −u∗ )−2M F u(k) −F(u∗ )−J (u∗ )(u(k) −u∗ ) ≤ 2M [J (u∗ )(u(k) − u∗ ) − F u(k) − F(u∗ ) − J (u∗ )(u(k) − u∗ )] ≤ 2M F u(k) . Because F(u(k) ) = 2f (u(k) ) → 0, we have u(k) → u∗ . Now, we complete the discussion for local convergence of the ASPIN method. 4. Convergence rate of the ASPIN method. In this section, we discuss the convergence rate of the ASPIN method. We will show that the ASPIN method is quadratically convergent under suitable conditions. The following theorem shows that under suitable assumptions, λk = 1 is acceptable for all k sufficiently large. Theorem 4.1. Assume that f (u) is twice continuously differentiable in N (u∗ , δ), and there exists γ > 0 such that for any v, w ∈ N (u∗ , δ), ∇f (v) − ∇f (w) ≤ γv − w, ∇2 f (v) − ∇2 f (w) ≤ γv − w. Let {u(k) } be the sequence generated by the ASPIN method such that u(k) → u∗ . If ηk → 0, then u(k+1) = u(k) + p(k) for all sufficiently large k. Proof. Because u(k) → u∗ , without loss of generality, we may assume that (k) {u } ⊂ N (u∗ , 2δ ). Thus, it follows from Proposition 3.1 that (4.1)
1 + ηmax (k) F u(k) + J u(k) p(k) ≤ F u . 2
ON CONVERGENCE OF THE ASPIN METHOD
1863
Therefore, by (I4 ) and (4.1),
(4.2)
−1 (k) (k) p + F u(k) − F u(k) ] [J u p(k) = J u(k) −1 ≤ J u(k) · [J u(k) p(k) + F u(k) + F u(k) ] 1 + ηmax + 1 F u(k) ≤ 2M 2 = M (3 + ηmax )F u(k) .
Because F(u(k) ) → 0, the above inequality shows that p(k) → 0 (k → ∞). By (I4 ), we have T −1 −1 (k) 1 ∇f u(k) = J u(k) F u(k) ≥ J u(k) F u(k) ; ≥ F u 2M thus, Lemma 3.2 in connection with (I3 ), (I4 ), and (4.1) shows that T |∇f u(k) p(k) | 1 − 1+η2max 1 · ∇f u(k) ≥ 1+ηmax · (k) (k) p ) κ(J u 1+ 2 1 − ηmax 1 1 ≥ · · F u(k) 2 3 + ηmax 4M 2M 1 − ηmax = F u(k) , 8M 3 (3 + ηmax ) or we have (4.3)
T 8M 3 (3 + ηmax ) F u(k) p(k) ≤ · |∇f u(k) p(k) | 1 − ηmax T ≡ a|∇f u(k) p(k) |,
where a=
8M 3 (3 + ηmax ) . 1 − ηmax
From (4.2) and (4.3), we obtain p(k) 2 ≤ M (3 + ηmax )F u(k) p(k) T 8M 4 (3 + ηmax )2 ≤ · |∇f u(k) p(k) | 1 − ηmax T ≡ b|∇f u(k) p(k) |, where b=
8M 4 (3 + ηmax )2 . 1 − ηmax
1864
HENG-BIN AN
Next we show that λk = 1 is acceptable for all k sufficiently large. First note that if F = (F1 , F2 , . . . , Fn )T , then ∇2 f (u) = J (u)T J (u) +
n
Fi (u)∇2 Fi (u)
i=1
≡ J (u) J (u) + S(u). T
Since u(k) → u∗ and F(u∗ ) = 0, it follows that S(u(k) ) → 0. For each k, by the mean value theorem, there exists a ξ (k) on the line segment between u(k) and u(k) + p(k) such that 1 T T 1 f u(k) + p(k) − f u(k) − ∇f u(k) p(k) = ∇f u(k) + ∇2 f ξ (k) p(k) p(k) . 2 2 This gives (k) f u + p(k) − f u(k) − 1 ∇f u(k) T p(k) 2 1 T = ∇f u(k) + ∇2 f (ξ (k) )p(k) p(k) 2 T 1 = ∇f u(k) + ∇2 f u(k) p(k) p(k) + (p(k) )T ∇2 f (ξ (k) ) − ∇2 f u(k) p(k) 2
2 T (k)
1 ≤ J u(k) F u + J u(k) p(k) · p(k) + S u(k) + γ p(k) p(k)
2
2
1 ≤ ηk J u(k) · F u(k) · p(k) + S u(k) + γ p(k) p(k)
2
T 1 ≤ − 2aM ηk + b S u(k) + γ p(k) ∇f u(k) p(k) 2 T 1 ≡ − k ∇f u(k) p(k) ; 2 therefore, (4.4)
1 T f (u(k) + p(k) ) − f u(k) ≤ (1 − k )∇f u(k) p(k) . 2
Since ηk , S(u(k) ), and p(k) all converge to zero, it follows that k → 0. Thus, for all k sufficiently large, we have k
0 such that ∇f (v) − ∇f (w) ≤ γv − w, ∇2 f (v) − ∇2 f (w) ≤ γv − w for any v, w ∈ N (u∗ , δ). If {u(k) } is a sequence generated by the ASPIN method such that u(k) → u∗ and u(k+1) = u(k) + p(k) for all sufficiently large k, then (i) u(k) → u∗ superlinearly if and only if rk = o(F u(k) ), where rk = F u(k) + B u(k) p(k) ; (ii) u(k) → u∗ quadratically if and only if rk = O(F u(k) 2 ).
ON CONVERGENCE OF THE ASPIN METHOD
1867
Proof. Since u(k+1) = u(k) + p(k) for all sufficiently large k, without loss of generality, we assume that u(k+1) = u(k) + p(k) for all k in the following argument. Assume that u(k) → u∗ superlinearly. Since rk = F u(k) + B u(k) p(k) = F u(k) − F(u∗ ) − J (u∗ ) u(k) − u∗ − B u(k) − J (u∗ ) u(k) − u∗ + B u(k) u(k+1) − u∗ , (4.10) by Lemma 3.2.10 in [20], continuity of B(u) at u∗ and the assumption that u(k) → u∗ superlinearly, we have rk ≤ o(u(k) − u∗ ) + o(1)u(k) − u∗ + o(u(k) − u∗ ). Therefore, by Lemma 4.3, it follows that rk = o(u(k) − u∗ ) = o(F u(k) ). Conversely, assume that rk = o(F u(k) ). Since
(4.11)
u(k+1) − u∗ = (u(k) − u∗ ) + p(k) −1 = (u(k) − u∗ ) + B u(k) [rk − F u(k) ] −1 (k) B u − J (u∗ ) (u(k) − u∗ ) = B u(k) + rk − F u(k) − F(u∗ ) − J (u∗ )(u(k) − u∗ ) ,
thus, by (I2 ), the continuity of B(u) at u∗ , the assumption that rk = o(F u(k) ), and Lemma 3.2.10 in [20], we have u(k+1) − u∗ ≤ 2M [o(1)u(k) − u∗ + o(F u(k) ) + o(u(k) − u∗ )]. Therefore, by Lemma 4.3, u(k+1) − u∗ = o(u(k) − u∗ ) + o(F u(k) ) = o(u(k) − u∗ ). Since F is twice continuously differentiable in N (u∗ , δ), Lemma 4.2 shows that both J (u) and B(u) are Lipschitz continuous in a neighborhood V ⊂ N (u∗ , 2δ ) of u∗ . Thus, there exists L > 0 such that for any u ∈ V , J (u) − J (u∗ ) ≤ Lu − u∗
(4.12) and
B(u) − B(u∗ ) ≤ Lu − u∗ .
(4.13)
By (4.12) and Lemma 3.2.12 in [20], (4.14)
F(u) − F(u∗ ) − J (u∗ )(u − u∗ ) ≤
L u − u∗ 2 2
∀ u ∈ V.
In a similar way, it follows from (4.10), (4.11), (4.13), (4.14), and Lemma 4.3 that u(k) → u∗ quadratically if and only if rk = O(F u(k) 2 ). This concludes the proof.
1868
HENG-BIN AN
From Theorems 4.1 and 4.4, we can obtain the following result. Corollary 4.5. Assume that both F (u) and f (u) are twice continuously differentiable in N (u∗ , δ), and there exists γ > 0 such that for any v, w ∈ N (u∗ , δ), ∇f (v) − ∇f (w) ≤ γv − w, ∇2 f (v) − ∇2 f (w) ≤ γv − w. Let {u(k) } be a sequence generated by the ASPIN method such that u(k) → u∗ . Then (i) u(k) → u∗ superlinearly if ηk → 0; (ii) u(k) → u∗ quadratically if ηk = O(F u(k) ). Corollary 4.5 reflects how the forcing term influences the convergence rate of the ASPIN method. This result is similar to Corollary 3.5 in [8] for the inexact Newton method. In particular, by Corollary 4.5, we can determine the convergence rate of the ASPIN method by choosing proper forcing terms. 5. Conclusion. The inexact Newton method is one of the effective tools for solving large sparse systems of nonlinear equations. By using nonlinear additive Schwarz preconditioning technique, Cai and Keyes [6] introduced the ASPIN method. This method is very effective for solving some nonlinear problems with strong nonbalanced nonlinearities. However, the convergence of the ASPIN method is not discussed by them or others. In this paper, we discussed the convergence property of the ASPIN method and thus we provided a theoretical support for the ASPIN method. The convergence result is local since the design of the ASPIN only concerns the local properties of the original function. Furthermore, we discussed the convergence rate for the ASPIN method, and the result shows that the convergence rate of the ASPIN method is similar to that of the inexact Newton method. Thus, we can obtain the desired convergence rate by choosing proper forcing terms. Appendix. We give a simple example to show that our main assumptions, inequalities (I1 )–(I6 ) in section 3, are not so strict. Consider the nonlinear equations 2u1 − u2 + λeu1 − λ = 0, −u1 + 2u2 + λeu2 − λ = 0, where λ > 0. It is obvious that u∗ = (0, 0)T is a solution of the system. Let S = {1, 2},
S1 = {1},
S2 = {2}
and T1 (u) =
T11 (u) , 0
T2 (u) =
0 , T22 (u)
where T11 , T22 : R2 → R. By FSi (u − Ti (u)) = 0, we have 2(ui − Tii (u)) − u3−i + λeui −Tii (u) − λ = 0,
(A.1)
i = 1, 2.
It is easy to see that Tii (u) are continuous functions. Furthermore, we can obtain T1 (u)
=
−1 1 − 2 + λeu1 −T11 (u) , 0 0
T2 (u)
=
0
− 2 + λeu2 −T22
(u) −1
0 . 1
ON CONVERGENCE OF THE ASPIN METHOD
Thus,
1
J (u) =
−1 − 2 + λeu2 −T22 (u)
Besides, it is easy to see that B(u) =
1 −1 − (2 + λeu2 )
Since B(u∗ ) = J (u∗ ) =
1869
−1 − 2 + λeu1 −T11 (u) . 1 −1
− (2 + λeu1 ) 1
1 −(2 + λ)−1
.
−(2 + λ)−1 , 1
so M = max{J (u∗ ), J (u∗ )−1 } =
2+λ . 1+λ
Now we will choose some proper δ such that inequalities (I1 )–(I6 ) hold. Let δ ∈ (0, 12 ) and assume that u = (u1 , u2 )T ∈ N (u∗ , δ). Note that if 1 a A= , b 1 then A ≤ 1 + max{|a|, |b|}. Therefore, B(u) ≤ 1 + max (2 + λeui )−1 < 1 + (2 + λe−δ )−1 < 2M. 1≤i≤2
In addition, since B(u)−1 =
1 u 1 1 − (2 + λe )−1 (2 + λeu2 )−1
1 (2 + λeu2 )−1
(2 + λeu1 )−1 , 1
one can easily verify that B(u)−1 ≤
2 + λe−δ < 2M. 1 + λe−δ
Because u ∈ N (u∗ , δ), by (A.1), we have |ui − Tii (u)|
1 + x,
x ∈ (0, 1),
1 and ex < 1 + x, 2
x ∈ (−1, 0),
we have |ui − Tii (u)| ≤
2 |u3−i |, 4+λ
i = 1, 2.
Therefore,
T11 (u) − u1 + (2 + λ)−1 u2
F(u) − F(u∗ ) − J (u∗ )(u − u∗ ) =
T22 (u) − u2 + (2 + λ)−1 u1
|T11 (u) − u1 | + (2 + λ)−1 |u2 |
≤
|T22 (u) − u2 | + (2 + λ)−1 |u1 |
3 ≤ u < 2M u − u∗ . 2+λ Thus, inequality (I6 ) holds. Summing up the above discussion, we know that if 1 (1 + λ)2 (1 − ηmax ) , , δ ≤ min 2 3eλ(2 + λ)(1 + ηmax ) then for any u ∈ N (u∗ , δ), inequalities (I1 )–(I6 ) hold. For this example, other assumptions in the paper can also be checked to be true. Acknowledgments. The author thanks the anonymous referees and Carol S. Woodward for many valuable suggestions. REFERENCES [1] S. Bellavia and B. Morini, A globally convergent Newton-GMRES subspace method for systems of nonlinear equations, SIAM J. Sci. Comput., 23 (2001), pp. 940–960. [2] S. Bellavia, M. Macconi, and B. Morini, A hybrid Newton-GMRES method for solving nonlinear equations, in Numerical Analysis and Its Applications, Lecture Notes in Comput. Sci. 1988, L. Vulkov, J. Wasniewski, and P. Yalamov, eds., Springer-Verlag, Berlin, 2000, pp. 68–75. [3] P. N. Brown and Y. Saad, Hybrid Krylov methods for nonlinear systems of equations, SIAM J. Sci. Statist. Comput., 11 (1990), pp. 450–481. [4] P. N. Brown and Y. Saad, Convergence theory of nonlinear Newton-Krylov algorithms, SIAM J. Optim., 4 (1994), pp. 297–330. [5] X.-C. Cai, M. Dryja, and M. Sarkis, Restricted additive Schwarz preconditioners with harmonic overlap for symmetric positive definite linear systems, SIAM J. Numer. Anal., 41 (2003), pp. 1209–1231. [6] X.-C. Cai and D. E. Keyes, Nonlinearly preconditioned inexact Newton algorithms, SIAM J. Sci. Comput., 24 (2002), pp. 183–200. [7] X.-C. Cai, D. E. Keyes, and L. Marcinkowski, Nonlinear additive Schwarz preconditioners and applications in computational fluid dynamics, Internat. J. Numer. Methods Fluid Mech., 40 (2002), pp. 1463–1470.
ON CONVERGENCE OF THE ASPIN METHOD
1871
[8] R. S. Dembo, S. C. Eisenstat, and T. Steihaug, Inexact Newton methods, SIAM J. Numer. Anal., 19 (1982), pp. 400–408. [9] J. E. Dennis and R. B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice–Hall, Englewood Cliffs, NJ, 1983. [10] M. Dryja and W. Hackbusch, On the nonlinear domain decomposition method, BIT, 37 (1997), pp. 296–311. [11] M. Dryja and O. B. Widlund, Domain decomposition algorithms with small overlap, SIAM J. Sci. Comput., 15 (1994), pp. 604–620. [12] S. C. Eisenstat and H. F. Walker, Globally convergent inexact Newton methods, SIAM J. Optim., 4 (1994), pp. 393–422. [13] S. C. Eisenstat and H. F. Walker, Choosing the forcing term in an inexact Newton method, SIAM J. Sci. Comput., 17 (1996), pp. 16–32. [14] D. R. Fokkema, G. L. G. Sleijpen, and H. A. van der Vorst, Accelerated inexact Newton schemes for large systems of nonlinear equations, SIAM J. Sci. Comput., 19 (1998), pp. 657–674. [15] I. E. Kaporin and O. Axelsson, On a class of nonlinear equation solvers based on the residual norm reduction over a sequence of affine subspaces, SIAM J. Sci. Comput., 16 (1995), pp. 228–249. [16] D. A. Knoll and D. E. Keyes, Jacobian-free Newton-Krylov methods: A survey of approaches and applications, J. Comput. Phys., 193 (2004), pp. 357–397. [17] S. H. Lui, Nonlinearly preconditioned Newton’s method, in Domain Decomposition Methods in Science and Engineering, Fourteenth International Conference on Domain Decomposition Methods, Cocoyoc, Mexico, I. Herrera, D. E. Keyes, O. B. Widlund, and R. Yates, eds., National Autonomous University of Mexico (UNAM), Mexico City, Mexico, 2003, pp. 95– 105. [18] L. Luksan, Inexact trust region method for large sparse systems of nonlinear equations, J. Optim. Theory Appl., 81 (1994), pp. 569–591. [19] M. Pernice and H. F. Walker, NITSOL: A Newton iterative solver for nonlinear systems, SIAM J. Sci. Comput., 19 (1998), pp. 302–318. [20] J. M. Ortega and W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Classics Appl. Math. 30, SIAM, Philadelphia, PA, 2000. [21] Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Comput., 7 (1986), pp. 856–869.
