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

Qiujiang Jin Tongzheng Ren Nhat Ho Aryan Mokhtari Department of Electrical and
Computer Engineering, UT Austin<EMAIL_ADDRESS>of
Computer Science, UT Austin<EMAIL_ADDRESS>of Statistics
and Data Sciences, UT Austin<EMAIL_ADDRESS>of Electrical
and Computer Engineering, UT Austin<EMAIL_ADDRESS>
###### 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. While GD has optimal
statistical and computational complexities for estimating the true parameter
under the high signal-to-noise ratio (SNR) regime of the GLMs, it has sub-
optimal complexities when the SNR is low, namely, the iterates of GD require
polynomial number of iterations to reach the final statistical radius. The
slow convergence of GD for the low SNR case is mainly due to the local
convexity of the least-square loss functions of the GLMs. Even though Newton’s
method can be used to resolve the flat curvature of the loss functions, its
computational cost is prohibitive in high-dimensional settings due to its per
iteration cost of $\mathcal{O}(d^{3})$. To address the shortcomings of GD and
Newton’s method, we propose to use Broyden–Fletcher–Goldfarb–Shanno (BFGS)
quasi-Newton method to solve parameter estimation of the GLMs, which has a per
iteration cost of $\mathcal{O}(d^{2})$. On the optimization side, when the SNR
is low, we demonstrate that 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. On
the statistical side, we prove that the iterates of BFGS reach the final
statistical radius of the low SNR GLMs after a logarithmic number of
iterations, which is much lower than the polynomial number of iterations of
GD. We also present numerical experiments that match our theoretical findings.

##  1 Introduction

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 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

$\min_{\theta\in\mathbb{R}^{d}}\mathcal{L}_{n}(\theta):=\frac{1}{n}\sum_{i=1}^{n}\ell(\theta;(X_{i},Y_{i})),$
(1)

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 approximate of the ground-truth solution
$\theta^{*}$, where $\theta^{*}$ is also a minimizer of the population loss
defined as

$\min_{\theta\in\mathbb{R}^{d}}\mathcal{L}(\theta):=\mathbb{E}\left[\ell(\theta;(X,Y))\right].$
(2)

If one can solve the empirical risk efficiently, the output model could be
close to $\theta^{*}$, when $n$ is sufficiently large. There are several works
on studying 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)
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 and its corresponding population risk minimization
problem could change substantially. More precisely, if we are in the setting
that $\|\theta^{*}\|$ is sufficiently large, which we refer to this case as
the high signal-to-noise ratio (SNR), the underlying population loss of the
problem of interest is locally strongly convex and smooth. On the other hand,
when we are in the regime that $\|\theta^{*}\|$ is close to zero, denoted by
the low SNR regime, then 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 final statistical radius 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 to the
final statistical radius and thus requires polynomial number of iterations in
terms of the sample size. To resolve this issue, in (Ren et al., 2022) the
authors recommended the use of GD with Polyak step size to accelerate the
convergence of GD in the low SNR setting, and they showed that the number of
iterations becomes logarithmic function of the sample size. However, as this
method is still a first-order method, its overall complexity scales linearly
by the condition number of the problem which depends on the condition number
of the feature vectors covariance as well as the norm $\|\theta^{*}\|$.
Moreover, implementation of Polyak step size requires the knowledge of optimal
objective function value. Even we can approximate the optimal value
sufficiently well, the Polyak iterates still have instability during training
due to the influence of noise in the models.

An alternative approach is using the update of Newton’s method to handle the
ill-conditioning issue in the low SNR case as well as eliminating the need to
estimate the optimal function value. 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}$, when the
sample size is infinite. 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 the following question: Can we find a method that performs well
in both high and low SNR settings at a reasonable per iteration computational
cost?

Contributions. In this paper, we show that the answer to the above question is
positive and the BFGS method is capable of achieving these goals. BFGS is a
quasi-Newton method that approximates objective function curvature using
gradient information and its per iteration cost is $\mathcal{O}(d^{2})$. It is
well-known that it enjoys from a superlinear convergence rate that is
independent of condition number in strongly convex and smooth settings, and
hence, in high SNR setting it performs well. In this paper, we also show that
in the low-SNR settings it improves the iteration complexity of GD. Next, we
formally summarize our contributions.

  * •

Infinite sample, low SNR: For the infinite sample case where we are minimizing
the population loss, when we are in the low SNR case, the iterates of BFGS
converge to the ground truth $\theta^{*}$ at an exponentially fast rate that
is independent 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 not strictly or strongly convex.

  * •

Finite sample, low SNR: Under the low SNR regime of the generalized linear
model where we solve the ERM problem and we have finite sample, we demonstrate
that the BFGS iterates converge to the final statistical radius
$\mathcal{O}(1/n^{1/(2p)})$ 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 the similar statistical radius.
Furthermore, while the iteration complexity of the BFGS is comparable to the
logarithmic number of iterations of the Polyak step size GD, the BFGS removes
the dependency of the convergence of the iterates on the condition number of
the problem as well as the need to estimate the optimal function value.

  * •

Experiments: We conduct numerical experiments on both infinite and finite
sample loss objective function to compare the performance of GD (with constant
stepsize and Polyak stepsize), Newton’s method, and BFGS. These empirical
results are consistent with our theoretical findings and show the advantages
of BFGS in the low SNR regime.

Outline. In Section 2, we first briefly recap the BFGS quasi-Newton method
that we study in this paper. Then, in Section 3, we introduce the generalized
linear model and discuss low SNR and high SNR settings. In Section 4, we
present the linear convergence rate of BFGS method and Newton’s method applied
to the population loss of the generalized linear model defined in (6). In
Section 5, we extend the BGS quasi-Newton method to the empirical loss
function of the generalized linear model. In Section 6, we conduct some
numerical experiments to empirically verify our theoretical results. Finally,
in Section 7, we close the paper with some concluding remarks.

##  2 BFGS algorithm

In this section, we briefly review the basics of the BFGS quasi-Newton method,
which is the main algorithm we analyze in this paper. Consider the case that
we aim to minimize a differentiable convex function
$f:\mathbb{R}^{d}\to\mathbb{R}$ with optimal solution $\hat{\theta}$. The
iterative method is defined as

$\theta_{k+1}=\theta_{k}-\eta_{k}H_{k}\nabla{f(\theta_{k})},\qquad\forall
k\geq 0,$ (3)

where $H_{k}\in\mathbb{R}^{d\times d}$ is the matrix and $\eta_{k}$ is the
step size. If we set $H_{k}$ as the identity matrix $I$, we recover the
gradient descent whose computational cost per iteration is $\mathcal{O}(d)$.
If we set $H_{k}$ as $\nabla^{2}{f(\theta_{k})}^{-1}$, we obtain the update of
Newton’s method with $\mathcal{O}(d^{3})$ computational cost per iteration.
Quasi-Newton methods serve as a middle ground between first-order and the
second-order methods. The main idea of quasi-Newton methods is to come up with
a Hessian inverse approximation matrix $H_{k}$ that is close to the exact
Hessian inverse $\nabla^{2}f(\theta_{k})^{-1}$ using only first-order
information. 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

$H_{k}=\left(I-\frac{s_{k-1}u_{k-1}^{\top}}{s_{k-1}^{\top}u_{k-1}}\right)H_{k-1}\left(I-\frac{u_{k-1}s_{k-1}^{\top}}{s_{k-1}^{\top}u_{k-1}}\right)+\frac{s_{k-1}s_{k-1}^{\top}}{s_{k-1}^{\top}u_{k-1}},\qquad\forall
k\geq 1,$ (4)

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

$s_{k-1}:=\theta_{k}-\theta_{k-1},\qquad
u_{k-1}:=\nabla{f(\theta_{k})}-\nabla{f(\theta_{k-1})},\qquad\forall k\geq 1,$
(5)

respectively. The main logic behind the update in (4) is to ensure that the
Hessian inverse approximation matrix $H_{k}$ satisfies the secant condition
$H_{k}u_{k-1}=s_{k-1}$. We can observe that the update in (4) only requires
the matrix-vector multiplication. Therefore, the computational cost per
iteration of BFGS method is $\mathcal{O}(d^{2})$, where $d$ is the dimension.

