Introduction

In this paper, we consider the unconstrained minimization problem min

where f : R d → R is strongly convex and twice continuously differentiable. We focus on the non-asymptotic global convergence properties of quasi-Newton methods for solving Problem (1). The core idea behind quasi-Newton methods is to mimic the update of Newton's method using only first-order information, i.e., the gradients of f . Specifically, the update rule at the k-th iteration is

where η k is the step size and B k ∈ R d×d is a matrix constructed from the gradients of f to approximate the Hessian ∇ 2 f (x k ). Various quasi-Newton methods have been developed, each distinguished by its strategy for constructing the Hessian approximation B k and its inverse. The key methods among them are the Davidon-Fletcher-Powell (DFP) method [11,17], the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [5,15,19,39], the Symmetric Rank-One (SR1) method [10,26], and the Broyden method [3]. Notably, these quasi-Newton methods directly maintain and update the inverse matrix B -1 k using a constant number of matrix-vector multiplications. This results in a computational cost of O(d 2 ) per iteration and thus makes quasi-Newton methods more efficient than Newton's method, which involves computing the Hessian and solving a linear system that could incur a computational cost of O(d 3 ) per iteration.

Compared to other first-order methods, such as gradient descent and accelerated gradient descent, the primary advantage of quasi-Newton methods is their ability to achieve Q-superlinear convergence, i.e., lim

where x * ∈ R d denotes the optimal solution of Problem (1). Specifically, [6] and [13] have established that both DFP and BFGS converge Q-superlinearly with unit step size η k = 1, where the initial point x 0 is required to be within a local neighborhood of the optimal solution x * . Later, it has also been extended to various settings [1,12,18,21,28,31,42,45]. However, these local convergence results are all asymptotic and fail to provide an explicit convergence rate after a finite number of iterations.

Recently, there has been progress regarding non-asymptotic local convergence analysis of quasi-Newton methods. The authors of [38] showed that, if the initial point x 0 is in a local neighborhood of the optimal solution x * and the initial Hessian approximation matrix B 0 is initialized as L I , then BFGS with unit step size attains a local superlinear convergence rate of the form ( d L μk ) k , where d is the problem's dimension, L is the Lipschitz parameter of the gradient, and μ is the strong convexity parameter. Later in [37], the local convergence rate of BFGS was improved to ( d log (L/μ) k ) k under similar initial conditions. Similar local superlinear convergence analysis has also been established for the SR1 method [43]. In a concurrent work [25], it was demonstrated that, if x 0 is in a local neighborhood of the optimal solution x * and B 0 is sufficiently close to the exact Hessian at the optimal solution (or selected as the exact Hessian at x 0 ), then BFGS with unit step size achieves a local superlinear rate of (1/k) k/2 , which is independent of the dimension d and the condition number L/μ. While these non-asymptotic results successfully characterize an explicit superlinear rate, they rely heavily on local analysis, requiring the initial point to be sufficiently close to the optimal solution x * and imposing conditions on the step size and the initial Hessian approximation matrix B 0 . Consequently, these results cannot be directly extended to a global convergence guarantee. We discuss this issue in detail in Sect. 6.

To guarantee global convergence, quasi-Newton methods must be combined with line search or trust-region techniques. The first global result for quasi-Newton methods was derived by Powell in [34], where it was established that DFP with exact line search converges globally and Q-superlinearly. Later, Dixon [14] proved that all quasi-Newton methods from the convex Broyden's class generate the same iterates using exact line search, thus extending Powell's result to the convex Broyden's class including BFGS. In order to relax the exact line search condition, the work in [35] considered BFGS using inexact line search based on Wolfe conditions and showed that it retains global superlinear convergence. This result was later extended in [9] to the convex Broyden class except for DFP. Moreover, [7,10,26] showed that the SR1 method with trust-region techniques achieves global and superlinear convergence.

However, all these results lack an explicit global convergence rate; they only provide asymptotic convergence guarantees and fail to characterize the explicit global convergence rate of classical quasi-Newton methods. The only exception is a recent work in [27], where the authors also studied the global convergence rate of BFGS with exact line search. Specifically, it was shown that BFGS attains a global linear rate of (1 -2μ 3 L 3 (1

Lk ) -1 ) k , where Tr(•) denotes the trace of a matrix. We note that after k = O(d) iterations, their linear rate approaches the rate of (1 -2μ 3 L 3 ) k , which is substantially slower than gradient descent-type methods. More importantly, their study does not extend to demonstrating any superlinear convergence rate and fails to fully characterize the behavior of BFGS.

The discussions above reveal a major gap in classical quasi-Newton methods: the lack of an explicit global convergence rate characterization. Contributions. In this paper, we present the first results that contain explicit nonasymptotic global linear and superlinear convergence rates for the BFGS method with exact line search. Note that due to the equivalence result by Dixon [14], our results also hold for other quasi-Newton methods in the convex Broyden class with exact line search. At a high level, our convergence analysis sharpens the potential functionbased framework first introduced in [8], leading to a unifying framework for proving both the global linear convergence rates and the superlinear convergence rates. Our convergence results are global as they hold for any initial point x 0 ∈ R d and any initial Hessian approximation matrix B 0 that is symmetric positive definite. Specifically, our analysis divides the convergence process into three phases, characterized by different convergence rates: (i) First linear phase: We show that

Here, B0 = 1 L B 0 is the scaled initial Hessian approximation matrix, Ψ (•) is a potential function defined later in (18), κ = L μ denotes the condition number, and

is based on the initial optimality gap with M as the Hessian's Lipschitz parameter. In particular, when k ≥ Ψ ( B0 ), this leads to a linear rate of

κ}, the algorithm attains an improved linear rate matching that of standard gradient descent:

is the normalized initial Hessian approximation matrix.

To make our convergence rates easily interpretable, we focus on the global linear and superlinear convergence rates of the special case where B 0 = α I for a given scalar α > 0. We also study the practical initialization, B 0 = cI , where c = s y

