Statistical and Computational Complexities of BFGS Quasi-Newton Method for Generalized Linear Models
Authors: Qiujiang Jin
Tongzheng Ren
Nhat Ho
Aryan Mokhtari
Venue: Accepted by TMLR ()

Abstract: The gradient descent (GD) method has been used widely to solve parameter estimation in generalized linear models (GLMs), a generalization of linear models when the link function can be non-linear. In GLMs with a polynomial link function, it has been shown that in the high signal-to-noise ratio (SNR) regime, due to the problem's strong convexity and smoothness, GD converges linearly and reaches the final desired accuracy in a logarithmic number of iterations. In contrast, in the low SNR setting, where the problem becomes locally convex, GD converges at a slower rate and requires a polynomial number of iterations to reach the desired accuracy. Even though Newton's method can be used to resolve the flat curvature of the loss functions in the low SNR case, its computational cost is prohibitive in high-dimensional settings as it is $\mathcal{O}(d^3)$, where $d$ the is the problem dimension. To address the shortcomings of GD and Newton's method, we propose the use of the BFGS quasi-Newton method to solve parameter estimation of the GLMs, which has a per iteration cost of $\mathcal{O}(d^2)$. When the SNR is low, for GLMs with a polynomial link function of degree $p$, we demonstrate that the iterates of BFGS converge linearly to the optimal solution of the population least-square loss function, and the contraction coefficient of the BFGS algorithm is comparable to that of Newton's method. Moreover, the contraction factor of the linear rate is independent of problem parameters and only depends on the degree of the link function $p$. Also, for the empirical loss with $n$ samples,  we prove that in the low SNR setting of GLMs with a polynomial link function of degree $p$, the iterates of BFGS reach a final statistical radius of $\mathcal{O}((d/n)^{\frac{1}{2p+2}})$  after at most $\log(n/d)$ iterations. This complexity is significantly less than the number required for GD, which scales polynomially with $(n/d)$.

Published in Transactions on Machine Learning Research (April/2024)
Statistical and Computational Complexities of BFGS Quasi- Newton Method for Generalized Linear Models
Qiujiang Jin
giaziong@Boustin.uteee.cmm wtw.:ssssssituttssed..
Deportment of Electricnl and Computer Engineering University of Tezas at Austin
Tongzheng Ren Department of Computer Science University of Texas ot Austin
Congzheng@uuteras edu
Nhat Ho
minknhhat@uie Satianaaauuueeas. edu
Department of Statistics and Data Sciences University of Teras at Austin
Aryan Mokhtari
mobhtari@austin. utw.as.saacccccmm... edu
Deportment of Electricol and Computer Engineering Universily of Tezus al Austin
Reviewed on OpenReview: https: // openreview. net/ form?? id= PIL3YYYa
Abstract
The gradient descent (GD) method has been used widely to solve parameter estimation in generalized linear models (GLMs), a generalization of linear models when the link function can be non-linear. In GLMs with polynomial link function, it has been shown that in the high signal-to-noise ratio (SNR) regime, due to the problem's 's strong convexily and smoothness, GD converges linearly and reaches the final desired accuracy in logarithmic mumber of iterations. In contrast, in the low SNR setting, where the problem becomes locally convex, GD converges at a slower rate and requires a polynomial number itcrations to reach the desired accuracy. Even though Newton's method can be used to resolve the flat curvature of the loss functions in the low SNR case, its computational cost is prohibitive in high-dimensional settings it is O(d), where d the is the problem dimension. To address the shortcomings of GD and Newton's 's method, quasi-Newton method to solve parameter estimation of the GLMs, which has a per iteration cost of O(P). When the SNR is low, for GLMs with a polynomial link function of degree we demonstrate that the iterates of BFGS converge linearly to the optimal solution the population least-square loss function, and the contraction coefficient. of the BFGS algorithm comparable that of Newton's 's method. Moreover, the contraction factor of the lincar rate independent of problem parameters and only depends on the degree of the link function p. Also, for the empirical loss with samples, we prove that in the low SNR setting of GLMs with polynomial link function of degree p. the iterates of BFGS reach final statistical radius of ((d/n)) कイंड) after at most. log(n/d) iterations. This complexity is significantly less than the number required for GD, which scales polynomially with (n/d).
propose the use of the BFGS
1 Introduction
In supervised machine learning, we are given a set of corresponding labels Yı,... Yn, that are drawn from some unknown distribution and our goal is to train a model that maps the feature vectors to their corresponding labels. We assume that the data is generated
n independent samples denoted by NI,..., X, with

---

