Exploiting Local Convergence of Quasi-Newton Methods Globally: Adaptive Sample Size Approach

Abstract:
In this paper, we study the application of quasi-Newton methods for solving empirical risk minimization (ERM) problems defined over a large dataset. Traditional deterministic and stochastic quasi-Newton methods can be executed to solve such problems; however, it is known that their global convergence rate may not be better than first-order methods, and their local superlinear convergence only appears towards the end of the learning process. In this paper, we use an adaptive sample size scheme that exploits the superlinear convergence of quasi-Newton methods globally and throughout the entire learning process. The main idea of the proposed adaptive sample size algorithms is to start with a small subset of data points and solve their corresponding ERM problem within its statistical accuracy, and then enlarge the sample size geometrically and use the optimal solution of the problem corresponding to the smaller set as an initial point for solving the subsequent ERM problem with more samples. We show that if the initial sample size is sufficiently large and we use quasi-Newton methods to solve each subproblem, the subproblems can be solved superlinearly fast (after at most three iterations), as we guarantee that the iterates always stay within a neighborhood that quasi-Newton methods converge superlinearly. Numerical experiments on various datasets confirm our theoretical results and demonstrate the computational advantages of our method.

text: In this paper, we study the application of quasi-Newton methods for solving empirical risk minimization (ERM) problems defined over a large dataset. Traditional deterministic and stochastic quasi-Newton methods can be executed to solve such problems; however, it is known that their global convergence rate may not be better than first-order methods, and their local superlinear convergence only appears towards the end of the learning process. In this paper, we use an adaptive sample size scheme that exploits the superlinear convergence of quasi-Newton methods globally and throughout the entire learning process. The main idea of the proposed adaptive sample size algorithms is to start with a small subset of data points and solve their corresponding ERM problem within its statistical accuracy, and then enlarge the sample size geometrically and use the optimal solution of the problem corresponding to the smaller set as an initial point for solving the subsequent ERM problem with more samples. We show that if the initial sample size is sufficiently large and we use quasi-Newton methods to solve each subproblem, the subproblems can be solved superlinearly fast (after at most three iterations), as we guarantee that the iterates always stay within a neighborhood that quasi-Newton methods converge superlinearly. Numerical experiments on various datasets confirm our theoretical results and demonstrate the computational advantages of our method.