, and x 1 , x 2 being two randomly chosen vectors. We further consider B 0 = L I and B 0 = μI as two specific cases. The global convergence results with these initializations are summarized in Table 1. Our analysis reveals a trade-off between the linear and the superlinear rates, depending on the choice of the initial matrix B 0 . Specifically, while both initializations lead to the same linear convergence rates, initiating with B 0 = L I allows the algorithm to reach this rate d log κ iterations earlier than with B 0 = μI . On the other hand, for the superlinear convergence phase, the difference between B 0 = L I and B 0 = μI essentially boils down to comparing dκ against (1 + 4C 0 )d log κ. Thus, when C 0 κ, initializing with B 0 = μI enables an earlier transition to the superlinear convergence compared

to B 0 = L I , as well as a faster superlinear convergence rate. As we shall see in Sect. 7, our experiments also demonstrate this trade-off.

Additional related work.

In addition to the standard quasi-Newton methods such as BFGS, the superlinear convergence of other variants of quasi-Newton methods has also been studied in the literature. The greedy variants of quasi-Newton methods were first introduced in [36] and developed in subsequent works [22,29,30]. Instead of using the difference of successive iterates to update the Hessian approximation matrix, the key idea is to greedily select basis vectors to maximize a certain measure of progress. In [36], greedy BFGS is shown to achieve a local superlinear convergence rate of (dκ(1

) k and the superlinear convergence phase begins after dκ ln (dκ) iterations. Similar superlinear convergence rates are extended to other greedy quasi-Newton updates in [22,29,30]. However, we note that their results are all local and require the initial point to be sufficiently close to the optimal solution x * . Recently, along a different line of work, the authors in [23,24] proposed quasi-Newton-type methods based on the hybrid proximal extragradient framework [32,40] and studied their global convergence rates. Specifically, it was shown that the quasi-Newton proximal extragradient method in [23] achieves a global linear convergence rate of (1 -1/κ) k and a global superlinear rate of the form (1 + k/O(κ 2 d)) -k . However, these methods are distinct from the classical quasi-Newton methods such as BFGS analyzed in this paper, since they formulate the update of the Hessian approximation matrices B k as an online convex optimization problem and follow an online learning algorithm to update B k . Outline. In Sect. 2, we provide an overview of the BFGS method with exact line search, outline our assumptions, and introduce some preliminary lemmas for the exact line search scheme. Sect. 3 presents our general analytical framework, which is employed to establish global linear and superlinear convergence results for the BFGS method, along with the intermediate results for the update of quasi-Newton methods. In Sect. 4, we establish the global linear convergence rate of BFGS using exact line search that applies to any choices of B 0 and x 0 and we consider specific initializations with B 0 = α I , B 0 = cI , B 0 = L I and B 0 = μI . Building on the linear convergence rates, Sect. 5 details our global superlinear convergence results. In Sect. 6, we compare our analytical framework to both classical asymptotic analysis and recent local non-asymptotic analysis of BFGS. Sect. 7 displays our numerical experiments that corroborate our theoretical findings. Finally, we finish the paper by presenting some concluding remarks in Sect. 8. Notation. We use • to denote the ℓ 2 -norm of a vector or the spectral norm of a matrix. We denote S d + and S d ++ as the set of symmetric positive semidefinite and symmetric positive definite matrices with dimension d × d, respectively. Given two symmetric matrices A and B, we denote A B if and only if B -A is positive semidefinite. Given a matrix A, we use Tr(A) and Det(A) to denote its trace and determinant, respectively.

Preliminaries

In this section, we first outline the assumptions, notations, and lemmas essential for our convergence proof. Following this, we explore the general framework of quasi-Newton methods incorporating exact line search and provide an overview of the principal concepts underpinning the update mechanism in the convex Broyden's class of quasi-Newton methods, which encompasses both the BFGS and DFP algorithms.

Assumptions

To begin with, we state our assumptions on the objective functions f .

Assumption 1

The objective function f is strongly convex with parameter μ > 0, i.e., ∇ f (x) -∇ f (y) ≥ μ xy , for any x, y ∈ R d .

Assumption 2

The objective function gradient ∇ f is Lipschitz continuous with parameter L > 0, i.e., ∇ f (x) -∇ f (y) ≤ L xy for any x, y ∈ R d .

Both Assumptions 1 and 2 are standard in the convergence analysis of first-order methods. Moreover, since f is twice differentiable, they imply that μI ∇ 2 f (x) L I for any x ∈ R d . Additionally, the condition number of f is defined as κ := L μ . We also remark that Assumptions 1 and 2 are sufficient in our analysis to prove a global linear convergence rate of BFGS with exact line search. In order to achieve a superlinear convergence rate, we need to impose an additional assumption on the Hessian of the function f , stated below.

Assumption 3

The objective function Hessian ∇ 2 f is Lipschitz continuous along the direction of optimal solution x * with parameter M > 0, i.e.,

Assumption 3 is commonly used in the analysis of quasi-Newton methods, such as in [8], as it provides a necessary smoothness condition for the Hessian of the objective function. Importantly, we do not require Hessian smoothness for arbitrary points x, y ∈ R d ; rather, we impose a weaker condition, assuming that the Hessian is Lipschitz continuous only along the direction of the optimal solution x * .

Quasi-Newton methods with exact line search

Next, we briefly review the template for updating quasi-Newton matrices, focusing specifically on the DFP and BFGS algorithms. Specifically, at the k-th iteration, the update in (2) can be equivalently written as

where

Here, η k ≥ 0 represents the step size, and B k ∈ R d×d is the Hessian approximation matrix. Replacing B k with the exact Hessian ∇ 2 f (x k ) turns the update into the classical Newton's method. Quasi-Newton methods aim to approximate the Hessian with firstorder information, typically adhering to a secant condition and a least-change property.

To elaborate, we define the variable difference s k and gradient difference y k as

The secant condition requires that B k+1 satisfies y k = B k+1 s k , ensuring the gradient consistency between the quadratic model

) and f at x k and x k+1 ; that is, [33,Chapter 6]). That said, the secant condition does not uniquely define B k+1 . Thus, we impose a least-change property to ensure that B k+1 , satisfying the secant condition, is closest to B k in a specific proximity measure. Various proximity measures have been proposed in the literature [16,19,20] and here we follow the variational characterization in [16]. Specifically, for any symmetric positive definite matrix A ∈ S d ++ , define the negative log-determinant function Φ(A) =log Det(A) and define the Bregman divergence generated by Φ by

Note that the Bregman divergence can be regarded as a measure of proximity between two positive definite matrices, and D Φ (A, B) = 0 if and only if A = B. For the BFGS update, it was shown in [16] that B k+1 is given as the unique solution of the minimization problem:

which admits the following explicit update rule:

Moreover, if we define H k := B -1 k as the inverse of the Hessian approximation matrix, it follows from the Sherman-Morrison-Woodbury formula that

The DFP update rule can be regarded as the dual of BFGS, where the roles of the Hessian approximation matrix B k+1 and its inverse H k+1 are exchanged. Specifically, the DFP update rules are given by

Both BFGS and DFP belong to a more general class of quasi-Newton methods, known as the convex Broyden's class [4]. In this class, the Hessian approximation matrix B k+1 is defined as

where φ k ∈ [0, 1] for any k ≥ 0. Accordingly, there exists ψ k ∈ [0, 1] such that the Hessian inverse approximation matrix H k+1 is given by

The convex Broyden's class exhibits a crucial property: if the initial Hessian approximation matrix B 0 is symmetric positive definite and the objective function f is strictly convex, then all subsequent B k matrices produced by this class maintain symmetric positive definiteness (see [33]).

To guarantee the global convergence of quasi-Newton methods in (4), it is necessary to employ a line search scheme to select the step size η k . In this paper, our primary focus is on the exact line search step size, where we aim to minimize the objective function along the search direction d k . Specifically,

Remarkably, it was shown in [14] that, when employing the exact line search scheme, the convex Broyden's class of quasi-Newton methods produce identical iterates given that the initial point x 0 and the initial matrix B 0 are the same. Thus, in the remainder of the paper, we focus on the BFGS update in (7) as all results hold for other algorithms in the convex Broyden family. Finally, we introduce some intermediate results related to the exact line search step size, as defined in (9). These standard results (see, e.g., [34]) are essential for the forthcoming demonstration of the convergence rate of the quasi-Newton method.

Lemma 1 Consider the standard quasi-Newton method in (4) with the exact line search specified in (9). The following results hold for any k ≥ 0:

Convergence analysis framework

In this section, we introduce our theoretical framework for establishing the global convergence rates of the BFGS algorithm with exact line search. Our framework builds on two key propositions. In Proposition 1, we characterize the amount of function value decrease in one iteration in terms of the angle θ k between the steepest descent direction -g k and the search direction d k given in (4). Subsequently, Proposition 2 presents a potential function for the BFGS update, which leads to a lower bound on cos(θ k ).

To formally begin the analysis, we first present a weighted version of key vectors and matrices as introduced from the previous work [13]. Specifically, given a weight matrix P ∈ S d ++ (also referred to as a transformation matrix), we define the weighted gradient ĝk , the weighted gradient difference ŷk , and the weighted iterate difference ŝk as

Similarly, we define the weighted Hessian approximation matrix Bk as

Note that the weight matrix P can be chosen as any positive definite matrix, and its choice will be evident from the context. In particular, as we shall see later, we use P = L I in Sect. 4 to prove the global linear convergence rate, and use P = ∇ 2 f (x * ) in Sect. 5 to prove the global superlinear convergence rate. Moreover, since the above weighting procedure amounts to a change of the coordinate system, the weighted versions of the vectors and matrices defined in (10) and (11) retain the same algebraic relations as their original forms. In particular, the weighted Hessian approximation matrices generated by the BFGS algorithm follow the subsequent update rule:

Before introducing our first key proposition, we define a quantity θk by

which is the angle between the weighted steepest descent direction -ĝk and the weighted iterate difference ŝk . It is well-known that the convergence of quasi-Newton methods can be established by monitoring the behavior of cos( θk ). We next quantify the link between functional value decrease and cos( θk ).

Proposition 1 Let {x k } k≥0 be the iterates generated by the BFGS method with exact line search. Given a weight matrix P ∈ S d ++ , recall the weighted vectors and matrices defined in (10) and (11). For any k ≥ 0, we have

where we define

As a corollary, we have that for any k ≥ 1,

Proof First, we use the definition of αk in (15) to write

Moreover, note that we have -ĝ k ŝk = ŷ k ŝk by Lemma 1(b). Hence, using the definition of θk in ( 13) and the definition of mk in (15), it follows that

Furthermore, we have ĝk

) from the definition of qk in (15). Thus, the equality in (17) can be rewritten as

By rearranging the term in the above equality, we obtain (14). To prove the inequality in (16), note that for any k ≥ 1, we have

where the last equality is due to (14). Notice that the term 1 -αi qi mi cos 2 ( θi ) are nonnegative for any i ≥ 0. Thus, by applying the inequality of arithmetic and geometric means twice, we obtain that

This completes the proof.

Remark 1

We note that similar results relating f (x k )f (x k+1 ) to cos 2 ( θk ) have appeared in prior work such as [9, Lemma 4.2] and [8], though they are used in the analysis of quasi-Newton methods with inexact line search. Compared with these prior results, Proposition 1 is more general in the sense that we consider the weighted iterates using a general weight matrix P. This flexibility enables us to obtain tighter bounds and, more importantly, to obtain a global superlinear convergence rate under the same framework (see Sect. 5). Another subtle yet important difference is that previous works typically upper bound the term mk by L prematurely, leading to the worst dependence on the condition number κ. Instead, we keep mk in ( 14) as is and lower bound the term cos 2 ( θk )/ mk together, as later shown in Proposition 2.

Remark 2 without loss of generality, we assume that x k = x * for any k ≥ 0 throughout the paper, meaning the exact optimal solution x * is never reached. Consequently, this ensures that ĝk = 0, η k > 0, and ŝk = 0 for any k ≥ 0. As a result, we have

since Bk is symmetric positive definite. Moreover, under this assumption, the definitions in equation ( 15) are valid since all the denominators-

and ŝk 2 -are nonzero. Similarly, mk = ŷ k ŝk ŝk 2 is well-defined and nonzero, as

Proposition 1 shows that BFGS's convergence rate hinges on four quantities: αk , qk , mk , and cos( θk ). Note that αk and qk can be bounded using Assumptions 1, 2 and 3, independent of the quasi-Newton update, with details deferred to Sect. 3.1. The focus here is to establish a lower bound for cos 2 ( θk )/ mk . This involves analyzing the dynamics of the Hessian approximation matrices {B k } k≥0 through their trace and determinant, leveraging the following potential function from [8] that integrates both:

Given (6), Ψ (A) can be regarded as the Bregman divergence generated by Φ(A) = log det(A) between the matrix A and the identity matrix I . In particular, Ψ (A) ≥ 0 and also we have Ψ (A) = 0 if and only if A = I . Now we are ready to state Proposition 2, which is a classical result in the quasi-Newton literature (e.g., see [33,Section 6.4]). For completeness, we provide its proof in Appendix A.

Proposition 2 Given a weight matrix P ∈ S d ++ , recall the weighted vectors and matrices defined in (10) and (11). Let { Bk } k≥0 be the weighted Hessian approximation matrices generated by the BFGS update in (12). Then we have

where mk and θk are defined in (15). As a corollary, we have for any k ≥ 1,

Taking exponentiation of both sides in (20), Proposition 2 provides a lower bound for the product k-1

and Ψ ( B0 ). We will use Assumptions 1, 2 and 3 to bound the term ŷk 2 ŝ k ŷk for any k ≥ 0, as shown in Lemma 5 of Sect. 3.1. Moreover, the second term Ψ ( B0 ) depends on our choice of the initial Hessian approximation matrix B 0 . Specifically, we will consider two different initializations: (i) B 0 = L I ; (ii) B 0 = μI . As we shall discuss in the upcoming sections, these two choices result in different bounds and thus lead to a trade-off between the initial linear convergence rate and the final superlinear convergence rate.

Having outlined our key propositions, Sects. 4 and 5 will merge Proposition 1 and Proposition 2 to demonstrate that BFGS achieves global non-asymptotic linear and superlinear convergence rates, respectively. Our approach involves selecting an appropriate weight matrix P and bounding the quantities in (16) to derive the overall convergence rate. Specifically, we set P = L I for global linear convergence and P = ∇ 2 f (x * ) for superlinear convergence. The following section presents intermediate lemmas that will be used to establish these convergence bounds.

Intermediate lemmas

Next, we provide some intermediate results that lower bound the quantities αk and qk defined in (15) and the term ŷk 2 ŝ k ŷk appearing in (19). To do so, we first define the average Hessian matrices J k and G k as

These two matrices play an important role in our analysis, since the fundamental theorem of calculus implies that y k = J k s k and g k = G k (x kx * ) for any k ≥ 0. We also define the weighted average Hessian matrix Ĵk = P -1 2 J k P -1 2 for the given weight matrix P ∈ S d ++ . Moreover, we define a quantity C k that depends on the function value at the iterate x k :

where M is the Lipschitz constant of the Hessian in Assumption 3 and μ is the strong convexity parameter in Assumption 1. Given these definitions, in the following lemma, we characterize the relationship between different matrices that appear in our convergence analysis.

Lemma 2 Suppose Assumptions 1, 2, and 3 hold, and recall the definitions of the matrices J k in (21), G k in (22), and the quantity C k in (23). Then, the following statements hold:

(a) For any k ≥ 0, we have that

(b) For any k ≥ 0, we have that

(c) For any k ≥ 0 and any τ ∈ [0, 1], we have that

(d) For any k ≥ 0 and τ ∈ [0, 1], we have that

Proof Please check Appendix B.

After establishing Lemma 2, in the following three lemmas, we will provide bounds on the quantities αk , qk and

in (15). Suppose Assumptions 1, 2, and 3 hold. Then, for any k ≥ 0, we have

Proof To begin with, note that we have -ĝ k ŝk = ŷ k ŝk due to Lemma 1(b) and ŷ k ŝk = y k s k for any choice of the weight matrix P. Thus, αk can be equivalently

. We first prove the first bound in (28). By Assumptions 1 and 2, the function f is μ-strongly convex and its gradient is L-Lipschitz. Then for any x, y ∈ R d , it holds that

This is also known as the interpolation inequality; see, e.g., [41,Theorem 4]. By setting x = x k , y = x k+1 in (29) and recalling that

and g k+1 = ∇ f (x k+1 ), we obtain that

Moreover, Lemma 1 shows that g k+1 s k = 0 due to exact line search. Thus, we can simplify the above inequality as

where we used Young's inequality in the second inequality and the fact that s k y k ≤ s k y k due to Cauchy-Schwartz inequality in the third inequality. Hence, we conclude that αk

κ . Now we proceed to establish the second lower bound on αk . Given Taylor's theorem, there exists

where we used g k+1 s k = 0. Moreover, based on (26) in Lemma 2(c), we have

Hence, we obtain that

By combining the inequalities in (30) and (31), the main claim follows. (15). Suppose Assumptions 1, 2, and 3 hold. Then we have the following results:

Lemma 4 Recall the definition qk

(b) If we choose P

Proof We first prove (a). When P = L I , we have qk = [2,Section 9.1.2]). Hence, we conclude that qk ≥ 2μ/L = 2/κ.

Next, we prove (b). When P = ∇ 2 f (x * ), we have ĝk

By applying Taylor's theorem with Lagrange remainder, there exists τk ∈ [0, 1] such that

where we used the fact that ∇ f (x * ) = 0 in the last equality. Moreover, by the fundamental theorem of calculus, we have

where we use the definition of G k in (22). Since ∇ f (x * ) = 0 and we denote g k = ∇ f (x k ), this further implies that

Combining (32) and (33) leads to

Based on (27) in Lemma 2(d), we have

Moreover, it follows from (25

Combining (35) and (36), we obtain that

and hence

By using (34) and the fact that ĝk

and the claim follows.

Lemma 5 Suppose Assumptions 1, 2, and 3 hold. Then we have

As a corollary, we have the following results:

(a) If we choose P = L I , then ŷk 2

Proof Note that by the fundamental theorem of calculus, we have y k = J k s k , which implies that ŷk = Ĵk ŝk . Hence, we can bound

Hence, if P = L I , then Ĵk = 1 L J k ≤ 1 by Assumption 2, which proves the result in (a). Moreover, if P = ∇ 2 f (x * ), then 24) in Lemma 2(a), which proves the result in (b).

Global linear convergence rates

In this section, we establish the explicit global linear convergence rates for the BFGS method using an exact line search step size, marking one of the first non-asymptotic global linear convergence analyses of BFGS with a line search scheme. The subsequent global superlinear convergence analyses are established based on these linear rates. Specifically, we combine the fundamental inequality ( 16) from Proposition 1 with lower bounds of the terms αk , qk , and cos 2 ( θk )/ mk from Lemma 3, 4, 5 and Proposition 2 to prove all the global linear convergence rates. In this section, we set the weight matrix P as P = L I and we define the weighted matrix Bk as:

In the following lemma, we prove the global linear convergence rate of the BFGS method for any choice of B 0 ∈ S d ++ .

Lemma 6 Let {x k } k≥0 be the iterates generated by the BFGS method with exact line search and suppose that Assumptions 1 and 2 hold. For any initial point x 0 ∈ R d and any initial Hessian approximation matrix B 0 ∈ S d ++ , we have the following global linear convergence rate for any k ≥ 1,

Proof Our starting point is applying Proposition 1 with the weight matrix P chosen as P = L I . Specifically, (16) shows that to obtain a convergence rate, it suffices to prove a lower bound on k-1 i=0 αi qi mi cos 2 ( θi ). It follows from Lemma 3 that αk =

for any k ≥ 0. Moreover, by applying Lemma 4 with P = L I , we obtain that qk = ĝk

κ for any k ≥ 0. Furthermore, applying Proposition 2 with P = L I , it follows from (20) that

where in the last inequality we used ŷi 2 ŝ i ŷi ≤ 1 by Lemma 5 with P = L I . Taking exponentiation of both sides, this further implies that

Combining all the pieces above, we get

Thus, it follows from Proposition 1 that

This completes the proof.

Notice that this result holds without the Hessian Lipschitz continuity assumption. In the next lemma, we present another version of the global linear convergence analysis with the additional assumption the Hessian of f is M-Lipschitz. We show that the BFGS method with exact line search will eventually reach a global linear convergence rate of (1 -1/O(κ)) k , which is the same as the gradient descent method.

Lemma 7 Let {x k } k≥0 be the iterates generated by the BFGS method with exact line search and suppose that Assumptions 1, 2 and 3 hold. For any initial point x 0 ∈ R d and any initial Hessian approximation matrix B 0 ∈ S d ++ , we have the following global linear convergence rate for any k ≥ 1,

