Broyden Updating, the Good and the Bad!

Andreas Griewank

Abstract.

2010 Mathematics Subject Classification: 65H10, 49M99, 65F30 Keywords and Phrases: Quasi-Newton, secant condition, least change, bounded deterioration, superlinear convergence

So far so good! We had an updating procedure (the ’full’ secant method) that seemed to work provided that certain conditions of linear independence were satisfied, but the problem was that it did not work very well. In fact it proved to be quite numerically unstable. Charles Broyden in On the discovery of the ‘good Broyden’ method [6].

The idea of secant updating

As Joanna Maria Papakonstantinou recounted in her comprehensive historical survey [29], regula falsi and other variants of the secant method for solving one equation in one variable go back to the Babylonian and Egyptian civilizations nearly 4000 years ago. They may be viewed just as a poor man’s version of what is now known as Newton’s method, though we should also credit Al Tusi [20]. During antiquity the very concept of derivatives was in all likelihood unknown, and in modern times the evaluation (and in the multivariate case also factorization) of Jacobian matrices is frequently considered too tedious and computationally expensive.

The latter difficulty was certainly the concern of Charles Broyden in the sixties, when he tried to solve nonlinear systems that arose from the discretization of nonlinear reactor models for the English Electric Company in Leicester [6]. Now we know that, due to diffusion, the resulting system of ODEs must have been rather stiff, but that property was only identified and analyzed a few years later by Dahlquist. Nevertheless, Broyden and his colleagues already used some implicit time integration schemes, which required solving sequences of slightly perturbed nonlinear algebraic systems F (x) = 0 for F : ℝ^{n} ↦ ℝ^{n}.

Broyden noted that one could avoid the effort of repeatedly evaluating and factoring the system Jacobian by exploiting secant information, i.e., function value differences

y_{i} ≡ F_{i} − F_{i−1} with F_{j} ≡ F (x_{j}) for j ≤ i

Here, x_{i} ∈ ℝ^{n} denotes the current iterate and x_{j}, for j < i, disctinct points at which F has been evaluated previously. With s_{i} ≡ x_{i} − x_{i−1} the new approximation B_{i} to the Jacobian F′(x_{i}) ∈ ℝ^{n×n}

B_{i}s_{i} = y_{i} = F′(x_{i})s_{i} + o(∥s_{i}∥)  (1)

The first order Taylor expansion on the right is valid if F has a Jacobian F′(x) ∈ ℝ^{n×n} that varies continuously in x. We will tacitly make this assumption throughout so that F ∈ C^{1}(D) on some open convex domain D ⊂ ℝ^{n} containing all evaluation points of interest.

In the univariate case of n = 1, one can divide by s_{i} to obtain B_{i} = y_{i}/s_{i} ≈ F′(x_{i}) uniquely. In the multivariate case, the secant condition merely imposes n conditions on the n^{2} degrees of freedom in the new approximating Jacobian B_{i}. A natural idea is to remove the indeterminacy by simultaneously imposing earlier secant conditions B_{i}s_{j} = y_{j}, for j = i − n + 1 … i. The resulting matrix equation for B_{i} has a unique solution provided the n + 1 points x_{i−n+j}, for j = 0 … n, are in general position, i.e., do not belong to a proper affine subspace of ℝ^{n}. Theoretically, that happens with probability 1, but in practice the step vectors s_{j}, for j = i − n + 1 … i, are quite likely to be nearly linearly dependent, which leads to the observation of instability by Broyden cited above.

Rather than recomputing B_{i} from scratch, Broyden reasoned that the previous approximation B_{i−1} should be updated such that the current secant condition is satisfied, but B_{i}v = B_{i−1}v in all directions v ∈ ℝ^{n} orthogonal to s_{i}. As he found out ‘after a little bit of scratching around’, these conditions have the unique solution [2]

B_{i} = B_{i−1} + r_{i}s_{i}^{⊤}/ s_{i}^{⊤}s_{i}, with r_{i} ≡ y_{i} − B_{i−1}s_{i}  (2)

Here the outer product C_{i} ≡ r_{i}s_{i}^{⊤}/s_{i}^{⊤}s_{i} of the column vector r_{i} and the row vector s_{i}^{⊤} represent a rank one matrix. This formula became known as the good Broyden update, because it seemed to yield better numerical performance than the so-called bad formula (6) discussed below. For a recent review of quasi-Newton methods see the survey by J. M. Martinez [25].

Broyden stated that the fact that C_{i} = B_{i} − B_{i−1} turned out to be of rank one was pure serendipity. Even though he claimed ’When I was at University they did not teach matrices to physicists’, he realized right away that the low rank property could be used to reduce the linear algebra effort for computing the next quasi-Newton step