The main advantage of BFGS is its fast superlinear convergence rate under the
assumption of strict convexity, i.e.,
$\lim_{k\to\infty}{\|\theta_{k}-\hat{\theta}\|}/{\|\theta_{k-1}-\hat{\theta}\|}=0$.
For the past few decades, all the superlinear convergence results of quasi-
Newton method are asymptotic. Most recently, non-asymptotic superlinear
convergence rates of quasi-Newton method have been established in the works of
(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, in (Jin and
Mokhtari, 2020) the authors 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. In this paper, as mentioned later, we will face settings in which the
Hessian at the optimal solution could be singular, and, hence the above
convergence guarantees do not hold and hence we need to establish new
convergence guarantees for BFGS.

##  3 Generalized linear model with polynomial link function

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$, i.e.,

$\displaystyle Y_{i}=(X_{i}^{\top}\theta^{*})^{p}+\zeta_{i},$ (6)

where $\theta^{*}$ is a true but unknown parameter, $p\in\mathbb{N}$ is a
given power, and $\zeta_{1},\ldots,\zeta_{n}$ are independent noises with zero
mean and variance $\sigma^{2}$. Note that, 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 that the feature
vectors are such that $X\in\mathbb{R}^{d}\sim\mathcal{N}(0,\Sigma)$ where
$\Sigma\in\mathbb{R}^{d\times d}$ is a symmetric positive definite matrix.

Note that in the above setting, when $p=1$, the model in (6) is the standard
linear regression model, while when $p=2$, it corresponds to the phase
retrieval model (Fienup, 1982). In this paper, we focus on the settings that
$p\geq 2$ and $p$ is an integer. It can be easily verified that the ERM
problem corresponding to the above model for $p\geq 2$ with quadratic loss
$\ell$ is a non-convex function with respect to $\theta$ and finding a global
minimizer of that could be a challenging task. On the other hand, it is
locally strongly convex when we are in the high signal-to-noise ratio (SNR)
case. In the low (SNR) setting, however, the problem becomes convex and the
strong convexity condition does not hold.

To showcase this issue, let us focus first on the population loss, which is
the limit of the ERM when the sample size goes to infinity. Note that the
population loss in the considered generalized linear model (6) with a
quadratic loss function $\ell$ is given by:

$\min_{\theta\in\mathbb{R}^{d}}\mathcal{L}(\theta)=\mathbb{E}_{X,Y}\left[\left(Y-(X^{\top}\theta)^{p}\right)^{2}\right],$
(7)

which based on the assumptions on the generalized linear model setting and the
distribution of the noise can be simplified as

$\min_{\theta\in\mathbb{R}^{d}}\mathcal{L}(\theta)=\mathbb{E}_{X}\left[\left((X^{\top}\theta^{*})^{p}-(X^{\top}\theta)^{p}\right)^{2}\right]+\sigma^{2}.$
(8)

Indeed, the ground truth parameter $\theta^{*}$ is an optimal solution of the
above. Next, we discuss the structure of the objective function in (8) for low
and high SNR settings.

High signal-to-noise regime. In the setting that the ground truth parameter
has a relative large norm, i.e., $\|\theta^{*}\|/\sigma\geq C$ for some
universal constant $C>0$, the population loss function in (8) is locally
strongly convex and smooth around $\theta^{*}$. More precisely, when
$\|\theta-\theta^{*}\|$ is small, we have

$(X^{\top}\theta^{*})^{p}-(X^{\top}\theta)^{p}=p(X^{\top}\theta^{*})^{p-1}X^{\top}(\theta-\theta^{*})+o(\|\theta-\theta^{*}\|).$

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

$\mathcal{L}(\theta)=p^{2}(\theta-\theta^{*})^{\top}\mathbb{E}_{X}\left[X(X^{\top}\theta^{*})^{2p-2}X^{\top}\right](\theta-\theta^{*})+\sigma^{2}+o(\|\theta-\theta^{*}\|^{2}).$

Indeed, if $\|\theta^{*}\|\geq C$ the above objective function behaves as a
quadratic function that is smooth and strongly convex, assuming that
$o(\|\theta-\theta^{*}\|^{2})$ is negligible.

Since the population loss in the high SNR case is almost a strongly convex
smooth quadratic function, the iterates of gradient descent (GD) converge to
the solution at a linear rate and hence 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 the optimal solution and the
rate would be independent of $\kappa$, however the cost per iteration would be
$\mathcal{O}(d^{2})$. Moreover, in this case Newton’s method converges at a
quadratic rate, but its per iterations cost is $\mathcal{O}(d^{3})$.

Low signal-to-noise 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 close to zero, the strong convexity parameter approaches zero 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$. Considering the
underlying distribution of $X$, for such a low SNR case, the population loss
can be written as

$\mathcal{L}(\theta)=\mathbb{E}_{X}\left[(X^{\top}\theta)^{2p}\right]+\sigma^{2}=(2p-1)!!\|\Sigma^{1/2}\theta\|^{2p}+\sigma^{2}.$
(9)

Since we focus on $p\geq 2$ it can be verified that this objective function is
neither strongly convex nor smooth, even 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.

##  4 Convergence analysis for the low signal-to-noise ratio case

In this section, we focus on the convergence properties of the BFGS method for
solving the population loss in the case of low signal-to-noise 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. Note
that the loss function in (9) can be considered as a special case of the
following convex optimization problem:

$\min_{\theta\in\mathbb{R}^{d}}f(\theta)=\|A\theta-b\|^{q},$ (10)

where $A\in\mathbb{R}^{m\times d}$ is a matrix, $b\in\mathbb{R}^{m}$ is a
given vector, and $q$ is the integer satisfying that $q\geq 4$. We should note
that for $q\geq 4$, the considered objective is not strictly convex because
the Hessian matrix is singular when $A\theta=b$. Indeed, if we set $m=d$ and
further let $A$ to be $\Sigma^{1/2}$ and choose $b=A\theta^{*}=0$, then we
recover the problem in (9) for $q=2p$. Instead of solving (9), we focus on the
convergence analysis of BFGS quasi-Newton method for solving the convex
function in (10), as it is a more general function and our results are of
general interest from an optimization point of view. To the best of our
knowledge, there is no global convergence theory (without line-search) for the
BFGS method when the function is not strictly convex, and our analysis
provides the first result for such general setting. First, we present two
assumptions of the objective function (10).

######  Assumption 1.

We assume that there exists $\hat{\theta}\in\mathbb{R}^{d}$, such that
$b=A\hat{\theta}$. This is equivalent to that vector $b$ is in the range space
of matrix $A$.

This assumption implies that the problem is realizable, $\hat{\theta}$ is an
optimal solution, and the optimal function value is zero. Note that the above
assumption is indeed satisfied in our considered setting low SNR case in (9)
as we assume $\theta^{*}=0$ which implies $b=0$.

######  Assumption 2.

The matrix $A^{\top}A\in\mathbb{R}^{d\times d}$ is invertible. This is
equivalent to that matrix $A^{\top}A$ is symmetric positive definite, i.e.
$A^{\top}A\succ 0$.

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

######  Theorem 1.

Consider the update of BFGS in (3), (4) and (5). Suppose Assumptions 1 and 2
are satisfied, and the initial point $\theta_{0}$ is an arbitrary vector in
$\mathbb{R}^{d}$ and the initial Hessian inverse approximation matrix is
selected as $H_{0}=\nabla^{2}{f(\theta_{0})}^{-1}$. If the step size of BFGS
is selected as $\eta_{k}=1$ for all $k\geq 0$, then the iterates of BFGS
converge to the optimal solution $\hat{\theta}$ at a linear rate of

$\|\theta_{k}-\hat{\theta}\|\leq
r_{k-1}\|\theta_{k-1}-\hat{\theta}\|,\quad\forall k\geq 1,$ (11)

where the contraction factors $r_{k}\in[0,1)$ satisfy the following conditions

$r_{0}=\frac{q-2}{q-1},\qquad
r_{k}=\frac{1-r_{k-1}^{q-2}}{1-r_{k-1}^{q-1}},\quad\forall k\geq 1.$ (12)

###### Proof.

The proof is available in Appendix A.1. ∎

The above theorem 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 stepsize is fixed as $\eta_{k}=1$ for any $k\geq 0$.
Moreover, the initial point $\theta_{0}$ could be any vector and there is no
restriction on the distance between $\theta_{0}$ and optimal solution
$\hat{\theta}$. Most analyses of quasi-Newton methods require the initial
point $\theta_{0}$ to be in a local neighborhood of $\hat{\theta}$ to
guarantee the linear or superlinear convergence rate, without line-search.

Note that the cost of computing the initial Hessian inverse approximation
matrix $H_{0}=\nabla^{2}{f(\theta_{0})}^{-1}$ is $\mathcal{O}(d^{3})$, but
this is only required for the first iteration and it is not required for
$k\geq 1$ as for those iterates we update the Hessian inverse approximation
matrix $H_{k}$ based on the BFGS update in (4) at a computational cost of
$\mathcal{O}(d^{2})$. In fact, for the first iterate even an approximation of
$H_{0}=\nabla^{2}{f(\theta_{0})}^{-1}$, which could possibly require a cost
less than $\mathcal{O}(d^{3})$, is sufficient.

Note that the result in Theorem 1 does not specify the exact complexity of
BFGS for solving problem(10), as the contraction factors $r_{k}$ are not
explicitly given. 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$.

######  Theorem 2.

Consider the linear convergence factors $\\{r_{k}\\}_{k=0}^{\infty}$ defined
in (12) from Theorem 1. If $q\geq 4$, then the sequence
$\\{r_{k}\\}_{k=0}^{\infty}$ converges to $r_{*}\in(0,1)$ that is determined
by the equation

$r_{*}^{q-1}+r_{*}^{q-2}=1,$ (13)

and the rate of convergence is linear with a contract factor that is at most
${1}/{2}$, i.e.,

$|r_{k}-r_{*}|\leq\left(\frac{1}{2}\right)^{k}|r_{0}-r_{*}|,\qquad\forall
k\geq 0.$ (14)

###### Proof.

The proof is available in Appendix A.2. ∎

We plot the convergence of factors $\\{r_{k}\\}_{k=0}^{\infty}$ to the fixed
point $r_{*}$ in Figure 1 for $q=4$ and $q=100$. We observe from (a) and (b)
of Figure 1 that $r_{k}$ is very close to $r_{*}$ after only $5$ iterations.
Therefore, the linear convergence rate of the BFGS method is approximately
$r_{*}$ after only a few iterations. We also observe from (c) and (d) of
Figure 1 that the factors $\\{r_{k}\\}_{k=0}^{\infty}$ converge to the fixed
point $r_{*}$ with a linear rate upper bounded by $1/2$. Note that $r_{*}$ is
simply the solution of (13). These plots verify our results in Theorem 2.

(a) $q=4$.

(b) $q=100$.

(c) $q=4$.

(d) $q=100$.

Figure 1: Convergence of factors $\\{r_{k}\\}_{k=0}^{\infty}$ to $r_{*}$.

Theorem 2 presented that eventually, the iterations generated by BFGS converge
to the optimal solution at the linear rate $r_{*}$ determined by (13). More
specifically, the factors $\\{r_{k}\\}_{k=0}^{\infty}$ converge to the fixed
point $r_{*}$ with the linear rate ${1}/{2}$. Therefore, the linear
convergence factors $\\{r_{k}\\}_{k=0}^{\infty}$ and the fixed point $r_{*}$
are totally determined by the parameter $q$. They are all independent of the
dimension $d$ and the condition number $\kappa_{A}$ of the matrix $A$. Hence,
the performance of the BFGS method is not influenced by high-dimensional or
ill-conditioned problems. This is indeed different from the common convergence
theories of most optimization algorithms, whose performance is deteriorated
heavily under the circumstance of high dimension or ill conditioning. We
should add that this result is independently important from the 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.

###  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 characterize the global
linear convergence of Newton’s method with unit step size applied to the
objective function in (10).

######  Theorem 3.

Consider applying Newton’s method to optimization problem (10) and suppose
Assumptions 1 and 2 hold. Moreover, suppose the step size is $\eta_{k}=1$ for
any $k\geq 0$. Then, the iterates of Newton’s method converge to the optimal
solution $\hat{\theta}$ at a linear rate of

$\|\theta_{k}-\hat{\theta}\|=\frac{q-2}{q-1}\|\theta_{k-1}-\hat{\theta}\|,\quad\forall
k\geq 1.$ (15)

Moreover, this linear convergence rate $\frac{q-2}{q-1}$ is smaller than the
fixed point $r_{*}$ defined in (13) of the BFGS quasi-Newton method,
i.e.,$\frac{q-2}{q-1}<r_{*}$ for all $q\geq 4$.

###### Proof.

The proof is available in Appendix A.3. ∎

The convergence results of Newton’s method are also global without any
backtracking line search method, and the linear rate $\frac{q-2}{q-1}$ is
independent of dimension $d$ and condition number $\kappa_{A}$. Furthermore,
the condition $\frac{q-2}{q-1}<r_{*}$ implies that iterations of Newton’s
method converge faster than the BFGS method, but the gap is not substantial as
we illustrate in our numerical results. On the other hand, as we mentioned
before, the computational cost per iteration of Newton’s method is
$\mathcal{O}(d^{3})$ which is much worse than the $\mathcal{O}(d^{2})$ of
BFGS.

Moving back to our main problem, one important implication of the above
convergence results is that in the low signal-to-noise ratio setting the
iterates of BFGS converge linearly to the optimal solution of the population
least-square loss function, while the contraction coefficient of BFGS is
comparable to that of Newton’s method which is $(2p-2)/(2p-1)$. For example,
for $p=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 $r_{*}$ are $0.755,0.857,0.922,0.963$, respectively. We use
this intuition to establish a similar result for the finite sample case in the
following section.

##  5 Statistical rate of BFGS for solving the sample least-square loss

Thus far, we have demonstrated that the BFGS iterates converge linearly to the
true parameter $\theta^{*}$ for solving the population least-square loss
function $\mathcal{L}$ of the generalized linear models in equation (9). It
provides important insight into the behaviors of the BFGS iterates when the
sample size $n$ is infinite. In this section, we would like to study the
statistical behaviors of the BFGS iterates for solving the least-square loss
function $\mathcal{L}_{n}$, which is given by:

$\displaystyle\min_{\theta\in\mathbb{R}^{d}}\mathcal{L}_{n}(\theta):=\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-(X_{i}^{\top}\theta)^{p})^{2}.$
(16)