Moreover, when k ≥ (1

Proof We follow a similar argument as in the proof of Lemma 6 but with a different lower bound for αk . Specifically, by Lemma 3, we also have αk

. Combining this with qk ≥ 2/κ and (39) leads to

To begin with, recall the definition that

Since the objective function is non-increasing by Lemma 1, it holds that C i ≤ C 0 for any i ≥ 0. Thus, from (42) we have

Thus, by using Proposition 1 we obtain (40).

To prove the second claim in (41), we use the fact that 1 + x ≤ e x for any x ∈ R to get

Combining ( 42) and ( 43) leads to

Next, we prove an upper bound on

, which implies that k ≥ Ψ ( B0 ). Then (38) in Lemma 6 and (40) together imply that

where we used the fact that e -Ψ ( B0 ) k ≥ e -1 ≥ 1 3 . Moreover, we decompose the sum

For the first part, we have

For the second part, by the definition of C i , we have

where we used

x for all 0 ≤ x ≤ 1 in the last inequality. Combining both inequalities, we arrive at

Thus, when the number of iterations k exceeds (1

Together with Proposition 1, this proves the second claim in (41).

We summarize all the global linear convergence results from the above two lemmas in the following theorem.

Theorem 1 Let {x k } k≥0 be the iterates generated by the BFGS method with exact line search and suppose that Assumptions 1, 2 and 3 hold. For any initial point x 0 ∈ R d and any initial matrix B 0 ∈ S d ++ , we have the following global linear convergence rate for any k ≥ 1,

where B0 is defined in (37). When k ≥ Ψ ( B0 ), we have that

Moreover, when k ≥ (1

In Theorem 1, we present three distinct linear convergence rates during different phases of the BFGS algorithm with exact line search. Specifically, the linear rate in (46) is applicable from the first iteration, but the contraction factor depends on the quantity e -Ψ ( B0 )/k , which can be exponentially small and thus imply a slow convergence rate. However, this quantity will be bounded away from zero as the number of iterations k increases, resulting in an improved linear rate. In particular, for k ≥ Ψ ( B0 ), the quantity e -Ψ ( B0 )/k is bounded below by 1/3, leading to the second improved linear convergence rate in (47). Furthermore, as shown in Lemma 7, after an additional C 0 Ψ ( B0 ) + 3C 0 κ min{2(1 + C 0 ), 1 + √ κ} iterations, we achieve the last linear convergence rate in (48), which is comparable to that of gradient descent.

From the discussions above, we observe that the quantity Ψ ( B0 ) (recall that B0 = 1 L B 0 ) plays a critical role in determining the transitions between different linear convergence phases, and a smaller Ψ ( B0 ) implies fewer iterations required to reach each linear convergence phase. Thus, we consider a special case to simplify our bounds, where B 0 = α I and α > 0 is an arbitrary positive scalar. Given this simplification, we obtain

We apply this to Theorem 1 to establish the corresponding global linear rates, as stated in Corollary 1.

Corollary 1 Let {x k } k≥0 be the iterates generated by the BFGS method with exact line search and suppose that Assumptions 1, 2 and 3 hold. For any initial point x 0 ∈ R d and the initial Hessian approximation matrix B 0 = α I with α > 0, we have the following global convergence rate for any k ≥ 1,

Moreover, when k ≥ (1

The above corollary characterizes the behavior of BFGS with exact line search when the initial Hessian approximation is a scaled identity matrix. In the following paragraphs, we refine these results by analyzing the bounds for specific values of α.

First, we examine two extreme cases for initialization: α = L (where B 0 = L I ) and α = μ (where B 0 = μI ). The former corresponds to an upper bound on the eigenvalues of the Hessian, while the latter uses a lower bound. Note that in the case where

= 0 and we achieve the global linear convergence rate in (50) for any k ≥ 1. Moreover, when k ≥ 3C 0 κ min{2(1 + C 0 ), (1 + √ κ)}, we reach the second linear convergence rate in (51). In the case where

We have the following global convergence rate for any k ≥ 1,

When k ≥ d log κ, we achieve the global linear convergence rate in (50). Moreover, when k ≥ (1

we reach the second linear convergence rate in (51). Comparing the above results, we observe that BFGS with B 0 = μI requires additional d log κ iterations to achieve a similar linear rate as in the first case. However, as we present in the next section, the choice of the initial Hessian approximation matrix B 0 = μI could achieve a superlinear rate faster. This trade-off between the linear and superlinear convergence phase is the fundamental consequence of different choices of the initial Hessian approximation matrix in our convergence analysis.

While the special cases discussed above are valuable for theoretical comparison, they may not be practical for selecting the initial Hessian approximation, as the constants μ and L are often unknown. A more practical choice, which can be easily computed, is to set B 0 = cI , where c is determined based on gradient and variable differences between two randomly selected points. Specifically, c is given by c = s y s 2 where s = x 2x 1 and y = ∇ f (x 2 ) -∇ f (x 1 ), with x 1 and x 2 being two randomly chosen vectors. This choice ensures that c ∈ [μ, L]. For this initialization of B 0 , we can establish the bound d c L -1 + log L c ≤ d log κ. Applying this upper bound to our linear convergence result in Corollary 1, we obtain that when k ≥ d log κ, we have the linear convergence rate in (50). When k ≥ (1

we have the linear rate in (51).

Global superlinear convergence rates

In this section, we establish the non-asymptotic global superlinear convergence rate of BFGS with exact line search, employing a similar approach to the global linear convergence rate analysis from the previous section. We utilize the framework from Proposition 1 and integrate the lower bounds from Lemmas 3, 4, 5, and Proposition 2.

The key distinction lies in the choice of the weight matrix: instead of P = L I used in the linear convergence analysis, we opt for P = ∇ 2 f (x * ) for the global superlinear convergence proof. We define the weighted matrix Bk as:

In the following proposition, we first provide a general global convergence bound with an arbitrary initial Hessian approximation matrix B 0 ∈ S d ++ . All the global superlinear convergence rates are based on the following proposition.

Proposition 3 Let {x k } k≥0 be the iterates generated by the BFGS method with exact line search and suppose that Assumptions 1, 2 and 3 hold. Recall the definition of C k in (23) and Ψ (•) in (18). For any initial point x 0 ∈ R d and any initial Hessian approximation matrix B 0 ∈ S d ++ , the following result holds for any k ≥ 1,

Proof Recall that we choose the weight matrix as P = ∇ 2 f (x * ) throughout the proof. From Lemma 3 and Lemma 4(b), we have αk ≥

Hence, using the inequality 1 + x ≤ e x for any x ≥ 0, it follows that

Moreover, by using the inequality (20) in Proposition 2 with P = ∇ 2 f (x * ), we obtain that

where in the last inequality we used the fact that ŷi 2

This further implies that

Combining ( 55), (56), and ( 16) from Proposition 1, we prove that

where the last inequality is due to the fact that 1e -x ≤ x for any x.

The above global result shows that the error after k iterations for the BFGS update with exact line search depends on the potential function of the weighted initial Hessian approximation matrix B0 , i.e., Ψ ( B0 ), and the sum of weighted function value optimality gap, i.e., k-1 i=0 C i . This result forms the foundation of our superlinear result, since if we can demonstrate that the sum k-1 i=0 C i is bounded above, it leads to a superlinear rate of the form O((1/k) k ).

Having established the non-asymptotic global linear convergence rate of BFGS in the previous section, we can leverage it to show that the sum k-1 i=0 C i is uniformly bounded above, allowing us to establish an explicit upper bound for this finite sum. In the following theorem, we apply the linear convergence results from Sect. 4 to prove the non-asymptotic global superlinear convergence rates of BFGS with exact line search for any initial Hessian approximation matrix B 0 ∈ S d ++ .

Theorem 2 Let {x k } k≥0 be the iterates generated by the BFGS method with exact line search and suppose that Assumptions 1, 2 and 3 hold. For any initial point x 0 ∈ R d and any initial Hessian approximation matrix B 0 ∈ S d ++ , we have the following superlinear convergence rate,

) where B0 and B0 are defined in (37) and (53).

Proof From (45) in Lemma 7, we know that for k ≥ 1,

Leveraging ( 58) and (54) in Lemma 3, we prove that for k ≥ 1,

and the proof is complete. This result indicates that BFGS with exact line search achieves a superlinear convergence rate when the number of iterations satisfies the condition k ≥ Ψ ( B0

The initial matrix B 0 critically influences the required iterations to attain this rate, as it appears in the numerator of the upper bound through

Thus, different choices of B 0 yield different values for Ψ ( B0 ) + 4C 0 Ψ ( B0 ), affecting the number of iterations required for superlinear convergence. Indeed, one can try to optimize the choice of B 0 to make the expression Ψ ( B0 ) + 4C 0 Ψ ( B0 ) as small as possible.

Now we consider the special case where B 0 = α I with α > 0 as any positive constant. For this case, we have

from Assumptions 1 and 2. Applying these bounds in Theorem 2, we obtain the following corollary.

Corollary 2 Let {x k } k≥0 be the iterates generated by the BFGS method with exact line search and suppose that Assumptions 1, 2 and 3 hold. For any initial point x 0 ∈ R d and the initial Hessian approximation matrix B 0 = α I with α > 0, we have the following superlinear convergence rate,

The above corollary characterizes the non-asymptotic superlinear convergence rate of BFGS with exact line search when the initial Hessian approximation is an identity matrix multiplied by a constant α. Similar to the linear convergence analysis, we present the superlinear convergence rates for specific values of α in the following paragraphs.

Hence, we obtain the superlinear convergence rate

Similarly, when α = μ (B 0 = μI ), we have

This leads to the superlinear convergence rate

As shown in the above two results, choosing B 0 = L I minimizes Ψ ( B0 ), resulting in Ψ ( B0 ) = 0. However, Ψ ( B0 ) in this case could be as large as dκ. On the other hand, setting B 0 = μI yields a more favorable upper bound, ensuring that both Ψ ( B0 ) and Ψ ( B0 ) are bounded by d log κ. Hence, initializing the Hessian approximation with B 0 = μI instead of B 0 = L I could result in fewer iterations to reach the superlinear convergence phase. Generally, during the initial linear convergence stage, the iterates generated by the BFGS method with B 0 = L I outperform those with B 0 = μI , due to a faster linear convergence speed. However, the BFGS method with B 0 = μI transitions to the ultimate superlinear convergence phase in fewer iterations compared to B 0 = L I . This phenomenon has also been observed in our experiments in Sect. 7.

As in the linear convergence analysis, we also consider the practical initial Hessian approximation: B 0 = cI , where c is s y s 2 , with s = x 2 -x 1 , y = ∇ f (x 2 )-∇ f (x 1 ), and x 1 , x 2 as two random vectors. For this choice of B 0 , we can derive the following upper bounds:

Applying these values of Ψ ( B0 ) and Ψ ( B0 ) to our superlinear convergence result in Corollary 2, we can obtain the following convergence guarantees for B 0 = cI :

(60) While all of our presented results are global and do not impose any initial condition on x 0 , in the following remark, we present a potential local result when B 0 = μI . Remark 3 Consider the scenario where BFGS starts at a point x 0 near the optimal solution x * such that the initial error condition

In this case, we can establish that (1

. Thus, when B 0 = μI , we obtain the local superlinear convergence rate of O( d log κ k ) k , which aligns with the local convergence result in [37]. It is noteworthy that the local result in [37] relied on a unit step size, while our local side result is derived using exact line search.

Discussions

Comparison with local non-asymptotic analysis. In this section, we discuss the recent non-asymptotic local convergence results for BFGS and DFP in [25,37,38] and explain why these results cannot be easily extended to achieve global complexity bounds.

To begin with, note that these results are crucially based on local analysis and only apply when the iterates are close to the optimal solution x * and the step size η k is set to 1 in this local region. Therefore, to extend their results into a global convergence guarantee, one plausible strategy is to employ a line search scheme to ensure global convergence, and then switch to the local analysis when the iterates enter the region of local convergence. However, this approach faces several challenges.

First, it remains unclear how to explicitly upper bound the number of iterations until the line search subroutine accepts the unit step size η k = 1. Moreover, assume that the iterates enter the region of local convergence after k 0 iterations and we have η k = 1 for all k ≥ k 0 . Even then, there is no guarantee that the Hessian approximation matrix B k 0 will satisfy the necessary conditions required for the local analysis in [25,37,38]. Specifically, for the analysis in [25] to hold, B k 0 must be sufficiently close to the exact Hessian matrix, which is not satisfied in general. Regarding [37,38], we note that their analyses depend on the condition number of B k 0 , which could be exponentially large and thus render the superlinear rate meaningless. To be more concrete, inspecting the proofs in [37,Lemma 5.4] and [38,Theorem 4.2] reveals that the superlinear convergence rate occurs when k = Ω(Ψ ( B-1 k 0 )) and k = Ω(Ψ ( Bk 0 )), respectively, where

with J k 0 defined in (21) and Ψ (•) is the potential function defined in (18). Consequently, it is essential to establish bounds for the smallest and largest eigenvalues of Bk 0 . However, the current theory indicates (see e.g. [38,Theorem 4.1]) that e -2κ Mλ 0 I Bk 0 e 2κ Mλ 0 I , where λ 0 = (∇ 2 f (x 0 )) -1 2 ∇ f (x 0 ) denotes the initial Newton decrement. This suggests that without a sufficiently small λ 0 , the extreme eigenvalues of Bk 0 will be exponentially dependent on the condition number κ, leading to Ψ ( B-1 k 0 ), Ψ ( Bk 0 ) = Ω(de 2κ Mλ 0 ). Hence, a superlinear rate will be achieved only after Ω(de 2κ Mλ 0 ) iterations.

Our convergence framework also diverges significantly from the previous works [25,37,38] in terms of the proof strategy. Specifically, the approach in the aforementioned studies employs an induction argument to control the largest and smallest eigenvalues of the Hessian approximation matrix B k and prove a local linear convergence rate. In comparison, as presented in Sects. 4 and 5, we prove global linear and superlinear convergence rates without explicitly establishing upper or lower bounds on the eigenvalues of B k . This marks a notable departure from the local convergence analysis in [38], [37], and [25].

Comparison with global asymptotic analysis. As mentioned in Sect. 3, our convergence analysis framework resembles the approach taken in [8,9,35] for proving asymptotic linear convergence rates of classical quasi-Newton methods such as BFGS and DFP. While these works considered inexact line search schemes and thus are different from our exact line search setting, they used a similar inequality as (16) in Proposition 1 to express the convergence rate in terms of the angle θk . Moreover, the authors in [35] and [9] analyzed the traces and the determinants of the Hessian approximation matrices {B k } k≥0 separately to lower bound k-1 i=0 cos ( θi ). Later, this process was simplified in [8] by introducing the potential function Ψ (•) given in (18), combining the trace and determinant together as in our Proposition 2. However, since their main focus is on asymptotic convergence, we note that these previous works only demonstrate that ( k-1 i=0 cos ( θi )) 1/k is lower bounded by a constant, without giving an explicit form. Furthermore, our work builds upon previous analyses by incorporating a weight matrix P, while earlier works correspond to setting P = I . Another notable difference is that we keep the term mk and lower bound the term cos 2 ( θk )/ mk as shown in Proposition 2, whereas previous works relied on a looser bound for mk . These refinements enable us to provide a tighter linear convergence rate for the BFGS method.

On the other hand, in demonstrating superlinear convergence, our approach deviates significantly from that of [8,9,35]. Specifically, the previous works relied on the Dennis-Moré condition, i.e., lim k→∞

= 0, to establish asymptotic superlinear convergence. In comparison, we use the same framework outlined in Sect. 3 to establish both linear and superlinear convergence rates. The key distinction lies in the choice of the weight matrix P: we choose P = L I for showing linear convergence and P = ∇ 2 f (x * ) for showing superlinear convergence. Thus, we provide a unified framework for studying the global non-asymptotic convergence of BFGS.

Numerical experiments

In this section, we present our numerical experiments to corroborate our convergence rate guarantees, and in particular, we explore the difference between the convergence paths of BFGS under different initializations of B 0 . We further compare these variants of BFGS implementations with the gradient descent algorithm when deployed with exact line search. In our numerical experiments, all the step sizes used in BFGS with different B 0 and gradient descent are computed by the exact line search condition defined in (9). Specifically, we use MATLAB's "fminsearch" function from its optimization toolbox to determine the exact line search step size for all algorithms. In our experiments, all initial points are chosen as random vectors in the corresponding Euclidean vector spaces.

We focus on a hard cubic objective function defined in [44, Sect. 5], i.e.,

and g : R → R is defined as

Fig. 1 Convergence rates of BFGS with different B 0 and gradient descent for solving the hard cubic objective function when condition number and dimension is varied where α, β, λ, Δ ∈ R are hyper-parameters and {v i } n i=1 are standard orthogonal unit vectors in R d . This hard cubic function is used to establish a lower bound for secondorder methods.

In Fig. 1, we compare the gradient descent method and BFGS with different initialization of B 0 : B 0 = L I , B 0 = μI , B 0 = 10L I , B 0 = 0.1μI , B 0 = √ LμI and B 0 = cI where c = s y s 2 . Here, s = x 2x 1 , y = ∇ f (x 2 ) -∇ f (x 1 ), and x 1 , x 2 as two randomly selected vectors. Note that c ∈ [μ, L] and the choice of B 0 = cI is the most commonly used initial Hessian approximation matrix in practice [33]. In (a), (b), and (c) of Fig. 1, we vary the problem's dimension while keeping the condition number as 1,000. Conversely, in (d), (e), and (f), we fix the problem's dimension as 600 and vary the condition number.

Several observations are in order.

-BFGS with B 0 = L I initially converges faster than BFGS with B 0 = μI in most plots, aligning with our theoretical findings that the linear convergence rate of BFGS with B 0 = L I surpasses that of B 0 = μI . -The transition to superlinear convergence for BFGS with B 0 = μI typically occurs around k ≈ d, as predicted by our theoretical analysis. Interestingly, this transition does not always coincide with the iterates approaching the solution's local neighborhood; in many cases, it occurs for BFGS with B 0 = μI even when its error is larger than that of gradient descent. -Although BFGS with B 0 = L I initially converges faster, its transition to superlinear convergence consistently occurs later than for B 0 = μI . Notably, for a fixed dimension d = 600, the transition to superlinear convergence for B 0 = L I occurs increasingly later as the problem condition number rises, an effect not observed for B 0 = μI . This phenomenon indicates that the superlinear rate for B 0 = L I is Additionally, to analyze the sensitivity of BFGS to different line-search schemes, we compare its performance when B 0 = cI under three distinct line-search strategies, as shown in Fig. 2.

The first approach is the Exact Line Search, which is the primary focus of our paper. It is implemented using MATLAB's "fminsearch" function.

The second approach is the Inexact Line Search, where the step size is determined by enforcing the well-known strong Wolfe conditions:

|∇ f (x t + η t d t ) d t | ≤ β|∇ f (x t ) d t |. (64) ŝ k Bk ŝk + ŷk 2 ŝ k ŷk , (65)

Taking the trace on both sides of the equation in (12) and using the fact that Tr(ab ) = a b for any vectors a and b, we obtain the equality in (65). For the proof of (66), we refer the reader to [38,Lemma 6.2] . Take the logarithm on both sides of the above equation, we obtain that log ŝ k ŷk ŝ k Bk ŝk = log Det( Bk+1 )log Det( Bk ).

Recall that mk = where the last inequality holds since xlog x + 1 ≥ 0 for any x > 0. Hence, (19) follows from the above inequality. Finally, the result in (20) follows from summing both sides of (19) where the last inequality holds since Ψ ( Bk ) ≥ 0 for any k ≥ 0.

B Proof of Lemma 2

(a) Recall that J k = 1 0 ∇ 2 f (x k + τ (x k+1x k ))dτ . Using the triangle inequality, we have

Moreover, it follows from Assumption 3 that ∇ 2 f (x * ) -∇ 2 f (x k + τ (x k+1x k )) ≤ M (1τ )(x *x k ) + τ (x *x k+1 ) for any τ ∈ [0, 1]. Thus, we can further apply the triangle inequality to obtain Since f is strongly convex, by Assumption 1 and f (x k+1 ) ≤ f (x k ), we have μ 2 x kx * 2 ≤ f (x k )f (x * ), which implies that x kx * ≤ √ 2( f (x k )f (x * ))/μ. Similarly, since f (x k+1 ) ≤ f (x k ), it also holds that

Moreover, notice that by Assumption 1, we also have J k μI and ∇ 2 f (x * ) μI . Hence, (67) implies that where we used the definition of C k in (23). By rearranging the terms, we obtain (24). (b) Recall that G k = 1 0 ∇ 2 f (x k + τ (x *x k ))dτ . Similar to the arguments in (a), we have

Moreover, notice that by Assumption 1 we also have G k μI and ∇ 2 f (x * ) μI . The rest follows similarly as in the proof of (a) and we prove (25). (c) For any τ ∈ [0, 1], we have

Together with (67), it follows from the triangle inequality that

Moreover, notice that by Assumption 1, we also have ∇ 2 f (x k + τ (x k+1x k )) μI and J k μI . The rest follows similarly as in the proof of (a) and we prove (26).

(d) For any τ ∈ [0, 1], we have that

Together with (68), it follows form the triangle inequality that

Moreover, notice that by Assumption 1, we also have ∇ 2 f (x k + τ (x *x k )) μI and G k μI . The rest follows similarly as in the proof of (a) and we prove (27).

C Approximate Exact Line Search Algorithm

The approximation exact line search is implemented using the bisection Algorithm 1 with as the approximation error. The key idea is to select η such that ∇ f (x t + ηd t ) d t ≈ 0, since the exact line search step size η exact satisfies the condition ∇ f (x t + η exact d t ) d t = 0. We begin with an initial step size of η = 1 and iteratively double it until ∇ f (x t + ηd t ) d t > 0. Once this condition is met, we apply the bisection algorithm, leveraging the sign of ∇ f (x t +ηd t ) d t to refine η. The bisection algorithm is well-suited for this task because the function h(η) = ∇ f (x t + ηd t ) d t is strictly increasing. This follows from the strong convexity of the objective function f , which ensures that h (η) = d t ∇ 2 f (x t + ηd t )d t > 0. Additionally, we note that h(0) = ∇ f (x t ) d t = -g t B -1 t g t < 0, since B t is symmetric positive definite, and that h(η exact ) = 0 with η exact > 0. In our experiments, we set the required accuracy for this scheme to be = 10 -8 and we observe that on average after 15 iterations the bisection method converges.