1. Introduction. The purpose of this paper is to analyze the convergence of the Energy Direct Inversion on the Iterative Subspace (EDIIS) algorithm [21]. EDIIS is a modification of Anderson acceleration [1] or the Direct Inversion on the Iterative Subspace (DIIS) method [34,37,22,21]. EDIIS relaxes the need for a sufficiently accurate initial iterate. EDIIS is the default solver for the self-consistent field (SCF) iteration in the widely used Gaussian [12] quantum chemistry code. We prove convergence from any starting point in a convex set in which the fixed point map is a contraction and then analyze local convergence. Our local convergence is an improvement of the result in [41] and applies to both EDIIS and Anderson acceleration.

We will begin this introductory section with a review of Anderson acceleration and some of the recent convergence results. We will then describe the EDIIS algorithm. In section 2 we prove our convergence results. Finally, in section 3 we report on a computation which both illustrates the theory and, as is also done in [21], shows how the convergence speeds for EDIIS and Anderson acceleration, while identical in theory, can differ significantly in practice.

Our notational convention is that vectors and vector-valued functions in R N are in bold. Scalars and elements of infinite dimensional spaces (e.g., integral operators and the functions acted upon by those operators) are in the usual italic math font.

Anderson acceleration [1] is an iterative method for fixed point problems of the form (1.1) u = G(u),

where u ∈ R N and G : R N → R N . The method was designed to accelerate Picard or fixed point iteration, i.e.,

(1.2) u k+1 = G(u k ).

Anderson acceleration was originally designed for integral equations and has been widely used in electronic structure computations (see [9] and many references since then) and is now very common in that field. Anderson acceleration is essentially the same as Pulay mixing [33,32], DIIS [34,37,22,21], and nonlinear GMRES [25,30,45,4]. Other applications include nuclear reactor design [42,16], stiff dynamics [13], hydrology [24], and fluid-structure interaction [10,15,23], where the method is called interface quasi-Newton. The analysis of Anderson acceleration is far from complete. In this paper we assume, as do all theoretical results about this algorithm, that the map G is a contraction. In practice, however, Anderson acceleration does very well for problems in which G is either definitely not a contraction [41] or not provably a contraction. The results here do not explain those cases.

Anderson acceleration was designed for a problem where Newton's method is not practical because obtaining approximate Jacobians or Jacobian-vector products is too costly. One should expect that Newton's method would perform better when derivative information can be had at a reasonable cost, and we have certainly found that to be the case in our own recent work [16]. Anderson iteration maintains a history of residuals

of size at most m + 1, where the depth m is an algorithmic parameter. When m is important, we will call the iteration Anderson(m). Anderson(0) is Picard iteration by definition.

The formal description in Algorithm 1.1 is most convenient for analysis and exposition, but not for implementation. We refer the reader to [44,41,7,38,43,39] for examples of efficient implementations.

The iteration uses the most recent m+1 residuals F(u j ) for k -m k ≤ j ≤ k where m k ≤ min(k, m). The key step in the iteration is solving the optimization problem

for the coefficients {α k j }. Any vector norm can be used in the optimization problem with no change in the convergence theory [41]. In particular, the optimization problem for the coefficients in either the 1 or the ∞ norm can be formulated as a linear programming problem [8]. The optimization problem is easier to solve if one uses the 2 norm, and that is standard practice. In this case optimization problem for the coefficients can be expressed as a linear least squares problem and solved very inexpensively. One way to do this is to solve the linear least squares problem (1.4) Minimize

The choice of m k is, in the original form, simply min(m, k). One can adapt m k as the iteration progresses to, for example, enforce well-conditioning of the linear least squares problem (1.4) [44,39].

One can also [11,35,34,44,31] show that Anderson acceleration is related to multisecant quasi-Newton methods or, in the case of linear problems, GMRES. None of these results lead to a convergence proof, even in the linear case, unless the available storage is large enough to allow GMRES to take a number of iterations equal to the dimension of the problem. The recent work of one of the authors and his students [41,40,39] contains the first convergence theory for Anderson acceleration as it is applied in practice.

1.1. Convergence theory. Theorem 1.1 is one of the convergence results from [41]. That paper also has results for several special cases. We assume that G is a contraction with contractivity constant c ∈ (0, 1) in a closed set

