Statistical and Computational Complexities of BFGS Quasi-Newton Method for Generalized Linear Models

In supervised machine learning, we are given a set of $n$ independent samples denoted by $X_{1},\dots,X_{n}$ with corresponding labels $Y_{1},\dots,Y_{n}$, 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 according to distribution $\mathcal{P}_{\theta^{*}}$ which is parameterized by a ground truth parameter $\theta^{*}$. Our goal as the learner is to find $\theta^{*}$ by solving the empirical risk minimization (ERM) problem defined as

where $\ell(\theta;(X_{i},Y_{i}))$ is the loss function that measures the error between the predicted output of $X_{i}$ using parameter $\theta$ and its true label $Y_{i}$. If we define $\theta_{n}^{*}$ as an optimal solution of the above optimization problem, i.e., $\theta_{n}^{*}\in\arg\min_{\theta\in\mathbb{R}^{d}}\mathcal{L}_{n}(\theta)$, it can be considered as an approximation of the ground-truth solution $\theta^{*}$, where $\theta^{*}$ is also a minimizer of the population loss defined as

If one can solve the empirical risk efficiently, the output model could be close to $\theta^{*}$, 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 $\theta$ (Balakrishnan et al., 2017; Ho et al., 2020; Loh and Wainwright, 2015; Agarwal et al., 2012; Yuan and Zhang, 2013; Dwivedi et al., 2020b; Hardt et al., 2016; Candes et al., 2011). However, when we move beyond linear models, the 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 a 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_{i}=(X_{i}^{\top}\theta^{*})^{p}+\zeta_{i}$, where $\zeta_{i}$ is an additive noise and $p\geq 2$ is an integer. Due to nonlinear structure of the generative model, even if we select a convex loss function $\ell$, the ERM problem denoted to the considered GLM could be non-convex with respect to $\theta$. Interestingly, depending on the norm of $\theta^{*}$, the curvature of the ERM problem and its corresponding population risk minimization problem could change substantially. More precisely, in the case that $\|\theta^{*}\|$ 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 $\|\theta^{*}\|$ is close to zero, denoted by the low SNR regime, the underlying problem is neither strongly convex nor smooth, and in fact, it is ill-conditioned.

These observations 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 sublinear rate and thus requires polynomial number of iterations 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 a method still remains a 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 $|\theta^{*}|$. 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 inaccuracies could potentially lead to a 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 distance-generating function can solve the population loss corresponding to the low SNR setting at a linear rate. However, the linear convergence rate of mirror descent, similar to GD with Polyak step size, also depends on the condition number of the problem, and hence could lead to a slow convergence rate.

A natural approach to handle the ill-conditioning issue in the low SNR case as 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 $\frac{2p-2}{2p-1}$ in the population case, where $p$ is the degree of the polynomial link function. Moreover, in the high SNR 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 $\mathcal{O}(d^{3})$ which is indeed more than the per iteration computational cost of GD that is $\mathcal{O}(d)$. These points lead to this question:

Is there a computationally-efficient 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 a quasi-Newton method that approximates the objective function curvature using gradient information and its per iteration cost is $\mathcal{O}(d^{2})$. It is well-known that it enjoys a 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 we minimize the population loss, we show that in the low SNR case the iterates of BFGS converge to the ground truth $\theta^{*}$ at an exponentially fast rate that is independent of all problem parameters except the power of link function $p$. 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 linear convergence of BFGS without line-search for a setting 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, we show that in the finite sample case, the BFGS iterates converge to the final statistical radius $\mathcal{O}(1/n^{1/(2p+2)})$ within the true parameter after a logarithmic number of iterations $\mathcal{O}(\log(n))$. It is substantially lower than the required number of iterations for fixed-step size GD, which is $\mathcal{O}(n^{(p-1)/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 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. Section 4 explores BFGS’s convergence in low SNR settings, highlighting its linear convergence rate, a marked improvement over gradient descent’s sublinear rate. This section also compares the convergence rates of BFGS and Newton’s method. Section 5 applies these insights to establish the convergence results of BFGS for the empirical loss $\mathcal{L}_{n}$ in the low SNR regime. Lastly, our numerical experiments are presented in Section 6.

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:\mathbb{R}^{d}\to\mathbb{R}$. The BFGS update is given by

where $\eta_{k}$ is the step size and $H_{k}\in\mathbb{R}^{d\times d}$ is a positive definite matrix that aims to approximate the true Hessian inverse $\nabla^{2}{f(\theta_{k})}^{-1}$. There are several approaches for approximating $H_{k}$ leading to different quasi-Newton methods, (Conn et al., 1991; Broyden, 1965; Broyden et al., 1973; Gay, 1979; Davidon, 1959; Fletcher and Powell, 1963; Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970; Nocedal, 1980; Liu and Nocedal, 1989), but in this paper, we focus on the celebrated BFGS method, in which $H_{k}$ is updated as

where the variable variation $s_{k}$ and gradient displacement $u_{k}$ are defined as

for any $k\geq 1$ respectively. The logic behind the update in (4) is to ensure that the Hessian inverse approximation $H_{k}$ satisfies the secant condition $H_{k}u_{k-1}=s_{k-1}$, while it stays close to the previous approximation matrix $H_{k-1}$. The update in (4) only requires matrix-vector multiplications, and hence, the computational cost per iteration of BFGS is $\mathcal{O}(d^{2})$.

The main advantage of BFGS is its fast superlinear convergence rate under the assumption of strict convexity,

where ${\theta}_{opt}$ 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 quasi-Newton methods (Rodomanov and Nesterov, 2021a, b, c; Jin and Mokhtari, 2020; Jin et al., 2022; Ye et al., 2021; Lin et al., 2021a, b). For instance, Jin and Mokhtari (2020) established a local superlinear convergence rate of $(1/\sqrt{k})^{k}$ for BFGS. However, all these superlinear convergence analyses 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 setting, as we discuss in the upcoming sections.

In this section, we formally present the generalized linear model (GLM) setting 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\in\mathbb{R}^{d}$ and their corresponding labels are denoted by $Y\in\mathbb{R}$. Suppose that we have access to $n$ sample points $(Y_{1},X_{1}),(Y_{2},X_{2}),\ldots,(Y_{n},X_{n})$ 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.,

where $\theta^{*}$ is a true but unknown parameter, $p\in\mathbb{N}$ is a given power, and $\zeta_{1},\ldots,\zeta_{n}$ are independent Gaussian noises with zero mean and variance $\sigma^{2}$. 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 as $X\sim\mathcal{N}(0,\Sigma)$ where $\Sigma\in\mathbb{R}^{d\times d}$ is a symmetric positive definite matrix. Here we focus on the settings that $p\in\mathbb{N}^{+}$ and $p\geq 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\geq 2$ also serves as the basis of the analysis on other popular statistical models. For example, as shown by Yi and Caramanis (2015); Wang et al. (2015); Xu et al. (2016); Balakrishnan et al. (2017); Daskalakis et al. (2017); Dwivedi et al. (2020a); Kwon et al. (2019); 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 $p=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. (2021) under the symmetric two-component location Gaussian mixture, there are essentially three regimes for estimating the true parameter $\theta^{*}$: Low signal-to-noise ratio (SNR) regime: $\|\theta^{*}\|/\sigma\leq C_{1}(d/n)^{1/4}$ where $d$ is the dimension, $n$ is the sample size, and $C_{1}$ is a universal constant; Middle SNR regime: $C_{1}(d/n)^{1/4}\!\leq\!\|\theta^{*}\|/\sigma\leq C$ where $C$ is a universal constant; High SNR regime: $\|\theta^{*}\|/\sigma\geq C$. The main idea is that with different $\theta^{*}$, 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\geq 2$:

• •

  (i) Low SNR regime: $\|\theta^{*}\|/\sigma\leq C_{1}(d/n)^{1/(2p)}$ where $d$ is the dimension, $n$ is the sample size, and $C_{1}$ is a universal constant;

• •

  (ii) Middle SNR regime: $C_{1}(d/n)^{1/(2p)}\!\leq\!\|\theta^{*}\|/\sigma\leq C$ where $C$ is a universal constant;

• •

  (iii) High SNR regime: $\|\theta^{*}\|/\sigma\geq 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 $\theta^{*}$ when $\theta^{*}$ approaches 0. 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:

To obtain insight about the landscape of the empirical loss function $\mathcal{L}_{n}$, a useful approach is to consider that function by its population version, which we refer to as population least-square loss function:

where the outer expectation is taken with respect to the data.

High SNR regime. In the setting that the ground truth parameter has a relatively large norm, i.e., $\|\theta^{*}\|\geq C$ for some constant $C>0$ that only depends on $\sigma$, the population loss in (8) is locally strongly convex and smooth around $\theta^{*}$. More precisely, when $\|\theta-\theta^{*}\|$ is small, using Taylor’s theorem we have

Hence, in a neighborhood of the optimal solution, the objective in (8) can be approximated as

Indeed, if $\|\theta^{*}\|\geq C\sigma$ the above function behaves as a quadratic function that is smooth and strongly convex, assuming that $o(\|\theta-\theta^{*}\|^{2})$ is negligible. As a result, the iterates of gradient descent (GD) converge to the solution at a linear rate and it requires $\kappa\log(1/\epsilon)$ to reach an $\epsilon$-accurate solution, where $\kappa$ depends on the conditioning of the covariance matrix $\Sigma$ and the norm of $\theta^{*}$. In this case, BFGS converges superlinearly to $\theta^{*}$ and the rate would be independent of $\kappa$, while the cost per iteration would be $\mathcal{O}(d^{2})$.

Low SNR regime. As mentioned above, in the high SNR case, GD has a fast linear rate. However, in the low SNR case where $\|\theta^{*}\|$ is small and $\|\theta^{*}\|\leq C_{1}(d/n)^{1/(2p)}$, the strong convexity parameter approaches zero when the sample size $n$ goes to infinity and the problem becomes ill-conditioned. In this case, we deal with a function that is only convex and its gradient is not Lipschitz continuous. To better elaborate on this point, let us focus on the case that $\theta^{*}=0$ as a special case of the low SNR setting. Considering the underlying distribution of $X$, which is $X\sim\mathcal{N}(0,\Sigma)$, for such a low SNR case, the population loss can be written as

Since we focus on $p\geq 2$ it can be verified that $\mathcal{L}(\theta)$ is not strongly convex in a neighborhood of the solution $\theta^{*}=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 are 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 setting of the generalized linear model, the eigenvalues of the Hessian matrix of the population least-square loss function approach 0 and their vanishing rates depend on some increasing function of $\|\theta^{*}\|$. The optimal statistical rate of the optimization algorithms, such as gradient descent algorithm, for solving $\theta^{*}$ depends strictly on the tightness of these vanishing rates in terms of $\|\theta^{*}\|$, which are non-trivial to obtain. In fact, to the best of our knowledge, there is no result on the convergence of iterative 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 a future work.

In this section, we focus on the convergence properties of BFGS for the population loss in the low SNR case introduced in (9). 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 to a more encompassing objective function, $f$. Detailed in the following expression, this function encompasses the population loss function $\mathcal{L}$ with $\theta^{*}=0$ as a specific example. This approach enables us to present our results in the most general setting possible. Specifically, we examine the function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$, defined as follows:

where $A\in\mathbb{R}^{m\times d}$ is a matrix, $b\in\mathbb{R}^{m}$ is a given vector, and $q\in\mathbb{Z}$ satisfies $q\geq 4$. We should note that for $q\geq 4$, the considered objective is not strictly convex because the Hessian is singular when $A\theta=b$. Indeed, if we set $m=d$ and further let $A$ be $\Sigma^{1/2}$ and choose $b=A\theta^{*}=0$, then we recover the problem in (9) for $q=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$ is also an integer.

Notice that the problem in (10), which serves as a 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 $\|A\theta-b\|^{2}$ with solution $(A^{\top}A)^{-1}A^{\top}b$ when $A^{\top}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 easier 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 properties 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 on the population loss which is a special case of (10) is closely related to the one for the empirical loss in (7).

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 we exploit to establish our result.

The above assumption is also satisfied for our considered setting as we assume that the covariance matrix for our input features is positive definite. Combining Assumptions 1 and 2, we conclude that $\hat{\theta}$ is the unique optimal solution of problem (10). Next, we state the convergence rate of BFGS for solving problem (10) under the disclosed assumptions. The proof of the next theorem is available in Appendix A.1.

Theorem 1 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 iterates generated by BFGS converge globally without any line search scheme and the step size is fixed as $\eta_{k}=1$ for any $k\geq 0$. Moreover, the initial point $\theta_{0}$ is arbitrary and there is no restriction on the distance between $\theta_{0}$ and the optimal solution $\hat{\theta}$.

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 a Newton’s update to solve the problem in (10), then the error vectors $(\theta_{k}-\hat{\theta})_{k\geq 0}$ would all be parallel. This property arises due to the fact that for the problem in (10) the Newton direction is $\nabla^{2}f(\theta_{k-1})^{-1}\nabla f(\theta_{k-1})=(\theta_{k-1}-\hat{\theta})-\frac{q-2}{q-1}(\theta_{k-1}-\hat{\theta})$. Consequently, the next error vector $\theta_{k-1}-\hat{\theta}$ computed by running one step of Newton would satisfy $\theta_{k}-\hat{\theta}=\frac{q-2}{q-1}(\theta_{k-1}-\hat{\theta})$. Therefore, the error vector at time $k$, denoted by $\theta_{k}-\hat{\theta}$, is parallel to the previous error vector $\theta_{k-1}-\hat{\theta}$, with the only difference being that its norm is reduced by a factor of $\frac{q-2}{q-1}$. Using an induction argument, it simply follows that all error vectors $(\theta_{k}-\hat{\theta})_{k\geq 0}$ remain parallel to each other while their norm contracts at a rate of $\frac{q-2}{q-1}$. This is a key point that is used in the result for Newton’s method in Theorem 3. In the convergence analysis of BFGS (stated in the proof of Theorem 1), we use a 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 $(\theta_{k}-\hat{\theta})_{k\geq 0}$ as parallel to each other. Additionally, we show that the norm is shrinking at a factor of $r_{k}<1$, which is always larger than $\frac{q-2}{q-1}$, yet remains independent of the problem’s condition number or dimensions. In the following theorem, we show that for $q\geq 4$, the linear rate contraction factors $\{r_{k}\}_{k=0}^{\infty}$ also converge linearly to a fixed point contraction factor $r_{*}$ determined by the parameter $q$. The proof is available in Appendix A.2.

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 $\{r_{k}\}_{k=0}^{\infty}$ converge to the fixed point $r_{*}$ at a linear rate with the contraction factor of ${1}/{2}$. Further, the linear convergence factors $\{r_{k}\}_{k=0}^{\infty}$ and their limit $r_{*}$ are all only determined by $q$, and they are independent of the dimension $d$ and the condition number $\kappa_{A}$ 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 linear 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_{k}\}_{k=0}^{\infty}$ to the fixed point $r_{*}$ in Figure 1 for $q=4$ and $q=100$. In plots (a) and (b), we observe that $r_{k}$ becomes close to $r_{*}$ after only $5$ iterations. Hence, the linear convergence rate of BFGS is approximately $r_{*}$ after only a few iterations. We further observe in plots (c) and (d) that the factors $\{r_{k}\}_{k=0}^{\infty}$ converge to the fixed point $r_{*}$ at a linear rate upper bounded by $1/2$. Note that $r_{*}$ is the solution of (13). These plots verify our results in Theorem 2.

Thus far, we have demonstrated that the BFGS iterates converge linearly to the true parameter $\theta^{*}$ when minimizing the population loss function $\mathcal{L}$ 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 $\mathcal{L}_{n}$ defined in (7) for the low SNR setting. The iterates of BFGS in this case follow the update rule

where $H_{k}^{n}$ is updated using the gradient information of the loss $\mathcal{L}_{n}$ by the BFGS rule.

We next show that the BFGS iterates (16) $\{\theta_{k}^{n}\}_{k\geq 0}$ 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 $\{\theta_{k}^{n}\}_{k\geq 0}$ generated based on the empirical loss and the iterates $\{\theta_{k}\}_{k\geq 0}$ generated according to the population loss. Assuming that they both start from the same initialization $\theta_{0}$, with the concentration of the gradient $\|\nabla\mathcal{L}_{n}(\theta)-\nabla\mathcal{L}(\theta)\|$ and the Hessian $\|\nabla^{2}\mathcal{L}_{n}(\theta)-\nabla^{2}\mathcal{L}(\theta)\|_{\mathrm{op}}$ from Mou et al. (2019); Ren et al. (2022b) 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 is available in Appendix A.5.

Our analysis hinges on linking the gradient and Hessian of the empirical loss to the population loss, subsequently demonstrating a convergence rate for the empirical loss based on the linear convergence established in Section 4 for the population loss. More details can be found in Appendix A.5. Theorem 4 shows that BFGS achieves an estimation accuracy of $\mathcal{O}((d/n)^{\frac{1}{2p+2}})$in $O(\log n)$ iterations, which is substantialy faster than GD that requires $O(n^{\frac{p-1}{p}})$ iteations to achieve the same guarantee as shown in (Ren et al., 2022a). A few comments about 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 $\mathcal{O}(n^{-1/(2p+2)})$ within the true parameter $\theta^{*}$ after $\mathcal{O}(\log(n))$ number of iterations. The statistical radius is slightly worse than the optimal statistical radius $O(n^{-1/(2p)})$. However, we conjecture that this is due to the proof technique and BFGS can still reach the optimal $O(n^{-1/(2p)})$ in practice. In our experiments, in the next section, we observe that when $d=4$ and $p=2$, the statistical radius of BFGS is closer to the optimal radius of $O(n^{-1/4})$ instead of $O(n^{-1/6})$ suggested by our analysis. We leave an improvement of the statistical analysis as the future work. On the other hand, the overall iteration complexity of BFGS, which is $\mathcal{O}(\log(n))$, is indeed better than the polynomial number of iterations of GD, which is at the order of $\mathcal{O}(n^{(p-1)/p})$ (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 $\mathcal{O}(\kappa\log(n))$ iterations (Corollary 1 in (Ren et al., 2022a)), where $\kappa$ is the condition number of the covariance matrix $\Sigma$. Note that while the iteration complexity of BFGS is comparable to that of 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 $\mathcal{L}_{n}$, which can be unstable in practice, and also removes the dependency on the condition number that appears in the complexity bound of GD with Polyak step size. A similar conclusion also holds for mirror descent with a distance-generating function $h(x)=\frac{1}{q+2}|x|^{q+2}+\frac{1}{2}|x|^{2}$, as its linear convergence rate has a contraction factor of $(1-1/\kappa)$, leading to an overall iteration complexity of $\mathcal{O}(\kappa\log(n))$.

Finally, the iteration complexity of BFGS is comparable to the $\mathcal{O}(\log(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 4 involve the minimum number of iterations, namely, this result holds for some $1\leq t\leq T$. It 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 it may sound limited, this issue can be handled via an early stopping scheme using the cross-validation approaches. We illustrate such early stopping of the BFGS iterates for the low SNR regime in Figure 3.

Our experiments are divided into two sections: the first focuses on iterative methods’ behavior on the population loss of 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 step size and distance generating function $h(x)=\frac{1}{q+2}\|x\|^{q+2}+\frac{1}{2}\|x\|^{2}$ applied to (10) which corresponds to the population loss. We choose different values of parameter $m$, dimension $d$ and the exponential parameter $q$ in (10). We generate a random matrix $A\in\mathbb{R}^{m\times d}$ and a random vector $\hat{\theta}\in\mathbb{R}^{d}$, and compute the vector $b=A\hat{\theta}\in\mathbb{R}^{d}$. The initial point $\theta_{0}\in\mathbb{R}^{d}$ is also generated randomly. To properly select the steps for GD with a fixed step size $\eta$, 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 $\|\theta_{k}-\hat{\theta}\|$ versus the number of iteration $k$ for different algorithms. All the values of different parameters $m$, $d$, $q$ and $\eta$ are mentioned in the caption of the plots.

In Figure 2, 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 $q$ 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.e., the dimension $d$ has no impact over the performance of BFGS and Newton’s method. Although the convergence speed of Newton’s method is faster than BFGS, their gap is not significantly large as we expected from our theoretical results in Section 4. We also observe that mirror descent outperforms GD, but its performance is not as good as the BFGS method.

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 $p=2$, namely, we consider the multivariate setting of the phase retrieval problem. The data is generated by first sampling the inputs according to $\{X_{i}\}_{i=1}^{n}\sim\mathcal{N}(0,\mathrm{diag}(\sigma_{1}^{2},\cdots,\sigma_{4}^{2}))$ where $\sigma_{k}=(0.5)^{k-1}$, and then generating their labels based on $Y_{i}=(X_{i}^{\top}\theta^{*})^{2}+\zeta_{i}$ where $\{\zeta_{i}\}_{i=1}^{n}$ are i.i.d. samples from $\mathcal{N}(0,0.01)$. In the low SNR regime, we set $\theta^{*}=0$, and in the high SNR regime we select $\theta^{*}$ uniformly at random from the unit sphere. Further, for GD, we set the step size as $\eta=0.1$, while for Newton’s method and BFGS, we use the unit step size $\eta=1$.

In plots (a) and (b) of Figure 3, we consider the setting that the sample size is $n=10^{4}$, 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 $\mathcal{L}_{n}$. Furthermore, for both Newton’s method and the BFGS algorithm, due to their instability, we also perform cross-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 method converge geometrically fast to the final statistical radius while those of the GD converge slowly to 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 to the optimal Polyak step size, while the faster convergence of BFGS and Newton’s method over GD is due to 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 $10^{2}$ to $10^{4}$ to empirically verify the statistical radius of these methods. As indicated in the plots of that figure, under the high SNR regime, the BFGS has statistical radius is $\mathcal{O}(n^{-1/2})$, while under the low SNR regime, its statistical radius becomes $\mathcal{O}(n^{-1/4})$.

To show that BFGS can also be applied to high dimension scenarios, we conduct additional experiments on the generalized linear model with input $d=50,100,500$ and the power of link function $p=2$. The inputs are generated by $\{X_{i}\}_{i=1}^{n}\sim\mathcal{N}(0,\mathrm{diag}(\sigma_{1}^{2},\cdot,\sigma_{d}^{2}))$ where $\sigma_{k}=(0.96)^{k-1}$, and the remaining setting and hyper-parameters are set identical to the low dimension scenarios. The results are shown in Figure 4 and 5. As the results show, the performance of BFGS in high dimensional scenarios are nearly identical to the low dimensional scenarios.

In this paper, we analyzed the convergence rates of 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 it requires a 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 linear in the dimension and scales as $\mathcal{O}(d^{2})$. One future research direction is to analyze some other iterative methods 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 is $\mathcal{O}(d)$.