A safeguard approach to detect stagnation of GMRES(m) with applications in Newton-Krylov methods

A safeguard approach to detect stagnation of GMRES(m) with applications in Newton-Krylov methods^*

Márcia A. Gomes-Ruggiero; Véra L. Rocha Lopes; Julia V. Toledo-Benavides

Department of Applied Mathematics, IMECC-UNICAMP, University of Campinas CP 6065, 13083-970 Campinas SP, Brazil, E-mails: [e-mail] / [e-mail] / [e-mail]

ABSTRACT

Restarting GMRES, a linear solver frequently used in numerical schemes, is known to suffer from stagnation. In this paper, a simple strategy is proposed to detect and avoid stagnation, without modifying the standard GMRES code. Numerical tests with the proposed modified GMRES(m) procedure for solving linear systems and also as part of an inexact Newton procedure, demonstrate the efficiency of this strategy.

Mathematical subject classification: 65H10, 65F10.

Key words: linear systems, restarting GMRES, inexact Newton method, nonlinear systems.

1 Introduction

The objective of this work is to improve the performance of the restarted GMRES [25]. A well known difficulty with the restart of GMRES, algorithm for solving Ax = b, A: n × n is that it can stagnate when the matrix A is not positive definite [24], in the sense that the residual sequence does not converge to zero within a reasonable time frame. Simoncini, [27], and Sturler, [31], modified the GMRES using spectral analysis. In Parks et al., [21], Ritz values are used to improve the performance of the restart. Morgan, [16], [17], also uses eigenvectors for improving convergence of restarted GMRES. He considers some vectors from the previous subspace, adding them to the new subspace in order to deflate eigenvalues. Convergence issues and stagnation are discussed by Simoncini, [29], and Zavorin, [36]. Van der Vorst and Vuik introduce a strategy to prevent stagnation in GMRES, by including LSQR steps in some phases of the process [33].

This work aims an early detection of stagnation; once detected stagnation, the initial residue of the next cycle is steered away from the stagnation zone by means of a simple hybrid safeguard which mostly involves, in addition, the current residue. The strategies proposed here have the following objectives:

• avoid stagnation;

• use previous information given by the GMRES, avoiding any modification in it. At most, the new program should ask for few information besides that usually provided by GMRES;

• take into account information from the previous cycles performed by the GMRES.

It is also important to analyze the GMRES as a solver for the linear systems generated by Inexact-Newton type methods for solving nonlinear systems of equations. Our interest in choosing these methods relies, basically, on their popularity among practitioners, and on our research interest in Newton-Krylov methods, [10, 32]. Since strategies that improve the efficiency of GMRES, are fundamental in a better performance of these methods, they became one of the main subjects of our research.

In Section 2 we describe briefly the GMRES and the restarted GMRES algorithms. We also study the effect of the GMRES cycle length on the decrease of the residual norm. In Section 3 we establish our stagnation criteria and describe hybrid safeguards which modify the GMRES method, obtaining a version called here GMRESH. In Section 4 we show that GMRESH is capable to reduce considerably the effect of stagnation at the resolution of some linear systems. We also compare GMRESH with the GMRESDR method proposed by Morgan, [17]. In Section 5 we discuss the implementation of GMRESH within the Newton-Krylov method and test its performance on a ray-tracing problem and on a set of boundary value problems. Some concluding remarks are given in Section 6.

2 GMRES(m)

The method GMRES was proposed in [25] for solving linear systems As = b, where A is a nonsingular n × n matrix (not necessarily symmetric) and b

ⁿ. If s₀ is a initial approximation for the solution and r₀ = b - As₀ is the corresponding residual vector, the Krylov subspace after l iterations of the GMRES will be:

At each iteration l of GMRES, a value s_l s₀ + _l is computed to minimize the residual vector, that is: r_l = ||b - As||₂. In what follows we always mean ||.||₂ whenever we use ||.||.

It is known that, computationally speaking, GMRES is more expensive than other Krylov subspace methods, such as Bi-CGSTAB, [14], QMR [24] for general square matrices, or LSQR [19], [20] for anti-symmetric matrices. Nevertheless, it is widely used for solving linear systems derived from the discretization of partial differential equations, since theoretically the 2-norm of the residual vector is minimized inside the Krylov subspace at each step.

We can describe GMRES as follows: given the subspace

_l and the initial approximation s₀, compute s_l, the approximate solution for As = b, where s_l s₀ + _l in such a way that r_l = b - As_l is orthogonal to A

_l.