heading: 1 Introduction
text: In various machine learning problems, we aim to find an optimal model that minimizes the expected loss with respect to the probability distribution that generates data points, which is often referred to as expected risk minimization. In realistic settings, the underlying probability distribution of data is unknown and we only have access to a finite number of samples (N samples) that are drawn independently according to the data distribution. As a result, we often settle for minimizing the sample average loss, which is also known as empirical risk minimization (ERM). The gap between expected risk and empirical risk which is known as generalization error is deeply studied in the statistical learning literature [1–5], and it is well-known that the gap between these two loss functions vanishes as the number of samples N in the ERM problem becomes larger.
Several variance reduced first-order methods have been developed to efficiently solve the ERM problem associated with a large dataset [6–13]. However, a common shortcoming of these first-order methods is that their performance heavily depends on the problem condition number κ := L/µ, where L is the risk function smoothness parameter and µ is the risk function strong convexity constant. As a result, these first-order methods suffer from slow convergence when the problem has a large condition number, which is often the case in datasets with thousands of features. In fact, the best-known bound
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
achieved by first-order methods belongs to Katyusha [11] which has a cost of O((N + √ κN) logN) to achieve the statistical accuracy of an ERM problem with N samples.
To resolve this issue, second-order methods such as Newton’s method and quasi-Newton methods are often used, as they improve curvature approximation by exploiting second-order information. Several efficient implementations of second-order methods for solving large-scale ERM problems have been proposed recently, including incremental Newton-type methods [14], sub-sampled Newton algorithms [15], and stochastic quasi-Newton methods [16–20]. However, these methods still face two major issues. First, their global convergence requires selection of a small step size which slows down the convergence, or implementation of a line-search scheme which is computationally prohibitive for large-scale ERM problems, since it requires multiple passes over data. The second disadvantage of these (deterministic and stochastic) second-order methods is that their local quadratic or superlinear convergence only appears in a small local neighborhood of the optimal solution. This often means that, when we solve an ERM problem, the fast superlinear or quadratic convergence rate of these methods happens when the iterates have already reached the statistical accuracy of the problem and the test accuracy would not improve, even if the ERM problem error (training error) decreases.
To address these two shortcomings of second-order methods for solving large-scale ERM problems, the adaptive sample size Newton (Ada Newton) method was proposed in [21]. In adaptive sample size methods, instead of solving the ERM corresponding to the full training set directly (as in deterministic methods) or using a subset of samples at each iteration (as in stochastic methods), we minimize a sequence of geometrically increasing ERM problems. Specifically, we solve the ERM problem corresponding to a subset of samples, e.g. m samples, within their statistical accuracy, then we use the resulting solution as an initial point for the next ERM problem with n = αm samples (α > 1), where the larger set with n samples contains the smaller set with m samples. The main principle of this idea is that all training samples are drawn from the same distribution, therefore the optimal values of the ERM problems with m and n samples are close to each other. The main result of Ada Newton in [21] is that if one chooses the expansion factor α properly, the solution for the problem with m samples is within the Newton quadratic convergence neighborhood of the next ERM problem with n = αm samples. Hence, each subproblem can be solved quickly by exploiting the quadratic convergence of Newton’s method at every step of the learning process. This approach is promising as the overall number of gradient and Hessian computations for Ada Newton is O(N), while it requires computing O(logN) Newton’s directions. Note that O notation only hides absolute constants and these complexities are independent of problem condition number κ. The main drawback of this approach, however, is the requirement of O(logN) Newton’s directions computation, where each could cost O(d3) arithmetic operations, where d is the dimension. It also requires O(N) Hessian evaluations which requires O(Nd2) arithmetic operations, and this could be computationally costly. A natural approach to reduce the computational cost of Newton’s method is the use of quasi-Newton algorithms, in which, instead of exact computation of the Hessian and its inverse required for Newton update, we approximate the objective function Hessian and its inverse only by computing first-order information (gradients). Since this approximation only requires matrix-vector product computations at each iteration, the computational complexity of quasi-Newton methods is reduced. It is known that quasi-Newton methods such as DFP and BFGS converge superlinearly in a neighborhood close to the optimal solution, i.e., the iterates {wk} converge to the optimal solution w∗ at a rate of limk→+∞ ‖wk+1−w∗‖ ‖wk−w∗‖ = 0, where k is the iteration number. To exploit quasi-Newton methods for the described adaptive sample size scheme, explicit superlinear convergence rate of these methods is needed to characterize the number of iterations required for solving each ERM problem. The non-asymptotic convergence rate of quasi-Newton methods was unknown until very recently, when a set of concurrent papers [22, 23] and [24] showed that in a local neighborhood of the optimal solution, the iterates of DFP and BFGS converge to the optimal solution at a rate of (1/k)k/2.
In this paper, we exploit the finite-time analysis of the quasi-Newton methods studied in [24] to develop an adaptive sample size quasi-Newton (AdaQN) method that reaches the statistical accuracy
of the ERM corresponding to the full training set after O(N) gradient computations and O(logN) matrix-vector product evaluations and only 1 matrix inversion. In Table 1, we formally compare the computation cost of our proposed AdaQN algorithm and the Ada Newton method in [21]. As shown in Table 1, the total number of gradient evaluations and matrix-vector product computations for AdaQN are the same as the ones for Ada Newton, while the number of Hessian computations for AdaQN and Ada Newton arem0 andO(N), respectively, wherem0 is the size of initial training set in the adaptive sample size scheme, and it is substantially smaller than N , i.e., m0 N . In our main theorem we show that m0 is lower bounded by Ω ( max { d, κ2s log d }) , where s is a parameter determined by the loss function and the structures of the training data (see details in Remark 2). Hence, in this paper, we focus on the regime that the total number of samples satisfies N max { d, κ2s log d } . However, our numerical experiments in Section 5 show that AdaQN often achieves small risk for m0 much smaller than this lower bound and the effectiveness of AdaQN goes beyond this regime. Also, the number of times that Ada Newton requires to compute the inverse of a square matrix (or solving a linear system) is O(logN), while to implement AdaQN we only require one matrix inversion. We should also add that the implementation of Ada Newton requires a backtracking mechanism for the choice of growth factor α. This backtracking scheme is needed to ensure that the size of training set does not grow too fast, and the iterates stay within the quadratic convergence neighborhood Newton’s method for next ERM problem. Unlike Ada Newton, our proposed AdaQN method does not require such backtracking scheme, and as long as the size of initial training set m0 is sufficiently large, we can double the size of training set at the end of each phase, i.e., we set α = 2.
It is worth noting that in [25], the authors study a similar approach to our proposed adaptive sample size quasi-Newton method called the batch-expansion training (BET) method. The BET method differs from our proposed framework in terms of algorithm development and convergence analysis. Specifically, the BET method uses a similar idea of geometrically increasing the size of the training set, while using the limited-memory BFGS (L-BFGS) method for solving the Empirical Risk Minimization (ERM) subproblems. The authors leverage the linear convergence rate of L-BFGS to characterize the overall complexity of the BET algorithm. We should mention that since the linear convergence rate of L-BFGS is not provably better than first-order methods, the overall iteration complexity of BET is similar to the one for adaptive sample size first-order methods [26] which depends on the problem condition number. In contrast to [25], in this paper, we leverage the local superlinear convergence rate of the BFGS method to improve the overall complexity of adaptive sample size first-order methods, and hence, it also improves the complexity of the BET method.