To simplify the proof argument and gain the insight into the statistical
behaviors of the BFGS iterates, we specifically consider the univariate
setting, namely, $d=1$. We would like to remark that these statistical
behaviors of BFGS iterates still hold in multivariate settings, namely, $d\geq
2$ (See our experiments in Figure 3). When $d=1$, the BFGS iterates for
solving the sample loss function take the following form:

$\displaystyle\theta_{k+1}^{n}$
$\displaystyle=\theta_{k}^{n}-\eta_{k}\frac{\theta_{k}^{n}-\theta_{k-1}^{n}}{\nabla\mathcal{L}_{n}(\theta_{k}^{n})-\nabla\mathcal{L}_{n}(\theta_{k-1}^{n})}\nabla\mathcal{L}_{n}(\theta_{k}^{n}).$
(17)

Throughout this section, we consider the step size $\eta_{k}=1$ in the updates
(17) of BFGS. We show that the BFGS iterates (17) $\\{\theta_{k}^{n}\\}_{k\geq
0}$ converge to the final statistical radius after a logarithmic number of
iterations under the low SNR regimes of the generalized linear models.

######  Theorem 4.

Given the low SNR regime of the generalized linear model (6) namely,
$\theta^{*}=0$, when $d=1$ and the initialization is
$\theta_{n}^{0}\in\mathbb{B}(\theta^{*},\rho)$ for some $\rho>0$, there exist
universal constants $C_{1},C_{2},C_{3}$ such that $n\geq
C_{1}^{\prime}\log^{2p}(n/\delta)$, and $k\geq C_{2}\log(n/\log(1/\delta))$,
we have

$\min_{1\leq t\leq k}|\theta_{t}^{n}-\theta^{*}|\leq
C_{3}\left(\frac{\log^{2p}(n/\delta)}{n}\right)^{\frac{1}{2p}}.$ (18)

###### Proof.

The proof is available in Appendix A.4. ∎

A few comments with the results of Theorem 4 are in order.

Comparing to GD, GD with Polyak step size , and Newton’s method: The result of
Theorem 4 indicates that under the low SNR regime, namely, when
$\theta^{*}=0$, the BFGS iterates reach the final statistical radius
$\mathcal{O}(n^{-1/(2p)})$ within the true parameter $\theta^{*}$ after
$\mathcal{O}(\log(n))$ number of iterations. This complexity indeed is 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)). It is also
comparable to the logarithmic number of iterations of GD with Polyak step size
which requires $\mathcal{O}(\kappa\log(n))$ iterations (Corollary 1 in (Ren et
al., 2022)), where $\kappa$ is the condition number of the covriance matrix
$\Sigma$, and the $\mathcal{O}(\log(n))$ of Newton’s method (Corollary 3 in
(Ho et al., 2020)) to reach the similar statistical radius.

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. Finally, though the BFGS algorithm and the Newton’s
method have similar computational complexity in terms of $n$, BFGS has lower
per iteration cost in comparison to Newton’s method.

On the minimum number of iterations: The results of BFGS in Theorem 4 involve
the minimum number of iterations, namely, these results only hold for some
$1\leq t\leq k$. It suggests that the BFGS iterates may diverge after they
reach the final statistical radius under each regime of the generalized linear
models. As highlighted in (Ho et al., 2020), such instability behavior of BFGS
is inherent to fast and unstable methods. While it may sound limited, the
minimum number of iterations can be overcome 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.

Generalization of the results to multivariate settings: While the results of
Theorem 4 are only established for the univariate setting, we remark that it
is mainly for the simplicity of the proof argument and of the BFGS iterates.
In the experiments of Figures 3, we run BFGS for the case that dimension is
$d=4$ and observe that both the statistical radius and the iteration
complexities of BFGS are still consistent with those in Theorem 4. We leave a
theoretical verification of these results in multivariate settings for the
future work.

##  6 Numerical experiments

(a) $d=10,q=4$.

(b) $d=10,q=10$.

(c) $d=10^{3},q=4$,

(d) $d=10^{3},q=10$.

Figure 2: Convergence rates of Newton’s method, BFGS, GD with constant step
size and GD with Polyak step size for different $d$ and $q$. In plot (a),
$m=100$ and $\eta=10^{-4}$. In plot (b), $m=100$ and $\eta=10^{-8}$. In plot
(c), $m=2000$ and $\eta=10^{-12}$. In plot (d), $m=2000$ and $\eta=10^{-15}$.

###  6.1 Numerical experiment for the population loss function

In this section, we compare the performance of Newton’s method, BFGS, GD with
constant step size, and GD with Polyak step size 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. The GD constant step
size $\eta$ is tuned by hand to achieve the best performance of GD on each
problem. We present the logarithmic scale of $\|\theta_{k}-\hat{\theta}\|$
versus the number of iteration $k$ for different algorithms. All the values of
different parameters $m$, $d$, $q$ and $\eta$ as well as the numerical results
of our experiments are shown in Figure 2.

We observe that GD with constant step converges very slowly since it can only
reach a 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 the BFGS method generate
iterations with linear convergence rates, and their linear convergence rates
are only affected by the parameter $q$. 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 the BFGS method, their gap is not significantly
large. We should add that these empirical results are consistent with the
theoretical results we obtained in Section 4.

###  6.2 Numerical experiment for the empirical loss function

(a) High SNR regime.

(b) Low SNR regime.

(c) High SNR regime.

(d) Low SNR regime.

Figure 3: Illustration of different methods with high SNR regime in (a) and
low SNR regime in (b). Illustration of the statistical radius of BFGS with
high SNR regime in (c) and low SNR regime in (d).

We now move to illustrate the statistical and computational complexities of
BFGS for parameter estimation of the generalized linear model. In our
experiments, we specifically consider the dimension to be $d=4$ and the power
of the link function to be $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}$, 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,1)$. 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. Furthermore, for the
GD, we choose the step size to be $\eta=0.1$, while for Newton’s method and
BFGS, we select the unit stepsize $\eta=1$.