s_{i+1} = −B_{i}^{−1}F_{i}

to O(n^{2}). That compares very favourably with the n^{3}/3 arithmetic operations needed for a dense LU factorization of the new Jacobian F′(x_{i}) to compute the Newton step −F′(x_{i})^{−1}F_{i}. If the previous step is given by s_{i} = −B_{i−1}^{−1}F_{i−1}, one can easily check that the secant error vector r_{i} defined in (2) is identical to the new residual, i.e., r_{i} = F_{i}, which we will use below.

Tacking on a sequence of rank one corrections to an initial guess B_{0}, and reducing the linear algebra effort in the process looks more like an engineering trick than an algorithmic device of mathematical interest. Yet after a few years and in close collaboration with his coauthors John Dennis and Jorge Moré, a beautiful theory of superlinear convergence theory emerged [7], which was later built upon by other researchers and extended to many update formulas. For a much larger class of methods named after Charles Broyden and his coauthors Abbaffy and Spedicato, see [1].

Least change interpretation

John Dennis credits Jorge Moré with a short argument showing that the good Broyden formula is a least change update. Specifically, if we endow the real space of n × n matrices A with the inner product

⟨A, B⟩ ≡ T r(A^{⊤}B) = T r(B^{⊤}A)

then the corresponding norm

∥A∥_{F} ≡ \sqrt{⟨A, A⟩} ≥ ∥A∥  (3)

is exactly the one introduced by Frobenius. It is bounded below by the consistent matrix norm ∥A∥ induced by the Euclidean vector norm ∥v∥ on ℝ^{n}. The affine variety

[y_{i}/s_{i}] ≡ {B ∈ ℝ^{n×n} : Bs_{i} = y_{i}}

has the n(n − 1) dimensional tangent space [0/s_{i}] and the n dimensional orthogonal complement

[0/s_{i}]^{⊥} ≡ {vs_{i}^{⊤} ∈ ℝ^{n×n} : v ∈ ℝ^{n}}

Hence, the smallest correction of B_{i−1} to obtain an element of [y_{i}/s_{i}] is given by the correction

C_{i} = r_{i}s_{i}^{⊤}/s_{i}^{⊤}s_{i} ∈ [r_{i}/s_{i}] ∩ [0/s_{i}]^{⊥}

For formal consistency we will set C_{i} = 0 if s_{i} = 0 = y_{i}, which may happen for all i ≥ j if we have finite termination, i.e., reach an iterate x_{j} with F_{j} = 0.

The geometry is displayed below and yields for any other element A_{i} ∈ [y_{i}/s_{i}] by Pythagoras

∥B_{i−1} − A_{i}∥_{F}^{2} − ∥B_{i} − A_{i}∥_{F}^{2} = ∥C_{i}∥_{F}^{2}

In particular, we have the nondeterioration property

∥B_{i} − A_{i}∥_{F} ≤ ∥B_{i−1} − A_{i}∥_{F}

This to hold for all A_{i} ∈ [y_{i}/s_{i}] is in fact equivalent to the least change property of the update. Broyden stated this property apparently for the first time in his survey paper [4], which he rarely cited afterwards. Moreover, nondeterioration can be equivalently stated in the operator norm as

∥B_{i} − A_{i}∥ ≤ ∥B_{i−1} − A_{i}∥  (4)

which makes sense even on an infinite dimensional Hilbert space where ∥ · ∥_{F} is undefined.

Sequential properties in the affine case

So far we have described the single least change update C_{i} = r_{i}s_{i}^{⊤}/s_{i}^{⊤}s_{i}, but the key question is of course how a sequence of them compound with each other. One can easily check that B_{i+1} = B_{i} + C_{i+1} = B_{i−1} + C_{i} + C_{i+1} satisfies the previous secant condition B_{i+1}s_{i} = y_{i} only if s_{i} and s_{i+1} are orthogonal so that C_{i+1}s_{i} = 0. In fact, exactly satisfying all n previous secant conditions is not even desirable, because that would lead back to the classical multivariate secant method, which was found to be rather unstable by Broyden and others. However, successive updates do not completely undo each other and thus eventually lead to good predictions B_{i−1}s_{i} ≈ y_{i}.

Now we will briskly walk through the principal arguments for the case when F is affine on a finite dimensional Euclidean space. Later we will discuss whether and how the resulting relations extend to nonlinear systems and infinite dimensional Hilbert spaces. Suppose for a moment that our equation is in fact affine so that

F (x) = A x + b with A ∈ ℝ^{n×n} and b ∈ ℝ^{n}.