heading: 2 Problem formulation
text: In this paper, we focus on minimizing large-scale empirical risk minimization (ERM) problems. Consider the loss function f : Rd ×Z → R that takes inputs as the decision variable w ∈ Rd and random variable Z ∈ Z with a probability distribution P . The expected risk minimization problem aims at minimizing the expected loss over all possible choices of Z,
w∗ := arg min w R(w) = arg min w EZ∼P [f(w;Z)]. (1)
Note that in this expression we defined the expected risk R : Rd → R as R(w) := EZ∼P [f(w;Z)] and w∗ as an optimal solution of the expected risk minimization. In this paper we focus on loss functions f that are strongly convex with respect to w, hence the optimal solution w∗ is unique.
The probability distribution P is often unknown, and we only have access to a finite number of realizations of the random variable Z, which we refer to as our training set T = {z1, . . . , zN}. These N realizations are assumed to be drawn independently and according to P . Hence, instead of minimizing the expected risk in (1), we attempt to minimize the empirical risk corresponding to the training set T = {z1, . . . , zN}. To formalize this, let us first define Sn = {z1, . . . , zn} as a subset of the training set T which contains its first n elements. Without loss of generality we defined an ordering for the elements of the training set T . The empirical risk corresponding to the set Sn is defined as Rn(w) := 1n ∑n i=1 f(w; zi) and its optimal solution w ∗ n is given by
w∗n := arg min w Rn(w) = arg min w
1
n n∑ i=1 f(w; zi). (2)
Note that the ERM problem associated with the full training set T is a special case of (2) when n = N , and its unique optimal solution is denoted by w∗N .
The gap between the optimal values of empirical risk Rn and expected risk R is well-studied in the statistical learning literature and here we assume the following upper bound holds
E [|Rn(w∗n)−R(w∗)|] ≤ Vn, (3)
where Vn is a function of the sample size n that approaches zero as n becomes large. The expectation in (3) is over the set Sn. In this paper, we assume that all the expectations are with respect to the corresponding training set. Classic results have established bounds of the form Vn = O(1/ √ n) [1, 27] and other recent papers including [2–5] show that under stronger regularity conditions such as strong convexity, we have Vn = O(1/n), for n = Ω(κ2 log d). In this paper, we assume the requirements for the bound Vn = O(1/n) are satisfied. Since the gap between the optimal empirical and expected risks is always bounded above by the error Vn, there is no point to reducing the optimization error of minimizing Rn beyond the statistical accuracy Vn. In other words, if we obtain a solution ŵ such that E [Rn(ŵ)−Rn(w∗n)] = O(Vn), there would be no benefit in further minimizing Rn. Due to this observation, we say ŵ solves the ERM in (2) within its statistical accuracy if E [Rn(ŵ)−Rn(w∗n)] ≤ Vn. The ultimate goal is to efficiently find a solution wN that reaches the statistical accuracy of the full training set T , i.e., E [RN (wN )−RN (w∗N )] ≤ VN . Next, we state the notations and assumptions. Assumption 1. For all values of z, the loss function f(w; z) is twice differentiable and µ-strongly convex with respect to w and its gradient is Lipschitz continuous with parameter L > 0.
Assumption 2. For all values of z, the loss function f(w; z) is self-concordant with respect to w.
Assumptions 1 and 2 are customary for the analysis of quasi-Newton methods. These conditions also imply that the empirical risk Rn is also self-concordant, strongly convex with µ, and its gradient is Lipschitz continuous with L. Hence, the condition number of Rn is κ := L/µ.
Computational cost notation We report the overall computational cost in terms of these parameters: (i) τgrad and τHess which denote the cost of computing one gradient of size d and one Hessian of size d× d; (ii) τprod which indicates the cost of computing the product of a square matrix of size d × d with a vector of size d; and (iii) τinv which denotes the cost of computing the inverse of a square matrix of size d× d or the cost of solving a linear system with d variables and d equations.