Since s_l s₀ + _l we can write s_l = s₀ + δ, δ

_l; then r_l = b - As_l = r₀ - Aδ. We obtain δ in such a way that r_l is orthogonal to A

_l. Geometrically, Aδ is the orthogonal projection of r₀ in A

_l, as shown in Figure 1.

The dimension of the Krylov space keeps increasing, so the memory cost and complexity of the l-th GMRES step increase with l. A modified version called GMRES(m) is used in large scale problems. In this version, the GMRES proceeds in cycles of m iterations, see [14], [25]. Basically, the process begins with some vector s₀, and a fixed number m of iterations are performed. Then, a new cycle begins with s_m as initial approximation and r_m = b - As_m as initial residue. Note that at each cycle an m-dimensional Krylov subspace is generated fromthe initial residue, following the usual GMRES procedure.

Whereas the restarted policy is computationally more feasible, convergence cannot be guaranteed in general, and stagnation becomes possible [11], [24], [28], [31] and [36]. A rather expensive remedy would be to monitor the eigenvalues of the Hessenberg matrices generated during the GMRES, [28]. Other schemes, such as the one mentioned in [31], store some vectors created at the j-th cycle and use them at the (j + 1)-th cycle. We present a different strategy.

3 Stagnation

In this section we present a new strategy to generate an approximation to the new cycle that bypasses the stagnation of the method. We use the following notation: is the initial residue of the j-th cycle and is the residual vector at the end of this cycle.

Stagnation in GMRES(m) is usually described as slowness in the decrease of the consecutive residual norms, ||||, l = 1,2,3,..., m. However, a situation where and are roughly linearly dependent could also be classifiedas stagnant.

To prevent stagnation we need to control the cycles in such a way that it is possible to make a comparison between the norm of the last residue, ||||, and the norm of the initial residual vector of this cycle, ||||. Moreover, we need to guarantee that the basis generated for the Krylov subspace in the (j + 1)-th cycle is linearly independent with respect to the basis from the last cycles.

Firstly, we need to establish a criterion to detect stagnation, which can be based on the norm of the residual vectors. However, even in the case of a reasonable decrease of the norm of residues, if the angle between the initial and final residual vectors of one cycle j is close to zero, the Krylov subspace of the new cycle, (j + 1), can be similar to the previous one. This is because both subspaces began with vectors that are almost linearly dependent. It is obvious that in such case the progress of the process towards the solution will be very slow.

In Figure 1, we can observe that, given the projection property of GMRES, if the cosine of the angle between the initial and final residual vectors is closeto 1, then the norm of these residues are very close to each other so that there is an equivalence between the tests

and

where θ_j is the angle between the vectors and .

We must also consider the possibility of linear dependence between , l = 1, ..., j - 1 and , where is the last residue of the j-th cycle whereas is the initial residue of the l-th cycle. Nevertheless we decided to make the comparison only between and . Our stagnation criterion is based on simplicity and heuristics. Simplicity dictates a small number of angle comparisons; heuristics is justified as follows: comparison with is justified by the fact that, empirically, the largest descent occurs during the first GMRES cycles. Angle comparison with is based on the heuristic assumption that once stagnation is present it will always repeat itself; in other words, if the angle between and were small, stagnation would have already been discovered in the previous cycle.

If a linear dependence between and occurs, the Krylov subspace of the cycle (j + 1) would be close to the Krylov subspace generated in the first cycle. This would lead to the stagnation of the process. If this is the case, there is no equivalence between the tests (2) and (3), since and belong to different Krylov subspaces. Thus it is not possible to use the test with the norms of the residual vectors. Linear dependence can be detected testing the cosine of the angle between the residues and :

In the strategies proposed in this work, the analysis will be always done at the end of each cycle, to reduce a too big computational costs. Stagnation will be declared when:

In case of stagnation, another initial approximation must be chosen for the new cycle. This new approximation is obtained using information generated during the process. However, it needs to be constructed in such a way to guarantee a reduction in the norm of the residual vector. The new strategy that is proposed generates an approximation for the (j + 1)-th cycle using a hybrid scheme, based on the strategy proposed by Brezinski and Redivo-Zaglia in 1994, [3]. This hybrid scheme uses the approximations and which, in some sense, take into account the information generated by GMRES(m). In this way, we are trying to get out of a sequence which yields little decrease for the residuals. In what follows, we will briefly describe the scheme proposed in [3] for linear systems.