In (a) and (b) of Figure 3, we consider the sample size $n=10^{4}$ and 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 (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 optimality of step size of, while the faster convergence of BFGS and
Newton’s method over GD is due to their independence on the problem condition
number.

Finally, in (c) and (d) of Figure 3, we run the BFGS when the sample size is
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 $\mathcal{O}(n^{-1/2})$, while under
the low SNR regime, its statistical radius becomes $\mathcal{O}(n^{-1/4})$.
These empirical results are consistent with the theoretical results of the
BFGS in Theorem 4.

##  7 Conclusions

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 very
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 conducted numerical experiments for both infinite and finite sample loss
functions, and we showed that our empirical results are consistent with our
theoretical findings. Perhaps one limitation of the BFGS method is that its
computational cost is sill not linear in the dimension and scales as
$\mathcal{O}(d^{2})$. Therefore, one future research direction is to analyze
some other iterative methods such as the limited memory-BFGS (L-BFGS)
algorithm which may be able to achieve a fast linear convergence rate in the
low SRN setting, while its computational cost per iteration is
$\mathcal{O}(d)$.

## Acknowledgements

This research of Q. Jin and A. Mokhtari is supported in part by NSF Grants
2007668, 2019844, and 2112471, ARO Grant W911NF2110226, the Machine Learning
Lab (MLL) at UT Austin, and the Wireless Networking and Communications Group
(WNCG) Industrial Affiliates Program. NH ratefully acknowledges support from
NSF IFML 2019844 award and research gifts by UT Austin ML grant.

## References

  * Agarwal et al. (2012) 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\.
  * Balakrishnan et al. (2017) 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.
  * Broyden (1965) C. G. Broyden.  A class of methods for solving nonlinear simultaneous equations.  _Mathematics of computation_ , 19(92):577–593, 1965.
  * Broyden (1970) C. G. Broyden.  The convergence of single-rank quasi-Newton methods.  _Mathematics of Computation_ , 24(110):365–382, 1970.
  * Broyden et al. (1973) C. G. Broyden, J. E. D. Jr., Broyden, and J. J. More.  On the local and superlinear convergence of quasi-Newton methods.  _IMA J. Appl. Math_ , 12(3):223–245, June 1973\.
  * Candes et al. (2011) E. J. Candes, Y. Eldar, T. Strohmer, and V. Voroninski.  Phase retrieval via matrix completion, 2011.
  * Carroll et al. (1997) 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–489, 1997.
  * Conn et al. (1991) A. R. Conn, N. I. M. Gould, and P. L. Toint.  Convergence of quasi-Newton matrices generated by the symmetric rank one update.  _Mathematical programming_ , 50(1-3):177–195, 1991.
  * Davidon (1959) W. Davidon.  Variable metric method for minimization.  Technical report, Argonne National Lab., Lemont, Ill., 1959.
  * Dwivedi et al. (2020a) R. Dwivedi, N. Ho, K. Khamaru, M. J. Wainwright, M. I. Jordan, and B. Yu.  Sharp analysis of expectation-maximization for weakly identifiable models.  _AISTATS_ , 2020a.
  * Dwivedi et al. (2020b) R. Dwivedi, N. Ho, K. Khamaru, M. J. Wainwright, M. I. Jordan, and B. Yu.  Singularity, misspecification, and the convergence rate of EM.  _Annals of Statistics_ , 44:2726–2755, 2020b.
  * Fienup (1982) J. R. Fienup.  Phase retrieval algorithms: a comparison.  _Appl. Opt._ , 21(15):2758–2769, Aug 1982.  doi: 10.1364/AO.21.002758.  URL http://www.osapublishing.org/ao/abstract.cfm?URI=ao-21-15-2758.
  * Fletcher (1970) R. Fletcher.  A new approach to variable metric algorithms.  _The computer journal_ , 13(3):317–322, 1970\.
  * Fletcher and Powell (1963) R. Fletcher and M. J. Powell.  A rapidly convergent descent method for minimization.  _The computer journal_ , 6(2):163–168, 1963.
  * Gay (1979) D. M. Gay.  Some convergence properties of Broyden’s method.  _SIAM Journal on Numerical Analysis_ , 16(4):623–630, 1979.
  * Goebel and Kirk (1990) K. Goebel and W. A. Kirk.  _Topics in Metric Fixed Point Theory_.  Cambridge University Press, 1990.
  * Goldfarb (1970) D. Goldfarb.  A family of variable-metric methods derived by variational means.  _Mathematics of computation_ , 24(109):23–26, 1970.
  * Hardt et al. (2016) M. Hardt, B. Recht, and Y. Singer.  Train faster, generalize better: Stability of stochastic gradient descent.  In M. F. Balcan and K. Q. Weinberger, editors, _Proceedings of The 33rd International Conference on Machine Learning_ , volume 48 of _Proceedings of Machine Learning Research_ , pages 1225–1234, New York, New York, USA, 20–22 Jun 2016. PMLR.  URL http://proceedings.mlr.press/v48/hardt16.html.
  * Ho et al. (2020) N. Ho, K. Khamaru, R. Dwivedi, M. J. Wainwright, M. I. Jordan, and B. Yu.  Instability, computational efficiency and statistical accuracy.  _Arxiv Preprint Arxiv: 2005.11411_ , 2020.
  * Jin and Mokhtari (2020) Q. Jin and A. Mokhtari.  Non-asymptotic superlinear convergence of standard quasi-newton methods.  _arXiv preprint arXiv:2003.13607_ , 2020.
  * Jin et al. (2022) Q. Jin, A. Koppel, K. Rajawat, and A. Mokhtari.  Sharpened quasi-newton methods: Faster superlinear rate and larger local convergence neighborhood.  _The 39th International Conference on Machine Learning (ICML 2022)_ , 2022.
  * Lin et al. (2021a) D. Lin, H. Ye, and Z. Zhang.  Explicit superlinear convergence of broyden’s method in nonlinear equations.  _arXiv preprint arXiv:2109.01974_ , 2021a.
  * Lin et al. (2021b) D. Lin, H. Ye, and Z. Zhang.  Greedy and random quasi-newton methods with faster explicit superlinear convergence.  _Advances in Neural Information Processing Systems 34_ , 2021b.
  * Liu and Nocedal (1989) D. C. Liu and J. Nocedal.  On the limited memory BFGS method for large scale optimization.  _Mathematical programming_ , 45(1-3):503–528, 1989.
  * Loh and Wainwright (2015) P.-L. Loh and M. J. Wainwright.  Regularized M-estimators with nonconvexity: Statistical and algorithmic theory for local optima.  _Journal of Machine Learning Research_ , 16:559–616, 2015\.
  * Mou et al. (2019) W. Mou, N. Ho, M. J. Wainwright, P. Bartlett, and M. I. Jordan.  A diffusion process perspective on posterior contraction rates for parameters.  _arXiv preprint arXiv:1909.00966_ , 2019.
  * Netrapalli et al. (2015) P. Netrapalli, P. Jain, and S. Sanghavi.  Phase retrieval using alternating minimization.  _IEEE Transactions on Signal Processing_ , 63(18):4814–4826, 2015.  doi: 10.1109/TSP.2015.2448516.
  * Nocedal (1980) J. Nocedal.  Updating quasi-Newton matrices with limited storage.  _Mathematics of computation_ , 35(151):773–782, 1980.
  * Ren et al. (2022) T. Ren, F. Cui, A. Atsidakou, S. Sanghavi, and N. Ho.  Towards statistical and computational complexities of Polyak step size gradient descent.  _Artificial Intelligence and Statistics Conference_ , 2022.
  * Rodomanov and Nesterov (2021a) A. Rodomanov and Y. Nesterov.  Greedy quasi-newton methods with explicit superlinear convergence.  _SIAM Journal on Optimization_ , 31(1):785–811, 2021a.
  * Rodomanov and Nesterov (2021b) A. Rodomanov and Y. Nesterov.  Rates of superlinear convergence for classical quasi-newton methods.  _Mathematical Programming_ , pages 1–32, 2021b.
  * Rodomanov and Nesterov (2021c) A. Rodomanov and Y. Nesterov.  New results on superlinear convergence of classical quasi-newton methods.  _Journal of Optimization Theory and Applications_ , 188(3):744–769, 2021c.
  * Shanno (1970) D. F. Shanno.  Conditioning of quasi-Newton methods for function minimization.  _Mathematics of computation_ , 24(111):647–656, 1970.
  * Shechtman et al. (2015) Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev.  Phase retrieval with application to optical imaging: A contemporary overview.  _IEEE Signal Processing Magazine_ , 32(3):87–109, 2015.  doi: 10.1109/MSP.2014.2352673.
  * Ye et al. (2021) H. Ye, D. Lin, Z. Zhang, and X. Chang.  Explicit superlinear convergence rates of the sr1 algorithm.  _arXiv preprintarXiv:2105.07162_ , 2021.
  * Yuan and Zhang (2013) X.-T. Yuan and T. Zhang.  Truncated power method for sparse eigenvalue problems.  _Journal of Machine Learning Research_ , 14(Apr):899–925, 2013.

##  Appendix A Appendix

######  Lemma 1.

Consider the objective function in (10) satisfying Assumption 1 and 2. Then,
the inverse matrix of its Hessian $\nabla^{2}{f(\theta)}$ can be expressed as

$\nabla^{2}{f(\theta)}^{-1}=\frac{(A^{\top}A)^{-1}}{q\|A\theta-b\|^{q-2}}-\frac{(q-2)(\theta-\hat{\theta})(\theta-\hat{\theta})^{\top}}{q(q-1)\|A\theta-b\|^{q}}.$
(19)

###### Proof.

Notice that the Hessian of objective function (10) can be expressed as

$\nabla^{2}{f(\theta)}=q\|A\theta-b\|^{q-2}A^{\top}A+q(q-2)\|A\theta-b\|^{q-4}A^{\top}(A\theta-b)(A\theta-b)^{\top}A.$
(20)

We use the Sherman–Morrison formula. Suppose that $X\in\mathbb{R}^{d\times d}$
is an invertible matrix and $a,b\in\mathbb{R}^{d}$ are two vectors satisfying
that $1+b^{\top}X^{-1}a\neq 0$. Then, the matrix $X+ab^{\top}$ is invertible
and

$(X+ab^{\top})^{-1}=X^{-1}-\frac{X^{-1}ab^{\top}X^{-1}}{1+b^{\top}X^{-1}a}.$
(21)

Applying the Sherman–Morrison formula with $X=q\|A\theta-b\|^{q-2}A^{\top}A$,
$a=q(q-2)\|A\theta-b\|^{q-4}A^{\top}(A\theta-b)$ and $b=A^{\top}(A\theta-b)$.
Notice that $A^{\top}A$ is invertible, hence $X$ is invertible and

$\begin{split}&1+b^{\top}X^{-1}a\\\ =&\quad
1+(A\theta-b)^{\top}A\frac{(A^{\top}A)^{-1}}{q\|A\theta-b\|^{q-2}}q(q-2)\|A\theta-b\|^{q-4}A^{\top}(A\theta-b)\\\
=&\quad
1+(q-2)(A\theta-b)^{\top}A\frac{(A^{\top}A)^{-1}A^{\top}A(\theta-\hat{\theta})}{\|A\theta-b\|^{2}}\\\
=&\quad 1+(q-2)\frac{(A\theta-b)^{\top}(A\theta-b)}{{\|A\theta-b\|^{2}}}\\\
=&\quad q-1\neq 0.\qquad\text{(Since $q\geq 4$.)}\end{split}$ (22)

Therefore, we obtain that

$\begin{split}&\nabla^{2}{f(\theta)}^{-1}\\\
=&\quad\frac{(A^{\top}A)^{-1}}{q\|A\theta-b\|^{q-2}}-\frac{\frac{(A^{\top}A)^{-1}}{q\|A\theta-b\|^{q-2}}q(q-2)\|A\theta-b\|^{q-4}A^{\top}(A\theta-b)(A^{\top}(A\theta-b))^{\top}\frac{(A^{\top}A)^{-1}}{q\|A\theta-b\|^{q-2}}}{q-1}\\\
=&\quad\frac{(A^{\top}A)^{-1}}{q\|A\theta-b\|^{q-2}}-\frac{(q-2)}{q(q-1)\|A\theta-b\|^{q}}(A^{\top}A)^{-1}AA^{\top}(\theta-\hat{\theta})(\theta-\hat{\theta})^{\top}AA^{\top}(A^{\top}A)^{-1}\\\
=&\quad\frac{(A^{\top}A)^{-1}}{q\|A\theta-b\|^{q-2}}-\frac{(q-2)(\theta-\hat{\theta})(\theta-\hat{\theta})^{\top}}{q(q-1)\|A\theta-b\|^{q}}.\end{split}$
(23)

As a consequence, we obtain the conclusion of the lemma. ∎

######  Lemma 2.

Banach’s Fixed-Point Theorem.  Consider the differentiable function
$f:D\subset\mathbb{R}\to D\subset\mathbb{R}$. Suppose that there exists
$C\in(0,1)$ such that $|f^{\prime}(x)|\leq C$ for any $x\in D$. Now let
$x_{0}\in D$ be arbitrary and define the sequence $\\{x_{k}\\}_{k=0}^{\infty}$
as

$x_{k+1}=f(x_{k}),\qquad\forall k\geq 0.$ (24)

Then, the sequence $\\{x_{k}\\}_{k=0}^{\infty}$ converges to the unique fixed
point $x_{*}$ defined as

$x_{*}=f(x_{*}),$ (25)

with linear convergence rate of

$|x_{k}-x_{*}|\leq C^{k}|x_{0}-x_{*}|,\qquad\forall k\geq 0.$ (26)

###### Proof.

Check (Goebel and Kirk, 1990). ∎

###  A.1 Proof of Theorem 1

We use induction to prove the convergence results in (11) and (12). Note that
$b=A\hat{\theta}$ by Assumption 1 and the gradient and Hessian of the
objective function in (10) are explicitly given by

$\nabla{f(\theta)}=q\|A\theta-b\|^{q-2}A^{\top}(A\theta-b)=q\|A\theta-b\|^{q-2}A^{\top}A(\theta-\hat{\theta}),$
(27)
$\nabla^{2}{f(\theta)}=q\|A\theta-b\|^{q-2}A^{\top}A+q(q-2)\|A\theta-b\|^{q-4}A^{\top}(A\theta-b)(A\theta-b)^{\top}A.$
(28)

Applying Lemma 1, we can obtain that

$\nabla^{2}{f(\theta)}^{-1}=\frac{(A^{\top}A)^{-1}}{q\|A\theta-b\|^{q-2}}-\frac{(q-2)(\theta-\hat{\theta})(\theta-\hat{\theta})^{\top}}{q(q-1)\|A\theta-b\|^{q}}.$
(29)

First, we consider the initial iteration

$\theta_{1}=\theta_{0}-H_{0}\nabla{f(\theta_{0})}=\theta_{0}-\nabla{f(\theta_{0})}^{-1}\nabla{f(\theta_{0})},$
(30)
$\theta_{1}-\hat{\theta}=\theta_{0}-\hat{\theta}-\nabla{f(\theta_{0})}^{-1}\nabla{f(\theta_{0})}.$
(31)

Notice that $b=A\hat{\theta}$ by Assumption 1 and

$\begin{split}&\nabla{f(\theta_{0})}^{-1}\nabla{f(\theta_{0})}\\\
=&\quad[\frac{(A^{\top}A)^{-1}}{q\|A\theta_{0}-b\|^{q-2}}-\frac{(q-2)(\theta_{0}-\hat{\theta})(\theta_{0}-\hat{\theta})^{\top}}{q(q-1)\|A\theta_{0}-b\|^{q}}]q\|A\theta-b\|^{q-2}A^{\top}A(\theta_{0}-\hat{\theta})\\\
=&\quad\theta_{0}-\hat{\theta}-\frac{q-2}{q-1}\frac{(\theta_{0}-\hat{\theta})^{\top}A^{\top}A(\theta_{0}-\hat{\theta})}{\|A\theta_{0}-b\|^{2}}(\theta_{0}-\hat{\theta})\\\
=&\quad\theta_{0}-\hat{\theta}-\frac{q-2}{q-1}\frac{(A\theta_{0}-b)^{\top}(A\theta_{0}-b)}{\|A\theta_{0}-b\|^{2}}(\theta_{0}-\hat{\theta})\\\
=&\quad\theta_{0}-\hat{\theta}-\frac{q-2}{q-1}(\theta_{0}-\hat{\theta}).\end{split}$
(32)

Therefore, we obtain that

$\theta_{1}-\hat{\theta}=\theta_{0}-\hat{\theta}-\nabla{f(\theta_{0})}^{-1}\nabla{f(\theta_{0})}=\frac{q-2}{q-1}(\theta_{0}-\hat{\theta}).$
(33)

Condition (11) holds for $k=1$ with $r_{0}=\frac{q-2}{q-1}$. Now we assume
that condition (11) holds for $k=t$ where $t\geq 1$, i.e.,

$\theta_{t}-\hat{\theta}=r_{t-1}(\theta_{t-1}-\hat{\theta}).$ (34)

Considering the condition $b=A\hat{\theta}$ in Assumption 1 and the condition
in (34), we further have

$A\theta_{t}-b=A(\theta_{t}-\hat{\theta})=r_{t-1}A(\theta_{t-1}-\hat{\theta})=r_{t-1}(A\theta_{t-1}-b),$
(35)

which implies that

$\nabla{f(\theta_{t})}=qr_{t-1}^{q-1}\|A(\theta_{t-1}-\hat{\theta})\|^{q-2}A^{\top}A(\theta_{t-1}-\hat{\theta}).$
(36)

We further show that the variable displacement and gradient difference can be
written as

$s_{t-1}=\theta_{t}-\theta_{t-1}=\theta_{t}-\hat{\theta}-\theta_{t-1}+\hat{\theta}=(r_{t-1}-1)(\theta_{t-1}-\hat{\theta}),$
(37)

and

$u_{t-1}=\nabla{f(\theta_{t})}-\nabla{f(\theta_{t-1})}=q(r_{t-1}^{q-1}-1)\|A(\theta_{t-1}-\hat{\theta})\|^{q-2}A^{\top}A(\theta_{t-1}-\hat{\theta}).$
(38)

Considering these expressions, we can show that the rank-1 matrix in the
update of BFGS $u_{t-1}s_{t-1}^{\top}$ is given by

$u_{t-1}s_{t-1}^{\top}=q(r_{t-1}^{q-1}-1)(r_{t-1}-1)\|A(\theta_{t-1}-\hat{\theta})\|^{q-2}A^{\top}A(\theta_{t-1}-\hat{\theta})(\theta_{t-1}-\hat{\theta})^{\top},$
(39)

and the inner product $s_{t-1}^{\top}u_{t-1}$ can be written as

$\begin{split}s_{t-1}^{\top}u_{t-1}&=q(r_{t-1}^{q-1}-1)(r_{t-1}-1)\|A(\theta_{t-1}-\hat{\theta})\|^{q-2}(\theta_{t-1}-\hat{\theta})^{\top}A^{\top}A(\theta_{t-1}-\hat{\theta})\\\
&=q(r_{t-1}^{q-1}-1)(r_{t-1}-1)\|A(\theta_{t-1}-\hat{\theta})\|^{q}.\end{split}$
(40)

These two expressions allow us to simplify the matrix
$I-\frac{u_{t-1}s_{t-1}^{\top}}{s_{t-1}^{\top}u_{t-1}}$ in the update of BFGS
as

$I-\frac{u_{t-1}s_{t-1}^{\top}}{s_{t-1}^{\top}u_{t-1}}=I-\frac{A^{\top}A(\theta_{t-1}-\hat{\theta})(\theta_{t-1}-\hat{\theta})^{\top}}{\|A(\theta_{t-1}-\hat{\theta})\|^{2}}.$
(41)

An important property of the above matrix is that its null space is the set of
the vectors that are parallel to $u_{t-1}$. Considering the expression for
$u_{t-1}$, any vector that is parallel to
$A^{\top}A(\theta_{t-1}-\hat{\theta})$ is in the null space of the above
matrix. We can easily observe that the gradient defined in (36) satisfies this
condition and therefore

$\begin{split}&\left(I-\frac{u_{t-1}s_{t-1}^{\top}}{s_{t-1}^{\top}u_{t-1}}\right)\nabla{f(\theta_{t})}\\\
=&\quad
qr_{t-1}^{q-1}\|A(\theta_{t-1}-\hat{\theta})\|^{q-2}\left(I-\frac{A^{\top}A(\theta_{t-1}-\hat{\theta})(\theta_{t-1}-\hat{\theta})^{\top}}{\|A(\theta_{t-1}-\hat{\theta})\|^{2}}\right)A^{\top}A(\theta_{t-1}-\hat{\theta})\\\
=&\quad
qr_{t-1}^{q-1}\|A(\theta_{t-1}-\hat{\theta})\|^{q-2}\left(A^{\top}A(\theta_{t-1}-\hat{\theta})-\frac{A^{\top}A(\theta_{t-1}-\hat{\theta})\|A(\theta_{t-1}-\hat{\theta})\|^{2}}{\|A(\theta_{t-1}-\hat{\theta})\|^{2}}\right)\\\
=&\quad 0.\end{split}$ (42)

This important observation shows that if the condition in (34) holds, then the
BFGS descent direction $H_{t}\nabla{f(\theta_{t})}$ can be simplified as

$\begin{split}&H_{t}\nabla{f(\theta_{t})}\\\
=&\quad\left(I-\frac{s_{t-1}u_{t-1}^{\top}}{s_{t-1}^{\top}u_{t-1}}\right)H_{t-1}\left(I-\frac{u_{t-1}s_{t-1}^{\top}}{s_{t-1}^{\top}u_{t-1}}\right)\nabla{f(\theta_{t})}+\frac{s_{t-1}s_{t-1}^{\top}}{s_{t-1}^{\top}u_{t-1}}\nabla{f(\theta_{t})}\\\
=&\quad\frac{s_{t-1}s_{t-1}^{\top}}{s_{t-1}^{\top}u_{t-1}}\nabla{f(\theta_{t})}\\\
=&\quad\frac{(r_{t-1}-1)^{2}(\theta_{t-1}-\hat{\theta})(\theta_{t-1}-\hat{\theta})^{\top}}{q(r_{t-1}^{q-1}-1)(r_{t-1}-1)\|A(\theta_{t-1}-\hat{\theta})\|^{q}}qr_{t-1}^{q-1}\|A(\theta_{t-1}-\hat{\theta})\|^{q-2}A^{\top}A(\theta_{t-1}-\hat{\theta})\\\
=&\quad\frac{1-r_{t-1}}{1-r_{t-1}^{q-1}}r_{t-1}^{q-1}(\theta_{t-1}-\hat{\theta})\frac{\|A(\theta_{t-1}-\hat{\theta})\|^{q-2}(\theta_{t-1}-\hat{\theta})^{\top}A^{\top}A(\theta_{t-1}-\hat{\theta})}{\|A(\theta_{t-1}-\hat{\theta})\|^{q}}\\\
=&\quad\frac{1-r_{t-1}}{1-r_{t-1}^{q-1}}r_{t-1}^{q-1}(\theta_{t-1}-\hat{\theta}).\end{split}$
(43)

This simplification implies that for the new iterate $\theta_{t+1}$, we have

$\begin{split}\theta_{t+1}-\hat{\theta}&=\theta_{t}-H_{t}\nabla{f(\theta_{t})}-\hat{\theta}=\theta_{t}-\hat{\theta}-\frac{1-r_{t-1}}{1-r_{t-1}^{q-1}}r_{t-1}^{q-1}\frac{(\theta_{t}-\hat{\theta})}{r_{t-1}}\\\
&=\frac{1-r_{t-1}^{q-2}}{1-r_{t-1}^{q-1}}(\theta_{t}-\hat{\theta})=r_{t}(\theta_{t}-\hat{\theta}),\end{split}$
(44)

where

$r_{t}=\frac{1-r_{t-1}^{q-2}}{1-r_{t-1}^{q-1}}.$ (45)

Therefore, we prove that condition (11) holds for $k=t+1$. By induction, we
prove the linear convergence results in (11) and (12).

One property of this convergence results is that the error vectors
$\\{\theta_{k}-\hat{\theta}\\}_{k=0}^{\infty}$ are parallel to each other with
the same direction as shown in (11). This indicates that the iterations
$\\{\theta_{k}\\}_{k=0}^{\infty}$ converge to the optimal solution
$\hat{\theta}$ along the same straight line defined by
$\theta_{0}-\hat{\theta}$. Only the length of each vector
$\theta_{k}-\hat{\theta}$ reduces to zero with the linear convergence rates
$\\{r_{k}\\}_{k=0}^{\infty}$ specified in (12) and the direction remains all
the same.

###  A.2 Proof of Theorem 2

Recall that for $q\geq 4$, we have

$r_{0}=\frac{q-2}{q-1},\qquad
r_{k}=\frac{1-r_{k-1}^{q-2}}{1-r_{k-1}^{q-1}},\qquad\forall k\geq 1.$ (46)

Define the function $g(r)$ as

$g(r):=\frac{1-r^{q-2}}{1-r^{q-1}},\qquad r\in[0,1].$ (47)

Suppose that $r_{*}\in(0,1)$ satisfying that $r_{*}=g(r_{*})$, which is
equivalent to

$r_{*}^{q-1}+r_{*}^{q-2}=1.$ (48)

Notice that

$g^{\prime}(r)=\frac{(q-1)r^{q-2}-r^{2q-4}-(q-2)r^{q-3}}{(1-r^{q-1})^{2}},$
(49)

and

$\begin{split}&(q-1)r^{q-2}-r^{2q-4}-(q-2)r^{q-3}\\\ =&\quad
r^{q-3}[(q-1)(r-1)-(r^{q-1}-1)]\\\ =&\quad
r^{q-3}(r-1)(q-1-\frac{r^{q-1}-1}{r-1})\\\ =&\quad
r^{q-3}(r-1)(q-1-\sum_{i=0}^{q-2}r^{i}).\end{split}$ (50)

Since $r\in[0,1]$, we have that

$r^{q-3}\geq 0,\quad r-1\leq
0,\quad\sum_{i=0}^{q-2}r^{i}\leq\sum_{i=0}^{q-2}1=q-1.$ (51)

Therefore, we obtain that

$(q-1)r^{q-2}-r^{2q-4}-(q-2)r^{q-3}\leq 0,$ (52)

and

$|g^{\prime}(r)|=\frac{r^{2q-4}+(q-2)r^{q-3}-(q-1)r^{q-2}}{(1-r^{q-1})^{2}}.$
(53)

Our target is to prove that for any $r\in[0,1]$,

$|g^{\prime}(r)|\leq\frac{1}{2}.$ (54)

First, we present the plots of $|g^{\prime}(r)|$ for $r\in[0,1]$ with $4\leq
q\leq 11$ in Figure 4. We observe that for $4\leq q\leq 11$,
$|g^{\prime}(r)|\leq{1}/{2}$ always holds.

(a) $q=4$.

(b) $q=5$.

(c) $q=6$.

(d) $q=7$.

(e) $q=8$.

(f) $q=9$.

(g) $q=10$.

(h) $q=11$.

Figure 4: Plots of $|g^{\prime}(r)|$ with $r\in[0,1]$ for $4\leq q\leq 11$.

Next, we prove that for $q\geq 12$ and any $r\in[0,1]$, we have

$|g^{\prime}(r)|=\frac{(q-1)r^{q-2}-r^{2q-4}-(q-2)r^{q-3}}{(1-r^{q-1})^{2}}\leq\frac{1}{2},$
(55)

which is equivalent to

$r^{2q-2}-2r^{2q-4}-2r^{q-1}+2(q-1)r^{q-2}-2(q-2)r^{q-3}+1\geq 0,\qquad\forall
r\in[0,1].$ (56)

Define the function $h(r)$ as

$h(r):=r^{2q-2}-2r^{2q-4}-2r^{q-1}+2(q-1)r^{q-2}-2(q-2)r^{q-3}+1.$ (57)

We obtain that

$\frac{dh(r)}{dr}=2r^{q-4}h^{(1)}(r),$ (58)

where

$h^{(1)}(r):=(q-1)r^{q+1}-2(q-2)r^{q-1}-(q-1)r^{2}+(q-1)(q-2)r-(q-2)(q-3).$
(59)

Hence, we have that

$\frac{dh^{(1)}(r)}{dr}=(q-1)h^{(2)}(r),$ (60)

where

$h^{(2)}(r):=(q+1)r^{q}-2(q-2)r^{q-2}-2r+q-2.$ (61)

Therefore, we obtain that

$\frac{dh^{(2)}(r)}{dr}=h^{(3)}(r):=(q+1)qr^{q-1}-2(q-2)^{2}r^{q-3}-2,$ (62)

and

$\frac{dh^{(3)}(r)}{dr}=r^{q-4}h^{(4)}(r),$ (63)

where

$h^{(4)}(r):=q(q+1)(q-1)r^{2}-2(q-2)^{2}(q-3).$ (64)

Now we define the function $l(q)$ as

$\begin{split}l(q):=&\quad 2(q-2)^{2}(q-3)-q(q+1)(q-1)\\\ =&\quad
q^{3}-14q^{2}+33q-24\\\ =&\quad q^{2}(q-14)+33(q-1)+9.\end{split}$ (65)

We observe that for $q\geq 14$, we have $l(q)>0$ and we calculate that
$l(12)=84>0$ and $l(13)=236>0$. Hence, we obtain that $l(q)>0$ for all $q\geq
12$, which indicates that for all $r\in[0,1]$,

$r^{2}\leq 1<\frac{2(q-2)^{2}(q-3)}{q(q+1)(q-1)},$ (66)
$q(q+1)(q-1)r^{2}-2(q-2)^{2}(q-3)<0.$ (67)

Therefore, for all $r\in[0,1]$, $h^{(4)}(r)$ defined in (64) satisfies that
$h^{(4)}(r)<0$ and from (63) we know that $\frac{dh^{(3)}(r)}{dr}<0$. Hence,
$h^{(3)}(r)$ defined in (62) is decreasing function and
$h^{(3)}(r)<=h^{(3)}(0)=-2<0$. We know that
$\frac{dh^{(2)}(r)}{dr}=h^{(3)}(r)<0$, which implies that $h^{(2)}(r)$ defined
in (61) is decreasing function. So we have that $h^{(2)}(r)\geq
h^{(2)}(1)=1>0$. From (60) we know that $\frac{dh^{(1)}(r)}{dr}>0$ and
$h^{(1)}(r)$ defined in (59) is increasing function for $r\in[0,1]$. Hence, we
get that $h^{(1)}(r)\leq h^{(1)}(1)=0$ and from (58) we obtain that $h(r)$
defined in(57) is decreasing function for $r\in[0,1]$. Therefore, we have that
$h(r)\geq h(1)=0$ and condition (56) holds for all $r\in[0,1]$, which is
equivalent to $|g^{\prime}(r)|\leq{1}/{2}$.

In summary, we proved that for any $q\geq 12$, we have
$|g^{\prime}(r)|\leq{1}/{2}$. Combining this with the results from Figure 4,
we obtain that $|g^{\prime}(r)|\leq{1}/{2}$ holds for all $q\geq 4$. Applying
Banach’s Fixed-Point Theorem from Lemma 2, we prove the final conclusion (14).

###  A.3 Proof of Theorem 3

Notice that the gradient and the Hessian of the objective function (10) can be
expressed as

$\nabla{f(\theta)}=q\|A\theta-b\|^{q-2}A^{\top}(A\theta-b)=q\|A\theta-b\|^{q-2}A^{\top}A(\theta-\hat{\theta}),$
(68)
$\nabla^{2}{f(\theta)}=q\|A\theta-b\|^{q-2}A^{\top}A+q(q-2)\|A\theta-b\|^{q-4}A^{\top}(A\theta-b)(A\theta-b)^{\top}A.$
(69)

Applying Lemma 1, we can obtain that

$\nabla^{2}{f(\theta)}^{-1}=\frac{(A^{\top}A)^{-1}}{q\|A\theta-b\|^{q-2}}-\frac{(q-2)(\theta-\hat{\theta})(\theta-\hat{\theta})^{\top}}{q(q-1)\|A\theta-b\|^{q}}.$
(70)

Hence, we have that for any $k\geq 1$,

$\theta_{k}=\theta_{k-1}-\nabla{f(\theta_{k-1})}^{-1}\nabla{f(\theta_{k-1})},$
(71)
$\theta_{k}-\hat{\theta}=\theta_{k-1}-\hat{\theta}-\nabla{f(\theta_{k-1})}^{-1}\nabla{f(\theta_{k-1})}.$
(72)

Notice that $b=A\hat{\theta}$ by Assumption 1 and

$\begin{split}&\nabla{f(\theta_{k-1})}^{-1}\nabla{f(\theta_{k-1})}\\\
=&\quad[\frac{(A^{\top}A)^{-1}}{q\|A\theta_{k-1}-b\|^{q-2}}-\frac{(q-2)(\theta_{k-1}-\hat{\theta})(\theta_{k-1}-\hat{\theta})^{\top}}{q(q-1)\|A\theta_{k-1}-b\|^{q}}]q\|A\theta-b\|^{q-2}A^{\top}A(\theta_{k-1}-\hat{\theta})\\\
=&\quad\theta_{k-1}-\hat{\theta}-\frac{q-2}{q-1}\frac{(\theta_{0}-\hat{\theta})^{\top}A^{\top}A(\theta_{k-1}-\hat{\theta})}{\|A\theta_{k-1}-b\|^{2}}(\theta_{k-1}-\hat{\theta})\\\
=&\quad\theta_{k-1}-\hat{\theta}-\frac{q-2}{q-1}\frac{(A\theta_{k-1}-b)^{\top}(A\theta_{k-1}-b)}{\|A\theta_{k-1}-b\|^{2}}(\theta_{k-1}-\hat{\theta})\\\
=&\quad\theta_{k-1}-\hat{\theta}-\frac{q-2}{q-1}(\theta_{k-1}-\hat{\theta}).\end{split}$
(73)

Therefore, we prove the conclusion that for any $k\geq 1$,

$\theta_{k}-\hat{\theta}=\theta_{k-1}-\hat{\theta}-\nabla{f(\theta_{k-1})}^{-1}\nabla{f(\theta_{k-1})}=\frac{q-2}{q-1}(\theta_{k-1}-\hat{\theta}).$
(74)

We observe that the iterations generated by Newton’s method also satisfy the
parallel property, i.e., all vectors
$\\{\theta_{k}-\hat{\theta}\\}_{k=0}^{\infty}$ are parallel to each other with
the same direction.

###  A.4 Proof of Theorem 4

Recall that, we utilize the BFGS for solving the least-square loss function
$\mathcal{L}_{n}$ in equation (16), which is given by:

$\displaystyle\mathcal{L}_{n}(\theta):=\frac{1}{n}\sum_{i=1}^{n}\left(Y_{i}-(X_{i}^{\top}\theta)^{p}\right)^{2}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}Y_{i}^{2}-\frac{2}{n}\sum_{i=1}^{n}Y_{i}(X_{i}^{\top}\theta)^{p}+\frac{1}{n}\sum_{i=1}^{n}(X_{i}^{\top}\theta)^{2p},$
(75)