Then the secant conditions over all possible steps s_{i} = −B_{i−1}^{−1}F_{i−1} are satisfied by the exact Jacobian A ∈ [y_{i}/s_{i}] since y_{i} = A_{i}s_{i} by definition of F. Moreover, let us assume that A and all matrices B with ∥B −A∥ ≤ ∥B_{0} −A∥ have inverses with a uniform bound ∥B^{−1}∥ ≤ γ. This holds by the Banach Perturbation Lemma [27] for all B_{0} that are sufficiently close to a nonsingular A.

Then we can conclude, as Broyden did in [3], that all B_{i} are nonsingular and, consequently, all steps s_{i} = −B_{i−1}^{−1}F_{i−1} are well defined and bounded by ∥s_{i}∥ ≤ γ∥F_{i−1}∥. Repeatedly applying Pythagoras’ identity we obtain for any i the telescoping result that

∑_{j=1}^{i} ∥C_{j}∥_{F}^{2} = ∥B_{0} − A∥_{F}^{2} − ∥B_{i} − A∥_{F}^{2} ≤ ∥B_{0} − A∥_{F}^{2} .

Hence, we derive from C_{j}s_{j} = r_{j} and the fact that the Frobenius norm is stronger than the operator norm that

lim_{j} ∥C_{j}∥_{F} → 0 and lim_{j} ∥r_{j}∥/∥s_{j}∥ ≤ lim_{j} ∥C_{j}∥ = 0.  (5)

Whereas these limits remain valid in the nonlinear case considered below, they hold in a trivial way in the affine case considered so far. This follows from the amazing result of Burmeister and Gay [12] who proved that Broyden’s good method reaches the roots of affine equations exactly in at most 2 n steps. The proof appears a little like an algebraic fluke and there is nothing monotonic about the approach to the solution. Moreover, the restriction that the ball with radius ∥B_{0} − A∥ contains no singular matrix can be removed by some special updating steps or line-searches as, for example, suggested in [26], [17], and [23], also for the nonlinear case.

The glory: Q-superlinear convergence

The property ∥r_{j}∥/∥s_{j}∥ → 0 was introduced in [8] and is now generally known as the Dennis and Moré characterization of Q-superlinear convergence. The reason is that it implies, with our bound on the stepsize, that ∥r_{j}∥/∥F_{j−1}∥ ≤ γ^{−1}∥r_{j}∥/∥s_{j}∥ → 0 and thus

\frac{∥F_{i+1}∥}{∥F_{i}∥} → 0 ⟺ \frac{∥x_{i+1} − x_{∗}∥}{∥x_{i} − x_{∗}∥} → 0

The equivalence holds due to the assumed nonsingularity of A so that, in any pair of norms, the residual size ∥F (x)∥ is bounded by a multiple of the distance Charles Broyden and his fellow quasi-Newton musketeers, J. Dennis and J. Moré

∥x − x_{∗}∥ and vice versa. Correspondingly, the central concept of Q-superlinear convergence is completely invariant with respect to the choice of norms, a highly desirable property that is not shared by the weaker property of Q-linear convergence, where the ratio of successive residual norms ∥F (x_{j})∥ or solution distances ∥x_{i} − x_{∗}∥ is merely bounded away from 1.

Under certain initial assumptions Q-superlinear convergence is also achieved in the nonlinear case, and under a compactness condition even in infinite dimensional space. All this without any exact derivative information or condition that the sequence of steps be in some sense linearly independent.

Originally, it was widely believed that to ensure superlinear convergence one had to establish the consistency condition that the B_{i} converge to the true Jacobian F′(x_{∗}). In fact, these matrices need not converge at all, but, theoretically, may wander around F′(x_{∗}) in a spiral, with the correction norms ∥C_{j}∥ square summable but not summable. This means that the predicted increments B_{i−1}s_{i}/∥s_{i}∥ in the normalized directions s_{i}/∥s_{i}∥ cannot keep being substantially different from the actual increments y_{i}/∥s_{i}∥ because the s_{i}/∥s_{i}∥ belong to the unit sphere, which is compact in finite dimensions.

The seemingly counterintuitive nature of the superlinear convergence proof caused some consternation in the refereeing process for the seminal paper by Broyden, Dennis and Moré [7]. It eventually appeared in the IMA Journal of Applied Mathematics under the editorship of Mike Powell. Broyden had analyzed the affine case, John Dennis contributed the concept of bounded deterioration on nonlinear problems and Jorge Moré contributed the least change characterization w.r.t. the Frobenius norm leading to the proof of superlinear convergence. All this is not just for good Broyden, but for a large variety of unsymmetric and symmetric updates like BFGS, where the Frobenius norms must be weighted, which somewhat localizes and complicates the analysis.