Consider the linear system of equations As = b. Once the approximations and ŝ are known, and the corresponding residues, = b - A and = b - Aŝ are computed, the objective is to construct a new approximation as a linear combination of and ŝ, s = α + βŝ, with the aim of reducing the residual norm. As a simplifying tool in obtaining these parameters is to fix β = 1 - α, and then the corresponding residue r will be given by

r = b - As = α + (1- α).

Therefore our problem is reduced to find α , the least square solution for:

min_α ||r|| = min_α ||( - ) α + ||,

for which the optimal α is given by:

The new approximation will be s = α + (1 - α)ŝ, and the corresponding residue is r = α + (1 - α)^{^}r.

Let us go back to the solution for As = b by the GMRES(m). Let and denote the initial approximation of the whole process, and the last approximation obtained after performing the m GMRES iterations of the j-th cycle, respectively. The following safeguards are tested and the computation of the new approximation is done using the hybrid scheme, where corresponds to and ŝ to .

Strategy H:

if j ≠ 1:

Safeguard 1: test the angle θ_j between and :

if |cos(θ_j)| ~ 1, compute = α + (1 - α) , with a given by (6).

Otherwise,

Safeguard 2: test the angle θ_j,1 between and :

if |cos(θ_j,1)| ~ 1, then compute = α + (1 - α) , with a given by (6).

if j = 1:

Test the angle θ₁ between and :

if |cos(θ₁)| ~ 1, compute s_a as a random vector and = αs_a + (1 - α),

with a given by (6).

In the case j = 1, due to the orthogonality and minimization properties of GMRES, the vector calculated from is the same as the one encountered by the hybrid process. Thus, should be modified. We add that the corresponding residual vector r_a, related to s_a, is normalized so as to guarantee the monotone decrease of the residues.

In Figure 2 we depict the situation tested by Safeguard 2, when the decrease in the residual norm is sufficient, so that Safeguard 1 is not triggered. Thus, is necessarily much smaller than . If indeed the angle between them is small, or close to π, the residue formed by them will show a pronounced decrease in norm.

The hybrid scheme also presents the advantage of maintaining the minimization of the 2-norm of r, as is the purpose of GMRES(m). The important difference here is that the GMRES(m) solves this problem in the Krylov subspace generated in the last cycle of GMRES(m). In the hybrid scheme the minimization is carried out in the plane generated by vectors belonging to two different Krylov subspaces: the first Krylov subspace (₁) and the last Krylov subspace (_j). Therefore our scheme keeps some information about these two subspaces.

The hybrid vectors are used quite often. This idea is similar to the residual smoothing used in [26], but in that case, it is used as a stopping criterion for the Conjugate Gradient method; in [34] the authors compared the behavior of the smooth residue and the usual residue for the MRS, QMR and BCG methods. Hybrid preconditioners are used in [23] and [33]. We do not use the term ''hybrid GMRES'' since in the literature it is sometimes used in other contexts, such as in [18], [30]. So we are calling it GMRESH.

4 Numerical experiments with GMRES and GMRESH

We present two examples comparing GMRES and GMRESH for 3 × 3 matrices.

Our procedure in this paper is the following: a hybrid restart is triggered in the first 5 occurrences of cos(θ_j) > 0.8 or cos(θ_j,1) > 0.8 and in the next 5 occurrences of cos(θ_j) > 0.9 or cos(θ_j,1) > 0.9. In all other next iterations, the usual GMRES(m) is used. The point is that if GMRESH(m) shows persistent stagnation then further progress is not likely to occur.

Example 1. Zavorin [35] brings the linear system Ax = b:

and

For this example, the null vector was taken as the initial approximation. GMRES was applied with restarts at each two iterations and it was allowed 100 cycles. The sequence of GMRES(2) residues is constant and equal to 1, showing the stagnation of this method. For the GMRESH(2), the safeguards were triggered and the process converged in 19 iterations with precision 10^-4.

We also tested GMRES(2) and GMRESH(2) on 64 problems, which were created as follows: consider the systems As = b_ijl where b_ijl = b + v_ijl, v_ijl being a vector in ³ with entries in [-0.1, 0.1]. We used 4 points in each interval. Figure 3 shows the logarithmic relative error ||r_end||/|||| for each problem, where r_end is the last residue of the whole process. Although both methods stagnated in some cases, GMRESH(2) shows a clear improvement.

Example 2. Embree [8] gives the linear system Ax = b:

as an example where GMRES(1) converges in three iterations but for GMRES(2) the relative residue stagnates near 0.3. We considered 1681 linear systems, As = b_nμ, where

= b_nμ: = (n,μ, 1)^T, n, μ [-10, 10]