heading: 3 Algorithm
text: Quasi-Newton methods Before introducing our method, we first briefly recap the update of quasiNewton (QN) methods. Given the current iterate w, the QN update for the ERM problem in (2) is
w+ = w − ηH ∇Rn(w), (4)
where η > 0 is a step size and H ∈ Rd×d is a symmetric positive definite matrix approximating the Hessian inverse∇2Rn(w)−1. The main goal of quasi-Newton schemes is to ensure that the matrix H always stays close to ∇2Rn(w)−1. There are several different approaches for updating the Hessian approximation matrix H, but the two most-widely used updates are the DFP method defined as
H+ = H− Hyy >H
y>Hy +
ss> s>y , (5)
and the BFGS update defined as
H+ = ( I− sy >
s>y
) H ( I− ys >
s>y
) + ss>
s>y , (6)
where s := w+ − w is the variable variation and y := ∇Rn(w+) − ∇Rn(w) is the gradient variation. If we follow the updates in (5) or (6), then finding the new Hessian inverse approximation and consequently the new descent direction −H+∇Rn(w+) only require computing a few matrixvector multiplications. Considering this point and the fact that each step of BFGS or DFP requires n gradients evaluations, the computational cost of each step of BFGS and DFP is O(nτgrad + τprod).
Algorithm 1 AdaQN
Input: The initial sample size m0; The initial argument wm0 within the statistical accuracy of Rm0 ; The initial Hessian inverse Hm0 = ∇2Rm0(wm0)
−1; 1: Set n← m0; 2: while n ≤ N do 3: Set the initial argument ŵ← wn, the initial matrix Ĥ← Hm0 and iteration = 1; 4: Increase the sample size: n← min{2n,N}; 5: while iteration ≤ tn do 6: ŵ+ ← ŵ − Ĥ ∇Rn(ŵ); 7: s← ŵ+ − ŵ; 8: y← ∇Rn(ŵ+)−∇Rn(ŵ); 9: Ĥ+ = ( I− sy >
s>y
) Ĥ ( I− ys >
s>y
) + ss >
s>y ;
10: Set ŵ← ŵ+, Ĥ← Ĥ+ and iteration← iteration+ 1; 11: end while 12: Set wn ← ŵ; 13: end while

heading: 3.1 Adaptive sample size quasi-Newton algorithm (AdaQN)
text: BFGS, DFP or other quasi-Newton (QN) methods can be used to solve the ERM problem corresponding to the training set T withN samples, but (i) the cost per iteration would be ofO(Nτgrad+τprod), (ii) the superlinear convergence only appears towards the end of learning process when the iterates approach the optimal solution and statistical accuracy is already achieved, and (iii) they require a step size selection policy for their global convergence. To resolve the first issue and reduce the high computational cost of O(Nτgrad + τprod), one could use stochastic or incremental QN methods that only use a subset of samples at each iteration; however, stochastic or incremental QN methods (similar to deterministic QN algorithms) outperform first-order methods only when the iterates are close to the optimal solution and they also require line-search schemes for selecting the step size. Hence, none of these schemes is able to exploit fast superlinear convergence rate of QN methods throughout the entire training process, while always using a constant step size of 1.
Next, we introduce the adaptive sample size quasi-Newton (AdaQN) method that addresses these drawbacks. Specifically, (i) AdaQN does not require any line-search scheme and (ii) exploits the superlinear rate of QN methods throughout the entire training process. In a nutshell, AdaQN leverages the interplay between the statistical accuracy of ERM problems and superlinear convergence neighborhood of QN methods to solve ERM problems efficiently. It uses the fact that all samples are drawn from the same distribution, therefore the solution for ERM with m samples is not far from the solution of ERM with n samples, where n samples contain the m samples. Hence, the solution of the problem with less samples can be used as a warm-start for the problem with more samples. Note that there are two important points here that we should highlight. First, if m and n are chosen properly, the solution of the ERM problem Rm corresponding to the set Sm with m samples will be in the superlinear convergence neighborhood of the next ERM problem corresponding to the set Sn with n samples, where Sm ⊂ Sn. Hence, each subproblem can be solved with a few iterations of quasi-Newton methods (see Figure 5 in the appendix). Second, in each subproblem we only use a subset of samples, therefore the cost of running QN methods for solving subproblems with m samples (where m N ) is significantly less than the update of QN methods for the full training set. The steps of AdaQN are outlined in Algorithm 1. We start with a small subset of the full training set with m0 samples and solve its corresponding ERM problem within its statistical accuracy. The initial ERM problem can be solved using any iterative method and its cost will be negligible as it scales with m0 instead of N , where m0 N . In the main loop (Step 2-Step 12) we implement AdaQN. Specifically, we first use the solution from the previous round (or wm0 when n = m0) as the initial iterate, while we set the initial Hessian inverse approximation as Ĥ = Hm0 = ∇2Rm0(wm0) −1 (Step 3). Then, we double the size of the training set by adding more samples to the active training set (Step 4). In Steps 5-10 we run the BFGS update for minimizing the loss Rn, while we keep updating the iterates and Hessian inverse approximation. Once, the required condition for convergence specified in Step 5 is obtained, we output the iterate wn as the iterate that minimizes Rn within its
statistical accuracy. In Step 8 we used the update for BFGS, but one can simply replace this step with the update of the DFP method. As we ensure that iterates always stay within the neighborhood that BFGS (or DFP) converges superlinearly, the step size of these methods can be set as η = 1. We repeat this procedure until we reach the whole dataset T with N samples. In Algorithm 1, for the phase that we minimize Rn (Steps 5-10), our goal is to find a solution wn that satisfies E [Rn(wn)−Rn(w∗n)] ≤ Vn. We use tn to indicate the maximum number of QN updates to find such a solution. Our theoretical result (Theorem 1) suggests that due to fast convergence of QN methods, at most tn = 3 iterations required to solve each subproblem within its statistical accuracy. Remark 1. Note that a natural choice for the initial Hessian inverse approximation at each phase of AdaQN with n samples is the Hessian inverse at the initial iterate of that phase, i.e., setting the initial Hessian inverse approximation as Ĥ = ∇2Rn(wm)−1, where wm is the solution of the previous phase with m = n/2 samples. However, if we follow this initialization, then for each phase of AdaQN we need to compute n new Hessians to evaluate ∇2Rn(wm) and one matrix inversion to find ∇2Rn(wm)−1, which would increase the computational cost of AdaQN. To avoid this issue, for all values of n, we always use the initial Hessian inverse approximation matrix corresponding to the first phase with m0 samples and set Ĥ = Hm0 = ∇2Rm0(wm0)
−1 (Step 3). We show that even under this initialization the required condition for superlinear convergence of DFP or BFGS is always satisfied if the initial sample size m0 is large enough. Note that by following this scheme, we only need to compute m0 Hessians and a single matrix inversion to implement AdaQN. This is a significant gain compared to Ada Newton in [21], which requires computing O(N) Hessians and O(logN) matrix inversions.