More specifically, suppose one starts at x_{0} in the vicinity of a root x_{∗} ∈ F^{−1}(0) near which the Jacobian is nonsingular and Lipschitz continuous. Then the nondeterioration condition (4) becomes a bounded deterioration condition with A_{i} replaced by F′(x_{∗}) and a multiplicative factor 1 + O(∥x_{i} − x_{∗}∥) as well as an additive term O(∥x_{i} − x_{∗}∥) on the right-hand side. From that one can derive Q-linear convergence provided B_{0} is close enough to F′(x_{∗}), which, in turn, implies Q-superlinear convergence by the perturbed telescoping argument. More generally, we have the chain of implications

Bounded deterioration

⟹ Linear Convergence ⟹ Q-superlinear Convergence.

Actually, R-linear convergence is enough for the second implication. This modularization of the analysis is a very strong point of the Broyden-Dennis-Moré framework [7] and has allowed many other researchers to communicate and contribute in an economical fashion.

Bad Broyden by inverse least change

The BDM mechanism also applies to so-called inverse updates, especially Broyden’s second unsymmetric formula. It can be derived by applying the least change criterion to the approximating inverse Jacobian

H_{i} = B_{i}^{−1} with H_{i} y_{i} = s_{i}

The equation on the right is called the inverse secant condition, which must be satisfied by H_{i} if B_{i} = H_{i}^{−1} is to satisfy the direct secant condition (1). After exchanging s_{i} and y_{i} and applying the good Broyden formula to H_{i} one obtains the inverse update on the left, which corresponds to the direct update of B_{i} on the right

H_{i} = H_{i−1} + \frac{(s_{i} − H_{i−1}y_{i})y_{i}^{⊤}}{y_{i}^{⊤}y_{i}} ⟺ B_{i} = B_{i−1} + \frac{r_{i}y_{i}^{⊤}}{y_{i}^{⊤}s_{i}}  (6)

The correspondence between the two representations can be derived from the so-called Sherman–Morrison–Woodbury formula [13] for inverses of matrices subject to low rank perturbations.

Broyden suggested this formula as well, but apparently he and others had less favourable numerical experience, which lead to the moniker Bad Broyden update. It is not clear whether this judgement is justified, since the formula has at least two nice features. First, the inverse is always well defined, whereas the inverse of the good Broyden update can be seen to blow up if y_{i}^{⊤}B_{i−1}s_{i} = 0. Second, the bad Broyden update is invariant with respect to linear variable transformations in that applying it to the system F̃(x̃) ≡ F (T x̃) = 0 with det(T) ≠ 0 leads to a sequence of iterates x̃_{i} related to the original ones by x_{i} = T x̃_{i}, provided one initializes x̃_{0} = T^{−1}x_{0} and B̃_{0} = B_{0}T. The good Broyden formula, on the other hand, is dependent on the scaling of the variables via the Euclidean norm, but is independent of the scaling of the residuals, which strongly influences the bad Broyden formula. However, even for quasi-Newton methods based on the good Broyden update, the squared residual norm often enters through the back door, namely as merit function during a line-search. The resulting stabilized nonlinear equation solver is strongly affected by linear transformations on domain or range. In this brief survey we have only considered full step iterations and their local convergence properties.

Whether or not one should implement quasi-Newton methods by storing and manipulating the inverses H_{i} is a matter for debate. Originally, Broyden and his colleagues had apparently no qualms about this, but later it was widely recommended, e.g., by the Stanford school [14], that one should maintain a triangular factorization of the B_{i} for reasons of numerical stability. Now it transpires that the required numerical linear algebra games, e.g., chasing subdiagonal entries, are rather slow on modern computer architectures. In any case, the trend is to limited memory implementations for large scale applications, in view of which we will first try to study the influence of the variable number n on Broyden updating.

Estimating the R-order and efficiency index

One might fault the property of Q-superlinear convergence for being not sufficiently discriminating, because it can be established for all halfway sensible updating methods. In view of the limiting case of operator equations on Hilbert spaces to be considered later, one may wonder how the convergence rate of quasi-Newton methods depends on the dimension n. A finer measure of how fast a certain sequence x_{i} → x_{∗} convergences is the so-called R-order

ρ ≡ lim inf_{i} |log ∥x_{i} − x_{∗}∥|^{1/i}

The limit inferior on the right reduces to a proper limit when the sequence x_{i} → x_{∗} satisfies ∥x_{i} − x_{∗}∥ ∼ ∥x_{i−1} − x_{∗}∥^{ρ}. This is well known to hold with ρ = 2 for all iterations generated by Newton’s method from an x_{0} close to a regular root x_{∗}. Generally, the R-order [27] of a method is the infimum over ρ for all locally convergent sequences (x_{i})_{i=1...∞}.