and 10^-6 was taken as precision. Figure 4 shows the logarithm of the relative residual norms, (||r_end||/||||), ranging from 0 (white) to -8 (black). Actually, some of the GMRESH relative residues calculated were less than 10^-16. We can see that the GMRESH(2) relative residue is much smaller than the GMRES(2) relative residue in the vast majority of cases. Observing the size of the white region in both graphics, it is easy to conclude the better performance of GMRESH(2).

We applied the deflated GMRESDR method proposed by Morgan in [17] at the resolution of the linear systems (7) and (8). The results were obtained using the program GMRESDR.m, kindly sent to us by R. Morgan (Baylor University). The deflation of small eigenvalues can greatly improve the convergence of restarted GMRES. This method is denoted by GMRESDR(m, k), where m represents the maximum dimension of the subspace and k be the number of approximate eingenvectors (harmonic Ritz values) retained at a restart.

We observe that the matrices A of these systems are not ill-conditioned. At the resolution of system (7), GMRESDR(2,1) performed 69 inner iterations and obtained relative residual of 5.8544e-005. Meanwhile, GMRESH(2) performed 19 inner iterations as seen before. Considering the 1681 problems of type (8), we had difficulties in running GMRESDR(2,1) due to the ill conditioning introduced by the change of the eigenvalue done by this method.

In conclusion, GMRESH(2) did a better job in avoiding stagnation in these problems.

5 Numerical experiments

Here we apply the GMRES and GMRESH procedures, as linear solvers for the inexact Newton method with a nonmonotone line search, [1], for solving a nonlinear system

F(x) = 0, F: ⁿ → ⁿ.

In the inexact Newton methods, [5], the sequence x_k (called sequence of outer iterations) is generated by

x_k+1 = x_k + s_k;

s_k solves approximately the linear system J(x_k)s = -F(x_k), using this stopping criterion:

where J(.) is the Jacobian matrix, η_k (0,1] is the tolerance which is called the forcing term, [7].

In the line search procedure, it is needed the following parameters: σ (0 , 1), _min and _max such that 0 < _min < _max < 1 and the sequence {μ_k} such that μ_k > 0 for all k = 0,1,2,... and μ_k = μ < ∞.

Now we present the inexact Newton algorithm with the nonmonotone line search procedure. Let x₀

ⁿ be an arbitrary initial approximation to the solution for F(x) = 0. Given x_k

ⁿ, and the tolerance ε > 0, the steps to obtain a new iteration x_k+1 are the following:

Algorithm 1. (Inexact Newton method with nonmonotone line search):

While ||F(x_k)|| > ε, perform steps 1 to 4:

Step 1: Choose η_k.

Step 2: Find s_k such that ||F(x_k) + J(x_k)s_k|| < η_k||F(x_k)||;

Step 3: Take ξ = 1, compute x_aux = x_k + s_k and F(x_aux).

While

||F(x_aux)|| > [1 - ξσ]||F(x_k)|| + μ_k,

perform the steps 3.1 and 3.2:

step 3.1: compute ξ_new [_minξ, _maxξ];

step 3.2: set ξ = ξ_new and compute x_aux = x_k + ξs_k.

Step 4: Take ξ_k = ξ, compute x_k+1 = x_k+ ξ_ks_k and update k.

In the Step 1 we examine the following choices for the forcing term:

Constant (Cte): we chose η_k = 0.1;

EW1: η_k = (see Eisenstat and Walker [7]);

EW2: η_k = , γ [0,1], α (1,2] (see Eisenstat and Walker [7]);

GLT: η_k = [1/(k + 1)]^ρ cos²(ø_k) , ρ > 1 and -π/2 < ø_k < 0, (see Gomes-Ruggiero et al. [10]).

The vector s_k of the Step 2 is obtained by GMRES(m) and GMRESH(m).

The line search performed at Step 3 by Algorithm 1 follows the one proposed in [1] which is a nonmonotone strategy similar to the one introduced by Li and Fukushima, [12]. So, Algorithm 1 has global convergence [10]. Besides that, with the choices EW1, EW2 and GLT the convergence rate is superlinear, [7], [10].

5.1 Implementation features

We give now more details about the implementation of the algorithms. The implementation details can be found in [32], pages 26 and 49. All the tests were performed in a Centrino Duo 1.6GHz with 1 GB Ram computer, using the software MatLab 6.1.