heading: 4 Convergence analysis
text: In this section, we characterize the overall complexity of AdaQN. We only state the results for BFGS defined in (6), as the proof techniques and overall complexity bounds for adaptive sample size versions of DFP and BFGS are similar.
First we state an upper bound for the sub-optimality of the variable wm with respect to the empirical risk of Rn, given that wm has achieved the statistical accuracy of the previous empirical risk Rm. Proposition 1. Consider Sm and Sn such that Sm ⊂ Sn ⊂ T , where there are m samples in Sm and n samples in Sn and n ≥ m. Consider the corresponding empirical risk functions Rm and Rn defined based on Sm and Sn, respectively. Assume that wm solves the ERM problem of Rm within its statistical accuracy, i.e., E [Rm(wm)−Rm(w∗m)] ≤ Vm. Then we have
E [Rn(wm)−Rn(w∗n)] ≤ 3Vm. (7)
Proposition 1 characterizes the sub-optimality of the variable wm for the empirical risk Rn which plays a fundamental role in our analysis. We use this bound later to show wm is close enough to w∗n so that wm is in the superlinear convergence neighborhood of Rn.
Next, we require a bound on the number of iterations needed by BFGS to solve a subproblem, when the initial iterate is within the superlinear convergence neighborhood. We establish this bound by leveraging the result of [24] which provides a non-asymptotic superlinear convergence for BFGS. Proposition 2. Consider AdaQN in the phase that the active training set contains n samples. If Assumptions 1-2 hold and the initial iterate wm and Hessian approximation∇2Rm0(wm0) satisfy
‖∇2Rn(w∗n) 1 2 (wm −w∗n)‖ ≤
1
300 ,
‖∇2Rn(w∗n)− 1 2 [∇2Rm0(wm0)−∇2Rn(w∗n)]∇2Rn(w∗n)− 1 2 ‖F ≤
1 7 ,
(8)
then after tn iterations we achieve the output wn with the following convergence result Rn(wn)−Rn(w∗n) ≤ 1.1 ( 1
tn
)tn [Rn(wm)−Rn(w∗n)]. (9)
Proposition 2 shows that if the initial Hessian approximation error is small and the initial iterate is close to the optimal solution, then the iterates of BFGS with step size η = 1 converge to the optimal
solution at a superlinear rate of (1/k)k after k iterations. The inequalities in (8) identify the required conditions to ensure that wm is within the superlinear convergence neighborhood of BFGS for Rn.
In the next two propositions we show that if the initial sample size m0 is sufficiently large, both wm and∇2Rm0(wm0) satisfy the conditions in (8) in expectation which are required to achieve the local superlinear convergence.
Proposition 3. Consider Algorithm 1 for the phase that the active training set contains n samples with empirical risk Rn. Suppose Assumptions 1-2 hold, and further suppose we are given wm which is within the statistical accuracy of Rm, i.e, E [Rm(wm)−Rm(w∗m)] ≤ Vm, and n = 2m. If the sample size m is lower bounded by m = Ω ( κ2 log d ) , then
E [ ‖∇2Rn(w∗n) 1 2 (wm −w∗n)‖ ] ≤ 1
300 . (10)
Proposition 3 shows that if the training set size is sufficiently large, then the solution of the previous phase is within the BFGS superlinear convergence neighborhood of the current problem, and the first condition in (8) holds in expectation. This condition is indeed satisfied throughout the entire learning process, if it holds for the first training set with m0 samples, i.e., m0 = Ω ( κ2 log d ) .
Next we establish under what condition the Hessian approximation used in adaptive sample size method which is the Hessian evaluated with respect to m0 samples, i.e., ∇2Rm0(wm0), satisfies the second condition in (8) in expectation which is required for the superlinear convergence of BFGS.
Proposition 4. Consider AdaQN in the phase that the active training set contains n samples. If Assumptions 1-2 hold and the initial sample size m0 satisfies m0 = Ω ( max { d, κ2s log d }) , where
s is defined as s := supw,n ( E[‖∇2Rn(w)−∇2R(w)‖F ] E[‖∇2Rn(w)−∇2R(w)‖] )2 . Then for any n ≥ 2m0 we have
E [ ‖∇2Rn(w∗n)− 1 2 [∇2Rm0(wm0)−∇2Rn(w∗n)]∇2Rn(w∗n)− 1 2 ‖F ] ≤ 1
7 . (11)
Based on Proposition 4, if the size of initial training set satisfies m0 = Ω ( max { d, κ2s log d }) , then the second condition in (8) holds in expectation. By combining the results of Propositions 3 and 4, we obtain the required conditions for m0. Moreover, Proposition 2 quantifies the maximum number of iterations required to solve each ERM problem to its corresponding statistical accuracy. In the following theorem, we exploit these results to characterize the overall computational complexity of AdaQN to reach the statistical accuracy of the full training set with N samples.
Theorem 1. Consider AdaQN described in Algorithm 1 for the case that we use BFGS updates in (6). Suppose Assumptions 1-2 hold, and the initial sample size m0 is lower bounded by
m0 = Ω ( max { d, κ2s log d }) , (12)
where s = supw,n ( E[‖∇2Rn(w)−∇2R(w)‖F ] E[‖∇2Rn(w)−∇2R(w)‖] )2 . Then at the stage with n samples, AdaQN finds wn
within the statistical accuracy of Rn after at most tn = 3 iterations. Further, the computational cost of AdaQN to reach the statistical accuracy of the full training set T is
τinv +m0τHess + 6Nτgrad + 3 (1 + log(N/m0)) τprod. (13)
Theorem 1 states that if the initial sample size is sufficiently large, the number of required BFGS updates for solving each subproblem is at most 3 iterations. Further, it characterizes the overall computational cost of AdaQN. Remark 2. Since ‖ · ‖ ≤ ‖ · ‖F ≤ √ d‖ · ‖, parameter s in Theorem 1 belongs to the interval
[1, d]. Hence, in the worst case s = d and m0 = Ω ( κ2d log d ) . However, for many common classes of problems including linear regression and logistic regression and for training datasets with specific structures s could be O(1). In those cases the initial sample size is lower bounded by Ω ( max { d, κ2 log d }) ; see Section C in the appendix for details. In fact, our numerical experiments also verify this observation and show that the choice ofm0 could be much smaller than the worst-case bound of m0 = Ω ( κ2d log d ) ; see Section 5.
Remark 3. In the above complexity bound, we neglect the cost of finding the initial solution wm0 , as the cost of finding an approximate solution for the first ERM problem with m0 samples is negligible compared to the overall cost of the AdaQN. For instance, if one solves the first ERM problem using the Katyusha algorithm [11], the cost of the initial stage would be O ( τgrad(m0 + √ m0κm0) logm0
) (where κm0 is the condition number for ERM with m0 samples), which is indeed dominated by the term O (Nτgrad), when we are in the regime that N m0. Remark 4. If one uses Katyusha, which is an optimal first-order method, to solve an ERM with N samples, the overall gradient computation to achieve accuracy = VN would be O((N +√ Nκ) logN). The gradient computation cost of AdaQN considering the initialization step is O(((m0 + √ m0κ) logm0) + N). Indeed, if we are in the regime that the size of initial set m0 is sufficiently smaller than the size of full training set N , (i.e., N m0 logm0), then AdaQN gradient complexity scales as O(N +√m0κ logm0) which is smaller than O(N logN + √ Nκ logN) cost of Katyusha. However, AdaQN requires evaluating one additional matrix inversion, m0 Hessian evaluations and O(log (N/m0)) matrix-vector product computations.