The result of Burmeister and Gay implies 2 n step quadratic convergence of Broyden’s good method on smooth nonlinear equations. That corresponds to an R-order of \sqrt[2n]{2} = 1 + 1/(2n) + O(1/n^{2}). We may actually hope for just a little more by the following argument adapted from a rather early paper of Janina Jankowska [21]. Whenever a Jacobian approximation B_{i} is based solely on the function values F_{i−j} = F (x_{i−j}) , for j = 0 … n, its discrepancy to the Jacobian F′(x_{∗}) is likely to be of order O(∥x_{j−n} − x_{∗}∥). Here we have assumed that things are going well in that the distances ∥x_{i} − x_{∗}∥ decrease monotonically towards 0, so that the function value at the oldest iterate x_{i−n} contaminates B_{i} most. Then the usual analysis of Newton-like iterations [9] yields the proportionality relation

∥x_{i+1} − x_{∗}∥ ∼ ∥x_{i−n} − x_{∗}∥ ∥x_{i} − x_{∗}∥

The first term on the right represents the error in the approximating Jacobian B_{i} multiplied by the current residual F_{i} of order ∥x_{i} − x_{∗}∥. Substituting the ansatz ∥x_{i} − x_{∗}∥ ∼ c^{ρ^{i}} for some c ∈ (0, 1) into the recurrence and then taking the log base c one obtains immediately the relations

ρ^{i+1} ∼ ρ^{i−n} + ρ^{i} ⟹ 0 = P_{n}(ρ) ≡ ρ^{n+1} − 1 − ρ^{n}

Hence, we can conclude that the best R-order we may expect from Broyden updating is the unique positive root ρ_{n} of the polynomial P_{n}(ρ).

For n = 1, both Broyden updating methods reduce to the classical secant scheme, which is well known [27] to have the convergence order ρ_{1} = (1+\sqrt{5})/2. The larger n, the smaller ρ_{n}, and it was shown in [19] that asymptotically

P_{n}^{−1}(0) ∋ ρ_{n} ≈ 1 + ln(n)/n ≈ \sqrt[n]{n}

Here a_{n} ≈ b_{n} means that the ratio a_{n}/b_{n} tends to 1 as n goes to infinity. The second approximation means that we may hope for n step convergence of order n rather than just 2 n step convergence of order 2 as suggested by the result of Burmeister and Gay.

The first approximation implies that the efficiency index [28] in the sense of Ostrowski (namely the logarithm of the R-order divided by the evaluation cost and linear algebra effort per step) satisfies asymptotically

\frac{ln(ρ_{n})}{OP S(F) + O(n^{2})} ≈ \frac{ln(n)/n}{OP S(F) + O(n^{2})} ≥ \frac{ln(2)}{n OP S(F) + O(n^{3})}

The lower bound on the right-hand side represents Newton’s method with divided difference approximation of the Jacobian, and dense refactorization at each iteration. As we can see there is a chance for Broyden updating to yield an efficiency index that is ln(n)/ ln(2) = log_{2} n times larger than for Newton’s method under similar conditions.

This hope may not be in vain since it was shown in [19] that the R-order ρ_{n} is definitely achieved when the Jacobian is updated by the adjoint Broyden formula

B_{i} = B_{i−1} + r_{i}r_{i}^{⊤}(F′(x_{i}) − B_{i−1})/ r_{i}^{⊤}r_{i}

However, this rank-one-update is at least twice as expensive to implement since it involves the transposed product F′(x_{i})^{⊤}r_{i}, which can be evaluated in the reverse mode of Algorithmic Differentiation. The latter may be three times as expensive as pure function evaluation, so that the efficiency gain on Newton’s method can be bounded below by (log_{2} n)/4 = log_{16} n.

Whether or not simple Broyden updating itself achieves the optimal R-order ρ_{n} has apparently not yet been investigated carefully. To be fair, it should be noted that taking roughly n/ log(n) simplified Newton steps before reevaluating and refactorizing the Jacobian in the style of Shamanskiĭ [22], yields the convergence order near 1 + n/ log(n) for any such cycle and the corresponding effort is approximately [n OP S(F) + O(n^{3})][1 + 1/ log(n)]. The resulting efficiency index is asymptotically identical to the optimistic estimate for Broyden updating derived above.

Pushing n to infinity

While Broyden updating is well established in codes for small and medium scale problems, its usefulness for large dimensional problems is generally in doubt. The first author who applied and analyzed Broyden’s method to a control problem in Hilbert space was Ragnar Winther [31]. Formally, it is easy to extend the Broyden method to an operator equation y = F (x) = 0 between a pair of Hilbert spaces X and Y. One simply has to interpret transposition as taking the adjoint so that v^{⊤} represents a linear function in X = X^{∗} such that v^{⊤}w ≡ ⟨v, w⟩ yields the inner product. The Good Broyden Update is still uniquely characterized by the nondeterioration condition (4) in terms of the operator norm ∥ · ∥. This implies bounded nondeterioration in the nonlinear case and everything needed to derive local and linear convergence goes through.