for all u, v ∈ D. The contraction mapping theorem implies that G has a unique fixed point u * ∈ D. As is standard, we let e = uu * and make the assumption from [41] on the smoothness of G and the Anderson iteration coefficients. The convergence of the Picard iteration for a contraction map is q-linear [19] with q-factor c, i.e., e k ≤ c e k-1 .

We will show in this paper that Anderson acceleration is r-linear with r-factor c, which means

for some r > 0.

There is M α such that for all k ≥ 0

The differentiability assumption is needed in the analysis, but not in the formulation or implementation of the algorithm. Our convergence result in section 2.2 relaxes the assumption to continuous differentiability. Theorem 1.1 ([41]). Let Assumption 1.1 hold, and let c < r < 1. Then if u 0 is sufficiently close to u * , the Anderson iteration converges to u * . In fact, for all k ≥ 0,

The interpretation of this result is that if the initial data are sufficiently good, then the r-factor for Anderson iteration is no worse than the q-factor of Picard iteration as predicted by the contractivity constant c. While r-linear convergence is weaker than q-linear convergence, Anderson acceleration is often faster than Picard iteration in practice. The requirement that the initial iterate be near the solution is also meaningful in practice [36,46,47] and motivated the EDIIS algorithm [21], which is the subject of this paper.

Both Picard iteration and Anderson acceleration can perform better than the prediction (see section 3). In practice, Anderson acceleration is often significantly better than Picard iteration, but there is no theory that explains this under practical (i.e., very limited storage) operating conditions.

1.2. The EDIIS algorithm. Anderson acceleration performs poorly for some applications. One example is electronic structure computations for metallic systems where the HOMO-LUMO gap is small and a good initial iterate is difficult to obtain. In this case both Picard iteration and Anderson acceleration perform poorly [21]. In such cases one can sometimes use a small mixing parameter to ensure convergence, especially when the initial iterate is poor. However, a small mixing parameter may degrade the performance of the iteration-especially when near the solution. The role of the damping parameter β in Picard iteration is simple damping:

then [40] one obtains

which is how damping is done in Anderson acceleration [1].

One attempt to solve these problems for small systems is the EDIIS algorithm from [21]. In [21] the authors also formulated the fixed point problem to directly minimize energy (hence the name of the method), but that does not affect the convergence analysis in this paper.

EDIIS differs from Anderson acceleration by imposing a nonnegativity constraint on the coefficients. So, the optimization problem becomes (1.7) Minimize

In [21] the authors present an example where EDIIS does well and both Picard and Anderson acceleration fail and another example where Anderson acceleration is successful and EDIIS, while converging, does not perform as well. We present another such example in section 3. One reason EDIIS might perform worse than Anderson acceleration could be that the optimization problem (1.7) for EDIIS has a more restricted feasible set and therefore a larger optimal value.

2. Convergence results. Our global convergence is Theorem 2.1. The proof does not require differentiability, but the convergence speed estimate is very pessimistic with an r-factor of c 1/(m+1) . We follow the global theorem with a local theorem that shows how the convergence behavior becomes locally r-linear with r-factor c, improving on the local results in [41]. In fact,

Proof. The proof does not use the optimality properties of the coefficients and only requires that the iteration {u k } have the form

where m k ≤ m, α k j ≥ 0, and m k j=0 α k j = 1. We induct on k. Clearly (2.1) holds for both m k = 0, by definition, and k = 1, m k = 0 because the iteration in that case is a single Picard iteration (i.e., one step of Anderson(0)). Assume that the result holds for k ≤ K. Then (2.2) and

Theorem 2.1 implies that for any δ > 0 there is K such that all iterations {u k } k≥K are in the set

Hence, starting an Anderson acceleration iteration after sufficiently many EDIIS iterations will result in local convergence at the rate predicted by Theorem 1.1, which is better than (2.1) since r can be arbitrarily near c and does not depend on m. However, it is not clear how to decide when to restart. The main result in section 2. [41] by both weakening the assumptions and improving the r-factor.

We will assume that an iteration begins with a history that lies in B(δ) for δ sufficiently small. This history could be either from the EDIIS iteration or from the Anderson acceleration iteration itself. Hence the assumption covers not only EDIIS but also allows us to improve the convergence theory from [41]. We will show that the residuals converge r-linearly to zero with an r-factor of c. Formally our assumption is as follows.

Assumption 2.1. G is a continuously differentiable contraction on D ⊂ R N with contractivity constant c and u * is the unique fixed point of G in D.