• Line search procedure (Step 3 of Algorithm 1):

  for the parameter σ, we took 1.d-04;

  if the vector x_aux = x_k + ξs_k does not give an acceptable decrease in the value of the function, then we compute the new step size as ξ_new = 0.5ξ; for the sequence μ_k, we define:

  f tip(0) = ||F(x₀)||,

  f tip(k) = min{||F(x_k)||, f tip(k - 1)}, if k is a multiple of 3 and

  f tip(k) = f tip(k - 1), otherwise;

then, we set:

• The initial value and safeguards for η:

  for all the choices for η_k we set the initial value η₀ = 0.1. For the choices EW1 and EW2 of [7] and for the choice GLT, we take η_k = min{η_k, 0.1} if k < 3, and η_k = min{η_k, 0.01} if k > 3. We also take η_k = 0.1 when ø_k > 0. At the final iterations we have adopted the safeguard introduced in [22] which can be described as: since the linear model is F(x) ~ F(x_k) + J(x_k)s, at the final iterations, we can have: ||F(x_k+1)|| ~ ||F(x_k) + J(x_k)s_k|| < η_k||F(x_k)||. In this case it is important to set η_k such that η_k ||F(x_k)|| ~ e where e is the precision required for the nonlinear system. A safeguard which represents these ideas is: if η_k < 2ε then we set η_k = (0.8e)||F(x_k)||.

• Parameters choice for η_k:

  for the choice EW2 it was taken γ = 1 and α = 0.5(1 + ) and for the choice GLT it was taken ρ = 1.1.

• Stopping criterion:

  the process is finished successfully if ||F(x_k)|| < 10^-6 and k < 100.

• Restarts and the maximum number of iterations in GMRES(m):

  we fix the restarts at each m iterations, m = 30 or m = 50, allowing initially a maximum of 100 cycles (100m iterations). This maximum number is called here by maxit and it is adjusted during the process. This is the case, if the value of ||F|| increases, when F is computed at the solution obtained after an inner iteration, before doing the line search. Also, maxit is adjusted when the number of GMRES iterations have exceeded a certain value. Indeed, maxit is computed according to: if 1 < ||F(x_new)||/||F(x_old)|| < 100, then maxit is fixed as 50; if ||F(x_new)||/||F(x_old)|| > 100 or if the maximum number of GMRES iterations has been exceeded, then maxit is fixed as 30. We observe that this procedure is repeated just in two consecutive iterations. After that, the value of maxit is taken as 100 again.

• Strategy H:

  after each GMRES iteration, a possible stagnation is detected by the tests (3) and (4). Initially, we use 0.9 as a tolerance for the value of the cosine of θ_j or θ_j,1. If the hybrid scheme is triggered 5 times, we change this tolerance to 0.8. In the case of the hybrid process be triggered at the beginning of 10 cycles, we stop the test as in Section 4.

We present some numerical results obtained from the solution for a ray-tracing problem [13] and also from a set of boundary value problems.

5.2 A ray-tracing problem, [2]

We will present here a very simpliflied ray-tracing problem. First of all, let us represent the earth surface, as if it was two dimensional. Imagine that an acoustic wave transmiter S and a geophone receiver G, are located in two different points on the earth surface. Each acoustic wave will cross a certain number of layers, under the surface, and when coming back, it will be captured by the geophone. Such layers considered elastic, isotropycs and homogeneous, compose the structure of the undersurfaces of the earth.

The problem that we are interested in, consists on determining the trajectory of a ray, when crossing the earth undersurface [13].

The interfaces between two consecutive layers, are defined as functions of the horizontal coordinate x. These functions are continuous and smooth - we will not consider the intersections between two consecutive interfaces. The functions will be represented by z = f_i(x), i = 1, ... , m, where m corresponds to the number of interfaces and z denotes the vertical coordinate.

The number of times that a ray crosses the layer - downwards or upwards - in a layer between the interfaces f_i and f_i+1 will be denoted by (a_i), i = 1, ... , m - 1. Figure 5 is an example of this, where a₁ = 2 and a₂ = 1. The vector a = (a₁, ..., a_m-1) is called the signature of the ray. For uniqueness, we assume that all the downward crossings in a layer occur before those in the next layer, as in Figure 5.

Under these conditions, the rays must satisfy the Snell law, and so, they describe a straight line trajectory in each layer. Let the initial point S, the end point G on the earth surface, and the velocity of the ray in the j-th layer v_j, j - 1 , ... , m - 1, be given. The problem consists on finding the points X_k = (x_k, f_ik(x_k)), for k = 1 , 2 ,... ,n that satisfy the Snell law in each reflector i_k. We will take X₀ = S and X_n+1 = G. The dimension of the problem is given by