where $(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$:

$\displaystyle Y_{i}=(X_{i}^{\top}\theta^{*})^{p}+\zeta_{i},$

where $\zeta_{i}\sim\mathcal{N}(0,\sigma^{2})$. In this proof, we focus on the
low SNR regime of the generalized linear model, namely, $\theta^{*}=0$.

Optimal solution $\theta_{n}^{*}$: For the case of $d=1$, we study the optimal
solution $\theta_{n}^{*}$ of the least-square loss function $\mathcal{L}_{n}$.
First of all, directly solving the gradient of the least-square loss function
leads to either $\theta_{n}^{*}=0$ or the following form of $\theta_{n}^{*}$:

$\displaystyle(\theta_{n}^{*})^{p}=\frac{\sum_{i=1}^{n}Y_{i}X_{i}^{p}}{\sum_{i=1}^{n}X_{i}^{2p}}.$
(76)

To bound $\theta_{n}^{*}$, we only focus on bounding the later value of
$\theta_{n}^{*}$ in equation (76). Given the generative model of the data, we
have

$\displaystyle\sum_{i=1}^{n}Y_{i}X_{i}^{p}=\left(\sum_{i=1}^{n}X_{i}^{2p}\right)(\theta^{*})^{p}+\sum_{i=1}^{n}\zeta_{i}X_{i}^{p}.$

Therefore, we obtain that

$\displaystyle(\theta_{n}^{*})^{p}=(\theta^{*})^{p}+\frac{\sum_{i=1}^{n}\zeta_{i}X_{i}^{p}}{\sum_{i=1}^{n}X_{i}^{2p}}.$

We then consider $\frac{1}{n}\sum_{i=1}^{n}\zeta_{i}X_{i}^{p}$ and
$\frac{1}{n}\sum_{i=1}^{n}X_{i}^{2p}$ separately. For the term
$\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}X_{i}^{p}$, note that, for any even
integer $q$, we have