heading: 5 Numerical experiments
text: Next, we practically evaluate the performance of AdaQN to solve large-scale ERM problems. We consider a binary classification problem with l2 regularized logistic regression loss function, where µ > 0 is the regularization parameter. The logistic loss is convex, and, therefore, all functions in our experiments are µ-strongly convex. We normalize all data points with the unit norm so that the loss function gradient is Lipschitz continuous with L ≤ µ+ 1, hence condition number is κ ≤ 1 + 1/µ. Moreover, the logistic regression loss function is self-concordant. Thus, Assumptions 1-2 hold.
We compare AdaQN with adaptive sample size Newton method (Ada Newton) [21], standard BFGS quasi-Newton method, the L-BFGS quasi-Newton method [28], the stochastic quasi-Newton (SQN) method proposed in [29], and three stochastic first-order methods, including stochastic gradient descent (SGD), the Katyusha algorithm [11], and SAGA which is a variance reduced method [7]. In our experiments, we start with the initial point w0 = c ∗ ~1 where c > 0 is a tuned parameter and ~1 ∈ Rd is the one vector. First we conduct several iterations of the gradient descent method on the initial sub-problem with m0 samples until the initial condition ‖∇Rm0(wm0)‖ ≤ √ 2µVm0 is satisfied. Note that Rm0(wm0) − Rm0(w∗m0) ≤ 1 2µ‖∇Rm0(wm0)‖
2 implies that this initial condition guarantees wm0 is within the statistical accuracy of Rm0 . We use wm0 as the initial point of AdaQN and Ada Newton as presented in Algorithm 1. We use w0 = c ∗ ~1 as the initial point of all other algorithms. We include the cost of finding the proper initialization for AdaQN in our comparisons. We compare these methods over (i) MNIST dataset of handwritten digits [30], (ii) Epsilon dataset from PASCAL challenge 2008 [31], (iii) GISETTE handwritten digit classification dataset from the NIPS 2003 feature selection challenge [32] and (iv) Orange dataset of customer relationship management from KDD Cup 2009 [33].1 More details provided in Table 2.
In our experiments, we observe that even when the initial sample size m0 is smaller than the threshold in (12), AdaQN performs well and converges superlinearly in each subproblem. This is because (12) is a sufficient condition to guarantee our theoretical superlinear rate, and in practice smaller choices of m0 also work. For Ada Newton we use the same scheme descried in Algorithm 1 and replace the QN update with Newton’s method with step size 1. The step sizes of the standard BFGS method and the L-BFGS method are determined by the Wolfe condition [35] using the backtracking line search algorithm to guarantee they converge on the whole dataset. All hyper-parameters (initialization parameter c, step size, batch size, etc.) of BFGS, L-BFGS, stochastic quasi-Newton method, SGD, SAGA, and Katyusha have been tuned to achieve the best performance on each dataset.
The convergence results are shown in Figures 1-4 for the considered datasets. We report both training error, i.e., RN (w) − RN (w∗N ), and test error for all algorithms in terms of number of effective passes over dataset and in terms of runtime. In general AdaQN mostly outperforms the first-order optimization methods (SGD, Katyusha and SAGA). This is caused by the fact that our considered problems are highly ill-conditioned. For instance, for the results of Epsilon dataset which has κ ≈ 104, there is a substantial gap between the performance of AdaQN and first-order methods.
1We use LIBSVM [34] with license: https://www.csie.ntu.edu.tw/~cjlin/libsvm/COPYRIGHT.
We observe that AdaQN outperforms BFGS, L-BFGS and SQN in all considered settings. As discussed earlier, this observation is expected since AdaQN is the only quasi-Newton algorithm among these methods that exploits superlinear convergence of quasi-Newton methods throughout the entire training process. Moreover, AdaQN does not require a line search scheme, while to obtain the best performance of BFGS and L-BFGS, we used a line-search scheme which is often computationally costly.
Next we compare Ada Newton and AdaQN. As shown in Figure 1, for the MNIST dataset with low dimension (d = 784) and moderate number of samples (N ≈ 11000), AdaQN performs very similar to Ada Newton. However, for settings where either the problem dimension d or the number of samples N is very large AdaQN outperforms Ada Newton significantly in terms of runtime. For instance, in Figures 2 and 3, which correspond to GISETTE dataset with N = 6000 samples and d = 5000 features and Orange dataset with N = 40000 samples and d = 14000 features, Ada Newton and AdaQN have almost similar convergence paths when compared in terms of number of passes over data, but in terms of runtime AdaQN is substantially faster than Ada Newton. This gap
comes from the fact that in these two examples the problem dimension d is large and as a result the cost of inverting a matrix at each iteration, which is required for Ada Newton, is prohibitive.
For Epsilon dataset that has a relatively large number of samples N = 80, 000 and moderate dimension d = 2000, as shown in Figure 4, Ada Newton outperforms AdaQN in terms of number of passes over data in both training and test errors. However, their relative performance is reversed when we compare them in terms of runtime. This time the slow runtime convergence of Ada Newton is due to the fact that sample size N is very large. Note that Ada Newton requires 2N Hessian computations, while for AdaQN we only need m0 Hessian evaluations. Remark 5. Theorem 1 requires the size of the initial training set m0 to satisfy the condition in (12). In our experiments, however, we observe that m0 could be much smaller than the threshold in (12). We believe that the main reason for this mismatch is that the recent non-asymptotic convergence analysis of BFGS that we exploit in our analysis may not be tight, and there is room for improving this bound. Specifically, establishing a less strict condition on the initial Hessian approximation error of BFGS required for achieving a superlinear convergence rate could lead to a much smaller lower bound for the parameter. Another possible reason for the gap between our theory and experiments is that our analysis is a worst-case analysis and it studies a lower bound that works for any setting. There could be some hard instances or corner cases that our current experiments do not capture, and for these hard instances we might need the number of samples required by our theoretical study. Identifying hard instances of empirical risk minimization problems for which a larger initial sample size is required is a research direction that requires further investigation.