However, the least change characterization and its consequences cannot be extended, because there is no generalization of the Frobenius norm (3) and the underlying inner product to the space B(X, Y) of bounded linear operators. To see this, we simply have to note that, in n dimensions, the Frobenius norm of the identity operator is n, the sum of its eigenvalues. That sum would be infinite for the identity on l^{2}, the space of square summable sequences to which all separable Hilbert spaces are isomorphic. There is apparently also no other inner product norm on B(X, Y) that is at least as strong as the operator norm so that the implication (5) would work.

These are not just technical problems in extending the superlinear result, since X is infinite dimensional exactly when the unit ball and, equivalently, its boundary, the unit sphere, are not compact. That means one can keep generating unit directions s̄_{i} ≡ s_{i}/∥s_{i}∥ along which the current approximation B_{i} is quite wrong. Such an example with an orthogonal sequence of s_{i} was given by Griewank [18]. There, on an affine bicontinuous problem, Broyden’s method with full steps converges only linearly or not at all.

To derive the basic properties of Broyden’s method in Hilbert space we consider an affine equation 0 = F (x) ≡ A x − b with a bounded invertible operator A ∈ B(Y, X). Then we have the discrepancies

D_{i} = A^{−1}B_{i} − I ∈ B(X, Y) and E_{i} ≡ D_{i}^{⊤}D_{i} ∈ B(X)

where D_{i}^{⊤} ∈ B(Y, X) denotes the adjoint operator to D_{i} and we abbreviate B(X) ≡ B(X, X) as usual. By definition, E_{i} is selfadjoint and positive semidefinite. Now the Broyden good update can be rewritten as

D_{i+1} = D_{i} (I − s̄_{i}s̄_{i}^{⊤}) ⟹ E_{i+1} ≡ E_{i} − r̄_{i}r̄_{i}

with r̄_{i} ≡ A^{−1}r_{i}/∥s_{i}∥.

In the finite dimensional case one could show that the Frobenius norm of the D_{i} decreases monotonically. Now we see that the operators E_{i} are obtained from the E_{0} = D_{0}^{⊤}D_{0} by the consistent subtraction of rank-one terms. Hence, they have a selfadjoint semidefinite limit E_{∗}. This implies, by a generalization of the interlacing eigenvalue theorem, that the eigenvalues (λ_{j}(E_{i}))_{j=1...∞} of E_{i} are monotonically declining towards their limits (λ_{j}(E_{∗}))_{j=1...∞}. Correspondingly, we find for the singular values σ_{j}(D_{i}) = \sqrt{λ_{j}(E_{i})} of the D_{i} that

σ_{j}(D_{i+1}) ≤ σ_{j}(D_{i}) and σ_{j}(D_{i}) → \sqrt{λ_{j}(E_{∗})} for i → ∞

Similarly, it was proven by Fletcher that the BFGS update monotonically moves all eigenvalues of the symmetric discrepancy B_{∗}^{−1/2}B_{i}B_{∗}^{−1/2} − I between the Hessian B_{∗} and its approximations B_{i} towards zero. With regards to convergence speed it was shown in [18] for C^{1,1} operator equations that Broyden’s method yields locally

lim sup_{i→∞} ∥A^{−1}F_{i+1}∥/ ∥A^{−1}F_{i}∥ ≤ σ_{∞}(D_{0}) ≡ lim_{j→∞} σ_{j}(D_{0})

In other words, the Q-factor is bounded by the essential spectrum σ_{∞}(D_{0}) of the initial relative discrepancy D_{0} = A^{−1}B_{0} − I. Hence, we must have Q-superlinear convergence if D_{0} or, equivalently, just B_{0} − A is compact, an assumption that is of course trivial in finite dimensions. Equivalently, we can require the preconditioned discrepancy D_{0} to be compact or at least to have a small essential norm. Thus we can conclude that Broyden updating will yield reasonable convergence speed in Hilbert space if D_{0} is compact or has at least a small essential norm σ_{∞}(D_{0}) = σ_{∞}(D_{j}). It is well known that the essential norm is unaffected by modifications of finite rank. On the other hand, all singular values σ_{j}(D_{0}) > σ_{∞}(D_{0}) are effectively taken out as far as the final rate of convergence is concerned.

Limited memory and data sparse