$\displaystyle\mathbb{E}\left[(\zeta
X^{p})^{q}\right]=\mathbb{E}\left[\zeta^{q}X^{pq}\right]\leq(\sigma^{2}q)^{q/2}\cdot(pq)^{pq/2}\leq\left(\max\\{\sigma^{2},p\\}q\right)^{(p+1)q/2}.$

Invoking Lemma 2 in (Mou et al., 2019), we have that

$\displaystyle\mathbb{P}\left[\left|\frac{1}{n}\sum_{i=1}^{n}\zeta_{i}X_{i}^{p}\right|\geq(4\max\\{\sigma^{2},p\\})^{(p+1)/2}\sqrt{\frac{\log
4/\delta}{n}}+\left(\max\\{\sigma^{2},p\\}\log\frac{n}{\delta}\right)^{(p+1)/2}\frac{\log
4/\delta}{n}\right]\leq\delta.$

For the term $\frac{1}{n}\sum_{i=1}^{n}X_{i}^{2p}$, applying Lemma 5 in
(Dwivedi et al., 2020a), we have that

$\displaystyle\mathbb{P}\left[\left|\frac{1}{n}\sum_{i=1}^{n}X_{i}^{2p}-(2p-1)!!\right|\geq\frac{C_{p}\log^{p}(n/\delta)}{\sqrt{n}}\right]\leq\delta,$