The tangent vector to the k-th interface at the point X_k, will be called τ_k, and is defined by τ_k = (1 , (x_k). The Snell law is given by

where ||.|| denotes the Euclidean norm. Thus, using equation (10) for k = 1 ,... ,n we obtain a nonlinear system of equations with n equations and n unknowns, given by

with Φ:

ⁿ→ ⁿ, Φ(x) = (Φ₁(x), Φ₂(x), ..., Φ_n(x))^T and the functions Φ_k: ⁿ→ (k = 1, ..., n) defined by

For the tests presented in Table 1, we considered a ray-tracing problem with two horizontal layers and signature a = (500,500). Each ordered pair represents the result of GMRES (left) and GMRESH (right). Column iterex represents the total number of external iterations that were performed; iterin, the number of inner iterations, that is, the number of iterations performed by GMRES(m) in the whole process; the column Feval represents the total number of function evaluations needed during the whole process; finally, CPU time represents the CPU time in seconds. The stopping criterion was ||Φ|| < 10^-5.

We observe that in this example the second safeguard in Strategy H was never triggered, while the first safeguard decreased considerably the number of inner iterations.

5.3 Boundary value problems

The general formulation of the boundary value problems solved in this work is finding u: Ω = [ 0, 1] × [ 0, 1] → , such that, for λ ,

The real valued function h(λ, u), the different values for the parameter l and the function f define the different problems tested. All the problems were discretized using central differences on a grid with L inner points in each axis. The discretized system obtained has L² equations and variables. In these problems we work with grids of L = 63 and L = 127 inner points. We now make a brief description of the particular problems that were solved:

• Bratu problem: the function h is given by h(λ, u) = λ exp(u), and the function f(s,t) is constructed so that u^*(s,t) = 10st(1 - s)(1 - t) is the exact solution for the problem. When λ < 0, the problem is considered relatively easy; not surprisingly, the hybrid strategy has never been triggered in our tests. The problem is more difficult for λ > 0, [9];

• a convection-diffusion problem: in this problem, the function h is given by h(λ, u) = λu(u_s + u_t), where u_s and u_t denote the partial derivatives of the function u with respect to s and t, and again the function f(s,t) is defined so that u^*(s,t) = 10st(1 - s)(1 - t) is the exact solution for the problem. This is a problem considered difficult to solve [9], in particular for values of λ greater than 50;

• a third problem: P3. This problem appears in the book of Briggs, Henson and McCormick [4], page 105. In this case, h(λ, u) is given by h(λ, u) = λue^u and the function

Table 2 shows the results of Newton-GMRES and Newton- GMRESH applied to the above problems, with: L = 63, n = 3969, x₀ = (0, 0, ..., 0)^T and λ = 100. This value for λ was chosen by the occurrence of stagnation. Each ordered pair should be read as in Table 1 and iterex, iterin, feval and CPU time have the same meaning as before.

For the Bratu problem, GMRESH shows a considerable decrease in the number of inner iterations, thus effectively mitigating the stagnation problem. The line search was not triggered and iterex is slightly reduced. In the convection-diffusion problem, the number of inner iterations was reduced with GMRES in most of the cases. For this problem the Jacobian matrices are always ill-conditioned near the solution. The line search was activated several times, indicating stagnation coupled with insufficient decrease. For all the problems considered, iterex was practically the same (or slighly better) using GMRESH. Problem P3 can be considered an easy problem with no ocurrence of stagnation.

We show in Table 3, the performance of the Algorithm 1 with GLT choice for the forcing term at the resolution of Bratu problem with λ = 100. Our purpose is to compare the number of inner iterations performed by GMRES and GMRESH at each outer iteration. We also show how many times the Strategy H is triggered. We observe the ocurrence of stagnation at 4-th outer iteration. After this iteration, the strategy H was triggered with a very good performance since the number of inner iterations was extremely reduced and the convergence was accelerated.

We chose the Bratu and convection-diffusion problems to exploit other parameters, such as the grid size L = 63 and L = 127. With these choices the dimension of the nonlinear systems solved were 3969 and 16129, respectively. We also implemented GMRES(m) with m = 30 and m = 50 for L = 63 and m = 50 and m = 80 for L = 127. Table 4 shows the results for these systems. Observe that even with a very large dimension there was no difficulty for the Algorithm 1 to solve the corresponding nonlinear systems, in both cases: GMRES(m) and GMRESH(m). As before, for the Bratu problem, GMRESH(m) shows a considerable decrease in the number of inner iterations, thus effectively mitigating the stagnation problem.

We also tested GMRESDR(30,5) at the resolution of the linear system generated at the first outer iteration of Algorithm 1. It was slightly better than GMRESH(30) in solving Bratu problem but much better when solving convection-diffusion problem which has a very ill conditioned Jacobian matrix.

5.4 The performance profile of the methods

The performance profile, proposed by Dolan and Moré [6], is a useful tool to compare a set of algorithms used for solving a set of problems. As comparison measures, we can use, for instance, the number of iterations performed, the number of function evaluations, the CPU elapsed time, etc.

In this subsection we analyze the performance profile of the AlgorithmsNewton-GMRES and Newton-GMRESH when applied to solve the boundaryvalue problems. A total of 18 problems were tested: Bratu, convection-diffusion and P3 with the following values for λ: = 10, 25, 30, 50, 75 and 100. We compare the performance of this Algorithms using for η_k only the choices EW2 and GLT, in order to get an understandable figure. We used a grid with L = 63 inner points and m = 30 because with these parameters the average behaviour of both Algorithms is well represented. In the legends of Figure 6, Algorithm Newton-GMRESH is indicated by an H, after the name of the choice.

The measures used for comparison were the total numbers of inner and outer iterations. Note that in Figure 6 a higher curve means better performance. With respect to both measures, it is evident that the two Newton-GMRESH (NGH) outperformed the two Newton-GMRES (NG) versions. In particular, the percentage of problems solved with a small number of inner (resp. outer) iterations was roughly 43-53% against 30% (resp. 80-95% against 60-70%). It can be observed that NG versions solved some problems with a number of inner iterations 3 and 4 times greater than the minimum required to solve that problem with NGH versions. Thus the algorithm NG had a worse performance than the new version NGH. The strategies worked successfully.

6 Conclusions

In this work we presented an inexact Newton-like Algorithm with a nonmonotone line search, in which it was introduced an strategy to prevent stagnation of the linear solver GMRES. This strategy showed as advantages: (i) the simplicity of implementation, since it does not requires interfering in the inner procedure of the linear solver GMRES; (ii) it can be monitored at each iteration; (iii) from the numerical results obtained, we can conclude that this strategy is efficient, either in the test for detecting the stagnation of the inner solver, or with respect to the safeguards triggered in this case. This conclusion can be seen in the solution of the ray-tracing problem, showed in Table 1, as well as in the solution of a bunch of boundary value problems whose performance profile is showed in Figure 6; (iv) GMRESH(m) is more efficient at the resolution of no ill conditioning stagnated systems. Maybe we can use GMRESH(m) together with GMRESDR(m,k) to obtain a more efficient algorithm.

Acknowledgements. The authors are indebted to an anonymous referee for the careful reading of this paper and also for the suggestions.

Received: 22/VI/07. Accepted: 28/II/08. thanks

#729/07.

• [1] E.G. Birgin, N. Krejić and J.M. Martínez, Globally convergent inexact quasi-Newtonmethods for solving nonlinear systems Numer. Algorithms, 32 (2003), 249-260.
• [2] N. Bleistein, Mathematical Methods for Wave Phenomena Academic Press, (1984).
• [3] C. Brezinski and M. Redivo-Zaglia, Hybrid procedures for solving linear systems Numer. Math., 67 (1994), 1-19.
• [4] W.L. Briggs, V.E. Henson and S.F. McCormick, A Multigrid Tutorial Second Edition, SIAM, Philadelphia (1987).
• [5] R.S. Dembo, S.C. Eisenstat and T. Steihaug, Inexact Newton methods SIAM J. Numer. Anal., 19 (1982), 401-408.
• [6] E.D. Dolan and J.J. Moré, Benchmarking optimization software with performance profile Math. Program. ser. A91 (2002), 201-213.
• [7] S.C. Eisenstat and H.F. Walker, Choosing the forcing terms in inexact Newton method SIAM J. Sci. Comput., 17 (1996), 16-32.
• [8] M. Embree, The tortoise and the hare restart GMRES SIAM Review, 45 (2003), 259-266.
• [9] M.A. Gomes-Ruggiero, D.N. Kozakevich and J.M. Martínez, A numerical study on large-scale nonlinear solvers Comput. Math. Appl., 32 (1996), 1-13.
• [10] M.A. Gomes-Ruggiero, V.L. Rocha Lopes and J.V. Toledo-Benavides, A new choice for the forcing term and a global convergent inexact Newton method Ann. Oper. Res., 157 (2008), 193-205.
• [11] W. Joubert, On the convergence behavior of the restarted GMRES algorithm for solving nonsymmetric linear systems Numer. Linear Algebra Appl., 1 (1994), 427-447.
• [12] D-H. Li and M. Fukushima, Derivative-free line search and global convergence of Broyden-like method for nonlinear equations Optim. Methods Softw., 13 (2000), 181-201.
• [13] H.B. Keller and D.J. Perozzi, Fast seismic ray-tracing SIAM J. Appl. Math., 43 (1983), 981-992.
• [14] C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations SIAM, Philadelphia (1995).
• [15] L.F. Mendonça and V.L. Rocha Lopes, Uma análise comparativa do desempenho de métodos quase-Newton na resoluçăo de problemas em sísmica Tema-SBMAC, 5 (2004), 107-116.
• [16] R.B. Morgan, A restarted GMRES method augmented with eigenvectors SIAM J. Matrix Anal. Appl., 16 (1995), 1154-1171.
• [17] R.B. Morgan, GMRES with deflated restarting SIAM J. Sci. Comput., 24 (2002), 20-37.
• [18] N.M. Nachtigal, L. Reichel and L.N. Trefethen, A hybrid GMRES algorithm for nonsymmetric linear systems SIAM J. Matrix Anal. Appl., 13 (1992), 796-825.
• [19] C.C. Paige and M.A. Saunders, Algorithm 583 LSQR: Sparse linear equations and least squares problems ACM Trans. Math. Software, 8 (1982), 195-209.
• [20] C.C. Paige and M.A. Saunders, LSQR: An algorithm for sparse linear equations and sparse least squares ACM Trans. Math. Software, 8 (1982), 43-71.
• [21] M. Parks, E. Sturler, G. Mackey, D.D. Johnson and S. Maiti, Recycling Krylov subspaces for sequences of linear systems SIAM J. Sci. Comput., 28 (2006), 1651-1674.
• [22] M. Pernice and H. Walker, NITSOL: a Newton iterative solver for nonlinear systems SIAM J. Sci. Comput., 19 (1998), 302-318.
• [23] Y. Saad, A flexible inner-outer preconditioned GMRES algorithm SIAM J. Sci. Comput., 14 (1993), 461-469.
• [24] Y. Saad, Iterative Methods for Sparse Linear Sistems Second Edition, SIAM, Philadelphia (2003).
• [25] Y. Saad and M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems SIAM J. Sci. Comput., 7 (1986), 856-869.
• [26] W. Schonauer, Scientific Computing on Vector Computers North-Holland, Amsterdam (1987).
• [27] V. Simoncini, A new variant of restart GMRES(m) Numer. Linear Algebra Appl., 6 (1999), 61-77.
• [28] V. Simoncini, On the convergence of restarted Krylov subspace methods SIAM J. Matrix Anal. Appl., 22 (2000), 430-452.
• [29] V. Simoncini and D.B. Szyld, Theory of inexact Krylov subspace methods and applications to scientific computing SIAM J. Sci. Comput., 25 (2003), 454-477.
• [30] G. Starke and R.S. Varga, A hybrid Arnoldi-Faber iterative method for nonsymmetric systems of linear equations Numer. Math., 64 (1993), 213-240.
• [31] E. de Sturler, Truncation strategies for optimal Krylov subspace methods SIAM J. Numer. Anal., 36 (1999), 864-889.
• [32] J.V. Toledo-Benavides, A globally convergent Newton-GMRES method with a new choice for the forcing term and some strategies to improve GMRES(m) PhD Thesis, T/UNICAMP T575m, (2005).
• [33] H.A. Van der Vorst and C. Vuik, GMRESR: a family of nested GMRES methods Num. Linear Algebra Appl., 1 (1994), 369-386.
• [34] H.F. Walker and L. Zhou, Residual smoothing techniques for iterative methods SIAM J. Sci. Comput., 15 (1994), 297-312.
• [35] I. Zavorin, Spectral factorization of the Krylov matrix and convergence of GMRES UMIACS Technical Report CS-TR-4309, University of Maryland, College Park (2002).
• [36] I. Zavorin, D.P. O'Leary and H. Elman, Complete stagnation of GMRES Lin. Algebra Appl., 367 (2003), 165-183.

*

Supported by FAPESP, CNPq, PRONEX-Optimization.

Publication Dates

• Publication in this collection
  21 July 2008
• Date of issue
  2008

History

• Accepted
  28 Feb 2008
• Received
  22 June 2007