For symmetric problems the idea of limited memory approximations to the Hessian of the objective [24] has been a roaring success. In the unsymmetric case things are not so clear. Whereas in the unconstrained, quadratic optimization case conjugate gradients generates the same iterates as BFGS in an almost memoryless way, there is, according to a result of Faber and Manteuffel [11], no short recurrence for unsymmetric real problems. Correspondingly, the more or less generic iterative solver GMRES for linear problems requires 2 i vectors of storage for its first i iterations. The same appeared be true of Broyden’s method, where starting from a usually diagonal B_{0}, one could store the secant pairs (s_{j}, y_{j}) for j = 1 … i.

The same appeared be true for Broyden’s method in inverse form, where starting from an usually diagonal H_{0} = B_{0}^{−1} one could store the secant pairs (s_{j}, z_{j}) with z_{j} ≡ H_{j−1}y_{j} for j = 1 … i. Then the inverse Hessian approximations have the product representation

H_{i} = [H_{i−1} + \frac{(s_{i} − z_{i})s_{i}^{⊤}H_{i−1}}{s_{i}^{⊤}z_{i}} ] = ∏_{j=1}^{i} [I + \frac{(s_{j} − z_{j})s_{j}^{⊤}}{s_{j}^{⊤}z_{j}} ] H_{0}

Deuflhard et al. noticed in [10] that for the fullstep iteration successive s_{j} and s_{j+1} = −H_{j}F_{j} satisfy the relation s_{j+1} = (s_{j} − z_{j})∥s_{j}∥^{2}/s_{i}^{⊤}z_{j}. Hence, one only needs to store the s_{j} and one can then cheaply reconstruct the z_{j} for applying the inverse in product form to any vector v usually the current residual F_{i}. Hence the storage requirement is only i + O(1) vectors of length n up to the i-th iteration. In contrast the storage requirement for i iterations of Bad Broyden appears to be twice as large [10], so at least in that sense the derogatory naming convention is justified. In either case, one normally wishes to limit the number of vectors to be stored a priori and thus one has to develop strategies for identifying and discarding old information. This issue has been extensively studied for the limited memory BFGS method and for Broyden updating it has been the focus of a recent PhD thesis [30]. Usually one wishes to get rid of information from earlier iterates because nonlinearity may render it irrelevant or even misleading near the current iterates. On discretizations of infinite dimensional problems, one may wish to discard all corrections of a size close to the essential norm σ_{∞}(D_{0}), since no amount of updating can reduce that threshhold.

In good Broyden updating the correction made to any row of the approximating Jacobian is completely independent of what goes on in the other rows. In other words we are really updating the gradients ∇F_{k} of the component functions F_{k} independently. That shows immediately that one can easily use the method for approximating rectangular Jacobians F′(x) for F : ℝ^{n} ↦ ℝ^{m} with m independent of n. Also in updating the k−th row one can disregard all variables that have no impact on F_{k} so that the corresponding Jacobian entries are zero. The resulting sparse update is known as Schubert’s method [5]. The least change characterization now applies in the linear subspace of matrices with the appropriate sparsity pattern, and the whole BDM locally linear and Q-superlinear convergence goes through without any modification. However, since the update matrices C_{j} are now of high rank, there is no longer any advantage compared to Newton’s method with regards to the linear algebra effort per step.

On the other hand, large sparse Jacobians can often be evaluated exactly, possibly using algorithmic differentiation [16], at an entirely reasonable cost. In particular it was found that none of the constraint Jacobians in the optimization test collection CUTEr takes more than 18 times the effort of evaluating the vector functions of constraints themselves. Since the sparsity patterns also tend to be quite regular, no methods based on Broyden type updating [15] can here compete with methods based on exact derivatives values.

Whether or not that situation is really representative for problems from applications is not entirely clear.

In any case we have to count the inability to effectively exploit sparsity as part of the Bad about Broyden updating. Still, there is a lot of Good as well, for which we have to thank primarily Charles Broyden, who passed away last year at the age of 78 after an eventful life with various professional roles and countries of residence.

Acknowledgement. The author is indebted to Jorge Moré, Trond Steihaug, and other colleagues for discussions on the historical record.

References

[1] Joszef Abaffy, Charles Broyden, and Emilio Spedicato. A class of direct methods for linear systems. Numer. Math., 45(3):361–376, 1984.

[2] C. G. Broyden. A class of methods for solving nonlinear simultaneous equations. Math. Comp., 19:577–593, 1965.

[3] C. G. Broyden. The convergence of single-rank quasi-Newton methods. Math. Comp., 24:365–382, 1970.

[4] C. G. Broyden. Recent developments in solving nonlinear algebraic systems. In Numerical methods for nonlinear algebraic equations (Proc. Conf., Univ. Essex, Colchester, 1969), pages 61–73. Gordon and Breach, London, 1970.

[5] C. G. Broyden. The convergence of an algorithm for solving sparse nonlinear systems. Math. Comp., 25:285–294, 1971.