where $C_{p}$ is a universal constant depends only on $p$.

For the case $\theta^{*}=0$ which corresponds to the low signal-to-noise
regime, with simple algebra, we know with probability at least $1-2\delta$,

$\displaystyle|\theta_{n}^{*}|\leq\left(\frac{(4\max\\{\sigma^{2},p\\})^{(p+1)/2}\sqrt{\frac{\log
4/\delta}{n}}+\left(\max\\{\sigma^{2},p\\}\log\frac{n}{\delta}\right)^{(p+1)/2}\frac{\log
4/\delta}{n}}{(2p-1)!!-\frac{C_{p}\log^{p}(n/\delta)}{\sqrt{n}}}\right)^{1/p}.$

In summary, it indicates that there exists universal constant $C$ such that
with probability $1-\delta$, we have $|\theta_{n}^{*}|\leq
C\left(\frac{\log^{2p}(n/\delta)}{n}\right)^{\frac{1}{2p}}$.

#### Statistical analysis of the BFGS:

Without loss of generality, we assume that $\theta_{n}^{*}$ takes the value in
equation (76) (the proof when $\theta_{n}^{*}=0$ can be argued in the similar
fashion). When we run the BFGS for solving the least-square loss function
$\mathcal{L}_{n}$, its updates take the following form

$\displaystyle\theta_{k+1}^{n}$
$\displaystyle=\theta_{k}^{n}-\frac{\theta_{k}^{n}-\theta_{k-1}^{n}}{\nabla\mathcal{L}(\theta_{k}^{n})-\nabla\mathcal{L}(\theta_{k-1}^{n})}\nabla\mathcal{L}(\theta_{k}^{n})$
$\displaystyle=\theta_{k}^{n}-\frac{(\frac{1}{n}\sum_{i=1}^{n}X_{i}^{2p})(\theta_{k}^{n})^{2p-1}-(\frac{1}{n}\sum_{i=1}^{n}Y_{i}X_{i}^{p})(\theta_{k}^{n})^{p-1}}{(\frac{1}{n}\sum_{i=1}^{n}X_{i}^{2p})(\sum_{j=0}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j})-(\frac{1}{n}\sum_{i=1}^{n}Y_{i}X_{i}^{p})(\sum_{j=0}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j})}.$

Given the formulation of $\theta_{n}^{*}$ in equation (76), we can rewrite the
update of $\theta_{k+1}^{n}$ as follows:

$\displaystyle\theta_{k+1}^{n}$
$\displaystyle=\theta_{k}^{n}-\frac{(\theta_{k}^{n})^{2p-1}-(\theta_{k}^{n})^{p-1}(\theta_{n}^{*})^{p}}{(\sum_{j=0}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j})-(\sum_{j=0}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j})(\theta_{n}^{*})^{p}}.$
(77)