The iteration begins with

(2.4)

Finally, there is ĉ ∈ (c, 1) so that

Theorem 2.1 implies that Assumption 2.1 will hold after sufficiently many EDIIS iterations. In the theorem there is no history if m = 0 and in that case the iteration is a Picard iteration. While we are motivated by a local iteration from the EDIIS algorithm, the local theory does not require that the coefficients be nonnegative.

Assumption 2.1 weakens the ones in [41] in two ways. The first is that we no longer assume that G is Lipschitz continuously differentiable. The second is that we do not assume that the coefficients {α k j } come from any particular optimization problem-only that the linear combination of residuals has norm no larger than that of the most recent residual.

The idea of the analysis is that as the iteration converges, the upper bound for the r-factor will approach c and therefore the r-factor is no larger than c. In the case where there is no history, this fact was implicit in the results from [41]. Adding the history makes the bookkeeping more difficult, and the proof of Theorem 2.2 must account for that.

Theorem 2.2. Let Assumption 2.1 hold. Assume that there is M α such that

for all k ≥ 0. Then if δ is sufficiently small, the iteration given by (2.3) and (2.4) converges to u * and

Proof. Let 0 < < ĉ -c. We will show that for e 0 sufficiently small,

This will complete the proof since is arbitrary and we can restart the proof once we have m vectors in the history which are near enough to u * to reduce further. We induct on k. Define L = (c/ĉ) m . We will show that (2.9)

for all k. Our assumption on the history that F(u l ) ≤ ĉl F(u 0 ) implies that (2.9) holds for 0 ≤ k ≤ m. Now suppose that (2.9) holds for all 0 ≤ l ≤ k with k ≥ m. We will establish the bound for k + 1. The analysis has three steps. We first set δ small enough for the iteration to remain in D. We then derive an estimate for F(u k+1 ) and finally use that estimate to continue the induction.

Step 1, initialization of δ. Since G is continuous in D, there is a nondecreasing function ρ ∈ C[0, ∞) with ρ(0) = 0 so that (2.10)

for all u ∈ D. This implies that for all u ∈ D,

where ∆(e) ≤ ρ( e ) e .

Reduce δ if necessary so that (2.12)

This implies that (2.13)

for sufficiently small δ because (2.14)

Step 2, estimation of F(u k+1 ). We may write, for k ≥ m -1,

We will estimate the two parts of the sum

and

separately.

Using only contractivity of G and (2.4), we have (2.15)

We now estimate B k . Using (2.12), we have for all u ∈ D with

The final stage in the proof is to show that, reducing δ if needed, (2.17)

Recall that

We use (2.11) to obtain

Hence

We will estimate terms separately. First (2.18)

Finally, using (2.14) and (2.16), (2. 19)

Adding the two estimates (2.18) and (2.19) leads to (2.17).

Step 3, continuation of the induction. Combining (2.15), (2.17), (2.9), and the induction hypotheses, we have (2.20)

This implies (2.8), which in turn implies (2.7) because is arbitrary.

Theorem 2.2 and nonsingularity of F (u * ) also imply r-linear convergence of the errors with r-factor c. This extends and sharpens (1.6). Proof. We will use Lemma 5.2.1 from [19], which states that if u is sufficiently near u * and F (u * ) is nonsingular, then

which is (2.21).

3. Numerical example. We will use an example [41] to show how the actual performance of EDIIS and Anderson acceleration can differ, even though the theoretical limiting convergence estimates are identical. Another point of this section is that the solver for the optimization problem can significantly affect the results.

The results in [21] also illustrate this point. Our example is simple enough to directly compare the iteration histories for Picard iteration, EDIIS, and Anderson with the worst-case prediction given by the contractivity constant. We find that when Anderson acceleration performs well, as it does in this example, EDIIS offers no advantage. Moreover, the additional constraint on the optimization problem for the coefficients leads to slower convergence, exactly matching Picard iteration in this case.