[6] C. G. Broyden. On the discovery of the “good Broyden” method. Math. Program., 87(2, Ser. B):209–213, 2000. Studies in algorithmic optimization.

[7] C. G. Broyden, J. E. Jr. Dennis, and J. J. Moré. On the local and superlinear convergence of quasi-Newton methods. JIMA, 12:223–246, 1973.

[8] J. E. Dennis, Jr. and Jorge J. Moré. A characterization of superlinear convergence and its application to quasi-Newton methods. Math. Comp., 28:549–560, 1974.

[9] J. E. Jr. Dennis and R. B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Prentice-Hall, 1996.

[10] Peter Deuflhard, Roland Freund, and Artur Walter. Fast secant methods for the iterative solution of large nonsymmetric linear systems. In IMPACT of Computing in Science and Engineering, pages 244–276, 1990.

[11] Vance Faber and Thomas Manteuffel. Necessary and sufficient conditions for the existence of a conjugate gradient method. SIAM J. Numer. Anal., 21(2):352–362, 1984.

[12] D. M. Gay. Some convergence properties of Broyden’s method. SIAM J. Numer. Anal., 16:623–630, 1979.

[13] D. M. Gay and R. B. Schnabel. Solving systems of nonlinear equations by Broyden’s method with projected updates. In Nonlinear Programming 3, O. Mangasarian, R. Meyer and S. Robinson, eds., Academic Press, NY, pages 245–281, 1978.

[14] Philip E. Gill, Walter Murray, and Margaret H. Wright. Numerical linear algebra and optimization. Vol. 1. Addison-Wesley Publishing Company Advanced Book Program, Redwood City, CA, 1991.

[15] A. Griewank and A. Walther. On constrained optimization by adjoint based quasi-Newton methods. Opt. Meth. and Soft., 17:869–889, 2002.

[16] A. Griewank and A. Walther. Principles and Techniques of Algorithmic Differentiation, Second Edition. SIAM, 2008.

[17] Andreas Griewank. The “global” convergence of Broyden-like methods with a suitable line search. J. Austral. Math. Soc. Ser. B, 28(1):75–92, 1986.

[18] Andreas Griewank. The local convergence of Broyden-like methods on Lipschitzian problems in Hilbert spaces. SIAM J. Numer. Anal., 24(3):684– 705, 1987.

[19] Andreas Griewank, Sebastian Schlenkrich, and Andrea Walther. Optimal r-order of an adjoint Broyden method without the assumption of linearly independent steps. Optim. Methods Softw., 23(2):215–225, 2008.

[20] Hermann Hammer and Kerstin Dambach. Sharaf al-tusi, ein vorläufer von newton und leibnitz. Der mathematische und naturwissenschaftliche Unterricht, 55(8):485–489, 2002.

[21] Janina Jankowska. Theory of multivariate secant methods. SIAM J. Numer. Anal., 16(4):547–562, 1979.

[22] C. T. Kelley. A Shamanskiĭ-like acceleration scheme for nonlinear equations at singular roots. Math. Comp., 47(176):609–623, 1986.

[23] Dong-Hui Li and Masao Fukushima. A derivative-free line search and global convergence of Broyden-like method for nonlinear equations. Optim. Methods Softw., 13(3):181–201, 2000.

[24] Dong C. Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. Math. Programming, 45(3, (Ser. B)):503–528, 1989.

[25] José Mario Martínez. Practical quasi-Newton methods for solving nonlinear systems. J. Comput. Appl. Math., 124(1–2):97–121, 2000. Numerical analysis 2000, Vol. IV, Optimization and nonlinear equations.

[26] J. J. Moré and J. A. Trangenstein. On the global convergence of Broyden’s method. Math. Comp., 30(135):523–540, 1976.

[27] J. M. Ortega and W. C. Reinboldt. Iterative Solution of Nonlinear Equations in Several Variables. Academic Press, 2000.

[28] A. Ostrowski. Solution of Equations and Systems of Equations. Academic Press, New York, 1966.

[29] Joanna Maria Papakonstantinou. Historical Development of the BFGS Secant Method and Its Characterization Properties. PhD thesis, Rice University, Houston, 2009.

[30] Bart van de Rotten. A limited memory Broyden method to solve highdimensional systems of nonlinear equations. PhD thesis, Mathematisch Instituut, Universiteit Leiden, The Netherlands, 2003.

[31] Ragnar Winther. A numerical Galerkin method for a parabolic control problem. PhD thesis, Cornell University, 1977.

Andreas Griewank

Institut für Mathematik

Humboldt Universität zu Berlin Unter den Linden 6 10099 Berlin Germany

               griewank@mathematik.hu-berlin.de

316

Documenta Mathematica · Extra Volume ISMP (2012)