heading: 6 Discussion and future work
text: We proposed AdaQN algorithm that leverages the interplay between the superlinear convergence of quasi-Newton methods and statistical accuracy of ERM problems to efficiently minimize the ERM corresponding to a large dataset. We showed that if the initial sample size m0 is sufficiently large, then we can double the size of training set at each iteration and use BFGS method to solve ERM subproblems, while we ensure the solution from previous round stays within the superlinear convergence neighborhood of the next problem. As a result, each subproblem can be solved with at most three BFGS updates with step size 1. In comparison with Ada Newton, AdaQN has three major advantages. (I) Ada Newton requires a backtracking technique to determine the growth factor, while for AdaQN we can double the size of training set if the size of initial set is sufficiently large. (II) Ada Newton implementation requires computing the Hessian inverse at each iteration, while AdaQN only requires computing a single matrix inversion. (III) Ada Newton requires 2N Hessian computations overall, while AdaQN only needs m0 Hessian evaluations at the initialization step.
In this paper, we showed that using quasi-Newton methods we only require constant number of iterations that is independent of the condition number κ and the dimension d to solve each subproblem to its statistical accuracy. More specifically, as the superlinear convergence rate of (1/ √ t)t for BFGS is independent of the problem parameters, each subproblem can be solved with at most three iterations. On the other hand, any other second-order method that enjoys a local convergence rate (either linear or superlinear) that is independent of the problem parameters could be capable of solving each subproblem to its statistical accuracy with a finite constant number of iterations that is independent of the condition number κ and dimension d. Hence, studying the application of such second-order or quasi-Newton methods for the considered adaptive sample size scheme is an interesting direction of research that requires further investigation.
One limitation of our results is that AdaQN provably outperforms other benchmarks only in the regime that N max { d, κ2s log d } . This is due to the lower bound on the size of initial training set m0. Notice that in many large-scale learning problems where N is extremely large, this condition usually holds. We believe this limitation of our analysis could be resolved by tighter non-asymptotic analysis of BFGS method. Moreover, our guarantees only hold for convex settings and extension of our results to nonconvex settings for achieving first-order or second-order stationary points is a natural future direction.

heading: Acknowledgement
text: This research is supported in part by NSF Grant 2007668, ARO Grant W911NF2110226, the Machine Learning Laboratory at UT Austin, and the NSF AI Institute for Foundations of Machine Learning.

Review summaries:
This paper  - introduces a new quasi-Newton (QN) algorithm called AdaQN to solve the ERM problem. This method is adaptive since it increases the sample size after few (3) runs of BFGS on a small ERM subproblem.  - gives its corresponding convergence analysis using a Matrix Bernstein bound - explains that AdaQN perfoms well because of the superlinear convergence on each subproblem which is possible when $m_0$ is large enough (greater than $\mathcal{O} (\kappa^2 log (d))$ if we leave out $s$) - provides clear comparison in term of theoretical computational cost is made against AdaNewton in Table 1  - performs numerical experiments showing the efficiency of the methods on small and medium size problems for l2 regularized logistic regression

This paper introduces an adaptive sample size quasi-Newton (AdaQN) method that exploits the superlinear convergence of quasi-Newton methods globally and throughout the entire learning process. The authors also provide convergence analysis to verify the AdaQN method, showing less gradient computation, as well as computational cost compared with previous SOTA methods. The numerical experiments on various datasets also confirm their theoretical findings.