The optimization problem for EDIIS requires more care than the linear least squares problem one must solve for Anderson acceleration. The reason for this is that one cannot simply use a QR factorization to solve (1.4). Instead one must apply a more sophisticated iterative solver. The approach of [21] is a direct examination of the boundary of the feasible simplex, which is not practical for a depth much greater than m = 3. Since m is small in practice, expressing the optimization problem as a boundconstrained quadratic program is an efficient alternative. References [27,26] survey the literature on this topic. For example, a bound-constrained quadratic programming code such as the MINQ [29] code is a reasonable choice. However, this approach squares the condition number and can (and did in our testing) result in a singular or nearly singular KKT system and failure of the optimization code's internal linear solvers. The method of [6], while still squaring the condition number, is more robust and terminated without error for this example. The classic method from [14] uses an active set method and the QR factorization to avoid using the normal equations. The approach in [14] performed better in the example here, where the least squares coefficient matrix for the optimization problem is ill-conditioned [41].

The example is the midpoint rule discretization of the Chandrasekhar H-equation [5,3].

We seek a solution H * ∈ C[0, 1]. When the parameter ω is important we will write H * as a function H * (µ, ω) of both µ and ω.

The integral equation and its midpoint discretization share the properties that the fixed point map

is a contraction for 0 ≤ ω < 1, but not for ω = 1. The Fréchet derivative (and the Jacobian for the discrete case) is singular at the solution for ω = 1, which is a simple fold singularity [17,28].

In this section we will compare the performance of Picard iteration, Anderson acceleration, and EDIIS for the case ω = .5 on an N = 100 point mesh. We terminated the iteration when the residual had decreased by a factor of 10 -12 .

One interesting result from [41] is that Anderson(m) is more efficient than Newton's method for this example, even in the singular case. In the context of this paper it is important to note that Picard iteration converges faster than one would expect from estimating the contractivity parameter by the spectral radius of the Fréchet derivative of G at the solution, which is a lower bound for the operator norm of G. From [41] ρ(G (H

However [18,20,2], the solution is analytic in ω and Picard iteration exploits that property to obtain q-linear convergence with q-factor ≤ ρ(G (H * )) and much less for small ω. In fact, if

is the Taylor expansion of H * in ω, then the coefficient functions {a m (µ)} are nonnegative for 0 ≤ µ ≤ 1. Moreover, the series converges for ω = 1. Hence, if H k is the kth Picard iteration and H 0 ≡ 0, then for all k ≥ 0 and ω, µ ∈ [0, 1],

All of the above statements about the singularity at ω = 1, the spectral radius of the Fréchet derivative, and the performance of Picard iteration apply to the discrete problem

In (3.2) µ i = (i -1/2)/N is the ith quadrature node for the N point composite midpoint rule, the vector h * is the solution of the discrete problem h = G(h), G(h * ) i is the ith component of G(h * ), and the ith component of h * is h * i ≈ H * (µ i ). As noted above, the optimization problem (1.7) for EDIIS is harder than the one for Anderson acceleration, and the choice of solver can be important. We compare the method of [14], as implemented in the MATLAB lsqlin code with the "active-set" option, with the method from [6], as implemented with the "interior-point" option in lsqlin. The method of [6] uses the normal equations and did exhibit problems with ill-conditioning. The computations were done on an Apple Macintosh running Table 3.1 compares ρ(G (H * )) to the r-factors of the residuals for Anderson acceleration, Picard iteration, and EDIIS. We estimate the r-factors by

, where the final iteration upon termination is h k . Note that, as discussed above, the q-factor for Picard iteration is smaller than the spectral radius. Anderson acceleration also does better than the theory predicts and, in fact, is more efficient than Newton-GMRES [41].

EDIIS-I is the only one of the methods which is sensitive to the ill-conditioning of the optimization problem. We examined this sensitivity by solving the problem twice, once with no limit on the condition number and again by reducing m k if necessary to limit the condition number to 10 5 . This has no effect on EDIIS-A and slightly slows Anderson acceleration down. We show the residual histories in Figure 3.1, where one can clearly see the effect of limiting the condition number. As reported in [41], the optimization problem becomes more ill-conditioned as the iteration progresses. The figures show that the convergence of EDIIS-I degrades at the 6th iteration, but to a lesser degree when the condition number is limited. Note that the estimated r-factor seems to stabilize near the end and is, in the condition number limited case, back to Picard iteration for the final three iterations, albeit from a worse starting point.

The EDIIS algorithm was designed to improve the global convergence properties of the DIIS algorithm, which is also known as Anderson acceleration. We prove global convergence of the iteration and prove a local convergence result that applies to both EDIIS and Anderson acceleration and improves the results in [41]. We observe, as did the inventors of the method [21], that the unmodified version of Anderson acceleration can have better local convergence in practice.