Assume that $|\theta_{k}^{n}|\geq 2|\theta_{n}^{*}|$ for all $1\leq k\leq T$
where $T$ indicates the first iteration that the BFGS iterates reach the
statistical radius. Without loss of generality, we assume that the
initializations $\theta_{0}^{n}$ and $\theta_{1}^{n}$ of the BFGS satisfy
$\theta_{0}^{n}>\theta_{1}^{n}>2|\theta_{n}^{*}|$, namely, these
initializations are positive and lie above the statistical radius (otherwise,
the conclusion of the theorem trivially holds). Then, we will first
demonstrate by induction that $0<\theta_{k}^{n}<\theta_{k-1}^{n}$ for all
$1\leq k\leq T$. Indeed, the induction hypothesis holds for any $k\leq 1$
based on the conditions of the initializations; we now prove the hypothesis
for $k+1$. From equation (77), we have

$\displaystyle\theta_{k+1}^{n}=\theta_{k}^{n}\left(1-\frac{A_{n}}{B_{n}}\right),$

where we define

$\displaystyle A_{n}=$
$\displaystyle(\theta_{k}^{n})^{p-2}[(\theta_{k}^{n})^{p}-(\theta_{n}^{*})^{p}],$
$\displaystyle B_{n}=$
$\displaystyle\left(\sum_{j=0}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}\right)-\left(\sum_{j=0}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}\right)(\theta_{n}^{*})^{p}.$

Note that, as $\theta_{k}^{n}>0$, we can check that $A_{n}>0$. Furthermore,
direct computation shows that

$\displaystyle B_{n}$
$\displaystyle=\sum_{j=0}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}[(\theta_{k}^{n})^{p}-(\theta_{n}^{*})^{p}]+\sum_{j=p-1}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}$
$\displaystyle\geq\sum_{j=0}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}(2^{p}-1)|\theta_{n}^{*}|^{p}+\sum_{j=p-1}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}>0.$

Similarly, by taking the difference between $A_{n}$ and $B_{n}$, we have

$\displaystyle B_{n}-A_{n}$
$\displaystyle=\sum_{j=1}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}[(\theta_{k}^{n})^{p}-(\theta_{n}^{*})^{p}]+\sum_{j=p-1}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}$
$\displaystyle\geq\sum_{j=1}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}(2^{p}-1)|\theta_{n}^{*}|^{p}+\sum_{j=p-1}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}>0.$

Putting these results together, it indicates that $0<A_{n}/B_{n}<1$;
therefore, we have $0<\theta_{k+1}^{n}<\theta_{k}^{n}$, which implies that the
induction hypothesis holds for $k+1$. As a consequence, we have
$0<\theta_{k}^{n}<\theta_{k-1}^{n}$ for all $0\leq k\leq T$.

We now prove the following useful lower bound of $\theta_{k+1}^{n}$ in terms
of $\theta_{k}^{n}$ when $1\leq k\leq T-1$:

$\displaystyle\theta_{k+1}^{n}\geq\frac{p}{p+1}\theta_{k}^{n}.$ (78)

Indeed, the inequality (78) is equivalent to show that

$\displaystyle B_{n}-A_{n}\geq\frac{p}{p+1}B_{n}.$

In light of the previous calculations, this inequality is equivalent to

$\displaystyle\sum_{j=1}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}[(\theta_{k}^{n})^{p}-(\theta_{n}^{*})^{p}]+\sum_{j=p-1}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}$
$\displaystyle\geq\frac{p}{p+1}\left(\sum_{j=0}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}[(\theta_{k}^{n})^{p}-(\theta_{n}^{*})^{p}]+\sum_{j=p-1}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}\right),$

which is equivalent to

$\displaystyle\sum_{j=1}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}[(\theta_{k}^{n})^{p}-(\theta_{n}^{*})^{p}]+\sum_{j=p-1}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}\geq
p(\theta_{k}^{n})^{p-2}[(\theta_{k}^{n})^{p}-(\theta_{n}^{*})^{p}].$ (79)

It is clear that the RHS of equation (79) is less than or equal to
$p(\theta_{k}^{n})^{2p-2}$ while the RHS of equation (79) is at least
$\sum_{j=p-1}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}$, which is
greater than or equal to $p(\theta_{k}^{n})^{2p-2}$. Therefore, we obtain the
conclusion of the bound (79) and equivalently the bound (78).

Now, we would like to demonstrate the linear convergence rate of the BFGS
iterates $\\{\theta_{k}^{n}\\}$ to the optimal solution $\theta_{n}^{*}$:

$\displaystyle|\theta_{k+1}^{n}-\theta_{n}^{*}|\leq\left(1-\frac{\left(1-\frac{1}{p+1}\right)^{2p-2}}{2p^{2}-p}\right)|\theta_{k}^{n}-\theta_{n}^{*}|,$
(80)

for all $2\leq k\leq T-1$. As $\theta_{k}^{n}>C|\theta_{n}^{*}|$ for all
$2\leq k\leq T$, the result in equation (80) is equivalent to

$\displaystyle(\theta_{k+1}^{n}-\theta_{n}^{*})\leq\left(1-\frac{\left(1-\frac{1}{p+1}\right)^{2p-2}}{2p^{2}-p}\right)(\theta_{k}^{n}-\theta_{n}^{*}),$
(81)

for all $2\leq k\leq T-1$. From direct computation, we have

$\displaystyle\theta_{k+1}^{n}-\theta_{n}^{*}=(\theta_{k}^{n}-\theta_{n}^{*})$
$\displaystyle\cdot\left(1-\frac{(\theta_{k}^{n})^{p-1}(\sum_{j=0}^{p-1}(\theta_{k}^{n})^{p-1-j}(\theta_{n}^{*})^{j})}{(\sum_{j=0}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j})-(\sum_{j=0}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j})(\theta_{n}^{*})^{p}}\right).$

As $0<C|\theta_{n}^{*}|<\theta_{k}^{n}<\theta_{k-1}^{n}$, we have

$\displaystyle\left(\sum_{j=0}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}\right)-\left(\sum_{j=0}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}\right)(\theta_{n}^{*})^{p}$
$\displaystyle=$
$\displaystyle\sum_{j=0}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}[(\theta_{k}^{n})^{p}-(\theta_{n}^{*})^{p}]+\sum_{j=p-1}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}$
$\displaystyle=$
$\displaystyle\sum_{j=0}^{p-2}(\theta_{k}^{n})^{p-2-j}(\theta_{k-1}^{n})^{j}\left[(\theta_{k}^{n}-\theta_{n}^{*})\left(\sum_{l=0}^{p-1}(\theta_{k}^{n})^{l}(\theta_{n}^{*})^{p-1-l}\right)\right]+\sum_{j=p-1}^{2p-2}(\theta_{k}^{n})^{2p-2-j}(\theta_{k-1}^{n})^{j}$
$\displaystyle\leq$ $\displaystyle
p(p-1)(\theta_{k-1}^{n})^{2p-3}(\theta_{k}^{n}-\theta_{n}^{*})+p(\theta_{k-1}^{n})^{2p-2}$
$\displaystyle\leq$ $\displaystyle(2p^{2}-p)(\theta_{k-1}^{n})^{2p-2}.$

Meanwhile, when $2\leq k\leq T$ we have that

$\displaystyle(\theta_{k}^{n})^{p-1}\left(\sum_{j=0}^{p-1}(\theta_{k}^{n})^{p-1-j}(\theta_{n}^{*})^{j}\right)\geq(\theta_{k}^{n})^{2p-2}\geq\left(1-\frac{1}{p+1}\right)^{2p-2}(\theta_{k-1}^{n})^{2p-2}.$

Putting the above results together, as long as $2\leq k\leq T-1$ we find that

$\displaystyle(\theta_{k+1}^{n}-\theta_{n}^{*})\leq\left(1-\frac{\left(1-\frac{1}{p+1}\right)^{2p-2}}{2p^{2}-p}\right)(\theta_{k}^{n}-\theta_{n}^{*}),$

which shows the linear convergence of the BFGS iterates to the optimal
solution.

By repeating the inequalities (81), we obtain that

$\displaystyle|\theta_{T}^{n}-\theta_{n}^{*}|\leq\left(1-\frac{\left(1-\frac{1}{p+1}\right)^{2p-2}}{2p^{2}-p}\right)^{T-1}|\theta_{1}^{n}-\theta_{n}^{*}|.$

As long as we choose
$\left(1-\frac{\left(1-\frac{1}{p+1}\right)^{2p-2}}{2p^{2}-p}\right)^{T-1}(\theta_{1}^{n}-\theta_{n}^{*})=|\theta_{n}^{*}|$,
which is equivalently to
$T=1+\frac{\log(|\theta_{n}^{*}|/|\theta_{1}^{n}-\theta_{n}^{*}|)}{\log\left(1-\frac{\left(1-\frac{1}{p+1}\right)^{2p-2}}{2p^{2}-p}\right)}$,
then we obtain $|\theta_{T}^{n}-\theta_{n}^{*}|\leq|\theta_{n}^{*}|$. A direct
application of the triangle inequality, we have

$|\theta_{T}^{n}-\theta^{*}|\leq|\theta_{T}^{n}-\theta_{n}^{*}|+|\theta_{n}^{*}|\leq
2|\theta_{n}^{*}|\leq
2C\left(\frac{\log^{2p}(n/\delta)}{n}\right)^{\frac{1}{2p}}$

with probability $1-\delta$. As a consequence, we obtain the conclusion of the
theorem.

###  A.5 Additional Experiment Results

Figure 5: Illustration of the statistical radius of BFGS with error bar. Left:
Low SNR regime. Right: High SNR regime.

We additionally report the statistical radius with error bar, by sampling $40$
different dataset with the same generative model and report the median, the
$25\%$ percentile and $75\%$ percentile. For low SNR regime, BFGS still reach
the statistical radius within $\mathcal{O}(n^{-1/4})$, while for high SNR
regime, BFGS reach the statistical radius within $\mathcal{O}(n^{-1/2})$.