Published in Transactions on Machine Learning Research (April/202))
according to distribution Pa- which is parameterized by ground truth parameter 8*. Our goal as the learner is to find 0 by solving the empirical risk minimization (ERM) problem defined as
min (0) = Σ 10; Yi), OER
(1)
where E(0; (X,, Yi) is the loss function that measwres the error between the predicted output of X; using parameter 0 and its true label Y, If we define 8% as an optimal solution of the above optimization problem, i.e., 0% E arg min (O), it can be considered as an approximation of the ground-truth solution 0*, where 6* is also a minimizer of the population loss defined as
min (()) = E (e(0; (X. AR
(2)
If one can solve the empirical risk efficiently, the output model could be close to 0", when n is sufficiently large. Several works have studied the complexity of iterative methods for solving ERM or directly the population loss, for the case that the objective function is convex or strongly convex with respect to 0 (Balakrishnan
et al.
2017 Ho et al.] 2020 Loh & Wainwright 2015
Agarwal et al. 2012 Yuan & Zhang 2013 Dwivedi However, when move beyond linear models, the
et al.
2020b Hardt et al.) 2016 Candes et al. 2011
underlying loss becomes non-convex and therefore the behavior of iterative methods could substantially change, and it is not even clear if they can reach neighborhood of a global minimizer of the ERM problem.
The focus of this paper is on the generalized linear model (GLM) (Carroll et al. 1997 Netrapalli et al., 2015 Fienup 1982 Shechtman et al. 2015 Feiyan Tian 2021 where the labels and features are generated
according to a polynomial link function and we have Y; = (X; 8 + where G is an additive noise and p 2 2 is an integer. Due to nonlinear structure of the generative model, even if we select a convex loss function t, the ERM problem denoted to the considered GLM could be non-convex with respect to 0. Interestingly, depending on the norm of 0", the curvature of the ERM problem and its corresponding population risk minimization problem could change substantially. More precisely, in the case that ||0*]| is sufficiently large, which we refer to this case as the high signal-to-noise ratio (SNR) regime, the underlying population loss of the problem of interest is locally strongly convex and smooth. On the other hand, in the regime that |0*|| is close to zero, denoted by the low SNR regime, the underlying problem neither strongly convex nor smooth, and in fact, it is ill-conditioned.
These obscrvations lead to the conclusion that in the high SNR setting, due to strong convexity and smoothness of the underlying problem, gradient descent (GD) reaches the desired accuracy exponentially fast and overall it only requires logarithmic number of iterations. However, in the low SNR case, as the problem becomes locally convex, GD converges at a sublincar rate and thus requires polynomial number of itcrations in terms of the sample size.
To address this issue, Ren et al. (2022a) advocated the use of GD with the Polyak step size to improve GD's convergence in low SNR scenarios. hey demonstrated that the number of iterations could become logarithmic with respect to the sample size, when GD is deployed with the Polyak step size. However, such method still remains first-order algorithm and lacks any curvature approximation. Consequently, its overall complexity is directly proportional to the condition number of the problem. This, in turn, depends on both the condition number of the feature vectors' covariance and the norm [0*). As a result, in low SNR settings, the problem becomes ill-conditioned with a large condition number, and hence GD with the Polyak step size could be very slow. Furthermore, the implementation of the Polyak step size necessitates access to the optimal value of the objective function. As precise estimation of this optimal objective function value may not always be feasible, any inaceuracies could potentially lead to reduced convergence rate for GD employing the Polyak step size.
Another alternative is to use a different distance metric instead of an adaptive step size to improve the convergence of GD for the low SNR setting. More precisely, Lu et al. (2018) have shown that the mirror descent method with a proper eettacc-generating function can solve the population loss corresponding to the low SNR setting at a lincar rate. However, the linear convergence rate of mirror descent, similar GD with Polyak step size, also depends on the condition number of the problem, and hence could lead to a slow convergence rate.
2

---

Published in Transactions on Machine Learning Research (April/2024)
A natural approach to handle the ill-conditioning issue in the low SNR case well as eliminating the need to estimate the optimal function value is the use of Newton's method. As we show in this paper, this idea indeed addresses the issue of poor curvature of the problem and leads to an exponentially fast rate with contraction factor 2p-2 in the population case, where is the degree of the polynomial link function. Morcover, in the high N setting, Newton's method converges at a quadratic rate as the problem is strongly convex and smooth. Alas, these improvements come at the expense of increasing the computational complexity of each iteration to (dd)) which is indeed more than the per iteration computational cost of GD that is (d)). These points lead to this question:
Is there a computationally-eficeen method that performs well in both high and low SNR settings at a reasonable per iteration computational cost?
Contributions. In this paper, we address this question and show that the BFGS method is capable of achieving these goals. BFGS is quasi-Newton method that approximates the objective function curvature using gradient information and its per iteration cost is (df). It is well-known that it enjoys superlinear convergence rate that is independent of condition number in strongly convex and smooth settings, and hence, in the high SNR setting it outperforms GD. In the low SNR setting, where the Hessian at the optimal solution could be singular, we show that the BFGS method converges linearly and outperforms the sublinear convergence rate of GD. Next, we formally summarize our contributions.
Infinite sample, low SNR: For the infinite sample case where minimize the population loss, we show that in the low SNR case the iterates of BFGS converge to the ground truth 0* at exponentially fast rate that is independent of all problem parameters except the power of link function We further show that the linear convergence contraction coefficient of BFGS is comparable to that of Newton's method. This convergence result of BFGS is also of general interest as it provides the first global línear convergence of BFGS without line-scarch for a sctting that is neither strictly nor strongly convex.
Finite sample, low SNR: By leveraging the results developed for the population loss of the low SNR regime, show that in the finite sample case, the BFGS iterates converge to the final statistical radius ((1/2/((P+2)) within the true parameter after logarithmic number of iterations (log(n)). It is substantially lower than the required number of iterations for fixed-step size GD, which is ((((--)/p), to reach a similar statistical radius. This improvement is the direct outcome of the linear convergence of BFGS versus the sublinear convergence rate of GD in the low SNR, case. Further, while the iteration complexity of BFGS is comparable to the logarithmic number of iterations of GD with Polyak step size, we show that BFGS removes the dependency of the overall complexity on the condition number of the problem as well as the need to estimate the optimal function value.
Experiments:
We conduct numerical experiments for both infinite and finite sample cases to compare the performance of GD (with constant step size and Polyak step size), Newton's 's method and BFGS. The provided empirical results are consistent with our theoretical findings and show the advantages of BFGS in the low SNR regime.
Outline. In Section 2 we discuss the BFGS quasi-Newton method. Section 3 details three scenarios in Generalized Linear Models (GLMs): low, middle, and high SNR regimes, outlining the characteristics of the population loss in each. Scction explores BFGS's convergence in low SNR scttings, highlighting its linear convergence rate, marked improvement over gradient descent's sublinear rate. This section also compares the convergence rates of BFGS and Newton's 's method. Section 5 applies these insights to establish the convergence results of BFGS for the empirical loss Ln in the low SNR regime. Lastly, our numerical experiments are presented in Section
2 BFGS algorithm
In this section, we review the basics of the BFGS quasi-Newton method, which is the main algorithm we analyze. Consider the case that we aim to minimize a differentiable convex function f Rd - R. The BFGS
3

---

Published in Transactions on Machine Learning Research (April/202)
update is given by
oc+1 = 02 - 4k 0,
(3)
where Th is the step size and H2 E Rddd is a positive definite matrix that aims to approximate the true Hessian
inverse ? [²x(0))= There are several approaches for approxinating H2 leading to different quasi-Newton
methods, CConn al.,
1991
Broyden 1965
Broyden al. 1973 Gay 1979
Davidon,
1959 Fletcher &
Powell 1963 Broyden 1970
Fletcher
1970
Goldfarb 1970 Shanno 1970: Nocedal] 1980 Liu & Nocedal
1989), but in this paper, we locus on the celebrated BFGS method, in which H& is updated
1 "n-(1-
Sk- sk-14
(
uk-18kk1 1k-1
+ sk-1šk_1
Ha-1
Vk 2 1,
(4)
1 1k-1
where the variable variation Sy and gradient displacement ug are defined as
8k-1 = 02 0k-1,
Uk-1 == (f(as) (fx==-));
(5)
for any h 2 1 respectively. The logic behind the update in (4 is to ensure that the Hessian inverse approximation H3 satisfies the secant condition Ha = 8k-1, while it stays close to the previous approximation matrix Hę-1- The update in ( only requires matrix-vector multiplications, and hence, the computational cost per iteration of BFGS is (d3)).
The main advantage of BFGS is its fast superlinear convergence rate under the assumption of strict convexity,
0 8ppt 0,
lim
| Oppl ||
where Oppt is the optimal solution. Previous results on the superlinear convergence of quasi-Newton methods were all asymptotic, until the recent advancements on the non-asymptotic analysis of Newstn methods
[Rodomanov & Nesterov 2021a b Jin & Mokhtari] 2020 Jin et. al. 2022 Ye et al., 2021 Lin et al. 2021a bb.
For instance, Jin & Mokhtari 2020) established local superlinear convergence rate of (I/Vk)= for BFGS.
However, all these superlinear convergence analysos require the objective function to be smooth and strictly or strongly convex. Alas, these conditions are not satisfied in the low SNR setting, since the Hessian at the optimal solution could be singular, and hence the loss function is neither strongly convex nor strictly convex; we further discuss this point in Section 3 This observation implies that we need novel convergence analyses to study the behavior of BFGS in the low SNR sctting, as we discuss in the upcoming scctions.
3 Generalized linear model with polynomial link function
In this section, we formally present the generalized linear model (GLM) sctting that we consider in our paper, and discuss the low and high SNR settings and optimization challenges corresponding to these cases. Consider the case that the feature vectors are denoted by X E IR and their corresponding labels are denoted by Y R. Suppose that we have access to n sample points (Y,, X1), X2), (Yn, X.) that are i.i.d. samples from the following generalized linear model with polynomial link function of power p (Carroll et al. 1997). i.e.,
Y = (XT O
where 8* is a true but unknown parameter, p E N is a given power, and C1...., Sn are independent. Gaussian noises with zero mcan and variance d². The Gaussian assumption on the noise is for the simplicity of the discussion and similar results hold for the sub-Gaussian i.i.d. noise case. Furthermore, we assume the feature vectors are generated X ~ (0, Σ) where Σ E dddd is symmetric positive definite matrix. Here we focus on the settings that E N+ and p 2 2.
The above class of GLMs with polynomial link functions arise in several settings. When p= 1, the model in (6) is the standard linear regression model, and for the case that p = 2, the above setup corresponds to the phase retrieval model (Fienup) 1982. Shechtman et al. 2015 Candes et al.) 2011 Netrapalli et al. 2015)) which has found applications in optical imaging, x-ray tomography, and audio signal processing. Moreover, the analysis of GLMS with p 2 2 also serves as the basis of the analysis on other popular statistical models.
4

---

Published in Transactions on Machine Learning Research (April/2024)
For example, shown by Yi & Caramanis (2015); Wang et al. (2015): Xu et al. (2016): Balakrishnan et al.
[2017] Daskalakis et al. 201 : Dwivedi et al. (2020a) Kwon et. al.
(201)) : Dwivedi et al. (2020b); Kwon
et al. (2021, the local landscape of log-likelihood for Gaussian mixture models and mixture linear regression models are identical to GLMs for 2.
In the case that the polynomial link function parameter in the GLM model in 6 is p= 2, by adapting similar arguments from Kwon et al. (002]) under the symmetric two-component location Gaussian mixture, there are essentially three regimes for estimating the true parameter 8": Low hignal-to-noise ratio (SNR) regime: |/0 I Ci (d/n))// where d is the dimension, is the sample size, and C1 is universal constant; Middle SNR regime: Ci (dd/n) < ||0* s C' where C' is universal constant; High SNR regime: || 2 C. The main idea is that with different 0", the optimization landscape of the parameter estimation problem changes. By generalizing the insights from the case p = 2, we define the following regimes for any p 2:
(i) Low SNR regime: ||0^1// - (dd/m) where d is the dimension, n. is the sample size, and C1 is universal constant;
(ii) Middle SNR regime: Gidd 1// S C where C a universal constant;
(iii) High SNR regime: ||00|//0 C.
Note that, the rate (d/n)1/2p that we use to define the SNR regimes is from the statistical rate of estimating the true parameter 0 when 0* approaches O. Next, we provide insight into the landscape of the least-square loss function for each regime. In particular, given the GLM in (6), the sample least-square takes the following form:
min C,) 6E
==
1 Σ n2
(xiTx))))))
Y
(7)
To obtain insight about the landscape of the empirical loss function L., a useful approach is to consider that function by its population version, which we refer to as population least-square loss function:
min L()) == (0), BER"
(8)
where the outer expectation taken with respect to the data.
High SNR regime. In the setting that the ground truth parameter has relatively large norm, i.e., ||01| 2 C for some constant C >> 0 that only depends a, the population loss in (8) is locally strongly
convex and smooth around 8* More precisely, when ||0 - 1 is small, using Taylor's theorem we have
(x (xT 0)"
0')" p(X
(0-00)))))) 0*) o - 1I))).
Hence, in a neighborhood of the optimal solution, the objective in
be approximated as
((0) = E[(Y - = + -
- + o² = =p-(0-0") Ex (X(xTT0)300221 (0- -000000
(x [E[((X (1
-x
Indeed, if [ 2 Co the above function behaves as a quadratic function that smooth and strongly convex, assuming that o(110 - 0* (12) is negligible. As a result, the iterates of gradient descent (GD) converge to the solution at linear rate and it requires llog(1/c) to reach e-accurate solution, where depends the conditioning of the covariance matrix Σ and the norm of 6*, In this case, BFGS converges superlinearly to 0* and the rate would be independent of K, while the cost per itcration would be (dd))
Low SNR regime. As mentioned above, in the high SNR case, GD has a fast linear rate. However, in the low SNR. case where ||0*|| is small and ||0*|| < C1 C1(d/n)/(22) the strong convexity parameter approaches zero when the sample size n goes to infinity and the problem bocomes ill-conditioned. In this case, we deal with a function that only conver and its gradient is not Lipschilz contiruous. To better elaborate on this point, let us focus on the case that 0r = 0 a special case of the low SNR setting. Considering the underlying distribution of X, which is X (0 M), for such low SNR case, the population loss can be written
(()) Ex [(X [)pp + 0² (2p 1) + o²
5

---

Published in Transactions on Machine Learning Research (April/202))
Since we focus on p 2 2 it can be verified that (()) is not strongly convex in a neighborhood of the solution 0* = 0. For this class of functions, it is well-known that GD with constant step size would converge at a sublinear rate, and hence GD iterates require polynomial number of iterations to reach the final statistical radius. In the next section, we study the behavior of BFGS for solving the low SNR setting and showcase its advantage compared to GD.
Middle SNR regime. Different from the low and high SNR regimes, the middle SNR regime is generally harder to analyze as the landscapes of both the population and sample least-square loss functions complex. Adapting the insight from middle SNR regime of the symmetric two-component location Gaussian mixtures from Kwon et al. (2021), for the middle SNR sctting of the generalized linear model, the cigenvalues of the Hessian matrix of the population least-square loss function approach O and their vanishing rates depend on some increasing function of [| The optimal statistical rate of the optimization algorithms, such as gradient descent algorithm, for solving 8* depends strictly on the tightness of these vanishing rates in terms of ||@*|], which are non-trivial to obtain. In fact, to the best of our knowledge, there is no result on the convergence of iteralive methods (such as GD or its variants) for GLMs with a polynomial link function in the middle SNR. Hence, we leave the study of BFGS for this setting as future work.
4 Convergence analysis in the low SNR regime: Population loss
In this section, we focus on the convergence properties of BFGS for the population loss in the low SNR case introduced in 99. This analysis provides an intuition for the analysis of the finite sample case that we discuss in Section 5 as we expect these two loss functions to be close to each other when the number of
samples n is sufficiently large. In this section, instead of focusing on the population loss within the low SNR regime, as described in (9). we shift our attention a more encompassing objective function, f. Detailed in the following expression, this function encompasses the population loss function L with 0* 0 a specific example. This approach enables us to present our results in the most general setting possible. Specifically, we examine the function f: Rd + R, defined as follows:
min f(0) = ||.10 -61%, DER
(10)
where A E Rmnd is a matrix, b E RT is a given vector, and E Z satisfies 2 4. We should note that for
q 2 4, the considered objective is not strictly convex because the Hcssian is singular when AO = b. Indoed, if we set m = d and further let A be 11/2 and choose b b= A0 = O, then we recover the problem in ( for
g = 2p. Note that since we focus on the generalized linear model with the polynomial link function, which necessitates that p be an integer, the parameter q = 2p also an integer.
Notice that the problem in (10), which serves as surrogate for the finite sample problem that we plan to study in the next section, has the same solution set as the quadratic problem of minimizing ||A0 - bll with solution (A A when A T A is invertible. Given this point, one might suggest that instead of minimizing (10) we could directly solve the quadratic problem which is indeed much casier to solve. This point is valid, but the goal of this section is not to efficiently solve problem (10) itself. Our goal is to understand the convergence propertics of the BFGS method for solving the problem in (10) with the hope that it will provide some intuitions for the convergence analysis of BFGS for the empirical loss (7) in the low SNR regime. As we will discuss in Section 5 the convergence analysis of BFGS the population loss which is special case of (10) is closely related to the one for the empirical loss in 77
One remark is that population loss in g holds for the restrictive assumption of 0" =0, which is only special case of the general low SNR regime of TTo1 I (d(d/n)/((2)) Our ultimate goal is to enalyze the convergence behavior of the BFGS method applied to the empérical loss 7 of GLM problems in the low SNR regime. The errors between gradients and Hessians of the population loss with 6* = 0 and | S Ci(d/n) 1/(0)) are upper bounded by the corresponding statistical errors belween the population loss 8 and the empirical loss 7 in the low SNR regime, respectively. Therefore, the errors between iterations of applying BFGS to
Remark 1.
the population loss with 6* and ||0"1| S (dd are negligible compared to the statistical errors. Instead of directly analyzing BFGS for the population loss R in the general low SNR regime, studying the convergence properties of BFGS for the population loss 9 with 0* = 0 can equivalently lay foundations for
6

---

Published in Transactions on Machine Learning Research (April/2024)
the convergence analysis of BFGS for the empirical loss of this paper. The details con be found in Appendir A.4
7 in the low SNR regime, which is the main target
Our results in this section are of general interest from an optimization perspective, since there is no global convergence theory (without line-search) for the BFGS method in the literature when the objective function is not strictly convex. Our analysis provides one of the first results for such general settings. That said, it should be highlighted that our results do not hold for general convex problems, and we do utilize the specific structure of the considered convex problem to establish our results. The following examples illustrate some of the specific structures that exploit to establish our result.
Assumption 1. There erists E IRd, such that b b= Aό. In other words, b is in the range of matrix A.
This assumption implies that the problem in 10 is realizable, Ô an optimal solution, and the optimal function value zero. This assumption is satisfied in our considered low SNR setting in as 0 0 which implies b = O.
Assumption 2. The malrir A A E Rddd is invertâble. This is equivalent to A 0.
The above assumption is also satisfied for our considered setting as we assume that the covariance matrix for
our imput features is positive definite. Combining Assumptions and 2 we conclude that the unique optimal solution of problem (10). Next, we state the convergence rate of BFGS for solving problem under the disclosed assumptions. The proof of the next theorem is available in Appendix A.1
Theorem 1. Consider the BFGS method in 3 1-.5. Suppose Assumptions ] and hold, and the initial
Hessian inverse approrimation matrir is selected Ho ('6)) where 6g R is arbitrary initial point. If the step size is Ti 1 for all 2 0, then the iterates of BFGS converge the optimal solution Ā at linear rate of
102- 011 S TT-- 10-1 011
vk 1 1,
(11)
where the contraction factors Th [0, 1) satisfy
4 ro q rk
1 1- --rr--
4--2
Vk: > 1.
(12)
1 h-1
Theorem 0 shows that the iterates of BFGS converge globally at a linear rate to the optimal solution of 10 This result is of interest as it illustrates the ilerates generated by BFGS converge globally without any line search scheme and the step size is fixed as = 1 for any k 2 0. Moreover, the initial point dg is arbitrary and there is no restriction on the distance between 6g and the optimal solution 0.
We should highlight. the above linear convergence result and our convergence analysis both rely heavily on the distinct structure of problem (10) and may not hold for a general convex minimization problem. Specifically, it can be shown that if we had access to the exact Hessian and could perform Newton's update solve the problem in TOP, then the error vectors (0x- - 0)k20 would all be parallel. This properly arises due to the fact that for the problem in (10 the Newton direction is V (0x--))--V((00
= (0-11-0 0) 3 (0 0). Consequently, the next error vector 0x-1 9 computed by running one step of Newton would satisfy O2. - 0 - 2-2 (0x-1 - Ø). Therefore, the error vector at time k, denoted by 0, θ, is parallel to the previous error vector a - 0, with the only difference being that its norm is reduced by factor of . Using an induction argument, it simply follows that all error vectors (0p- - -)x20 remain parallel to each other while their norm contracts at rate of 2=?. This is key point that is used in the result for Newton's is method in Theorem 3 In the convergence analysis of BFGS (stated in the proof of Theorem T. we use similar argument. While we cannot guarantee that the Hessian approximation matrix in the BFGS method matches the exact Hessian, we demonstrate that it can inherit this property, maintaining all error vectors (0,- - 0))20 as parallel to each other. Additionally, we show that the norm is shrinking at a factor of Fа < 1, which is always larger than q-1+ yet remains independent of the problem's condition number or dimensions. In the following theorem, we show that for q 2 4, the linear rate contraction factors {r also converge lincarly to fixed point contraction factor r. determined by the parameter q. The proof is available in Appendix A.2
9-1
g-2
7

---

Published in Transactions on Machine Learning Research (April/202)) 2024)
2.04
0.7
12*
0.24
Pk T2
7
Pb
046
Number of Iteration
Number of Iteration
Number of Iteration k
Number of Iteration
(a) 44
(b) = 100.
(c) = 4.
(d) g = 100.
Figure Convergence of factors {ra]ã= to T..
Theorem 2. Consider the linear convergence factors {r defined in 12 from Theorem 4f 4 and Z, then the sequence (rk lk= converges to r. (0) that is determined by the equation
- +
(13)
and the rate of convergence linear with a contraction factor that at most 1/2, i.e.,
Ir - r. s (1/2)* po-
4k 2/20 0.
(14)
Based on Theorem 2 the iterates of BFGS eventually converge to the optimal solution at the linear rate of r. determined by 13 Specifically, the factors {ri] converge to the fixed point r. at a linear rate with the contraction factor of 1/2. Further, the linear convergence factors {r and their limit r. are all only determined by q, and they independent of the dimension d and the condition number кд of the matrix A. Hence, the performance of BFGS is not influenced by high-dimensional or ill-conditioned problems. This result is independently important from an optimization point of view, as it provides the first global lincar convergence of BFGS without line-search for a setting that is not strictly or strongly convex, and interestingly the constant of linear convergence is independent of dimension or condition number.
We illustrate the convergence of factors {r to the fixed point r. in Figure I for q= and q = 100. In plots (a) and (b), we observe that rk becomes close r. after only iterations. Hence, the linear convergence rate of BFGS is approximately r, after only few iterations. We further observe in plots (c) and (d) that the factors {Tx 12 o converge to the fixed point at linear rate upper bounded by 1/2. Note that is the
solution of 13 These plots verify our results in Theorem 2
Remark 2. The cost. of computing the inilial Hessian inverse approzimation Hg V²5(0)) is Odd), but this cost is only required for the first iteration, and it is not required for k 2 1 for those iterates update the Hessian inverse approzimation matrix Hu based on the update in (4) at cost of O(d²).
Remark 3. Although the vanilla. Gradient Descent (GD) method converges at sub-linear rate when applied to problem (10). it has been shown that other first-order methods, such as Mirror Descent with a distance- generating function h(x)) = -t/|=|ot2 + |x/|2 can solve problem 110 at a linear rate (refer to Lu et al. (2018)). However, the linear convergence rule of Mirror Descent is dependent on the condition wammbr k of the problem, characterized by the rate (1= In contrast, as demonstrated by Theorems 1 and 2 the linear convergence rate of BFGS is independent of the problem' 's condition number. Consequently, in settings with low SNR, which can be ill-conditioned, we anticipate a faster convergence rate for BFGS compared to Mirror Descent. We illustrate this point in our numerical experiments presented Figure 2
4.1 Comparison with Newton's method
Next, we compare the convergence results of BFGS for solving problem (10) with the one for Newton's method. The following theorem characterizes the global lincar convergence of Newton's 's method with unit step size applied to the objective function in 110
Theorem 3. Consider applying Newton's method optimization problem 10) and suppose Assumptions and 2 hold. Moreover, suppose the step size is ne for any k: 2 0. Then, the ilerates of Newton's method
8

---

Published in Transactions on Machine Learning Research (April/2024) 2024)
converge to the optimal solution at linear rate of
|10g 01 9 2
|10 -
Ak 2 1.
(15)
Moreover, this linear convergence rate 2- is smaller than the fixed point T. defined in 13 of the BFGS quasi-Newton method, i.e., 3 <r.< for all q
The proof is available in Appendix A.3 The convergence results of Newton's method are also global without any line search method, and the linear rate 2=2 is independent of dimension d and condition number KA- q-1
Furthermore, the condition q-1 < r. implies that iterates of Newton's method converge faster than BFGS, but the gap is not substantial as r. < 2g-3 On the other hand, the computational cost per iteration of Newton's method is Odd) which worse than the O(d2) of BFGS.
2q-2"
Moving back to our main problem, one important implication of the above convergence results is that in the low SNR setting the iterates of BFGS converge linearly to the optimal solution of the population loss function, while the contraction coefficient of BFGS is comparable to that of Newton's method which is (2p - 2) (2p 1). For instance, for = 2,3,5, 10, the linear rate contraction factor of Newton's method are 0.667, 0.8, 0.889, 0.947 and the approximate linear rate contraction factor of BFGS denoted by т, are 0.755, 0.857, 0.922, 0.963, respectively.
5 Convergence analysis in the low SNR regime: Finite sample setting
Thus far, we have demonstrated that the BFGS iterates converge linearly the true parameter /* when minimizing the population loss function of the GLM in (9) in the low SNR regime. In this section, we study the statistical behavior of the BFGS iterates for the finite sample case by leveraging the insights developed in the previous section about the convergence rate of BFGS in the infinite sample case, i.e., population loss. More precisely, we focus on the application of BFGS for solving the least-square loss function C,, defined in for the low SNR setting. The ilerates of BFGS in this case follow the update rule
0 0g - Mh Vn
(16)
where HT is updated using the gradient information of the loss L. by the BFGS rule.
We next show that the BFGS iterates (66 Jx20 converge to the final statistical radius within a logarithmic number of iterations under the low SNR regime of the GLMs. To prove this claim, we track the difference between the iterates {Øf)a2 generated based on the empirical loss and the iterates {0&) generated according to the population loss. Assuming that they both start from the same initialization tp. with the cancentration of the gradient (0) - VL())] and the Hessian [)) V² from Mou et al. (2019); Ren et al. (202bb) we control the deviation between these two sequences. Using this bound and the convergence results of the iterates generated based on the population loss discussed in the previous section, we prove the following result for the finite sample setting. The proof available in Appendix AA.5
Theorem 4. Consider the low SNR regime of the GLM in (6) namely, ||@*|| S Ci (d//))/(een) Apply the BFGS method to the empirical loss 7 with the initial point 0% where 0% E B(0", r) for some r > 0, the initial Hessian inverse opprocimation matrir as Ho = V²L..(03) (Og) where V²L, is posilive definite and step size T = 1. For any failure probability ó E (0,1), if the number of samples is n 2 Cz(dlog(d/6))² and the number of iterations satisfies T 2 C3 (oggg/((oog(////// (1)))), then with probability 1 8, we have
0 S C4 (dlog(1/6) 104 |1p|-
min tE[
(17)
n
where C1, C2, Cs, and C4 are constants independent of and d.
Our analysis hinges on linking the gradient and Hessian of the empirical loss to the population loss, subsequently demonstrating convergence rate for the empirical loss based on the linear convergence established in Section
9

---

Published in Transactions on Machine Learning Research (April/2024)
Nowton
BBGS
log (0x
Newton
log (0x
Newson
log| (0g
GG Polyuk
-GD constand
log 00
10-
Newton
GD consland
BFGS
10
GD conslant
BFGS
10*
Mirror Descent
10"
CD constant
BFGS
Polvak
GD Polvak
GD Polyak
Descent
Descent
Descent
Number of Iteration &
Number of Iteration &
Number Iteration
Number of Iteration
(a) = 10 --
(b) 10, 10.
(c) d 10 0 -0= 4.
(d) d 10", ,9 = 10.
Figure 2: Convergence of Newton's method, BFGS, CD with constant slep size, GD with Polyak step size and
mirror descent with constant step size for different values of and g. In plot (a), m 100 and r) = 10 4. In plot (b), m 100 and r 10-5 In plot (c), m 2000 and 1j = 10-12. In plot (d), m 2000 and 1) 10-13
for the population loss. More details can be found in Appendix A.5 Theorem shows that BFGS achieves an ostimation accuracy of (((d/n)=zi æगス jin O(log n) iterations, which is substantialy faster than GD that requires itcations to achiove the same guarantce as shown in (Ren et al. 2022a). A few comments about
(n p-1
Theorem 4 are in order.
Comparison to GD, GD with Polyak step size, and Newton's method: Theorem 4 indicates that under the low SNR regime, the BFGS iterates reach the final statistical radius O(n-1/(2x+2)) within the true parameter 0" after ((log(n) number of iterations. The statistical radius is slightly worse than the optimal statistical radius O(n 1/(2)))). However, we conjecture that this is due to the proof technique and BFGS can still reach the optimal Onn11(2p) practice. In our experiments, in the next scction, we observe that when d = 4 and p p= 2, the statistical radius of BFGS is closer to the optimal radius of O( instead of ((n-1/6) suggested by our analysis. We leave an improvement of the statistical analysis the future work. On the other hand, the overall iteration complexity of BFGS, which is (log(n)), is indeed better than the polynomial number of iterations of GD, which is at the order of (nn(--1/mr) (Corollary 3 in (Ho et al. 2020 ).
Moreover, the complexity of BFGS is better than the one for GD with Polyak step size which is O( log(n)) iterations (Corollary in (Ren et al. 2022a)), ), where k is the condition number of the covariance matrix Σ. Note that while the iteration complexity of BFGS comparable to that GD with Polyak step size in terms of the sample size, the BFGS overcomes the need to approximate the optimal value of the sample least-square loss L., which can be unstable in practice, and also removes the dependency on the condition mumber that appears in the complexity bound of GD with Polyak step size. A similar conclusion also holds for mirror descent with a eestance-generatinn
function h(x) = 442 |z| + = as its linear convergence rate has a contraction factor of (1 1/x), leading to overall itcration complexity of (og)))))
Finally, the ileration complexity of BFGS is comparable to the O(lo((n) of Newton's method (Corollary 3 in (Ho et al. 2020 while per iteration cost of BFGS is substantially lower than Newton's method.
On the minimum number of iterations: The results in Theorem involve the minimum number of iterations, namely, this result. holds for some S / S T. suggests that the BFGS iterates may diverge after they reach the final statistical radius. As highlighted in (Ho et al.) 2020), such instability behavior of BFGS is inherent to fast and unstable methods. While may sound limited, this issue can be handled via an early stopping scheme using the vols-validation approaches. We illustrate such early stopping of the BFGS iterates for the low SNR rogime in Figure 3
6 Numerical experiments
Our experiments are divided into two sections: the first. focuses on iterative methods' behavior on the population loss GLMs with polynomial link functions, and the second examines the finite sample setting.
Experiments for the population loss function. In this section, we compare the performance of Newton's method, BFGS, GD with constant step size, GD with Polyak step size, and and mirror descent with constant.
10

---

Published in Transactions on Machine Learning Research (April/2024)
Slope: -0 -05014 Error
Sloper -0.2847 Eiror
log (0x 00
Newton
GD constant
log O0x
constane
log 10,
log 10,
GD Polyak
CD Polyak
Number of Iteration &
Number of Iteration &
Log of Sample Size n
Log of Sample Size 12
(a) High SNR regime.
(b) Low SNR regime.
(c) High SNR regime.
(d) Low SNR regime.
Figure 3: Convergence of different methods when d 4 in high SNR (a) and low SNR (b) regimes. Illustration of the statistical radius of BFGS in high SNR (c) and low SNR (d) regimes.
||/=]]] applied to (10) which corresponds to the population loss. We choose different values of parameter m, dimension d and the exponential parameter in (10). We generate a random matrix A E Rmxdd and a random vector @ E Rd, and compute the vector b AB e R.. The initial point 0g R is also generated randomly. To properly select the steps for GD with a fixed step size 1), we employed a manual grid search across the following values [10-10, 10-9.... ..... 10-1,1] and selected the value that yielded the best performance for each specific problem. In our plots, we present the logarithm of error ||0g - O versus the number of iteration k for different algorithms. All the values of different parameters m, d, and 77 are mentioned in the caption of the plots.
step size and distance generating function h(x) = q42 2 ||==| 2
In Figure D we observe that GD with constant step converges slowly due to its sub-linear convergence rate. The performance of GD with Polyak step size is also poor when dimension is large or the parameter g is huge. This is due to the fact that as dimension increases the problem becomes more ill-conditioned and hence the linear convergence contraction factor approaches 1. We observe that both Newton's method and BFGS generate iterations with linear convergence rates, and their linear convergence rates are only affected by the parameter q. i.c., the dimension d has no impact over the performance of BFGS and Newton's method. Although the convergence speed of Newton's 's method is faster than BFGS, their gap is not significantly large as we expected from our theoretical results in Section GD, but its performance is not as good as the BFGS method.
| We also observe that mirror descent outperforms
Experiments for the empirical loss function. We next study the statistical and computational complexities of BFGS on the empirical loss. In our experiments, we first consider the case that d = 4 and the power of the link function is = 2, namely, we consider the multivariate setting of the phase retrieval problem. The data is generated by first sampling the inputs according to {X.]]= ~ [O0 diag(o), ,o)) where a = (0, and then generating their labels based on Y, = (X, 0²)2 + S. where (G)= are i.i.d. samples from (V 0.01). In the low SNR regime, we set H- 0, and in the high SNR regime we select H* uniformly random from the unit sphere. Further, for GD, set the step size as n = 0.1, while for Newton's method and BFGS, we use the unit step size 17 1.
In plots (a) and (b) of Figure 3 we consider the setting that the sample size is n = 10t and we run GD, GD with Polyak step size, BFGS, and Newton's method to find the optimal solution of the sample least-square loss C. Furthermore, for both Newton's method and the BFGS algorithm, due to their instability, we also perform voss-validation to choose their early stopping. In particular, we split the data into training and the test sets. The training set consists of 90% of the data while the test set has 10% of the data. The yellow points in plots (a) and (b) of Figure 3 show the iterates of BFGS and Newton, respectively, with the minimum validation loss. As we observe, under the low SNR regime, the iterates of GD with Polyak step size, BFGS and Newton's 's method converge geometrically fast to the final statistical radius while those of the GD converge slowly lo that radius. Under the high SNR regime, the iterates of all of these methods converge geometrically fast to the final statistical radius. The faster convergence of GD with Polyak step size over GD is due the optimal Polyak step size, while the faster convergence of BFGS and Newton's method over GD is due their independence on the problem condition number. Finally, in plots (c) and (d) of Figure 3 we run BFGS when the sample size is ranging from 102 to 10a to empirically verify the statistical radius of
11

---

Published in Transactions on Machine Learning Research (April/2024)
Slope: +0.5028 Beror
Slope: 40.2884 Error
Newton
log (0
constant
log O0
log i0.
log 0
GD Polyak
Newton
BFGS
CD constaut
Polyak
Number of Iteration &
Number of Iteration
Log of Sample Size n
Log of Sample Size 12
(a) High SNR regime.
(b) Low SNR regime.
(c) High SNR regime.
(d) Low SNR regime.
Figure 4: Convergence of different methods d = 50 in high SNR (a) and low SNR (b) regimes. Statistical radius of BFGS in high SNR (c) and low SNR (d) settings.
011
log Oz Ø1
Newton BFGS
Newton BRGS
log (0x
GD comslant Pulyak
log O2
GD comstant
log 0
Newlou BFGS
GD constauz:
GD constamt
CD Polyak
GD Polyak
28
66
30
20
30
20
10
16"
29
60 se 100
Number of Iteration k
Number of Iteration &
Number of Iteration &
Number of Iteration k
(a) High SNR (d 100).
(b) Low SNR (d 100).
(c) High SNR (d 500).
(d) Low SNR (d 500).
Figure 5: Convergence of different methods with d = 100 for high SNR regime are shown in (a) and low SNR regime in (b). Convergence of different methods with d = 500 for high SNR regime are shown in (c) and low SNR regime in (d).
these methods. As indicated in the plots of that figure, under the high SNR regime, the BFGS has statistical radius is (0n-112))))))) while under the low SNR regime, its statistical radius becomes O(n-1/4).
To show that BFGS can also be applied to high dimension scenarios, conduct additional experiments on the generalized linear model with input d = 50, 500 and the power of link function = 2. The inputs are generated by {X.] ~ NO(, diag(a?, o3)) where a = (0.96)*-1, and the remaining setting and hyper-parameters are set identical to the low dimension scenarios. The results are shown in Figure S and 5 As the results show, the performance BFGS in high dimensional scenarios are nearly identical to the low dimensional scenarios.
7 Conclusions
In this paper, we analyzed the convergence rales BFGS on both population and empirical loss functions of the generalized linear model in the low SNR regime. We showed that in this case, BFGS outperforms GD and performs similar to Newton's method in terms of iteration complexity, while requires lower per iteration computational complexity compared to Newton's method. We also provided experiments for both infinite and finite sample loss functions and showed that our empirical results are consistent with our theoretical findings. Perhaps one limitation of the BFGS method is that its computational cost is still not lincar in the dimension and scales as (dd)). One future research direction is to analyze some other iterative hods such as limited memory-BFGS (L-BFGS) which may be able to achieve a fast linear convergence rate in the low SNR setting, while its computational cost per iteration (dd).
Acknowledgements
The rescarch of Qiujiang Jin and Aryan Mokhtari is partially supported by NSF Award CCF-2007668 and the NSF AI Institute for Foundations of Machine Learning. Tongzheng Ren and Naht Ho acknowledge support from the NSF IFML, 201 9019844 and the NSF AT Institute for Foundations of Machine Learning.
12

---

Published in Transactions on Machine Learning Research (April/202))
References
A. Agarwal, S. Negahban, and M. J. Wainwright. Fast global convergence of gradient methods for high- dimensional statistical recovery. Annals of Statistics, 40(5):2452 2482, 2012.
S. Balakrishnan, M. J. Wainwright, and B. Yu. Statistical guarantees for the EM algorithm: From population to sample-based analysis. Annals of Statistics, 45:77-120, 2017.
C. G. Broyden, J. E. Dennis Jr., Broyden, and J. J. More. On the local and superlincar convergence of quasi-Newton methods. IMA J. Appl. Math, 12(3):223-22 22222,,, June 1973.
Charles G Broyden. A class of methods for solving nonlinear simultaneous equations. Mathematics of compulation, 19(9)):577 593, 1965.
Charles Broyden. The convergence of single-rank quasi-Newton methods. Mathematics of Computation, (110):365-8) 382, -8222, 1970.
Emmanuel J. Candes, Yonina Eldar, Thomas Strohmer, and Vlad Voroninski. Phase retricval via matrix completion, 2011.
R. J. Carroll, J. Fan, I. Gijbels, and M. P. Wand. Generalized partially linear single-index models. Journal of the American Statistical Association, 92:477-498 489, 1997.
Andrew R. Conn, Nicholas I. M. Gould, and Ph L Toint. Convergence of quasi-Newton matrices generated by the symmetric rank one update. Mathemalical programming, 50(1-33:117-199)) 1991.
C. Daskalakis, C. Tzamos, and M. Zampetakis. Ten steps of EM sullice for mixtures of two Gaussians. In Proceedings of the 2017 Conference on Learning Theory, 2017.
WC Davidon. Variable metric method for minimization. Technical report, Argonne National Lab., Lemont, MI.. 1959.
R. Dwivedi, N. Ho, K. Khamaru, M. J. Wainwright, M. I. Jordan, and B. Yu. Sharp analysis of expectation- maximization for weakly identifiable models. ATSTATS, 2020a.
Dwivedi, N. Ho, K. Khamaru, M. J. Wainwright, M. Jordan, and B. Yu. Singularity, misspecification, and the convergence rate of EM. Annals of Statistics, 44:2726 2755, 2020b.
Xiaoming Chen Feiyan Tian, Lei Liu. Generalized memory approximate message passing. https:/ /ttp:///rrv/o/oe org/aba aag/abs/// 221 20.00000000000010.1 06066 2021.
J. R. Fienup. Phase retrieval algorithms: a comparison. Appl. Opt., (15) Aug 1982. doi:
10.1364/AO 21.002758. URL http:///www..sappliihhiggg arg////bbstaaattttttttooo ot http:///.........aa.ccc.c.c.oooa
Roger Fletcher. A new approach to variable metric algorithms. The computer journal, 13(3):317-33) -22,, 1970.
Roger Fletcher and Michael JD Powell. A rapidly convergent descent method for minimization. The computer journal, 6(2):163-168) 168, 1963.
David M Gay. Some convergence properties of Broyden's method. SIAM Journal Numerical Analysis, 16 (4):623 630, 1979.
Kazimierz Goebel and W. A. Kirk. Topics in Metric Fized Point Theory. Cambridge University Press, 1990.
Donald Goldfarb. A family variable-metric methods derived by variational means. Mathematics of compulation, 24(109):23 26, 1970.
Moritz Hardt, Ben Recht, and Yoram Singer. Train faster, generalize better: Stability of stochastic gradient descent. In Maria Florina Balcan and Kilian Q. Weinberger (eds.), Proceedings of The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pp. 1225-1234, New York, New York, USA, 20 22 Jun 2016. PMLR. URL http: http:iinswsml....pesssseee html
mlr, wwwpppeesseeeeee v48 2ress//488harrt16 hardtt66
13

---

Published in Transactions on Machine Learning Research (April/2024)
N. Ho, K. Khamaru, R. Dwivedi, M. J. Wainwright, M. I. Jordan, and B. Yu. Instability, computational efficiency and statistical accuracy. Arriv Preprinl Arain: 2005. 11411, 2020.
Qiujiang Jin and Aryan Mokhtari. Non-asymptotic superlinear convergence of standard quasi-newton methods. arXiv preprint arXic:2003.138 13607, 2020.
Qiujiang Jin, Alec Koppel, Ketan Rajawat, and Aryan Mokhtari. Sharpened quasi-newton methods: Faster superlinear rate and larger local convergence neighborhood. The 39th International Conference on Machine Learning (ICML 2022), 2022.
J. Kwon, W. Qian, C. Caramanis, Y. Chen, and D. Damck. Global convergence of the EM algorithm for mixtures of two component linear regression. In Conference on Learning Theory (COLT), 2019.
J. Y. Kwon, N. Ho, and C. Caramanis. On the minimax optimality of the EM algorithm for learning two-component mixed linear regression. In AISTATS, 2021.
Dachao Lin, Haishan Ye, and Zhihua Zhang. Explicit superlinear convergence of broyden's method in nonlinear equations. arXiv preprint arXiw:2109.012 01974, 2021a.
Dachao Lin, Haishan Ye, and Zhihua Zhang. Greedy and random quasi-newton methods with faster explicit superlinear convergence. Advances in Neural Information Processing Systems 94, 2021b.
Dong C Liu and Jorge Nocedal. On the limited memory BEGS method for large scale optimization. Mathematical programming, 45(1-3):503-522 528, 1989.
Po-Ling Loh and Martin J Wainwright. Regularized M-estimators with nonconvexity: Statistical and
algorithmic theory for local optima. Journal of Machine Learning Research, 16:559 11659-616, 2015.
Haíhao Lu, Robert M Freund, and Yurii Nesterov. Relatively smooth convex optimization by first-order methods and applications. In SIAM Journal on Optimization, 28(1):333 354, 2018.
Wenlong Mou, Nhat Ho, Martin J Wainwright, Peter Bartlett, and Michael I Jordan. A diffussin process perspective on posterior contraction rates for parameters. arXiv preprint arXiv:1909. (9X:::1999,00666, 2019.
Praneeth Netrapalli, Prateek Jain, and Sujay Sanghavi. Phase retrieval using alternating minimization. 1EEE Transactions on Signal Processing, 63(18):48144488) 2015. doi: 001100//SPP2015.244448
Jorge Nocodal. Updating quasi-Newton matricos with limited storage. Mathematics of computation, 35(151): 773 782, 1980.
Tongzheng Ren, Fuheng Cui, Alexia Atsidakou, Sujay Sanghavi, and Nhat Ho. Towards statistical and computational complexities of Polyak step size gradient descent. Artificial Intelligence and Statistics Conference, 2022a.
Tongzheng Ren, Jiacheng Zhuo, Sujay Sanghavi, and Nhat Ho. Improving computational complexity in statistical models with second-order information. arXiv preprint arXin:2202.04219 2022b.
Anton Rodomanov and Yurii Nesterov. Greedy quasi-newton methods with explicit superlincar convergence. SIAM Journal on Oplimization, 31(1):785 811, 2021a.
Anton Rodomanov and Yurii Nesterov. Rates of superlinear convergence for classical quasi-newton methods. Mathematical Programming, pp. 1-32, 2021b.
Anton Rodomanov and Yurii Nesterov. New results on superlinear convergence of classical quasi-newton methods. Journal of Optimization Theory and Applications, 18843))74,,769, 2021c.
David F Shanno. Conditioning of quasi-Newton methods for function minimization. Mathematics of
compudation, 24(111):647-65) 1970.
14

---

Published in Transactions on Machine Learning Research (April/2024)
Yoav Shechtman, Yonina C. Eldar, Oren Cohen, Henry Nicholas Chapman, Jianwei Miao, and Mordechai Segev. Phase retrieval with application to optical imaging: A contemporary overview. IEEE Signal Processing Magazine, 32(3):87-11) 109, 2015. doi: 10.1099/ASP-221442555777
7. Wang, Q. Gu, Y. Ning, and H. Liu. High-dimensional expetation-maxiiiaation algorithm: Statistical optimization and asymptotic normality. In Advances in Neural Information Processing Systems 28, 2015.
J. Xu, D. Hsu, and A. Maleki. Global analysis of expectation maximization for mixtures of two Gaussians. In Advances in Neural Information Processing Systems 29, 2016.
Haishan Ye, Dachao Lin, Zhihua Zhang, and Xiangyu Chang. Explicit superlinear convergence rates of the srl algorithm. xiv SisprimlaaXix:210000 071612, 2021.
X. Yi and C. Caramanis. Regularized EM algorithms: A unified framework and statistical guarantees. In Advances in Newral Informalion Processing Systems 28, 2015.
Xiao-Tong Yuan and Tong Zhang. Truncated power method for sparse eigenvalue problems. Journal of Mochine Learning Research, 14(Apr):899 925, 2013.
15

---

Published in Transactions on Machine Learning Research (April/2024)
Appendix
A Proofs
Lemma 1. Consider the objective function in of its Hessian f(0) can be expressed as
10
satisfying Assumption and
Then, the inverse matrix
(A A)-1 -²5(0)-1= 9l1
( () 0)
(18)
b/
q(q 1) b11
Proof. Notice that the Hessian of objective function
10
be expressed as
f0) q/| -
4(4 - A) 6 AA" (A0 - () b) A.
(19)
We use the Sherman Morrison
Sherman-Morrison formula. Suppose that X Rd x is invertible matrix and a, E Id are two vectors satislying that 1 + b x # 0. Then, the matrix X + ab is invertible and
ab 1
(x + ab ) --==
Applying the Sherman-Morriset formula with x= A0 - - and b A (A0 b). Notice that A A is invertible, hence is invertible and
AT A, a= q-g A0 bl b)
1 x-1a
-T)0-)4AA A (A A)
1+ (A0 -
a|I bl19-2 044 - 2) |0 b/1 AT (A0 - b)
1+ (g (A0 6)T (AT AA A(O 0)
||A0 - -112
(A0 T)T (A0 b)
+ (g 2)
|0 b/2
= g 1 O.
(Since q 2 4.)
Therefore, we obtain that
V ((0)
(ATA)
A0 61 14 2)11 - 6115-1 AT (A0-
(ATA)) A)-
ATAO -100-b)) ()(AT((0 T))
(ATA)
1140-619-2
g||A0 b/10-2
9-1
A)-1
(4 2)
(20)
(AT A) -AA AAT (0 0) - T)T AAT(AA AA-
A|| 0110-2
4(9 1)] A0 bll0
(AT A)
2) (0 0)
A0 b1 q(q 1) A0 b1ld
As a consequence, we obtain the conclusion of the lemma.
Lemma 2. Banach's 's Fized-Point Theorem. Consider the differentiable function f DC R. Suppose that there exists C E 10,1) such that [f'(x) S C. for any x E D. Now let To E be arbitrary and define the sequence ( as
Øk+1 f(xx), Wk 2
(21)
Then, the sequence (+
converges the unique fixed point x. defined as
x= = fx.),
(22)
with linear convergence rate of
Vk 0. Ia x.1 S C 0+1200 120 x0/
(23)
Proof. Check Goebel & Kirk 1990
16

---

Published in Transactions on Machine Learning Research (April/202)
A.1 Proof of Theorem [
We use induction to prove the convergence results in and 12 Note that b b= AO by Assumption the gradient and Hessian of the objective function in
and
10 are explicitly given by
Vf(0) <l1 b (A0 b) all 6 A(D - 0), -0) - A + a(4 2) - b (A0 (A0 b) A.
(24)
(25)
Applying Lemma
we can obtain that
(A A)-1
(1 2) () 0) T
(26)
9|| b119-2
q(9 1) b11
First, we consider the initial iteration
81 8 Ho ((oo) 0 Vf(06)-\vf(())
01 0 0o 0 vf(0)) (V(()))
Notice that Aθ by Assumption
and
((0)) (VS
(ATA)-
- () - ()0 - 0)T
<ll - 21 A(00 - 0)
2ll A0o b110-2
<(g 1)ll A00 - b/10
9 2 (0u
-0)T A(Ou- 0) (00 0)
0o-
(27)
9 - 1
||A0 - b112
0 Q 2 (AOp
- TTTT))())))) (A0o (0o- 0)
0o
|| A0 b11
- 2
0o
(0o
0)-
q
Therefore, we obtain that
81 - 0 0u - v((o)) (V(o)) q 9-2 (0u 8))
Condition 11 holds for Ao = 1 with ro
9-2 q-1
Now we assume that condition
holds for k = / where
2 1, i.e.
A. 0 = r r 8).
(28)
Considering the condition b Aθ in Assumption and the condition
28 we further have
A0 = A(o; 0) = - - 0) = (A0 b),
which implics that
v)(00) = - 2 A(0,11 -)-
(29)
We further show that the variable displacement and gradient difference be written as
- 0--0-00-1 + 11)00-11 O), 8t_1 = 0 -
and
U4-1 = ((,)))) - = = - - AT ATT,((111 - 0). Considering these expressions, we can show that the rank-1 matrix the update BFGS 11118 given by
11-18-1= = ()(r-1 1) A(0, 0 - A(O -0)(0-11-0)) -
17

---

Published in Transactions on Machine Learning Research (April/202)
and the inner product si 14-1 can be written
T 4-1 = -(1p- ()(r=11 1)]1-(0:-1 -)
A(D(-1 -0)
= = - ([(22-11)]]A((())
These two expressions allow to simplify the matrix 1 " 1861 in the update of BFGS 81-1
A A(00-1 (0,-111 T)
Up 18-18-11
(30)
Uc-1
||A(0,-1 0)||²
An important property of the above matrix that its null space is the set of the vectors that are parallel to Ut-1- Considering the expression for ut-1, any vector that parallel to AA A(0,1 - 0) in the null space of the above matrix. We can casily observe that the gradient defined (29) satisfies this condition and therefore
1
16-18111
Vf((s)
81114--1
11-1
A 1 () -OT A A(0--1
0)
=
|A(0,-1 -0)12
(31)
AT A(0, 1
= -
A A(0(-1 1
|A(0,-1 -0)12
= 0.
This important observation shows that if the condition in 28 holds, then the BFGS descent direction H.V can be simplified
((b)(
81 11-1
H4-1 (1 (1-
(1-
Ut - 186-1
81-18 (((())
f(..))
Up-1
8 4
s 44-1
81- 81-1811 (0,)
4-1
(ri-1 - 1)3 () 1 T)
(32)
0019_ 1)(r:-1 A() 0)]10 -0)10----
00)11-227 ATA(0,-11 -0)10-2ATA((0-1 0)
r--1 / (0|-1 - 8) |1 A(O00 0)) -O)TAT A(O,-
1- [-]
|A 8)110
7- (0,-1 0)-
Tt- [
This simplificatinn implies that for the new iterate Oz41, we have
(0z 0)
8.11 θ =0 0, H.V)(())) - 0 0, -
T
g
r
Tt-1
(33)
1- 1 (0, 0) r0 0)
r[=1
where
1- 1- -
(34)
1 r
Therefore, we prove that condition results in 11 and 12
11 holds for l: + 1. By induction, we prove the linear convergence
18

---

Published in Transactions on Machine Learning Research (April/2024) 2024)
1.-
1.
08
Newtom
Newton
Newton
Newfon
Number of Iteration
Number of Iteration &
Number of Iteration
Number of Iseration
(a) 10 4.
(b) d 10, 10.
(c) d= 10", -0=4.
(d) d= 1dd=000gg== - 10.
Figure Values of cos 02 generated by Newton's method for different values of d and q.
8 8 8
BPGS
BFGS
BFGS
BFGS
Number Iteration
Number Iteration 3
Nomber of Recation
Number of Heration
(a) 10.4 =
(b) 10, 10.
(c) d= 10 .9
(d) d 10 .9
Figure 7: Values of cos 0g generated by the BFGS method for different values of and q.
One property of this convergence results is that the crror vectors (0x - 0 are parallel to each other with the same direction shown in [1].. This indicates that the iterations {0x] converge to the optimal solution i along the same straight line defined by 0p - ô. Only the length of each vector 0g - Ô reduces to zero with the lincar convergence rates (rx fm specified in (12 and the direction remains all the same. Notice that this is not common property for BFGS and Newton's 's method. This properly requires that the initial Hessian approximation matrix is Hg 2 f(06)= and the convex problem must be quadratic problem (T0). When the objective function is general convex function the initial Hessian approximation matrix is not Ho = 2 f(0o)- there is no guarantee that the error vectors {A,- 0) are parallel to each other. To better visualize this property, Newton's method and BFGS in Figures 0 and with population loss function defined in (10b. We observed that all values of cos 0g one, which indicates that all vectors (0x. 0 a parallel to each other.
have plotted cos 0g (0,- 0) (fg 0)/ 0|1|, O1 of k > 1 for both
A.2 Proof of Theorem
Recall that we have
2
Tk =
Vk 1
(35)
q
1
Consider that q 2 4 and define the function g(r) as
g(r) 1 - 2-2
re [0,
(36)
Suppose that r- E 1) satisfying that r, = g(r=), which is equivalent to
r + 2-2
(37)
Notice that
g((r) = (g - 1)rl-2 22g-4 (g - 2)r
(1 - r 111
19

---

Published in Transactions on Machine Learning Research (April/202)
(a) 0 = 4.
(b) q 5.
(c) = 6.
(d) y = 7.
04
04-
1900
(e) 8.
(I) 9.
(g) 10.
(h) g 11.
Figure Plots of [g'(rr] with r e [0, for 4 45
and
(g - 1)r9-2 - 220-4 - (g 9-3
p0-3(10- 1)(r 1) (1 1)]
3r 1)(1 78-1 1
9-2
,- 1)(1 1 -
Σnr). rl).
1=
Since [0, 1], we have that
Q-2
9-2
19-3 20, 5-100,
Σ
S
= q 1.
1=0
1=0
Therefore, we obtain that
(q 1)r0-2 2 220-4 - (q - 933 <
and
lo((r)) 22 (q 2)r - ( 1)r9-2
(1 - 1--1)2
Our target is to prove that for any [0.
(0((r)) < 2
First, we present the plots of lg' (r)) for r E [0, 1] with 4 < g S 11 in Figure
We observe that for 4 5 S 11,
lg((r) s 1/2 always holds.
Next, we prove that for q 2 12 and any r E [0, 1], we have
[g((r) = (q - 1)r9-2 2 22 - (q 2)
(1- (1- 19-1
2
20

---

Published in Transactions on Machine Learning Research (April/202)
which is equivalent to
20- 2r22 - + - 11)g-1)r--2-- - 2(4
2), +12 0, e
(38)
Define the function h(r) as
h(r) := 222-2 2222 - 2, + 2(g - 1)y9-2 - 2( - 2),
(39)
We obtain that
dh(r) 2, 1998-1100000
(40)
dr
where
(r) = (q 1)r
2(g 2)r9-1 69 1)r² + (q 1) 2)r (g -q 3).
(41)
Hence, we have that
dh (r) (g (r),
(42)
dr
where
(r) = (1 1)r9 - 2(q 2 2)r 2r + -2
(43)
Therefore, we obtain that
dh ( = (r) = (g 1)qr -2(1 2(g-2((22))) -2)2y0-33 -
(44)
dr
and
(r)
dr = (r-))(00(r),
(45)
where
(r) 7 1)(q 1)y² -2
2) 3).
(46)
Now we define the function ((q) as
I(4)= == 2(g () - a(a + 1)(g 1)
= q² 149² + 33g - 24
24 - 14) + -3(g 1) +
We observe that for 2 14, we have ((q) > 0 and we calculate that ((12) = 84 > 0 and ((13) 236 > 0. Hence, we obtain that ((q) > 0 for all 2 12, which indicates that for all r E [0. 1].
2(g 2(q-2))(( - 3)
,² < 1
q(a 1)(9-1))
a(g + 1)(q - 2) - -3) (3){3 < -
Therefore, for r [0, (r) defined in 46 satisfies that (r) < 0 and from 45 we know that
(r)3)() < 0. Hence, (1() (r))))) defined in 44 is decreasing function and (1) (r) <= (3) (0) 2 0. We know dr
that (r) (r) < 0, which implies that (r) defined in 43 is decreasing function. So we have dr that (((2(()) 2 (1) = > 0. From 42 we know that (r > 0 and (()(r)) defined in 41 is increasing
dr
function for E []. Hence, we get that (r) x (10)11) = 0 and from 440 we obtain that h(r) defined in is decreasing function for r E [0, 1]. Therefore, we have that h(r) > for all r E [0, 1). which is equivalent to (g'(r) < 1/2.
0 and condition
35 holds
In summary, we proved that for any q 2 12, we have [g'(r) S 1/2. Combining this with the results from Figure 8 we obtain that g (r() < 1/2 holds for all q 4. Applying Banach's Fixed-Point Theorem from
Lemma 2 we prove the final conclusion (4p.
21

---

Published in Transactions on Machine Learning Research (April/222]
A.3 Proof of Theorem 3
Notice that the gradient and the Hessian of the objective function
10 can be expressed as
Vf(0) 4| b1| - A (A0 b) gll - 6110- 2 A A(B 6).
((0 = q 2 + q(q - 2) A0 TT (AB () b) T A.
Applying Lemma we can obtain that
(ATA)-
2 ( 0) T
-11 b110-2
9/9 1)]1 - b11
Hence, we have that for any h 2 1,
02 02 1 - 10-00)1)~~
0g 8- ON 1 - 0 -
Notice that b == A0 by Assumption and
(A A)-1
(q (0 () 0)T
9(q 1)1 - b118 19||A0 - 21 100-2AA(0k--11
A(0k-1-0 1
A|1A0R b1|--2
0 0- 2 (0o - -0) AT A(0s-1 0) (0x -0) 0)
9 1
||A0 - -bl2
0 1 θ 9 2 (A0x-
1 b) (A0k(ARk)))))))) (0x-1 8)
1
| - b112
2
= 0k- - 9 1 (0k-1 - -0).
Therefore, we prove the conclusion that for any 2 1,
0x - 0 = 02-1 - 0 v2 f(0x-) 'V(00--1) 9-2 (0k-1 (0--1-0)
We observe that the iterations generated by Newton's method also satisfy the parallel property, ie., all vectors (0p- - _ are parallel to each other with the same direction.
Notice that the function h(r) = r2-1 + p9-2 is strictly increasing and h ( < 1, h(r.) 1 as well as
h1 3 > 1. Hence, we know that 2-2 < r. 2
2--2"
A.4 Elaboration of Remark 1
For the ease of presentation, we use Cp to denote any constant that is independent of d, n, and Cp can be varied case by case to simplify the proof. From Mou et al. (2019) and Lemma 6 of Ren et al. (002bb), we have that, as long as n = lod d/8)2)), we have the following two uniform concentration results holding with probability 1 d:
+ +r) r)--1 (V - VC(0)) =0,00
Valoog \11og(
sup DER(U' .))
(47)
(V)),
sup (c)(
Notice that we have
((0)= E((Y- ()))) = E[( + (X (XTP)
E[((XT 0')P - (XT0)))] + + -0 (XT0)))
=E[((XT0') 10-00 - ( +02
22

---

Published in Transactions on Machine Learning Research (April/2024)
where use the fact that Ç is independent of X and E[{] = 0, E[5²] = 0². Hence, we have that
VL(0) - 2pB[((X 0) - (x T (x)) O).
(48)
V=C(0) 2x'))= 2p(p - 1) 8) (x (x -0x
(49)
We denote _ the population loss function with respect the assumption of 0* 0. Therefore, we have =0 = B[ +0².
(50)
(C) 2x (0))
(51)
(2 = 2p(2p - 1)E () 20) xx
(52)
Hence, we have
Recall that X is a Gaussian or sub-Gaussian random variable with |||| बل < +oc For any 0, we alve that ||0| || $ ||0*1 |00
0*l and in the low SNR regime, have ||0"1| - ( मं Hence, have
sup (EB(O*.))
[V) (0)I sc,(1100 + + (algg///))
(53)
Similarly, we have that
-22 -CE|(X (0 (0-XXX000-xxX]11 "F x
sup (E
[V2 - v2c 5cppl +r)p- /logg1/6)/n
(54)
Leveraging 47
53
and
54
we obtain that
sup BEB(
(O) - VC' +))00~11)))--------11 /llg(//////////
(55)
(o)) ( (1) (olg+)/))))//og1////1)
sup (E
This explains the Remark of assumption 0" 0 in section The errors between gradients and Hessians of the population loss with 0 = 0 and ||0*|| < (d/nn) are upper bounded by the corresponding statistical errors between the population loss and the empirical loss (7) in the low SNR regime, respectively. Hence, _ and L can be treated as equivalent.
A.5 Proof of Theorem
In this section, we present the proof of Theorem We use L and C, to denote the population objective and empirical objective. Also, as discussed in the previous section, we use L to refer to the population loss when the optimal solution is zero 8* = O. For the ease of presentation, we use Cp to denote any constant that is independent of d, n, and Cp can be varied case by case to simplify the proof. To prove the result for the iterates generated by the finite sample loss function, we control the gap between the gradient of show that in the low SNR setting, we have
Co and Ln- As discussed previously in (5)) from A.4
sup DEB(Or.r)
[VL. VL
(C, <C,([100 + +) log (1g/)//)
(56)
sup (E3(0*.))
|V
1//
(57)
23

---

Published in Transactions on Machine Learning Research (April/2024)
We further can show that ||0-|| - |10 -1 for any 0. If this does not hold, then we have ||0 8 < ||0*11 N
(1 to by the definition of the low SNR regime, which indicates that has already achieved the optimal statistical radius. Therefore, the following bounds hold with probability of at least 1- 6,
sup 0E3(0",r) (0) (O)Il SC " VVlog(116)/m)
(58)
sup (EB(O'.))
- (0) SCpr //log(1/0///
(59)
In the following proof, we assume that all the results hold with probability of at least 1- Given the fact that the difference between the gradient and Hessian of the finite sample loss, denoted as C., and the population loss with the optimal value at zero, denoted by L is of the order of statistical accuracy, and considcring that the iterates generated by running BFGS on _ converge to the optimal solution at linear rate, show that the iterates generated by running BFGS the finile sample loss reach the statistical accuracy in a logarithmic number of iterations of the required statistical accuracy.
To simplify our prosentation, we use [0e) as the iterates generated by running BFGS on the population loss _ and {0T h220 as the iterates generated by running BFGS on the finite sample loss L.. We assume that 0g = 0% i.c., the first. iterates are the same for both population loss and empirical loss. Similarly, we denote the inverse Hessian approximation matrix of BFGS applied the loss L by {Hii)10 and the inverse Hessian approximation matrix of BFGS applied to the loss fn by (HF) We have that Ho = v² and Hy LzL, (08) Given the results in Theorems 1 and 2 we know that the iterates generated by BFGS on L converge to the optimal solution at a linear rate, i.e.,
0+41 0- rr 0*),
0cn O<nSn STa 1
Vt V2 0,
(60)
where {ra) are the linear convergence rates and rirh E (0, are the lower and upper bounds of the corresponding linear convergence rales.
More precisely, we show that if the total number of iterations T that we run BFGS on Cn is order of log of the inverse of statistical accuracy at most, i.c., T= (/lloogn/d)), then for any t < T we can show that the difference between the iterales generated by running BFGS on L and bn controlled and bounded above by
| - 0el S cell0--1 - 11-11-1 Vdlog
(61)
where = (eex((t)) Now given the fact that |10,-1 - 0*ll approaches zero at a linear rate, we will use induction to prove the main claim of Theorem 4
We assume that for any t < T, we have that
1 Cp
col10,
1.
(62)
Otherwise, our results in Theorem 4 simply follow. To prove the claim in 61 will use induction argument. Before doing that, in the upcoming sections prove the following intermediate results that will be used in the induction argument.
Lemma 3. We denote (A) and (min(A) the largest and smallest eigenvalue of the matriz A, we have
[VL < Cpoll 8- 11201 Amare (V² Cpt 10
0 1|20-2
AA E
(V) - (0) S Cn1 0'11(110 - 0- + 110 10 20- V0, 0' E RR
Amin (0) 2 C22l 20-2
Amin (V2 (ooo)) 2 Cpal10p- 10-2
where Cpo. Opi, Cp2 and Cps are constants that only depend on p and the last inequalities only holds for Bo-
Proof. Recall that from 51 and 52
we have
VL (0)
E
(00 1x),
F2C (:00))22(2pppp))) (
24

---

Published in Transactions on Machine Learning Research (April/2024)
Using Seechy-Schwazz inoquality, we get that
[ = (x]| 2) X1|| - |||||[[|| 20 xi|| - || 2--1 || |
I 2 10 + | 2- - || - 0 11) - Cpoll 8 1|
where we use the fact that ||0"1 || = ||0 11 for any 8. If this does not hold, then we have ||e 8' < 10 = (1()) 2u by the definition of the low SNR regime, which indicates that 0 has already achieved the optimal statistical radius. Similarly, we have that
Amax (V)C - (O)] 2p(2p- - 1)]l 0) - III - 2 1XI00 II S 2p(2p- - 1)B[ X | | 11]]00000 2p|22 < 2p(2p 1) [| || (110'' || + 110 - 0" 11))) - 1p ||0 0-1120-2,
-
where we use the same argument that |00|| S 10 0-11 Using Taylor's Theorem, we have that
(C) - (0)= = (-0 + (1 -1- (1- 0¹)
where T E [0, 1]. Therefore, for any 0, 0, have that
- - - 17 +)0 - 0' - Cp (1 T)0' 0- 2 0'l 0'1 - Ctt 10 0' [110 0-1 + 10 0 11
=
Notice that (0) are positive definite for any 0 and V Ln (0)) is positive definite, hence there exists Cp2 and C33 such that Amin (O) 2 C22 0 0-122p-2 for any 0 and Amin (V² (00)) 2 Cpa 50 0" 0-112--2
A.5.1 Induction base
Now we use the induction argumentt to prove the claim in 66 Since the initial iterate for running BFGS on C and L, are the same, we have Og = fn- Now for the iterates generated after the first step of BFGS, we can show that the error between the iterates of two losses
10 0:ll is bounded above by
- = 100
- (v² (0))
(63)
C. -200000-~~~~000~~~0>>00>>>>- + | (66) (0o) (061l
We observe that for invertible matrices A and B, we have (A-1 - - observation and the results in Lemma 3 and (59 we can show that
= A- 1-1BB-al A)B
Given this
I 2 S IIV (0.) 'IIV C. (0n) 42 --="()))
S 0 1-1 / /dlog(1/)))))
Applying this bound into 63) and using the results in Lemma the bounds 158 we obtain that
|00 - 011l S Cp - 111| 11 V (lo((1/6)/n
(64)
Notice that from 64/. we know that 61 holds for t. t 1. Hence, the base of induction is complete.
A.5.2 Induction hypothesis and step
Now we assume for any k < 1, 61 holds, i.e., we have 100 Oxll
-0'11--- /llgg(///1 00110 Consider
k + 1, we have that
|0 - 0-- = ||0|- - H, - 0, + H.VC (0.)l < 100 0.ll + VH, -H,-- (65)
Recall update rules for the Hessian approximation matrices:
" H [
n n/ (3
+ 8 " T
H = [
(u._ T sT 2/
(sp_1) 2/t
(sp_,) أاا
25

---

Published in Transactions on Machine Learning Research (April/2024)
Hr- = 1 ((-1(((11)
-1-1((-1)
+ 8 T
Hz-1 I
(u--1) 81-1
(8--1) T u/-1
(8+-1) 11-1
s.-1 00 - 00--1 - = (OT) - (o
Sc-1 = 01 - -8 4-1=VL - (0-1))
Moreover, based on 31 we have
I u6-18/-1
= 0.
(66)
-
which implies that H. (0) can be simplified as -ل1 VL Now we proceed to bound the gap ,,-1
between the BFGS descent direction on Cn and L which can be upper bounded by
$ $1-18|11 (
V - H.VL
= H (o'))
8-1
$
H;
[ (
uft1 (s"-1)
(0%(
||
St ( (0,)
(uf1) s T
1s"_1) I) u
(s ) u
8 1
(67)
Now we bound these two terms separately. The first term in 67 can be bounded above by
$ (up-1)
uy (s"[1)
H
I
L (833)
u T 8
(s
(8
sp_ (up
(
u 2 (8
|
V V (07)
(68)
(b)-)) T uu
[ 2 (8
S
| ill
VL (87)
(sp_)) ) 26-
|
a (ut ,)
S 1. Now using the fact that (I- 14-14 (B)) 0, we
where we use the fact that can further upper bound
u1-
||
(
(o by
" (sb_1) (L
-= |
11-18|_1
(09) I
(577 ) uP
(s) up
1 4111
2ef 1 (sp_,) ) uuf11
Cn
up 187 T (0.)
< /Ln (0.)) - <C.(ll)1l+ +
Ue-1
(69)
+ u 1 (B,)
m M-1ST_1 V (0,)
11144-1
$ 111 41
TVL (s)
ST (()) 111
ST V (0.)
< (0') - ( + |لاي|||
+ [|u6_
u-11l
8 Wr-1
Next we proceed to bound the second term in
67
The second term can be bounded as
s (s) V..
S1-18,11
VL (0:)
VL
(s] T ui=1
81-1 U
s (sp 1) T
$ S - v - (0.)
shh S T
$1-18T (0.)
V T. (8th
+
1 (O.)
(s]_) T u
8 11-1
11114-1
$ u1
(sp T V (8
ST V 0o (B2)
s 1 ||
(()
|
+
+|
(70)
(s7_ w/11
s 1 1
s | यج-
26

---

Published in Transactions on Machine Learning Research (April/2024)
Putting together the upper bounds in
68
69
and
70
into
67
obtain that
HH - -,,
VLn (6f)
sT_ (0:)
< 1 VL. VZ (11 11llll1-11
1|
(s%-1) up-1
11-14-1
(71)
+ [[|-@|-|lll]]_ u. 1|8--1 و |I)| أ (
Wz-1
In the following lemmas, establish upper bounds for the expressions in the right hand side of 71 Lemma 4. The norm of voriable variation for the finite sample loss 5t-1 and its difference form the variable variation for the population loss St-1 are bounded above by
18--1- st-ll (cr 01-1) -0'11' /ill-V/dlog's /dlog(1/8//m
Moreover, this result implies that
|| 1|0-1 0-11-
Proof. The first result simply follows from the fact that
|| - -s--1l S 1000- 0ell + |, - -1--1l
(ce|(0,-1 - 0'11' c--1 - "ell-="\//loo((/// S (a + (1 -0'11'~~ \llg////8)/
(ct + --1)ll0,-11 / log(l//))/
where we used the induction assumption 61 and the linear convergence results (60) for the iterates generated on the population loss. The second claim can be also proved following
| || S ||: || | s--1ll 1ll 110, 0- (0, - 0*) 100
0.l + 110 0--1l
$ - - ||0 10 + c - 1/0
0- 112 - /dlog(1/0)//
$ -1 0 + 1 1/rr 10 - 00 s -
C
where we used the induction assumplion 61 linear convergence results (60 and the assumption (62 that all0,-1- - /ololl//ooo(//)//1 Jr with sufficiently large Cp make the last inequality hold.
The gap between the population loss and its finite sample version evaluated at the current iterates are bounded above by
Lemma 5.
- (C, 1p
0- -0 Vllog(1/0)/n
Moreover, this result implies that
VVC. (0")) S Cp1 1-110011
Proof. We have that.
(L,,(0)) - < 1C0000/00000 C 1011---- Vdlog(1/0))
27

---

Published in Transactions on Machine Learning Research (April/2024)
where the first inequality is due to (58) the third inequality is due the induction hypothesis 61 and linear
convergence results in 60) and the forth inequality is due to assumption (62 and rt-1 S 1. We also have
(VV) V (0)) & 1, ( 0- || ||0, 0 1|
$ Cpel10--1 /d1g
0ell + 2 0- (]) 2p-2
Cpeille-11 0- " log (10))//////((ell0--1 0-1|"- / (llog( (10gg(////++++ 2r 2re-1110--1 0011)22-2
s -peal10-11 0'11'~0 log lot(l/0)/0((00-11--0 0*1 + 2 0-
S Sp41l00-1 8* 1 log((/6)/mm
where the first inequality is due to results Lemma 3 the second inequality is due to the induction hypothesis
61) and -100- 0*I S - +++ 0* the third inequality due the induction hypothesis 61 and
the linear convergence results (60) and the forth inequality is due to the assumption 62 and rz-1 1. The
first claim follows using the fact that
[VC, (O)I -
V (0.1)
$ - /dlog//)))
Given this result, the second claim simply follows from the fact that
[V,,( < (0)) v (0.1l
I Cp110 - 1- (Cp C,c₁) 10:1 0" |p- / (olo((//))// s Cpl10, 0-1100-1 Cp110.-1
27112-11
Cp
$ Cp 0-| 1-112-11 2p-1
110 0r 1-120-11 < Cp110: 0r 00122-11,
where we used results from Lemma the linear convergence rate in 60 and the assumption
62
that
a110,-1- /'loo C again.
Lemma 6. The norm of gradient variation for the finite sample loss ui-1 and its difference form the gradient variation for the population loss u are bounded above by
|u4-1 u-11l S (Op + Opec Cpct - -001--00---- 1-1pp-1 /log(118)))
Morcover, this result implies that
|u S ---100-1--01133- 0-1120-1
Proof. The first claim simply follows from the following bounds,
- u- (1 V
dlg(11/))/
Cp+ Cpe-1 Cpp-- - log 1/00//n
I (C, + 0004+-/000011100000/11111111/110
3 SC,, + + -
28

---

Published in Transactions on Machine Learning Research (April/2024)
where we used the results in Lemma and the linear convergence results (60) of BFGS on the population loss. The second claim simply follows from,
- - |1--1 u || + ||4--1|
$ (V - -
(
$ (Cp + (pa))l0 1:00|/0000001100--1000- 0'1'"' -dlgp= /dlog(1/0)/n (C, Cpc-1))l0-- - Vdlog(1/00
C,110; [--1[[[10, 00 + |0:-1 0 1))20=2
3 02110-11 - 0-1 1 + Cp 0 - 0r 11 1---1 Cp110, - -0" - - -0 ]1(11 - 0-1 + 1100011 - 0-11)20-2
Cp
< 12110-11 - 0 20
0,-1 0-1120-1
2p-1
22-
c1l(1 r2- (1/-1 8*) |1( 1-- 1110-1101111
10 + ||0 f 201))
Cp
-cp10 8*1 1-1 +
01-1
1-1120-1 + (1 - 11-1)(1 + T 20 11-1120-1
2p-1
ccp101-1- 2-120-1
where the third inequality is due to results in Lemina and Lemma 3
the forth inequality is due to assumption
62 and the fifth inequality is due to linear convergence results 60
Lemma 7. We have the following bounds:
(s (L, S V (0:) (Cp + Cper + C,e-1) 0 /
1( (ufhi) - - 4-1 & + Cpcc + --
Proof. The first claim holds since
- +4
$ (a + -10 00110-1-
+ 11001 1C, + 110
I (Cper + (:
(Cp + /
+ 1a0c=:)|000110001/////l0001//111 (C, Gpe) (Ca 10"1" 10'11" /dog(/ (1/8)/n
S (C, + +
where the second inequality is due to results in Lemma and Lemma 5 the third inequality is due to linear convergence rates (60) and the forth incquality is duc to rt-1 S 1. The second claim holds since
[ss;_)) (1 ST_ " || s -
111-1 u-1ll
+ (1444410)00000000000000010111110
10-1
+ - l 1M2- - Valog(////
(Cpa +
+ ||0, - - 61-1 + |
+ C 1-1
< (Gra +
+ (1 (1- (1-rr-))l00-11 -0,, -|
1 dlgg (188//n
< + 11 CCpt + Cptt - 10
< (Cp Cpa+ - (0ggg
where the second incquality due to results in Lemma and Lemma the third inequality is due linear
convergence rates 66
29

---

Published in Transactions on Machine Learning Research (April/2024)
Lemma 8. The following bounds hold:
-0-1 1s (1. - -
1s 151-14-11 2 Cpl10: - -0
[(sp_1) (u 2 Cpll0, - 0-112p
Proof. The first claim holds since
[ V. (0.) = 10, - 0 - = 100 - 0- - -0,-1 +00 VV) - (0, (0:00-0-0-0000000 -0)
T1-1
< 1
Tt-1 11 0- (II s
1
- 1110, -
0" 11,, 0r 1011000~~~00000000 S C,110, 0-1120
Fh
where we use the linear convergence results
60 and the results in Lemma
The second claim holds since
= (1- -
= =(0 - 0- - 01-1 -0- 0c-1)
(11
(0: 03) V (1- 0" 01-1 +0'
Tt-1
[1
(0, (1
+1- - 0')
Tt-1
2 (1
)1000 11 X Mmn (V2 (1 7)0,-1))
T=-1
1 110, -0- 0'113C + (1 (1-r)0;-1 - 1- 2012 2 Cp110, - 10 - 0*) + (1 () 0-)1
2 (1
2 0p 00 1(( + - (0: 0-)1| Cp(T+ +
1 - 2 - 0 21 > Cp 10 11-1
Tt-1
rh
where T E [0, 1] and we use the Taylor's Theorem, the linear convergence rates The last cliam holds since
60 and results Lemma 3
[<-_1) ui-1 Tup_ = (s|_1) 1 u/- 81-14--1 + + 2 1s 1 1(05_1) " "h[1
T 14-1
2 Cpl10 - 1-1200 (Cp + Ope + Cpct-)) - 8- / /log(1/5)/n
2 Cp110 0-122p (Cp Cpez) -
-110
1
110 01'1" d - 11 /dlog//)111/)))
2 C, - 8*| 1 ||0, 8* 11
||0 8" 1 2 C - -
where we use the second claim, the results from Lemma and the assumption 62
Lemma
We have the following bounds:
V (0c)
S Cp.
14-1
(sp_) T VLn (8 (op
ST (.)
5 (C + Cpet + -000
2
Proof. The first claim holds since
(0x)
T (
Cp 01-1 - 12
S C.
Uc-1
1t-1
Cp 0-11 8" 20
30

---

Published in Transactions on Machine Learning Research (April/2024)
where we use the results in Lemma g
The second claim holds since
VLn()) (L,(
8 VLLVL(00
(s/_1) ut-1
8, 1
(s%_1) ) VLn _ VC [0)]1s+] Up-1 + |s|
C (0.)lss-- (s%_2) u8-11 Tu6
1<07_ T "4 11-181111 1
(0p_1) ) (0T) 1 (0,)
|sT_ V CV
1s 11-111111 (8% T 16-111
+
(<(_1)] 1 Tup-1l
5/-_ 1 16-11
(s/-_) [T--) ub--1l
(Gp + Gpa Opct-1) )100)1101111 00 11' (log(11/)) 1/5)/n
(Cp + Cpet + Cp(t-1)l00t-1 1011 V (l0g()//)//
+ +Cp
Cp110,- - 0- 2
0p100,11 - 1-1
Cp + Oper + 1p0r-1 100-1 0- /dlo- (Cp + Oper + Opet ))10=1110110000 0- 1-P /dloog(1//)/1
S (Cp + Opec + C,a--) -
where we use the results from Lemma 7 and Lemma
Finally, we present an upper bound for the norm of the inverse Hessian approximation matrix He.
Lemma 10. The norm of the inverse Hessian approzimation matrix is upper bounded by
|| S 1-110-11- - -0-12-20.
Proof. Recall the update of H,",
H 8p (ufp1) TT) (uf-1) T s
H.a (I- up (s}-1))
+ 8 "
(s/ ) T 2
(s 1
(ااماب T
( 1)
With the property of I|
I 1 and ||
|| I 1, we have that
[w -i
( -,
s (ut T
-P-1
(sp_1) T)
+ +|| s (37
|
S ||
|||
(u. ) ss $/-1
(sh_1) T "
|| 112
||
(V' -1
< H
(sp_1) up-1
(5,7) (u,")
From results in Lemma and Lemma
we know that for any 0 0sist- 1-
| (w²)
110: 1-112
Cp1 1-1 S Cpll0, - -0'11-22 2p
Hence, using linear convergence results in 60
we have that for all 0 0sigi t 1-1, 1,
10 - 0-
11 0-1 = 1-10 0-
j=
110 1 2
r/
-+ + 101-+110-11111 - 0-|,
110: - 0-112-20 2p < <Th
(2p 2)(((1-0)) 2010-1-00100-1111111111111--1 -0012-227
Σ
10 - 1-112-20
(2p 2)((-)) 10,0 2-1 <
1
2p 2 10 -2112-20
h
1=0
h
31

---

Published in Transactions on Machine Learning Research (April/202))
Therefore, we obtain that
|
5 Amin
+
cr110 - 0-112-20
(u.))
L2 (0n))
1
(2p-2)((-11)
Cp
- 12 - 0r 22pp- + Cp
10 - 0'112-2p- 2
11:-1 - 0-112-20
2p
Cp
2p-2 1 r/
-p1101--1 - 8-112-20.
Hence, we have that
| & Cpl10,-2 -8-11--2
s Cp
2- T
1
|0,-1
1- < Cprh > 21001 0-112-20 < cc 10 - -1112-22.
With all the above lemmas, we can complete the induction hypothesis. Applying results in Lemma
5 6
and 10 into 71 we have that
H. V. (0 (.)]l
8
ST
(8)
< 1 Vl. - (O)Il + H1) [I]| 1- + (1=2-111)
(5%_1) [ u/-1
4
(C)
+ +[) 4- || +
+||1111111 || - 11-
8 4
$ 1-116111 87| 2- 20 (C, C2c) C 11 V /log/1//)//
+ (C,110:- 0- 2-12-200,,10---0 0-122--1 |10 0² () Opei + Cpe--1)l00-11 /olog//)
+ (C,110--1 0- 1|2-222 (2, + CGa, + -- (llooo//
+ (c++ 8- -
$ (Cpp + Cxe) ( +Cp Oper + -pe--))100-1-00111 10 /dlog (10)/1
dloog(1 (1/6)/n
+ (C,+ + CGp + C,0-1)10,,-1 8 -''''" /11og(////)
< (Cp + Cper + Cpc- )| - 0- 1
Hence, we prove that
VL, - H1, (1, < (C, + Cper -pct-1))0:-1100011 0000:00000000011>>//////0/1//11
dlog(1/0)/n.
(72)
Notice that by induction and linear convergence rates
60
we observe that
|10| - 0,11 - ce|10.-1 - 0- 1- -
10, - 8" 1-'1'-0 Vdlogg (lo0////)
ct 1-p
r(-
(73)
< Sar 1 10 - 0-1| -01o0//////oog(////// S -000 - 0r 01100/
p-
p /dlog(1/6)///
Therefore, leveraging
65
72
and
73
we have that
|10%11 0.+1l || < ||0," 8g
+ V.. - H.V
S ct - 0-1| ""1''>>//oog(//)//m (C,, + Cpes + Cpct- - (Vloog//))
$ (Cp Cpe Cre--))l(8
0
We define that
C+1 0+1 +
+ Cper + Cper-1-
Then, we have that
|10% - 0++1ll ct+1l10: -
With the standard recursion, know that aa (Cp)² for Cp large enough. Hence,
61
holds for 1.
32

---

Published in Transactions on Machine Learning Research (April/202)
A.5.3 Final conclusion
Therefore, using induction we proved that
661 holds for all 2
|0 - -
where q (C)) (exp(tt)) Notice that
|000 - 0- - 1000 0ell + 10: 0011 I1, - / + 110, - 1.
The optimal T with minimum |00 0*1 should satisfy that
-0100-001///0
Cp.
for which we obtain T = [dg(11//doo((1111 for some constant C that is independent of d and n. Therefore, we prove the final conclusion that
1(+1)
109 - 0*11 S C"(ddlog(1/8))/////rrrr
where C' is constant that is independent of d and n.
B
Additional experiments for the medium SNR regime
Here briefly illustrate the behavior of BFGS in Medium SNR regime. We consider the generalized linear model with d = 50, 100, 500 and p 2. The inputs are still generated by (X:) but 0* now is uniformly sampled from the sphere with radius n 1/6
The results are shown in Figure D We can see BFGS still converges fast, and the statistical radius of middle SNR regime lies between the High SNR and Low SNR. A rigorous characterization of the statistical radius of middle SNR regime will be left as future work.
Newton
Slopes -0.5556
BFCS
Heror
611
CD
log O0 011
log (0x
log(
GD constand
log 10
GD Potvak
Newton BPGS
Polyak
Number of Iteration &
Log of Sample Size
Number Iteration
Number of Iteration
50.
(b) d = 50.
(c) d 100.
(d) 500.
Figure 9: Convergence results and statistical results for medium SNR regime with 50 are shown in (a)
and (b). Convergence of different methods with d 100 and d = 500 for medium SNR regime are shown
(c) and (d).
33


Reviews:
[{"review_id": "zJPaSysVXs", "reviewer": "Reviewer_uJae", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths:**\n\n1. The authors approach a theoretically important problem of understanding convergence of BFGS algorithm on a GLM problem without strict convexity. The proof of the main results (Theorems 1 and 2) seems solid and easy to follow. \n\n2. The paper introduces the problem well explaining the background and the technical challenges. \n\n\n**Weaknesses:**\n\n1. The statement in the second sentence of the abstract seems contradictory. It says GD is optimal in high SNR regime, but it is later revealed that in the same regime BFGS is known to be faster, achieving super-linear convergence due to strong convexity. Of course, the iteration cost are different, but it should be clarified what is meant by optimality here. \n\n2. Although vanilla GD might converge sub-linearly on problem (10), it is known that mirror descent with a distance generating function $h(x) = \\frac{1}{q + 2} \\\\|x\\\\|^{q+2} +  \\frac{1}{ 2} \\\\|x\\\\|^{2}$ will converge linearly, see e.g., Lu et al. (2018). It would make sense to compare to such approach both in theory and experiments. Notice that the iteration complexity is merely $O(d)$ for mirror descent. \n\n3. The property of BFGS and the Newton method that the iterates $(\\theta_k - \\hat \\theta )_{k\\geq0}$ are always parallel to each other seems surprising to me (see the last paragraph of Appendix A.1). Is it a common property for quasi-newton methods? If not, it would make sense to elaborate more on this property in the main part rather than merely commenting in the proof. Some illustration in small experiment would help dispel the doubts. \n\n4. The proof of Theorem 2 in appendix A.2 is not complete. For $q \\leq 11$, only integer values are checked by the picture (Fig. 3). However, the property might fail for non-integer $q$. \n\n5. In Equations (70), (71), (72), (73) the inverse of the gradient is taken. I assume it's supposed to be a Hessian. In equations (104) and (105), the absolute value is missing for the last terms (are they positive?). \n\n6. Markers should be added to all pictures. \n\n7. I was not able to parse the proof in section A.4. Some revision of it would be helpful. The main problem is that it is 7 pages without any structure. E.g., splitting the proof into lemmas could be very helpful. \n\n8. On page 19, referring to a specific lemma in Ren et al. (2022b) would help.\n\n**Questions:**\n\n1. Is it possible to show global linear convergence of BFGS for a general convex function without the special structure (10)? \n\n2. Authors explain that middle SNR regime is technically more challenging, but is there any negative results for this setting? E.g., does the problem have spurious local minima or saddle points? I think it would make sense to draw connections with the literature on linear neural networks, e.g., Kawaguchi (2016). \n\n\n**Typos and other suggestions:**\n\n1. Page 26 \"Hign\". Page 18 \"inequation\". Page 19 \"simply the proof\"\n\n2. Consider removing the equation counter where it is not necessary. \n\nHaihao Lu, Robert M Freund, and Yurii Nesterov. Relatively smooth convex optimization by first- order methods, and applications. SIAM Journal on Optimization, 28(1):333–354, 2018.\n\nKenji Kawaguchi. Deep learning without poor local minima. Advances in neural information processing systems, 29, 2016.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The work studies statistical and computational complexities of the classical BFGS method applied to generalized linear models (GLM). In particular, the authors focus on low signal-to-noise ratio (SNR) regime, where the loss function is merely convex (not strictly convex). In this regime, the loss reduces to least squares problem with arbitrary power $q \\geq 4$, i.e., $f(\\theta) = \\\\| A \\theta - b\\\\|^p$, where $A^{\\top} A$ is invertible. Utilizing this special structure, the authors prove global linear convergence of BFGS method. Importantly, the linear rate is independent of the condition number, although only first-order information is utilized. The work also shows (faster) global convergence of Newton method for this problem and extends the result for BFGS to finite sample setting.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Addressing points 2., 3., 4. and 7. in the weaknesses above would be the most critical for acceptance.", "broader_impact_concerns": "---", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "j2uUfXeOTR", "reviewer": "Reviewer_nEq8", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This paper provides a nice analysis of the BFGS algorithm, which is usually hard to study, especially without line search. Although the proofs are quite computation-heavy, they are presented nicely enough to be understood and checked without too much effort. Thus, this article completely satistfies the criteria for a TMLR publication.\n\nMy main issue with the paper (outside of minor qualms in the presentation) is the bound on sample complexity that you obtain, which is in $n \\asymp (d\\log(d))^{2p}$. This seems extremely wasteful, since we can learn a polynomial of degree $p$ with many methods in $O(d^p)$ samples, and (since $x^p$ has information exponent 1) [this paper](https://www.jmlr.org/papers/volume22/20-1288/20-1288.pdf) implies that if we know the link function, $O(d)$ samples are enough. Further, since we need to compute gradients of $\\mathcal L_n$, I feel like the actual complexity of one BFGS iteration is $O(n \\cdot d^2)$, which makes the previous issue all the more salient.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies the time complexity of the BFGS algorithm in a specific generalized linear model. It shows that when the link function of the GLM is $x^p$, the algorithm converges to a radius of $n^{1/(2p)}$ around the true signal vector in $O(\\log(n))$ iterations, where $n$ is the number of samples.\n\nThe proof proceeds in two steps: first, it studies the BFGS algorithm for infinite sample complexity and zero signal, and shows a linear convergence rate. Then, for large enough sample complexity, it bounds the difference between the Hessian of the finite and infinite sample setting, and leverages those bounds to show convergence of the actual BFGS algorithm.\n\nThe authors also provide numerical experiments showing that the BFGS algorithm outperforms classical Gradient Descent, even with optimal step-sizes.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "All mentioned changes are mostly for improvement.\n\n- you should use `\\eqref`for referring to equations, and check some of the references (e.g. \"empirical loss equation 7\" would be better as \"empirical loss (7)\")\n- you mention extensively a notion of \"statistical radius\", that is never defined in the paper; it does not seem to be an actual lower bound for recovering $\\theta^\\star$, so what does it correspond to ?\n- do you have any lower bound for the complexity of learning a GLM with link function $x^p$, if we want to be independent from the condition number of $\\Sigma$ ?\n- top of page 5 : your obtained expression for $\\mathcal L(\\theta)$ seems convex independently from the norm of $\\theta^\\star$ (and its relation to $\\sigma$), so the explanation is a bit unconvincing.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "cRsrqhl0Km", "reviewer": "Reviewer_iVPs", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n**S1.** Suggested approach (BFGS method) is known for its good practical performance  for practical optimization problems, but until recently has been lacking good theoretical analysis. The idea of using this quasi-Newton method seems very reasonable and worth studying more in the machine learning context. The authors show that their approach improves upon Gradient Descent as its complexity does not depend on the condition number of the problem. Moreover per-iteration cost $\\mathcal{O}(d^2)$ is cheaper than Newton’s method $\\mathcal{O}(d^3)$.\n\n**S2.** The paper analyses performance of the method on both population (1) and empirical (2) loss with a focus on the low signal-to-noise ratio (SNR) regime. In the first (1) case a linear convergence rate $r_k$ for the iterates of the Algorithm is established. Later $r_k$ is analyzed more accurately and compared with a result for Newton’s method. For the second (2) case of empirical loss, the authors show that minimum distance between iterations and global solution also converges to a statistical radius after a certain number of iterations with high-probability.\n\n**S3.** The authors claim that this is the first global linear convergence analysis of BFGS without line-search for a setting that is neither strictly nor strongly convex. This may be interesting for a broader optimization community as the area of quasi-Newton method have recently attracted novel attention.\n \n**S4.** Overall the paper is well-written and structured, but the clarity and presentation can be improved by expanding and moving some of the explanations from the Appendix to the main part. More on this in the Weaknesses section later.\n\n**S5.** Numerical experiments are well-designed to support the main claims of the work. Several parameter settings are explored both in the main part and in the Appendix.\n\n# Weaknesses\n\n**W1.** Some of the claims made in the paper need to be justified (e.g. by providing a reference where a certain effect was observed)\n\n> the Polyak iterates still have instability during training due to the influence of noise in the models.\n\n> Since we focus on $p \\geq 2$ it can be verified that $L(\\theta)$ is not strongly convex in a neighborhood of the solution $\\theta^*= 0$.\n\n**W1.** **Presentation and clarity**\n\n- Separate notation for the loss $\\mathcal{L}$ and function $f$ to be optimized can be confusing.\n- Presentation of different settings and special problem cases is somewhat overloaded and lacks clarity\n- Second paragraph on page 6 is tricky to understand, especially the statements about error between gradients and Hessians. I believe that some of the Appendix details are needed in the main text.\n\n## Questions\n**Q1.** Why does the norm of the optimum $\\|\\theta^*\\|$ has to converge to zero for increasing the sample size $n$.\n\n**Q2.** Why high-probability analysis is needed for Theorem 4?\n\n**Q3.** I would like to ask for a more precise comparison to theoretical results of Newton method.\n\n**Q4.** How exactly the constant step size $\\theta$ is tuned “by hand” to achieve the best performance of GD on each problem? Why the “tuned” values are so small: {$10^{-4}, \\dots, 10^{-15}$}? Why in the experiments for empirical loss the step size is chosen as $\\eta = 0.1$?\n\n**Q5.** What are the main limitations of the proposed approach?\n\n### Minor\nAt the end of page 4 it is mentioned that $\\mathcal{L}(\\theta) := \\mathbb{E} \\mathcal{L}_n (\\theta)$ is\n> an approximation of that function by its population version\n\nIsn’t it usually the other way around that empirical loss approximates the population one with a finite sample Monte-Carlo estimate?\n\n**Typo** on page 7\n> problemequation", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper studies the problem of the generalized linear model (GLM), with a polynomial link function. Authors analyze statistical and computational complexities of quasi-Newton (BFGS) method with a main emphasis on the low signal-to-noise ratio (SNR) regime. This work shows theoretical and practical (experimental) benefits of BFGS over Gradient Descent and Newton’s method on synthetic problems.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**C1.** I would like to ask the authors to give a precise characterization of the problem being solved in terms of smoothness and strong convexity parameters. BFGS is introduced for the convex setting at the beginning of the paper. It is said that under Assumptions 1 and 2 problem 10 has a unique solution. How far then this problem is from the strongly convex setting? Does it satisfy Polyak-Łojasiewicz (PL) condition?\n\n**C2.** It would be helpful to add more details (maybe from the Appendix) at the beginning of page 5 to make the reasoning more rigorous and clear. It is unclear why the equations are split into multiple lines.\n\n**C3.** The paper would benefit from a more comprehensive overview of the technical innovations and analytical basis of this study, especially in relation to previous works on BFGS. A comparative analysis with prior literature, particularly focusing on the analysis and settings of BFGS, would provide valuable context and highlight the novel contributions of this research. Clarifying these points will help in distinguishing this study from existing works and underscore its originality and significance.", "broader_impact_concerns": "No for this work", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper investigates the convergence of the BFGS method in generalized linear problems. Unlike previous studies, this paper does not depend on the assumption of strong convexity. Instead, it leverages that the BFGS iterates remain within the linear span of the data vectors if the method is initialized in a specific way. This approach enables the derivation of novel convergence results, particularly in low signal-to-noise (SNR) settings.

This contribution makes an essential step towards a better understanding of quasi-Newton methods.
Decision: Accept (poster)