On the acceleration of the Barzilai-Borwein method

First, let us investigate the behavior of gradient method (8) with $\Psi(A)=I$ (i.e., the BB1 method). Particularly, we apply it to the non-random quadratic minimization problem proposed in de2014efficient , which has the form (4) with a diagonal matrix $A$ given by

and $b$ being a null vector. Here, $ncond=\log_{10}\kappa$ and $\kappa>0$ is the condition number of $A$. We set $n=10$, $\kappa=10^{3}$ and use $(10,10,\ldots,10)^{T}$ as the initial point. The iteration was stopped once the gradient norm is reduced by a factor of $10^{-6}$. Denote the $i$-th component of $g_{k}$ by $g_{k}^{(i)}$ and the indices of the components of $g_{k}$ with two largest magnitudes by $i_{1}$ and $i_{2}$, respectively. Then, the percentage of the magnitudes of the first two largest components to that of the whole gradient can be computed by

This $\Upsilon(g_{k})$ is plotted in Fig. 1 (left), where we can see that $\Upsilon(g_{k})\geq 0.8$ holds for more than half of the iterations (145 out of 224 total iterations). Hence, roughly speaking, the searches of the BB1 method are often dominated in some two-dimensional subspaces. The history of index $i_{1}$ against the iteration number is also plotted in Fig. 1 (right), where we can see that $|g_{k}^{(i_{1})}|$ corresponds more frequently to the largest eigenvalues $\lambda_{10}$ or $\lambda_{9}$. Since

and $1/\lambda_{n}\leq\alpha_{k}\leq 1/\lambda_{1}$, the history of $i_{1}$ in Fig. 1 (right) in fact indicates that, most stepsizes generated by the BB1 method are often much larger than $1/\lambda_{10}$ or $1/\lambda_{9}$. As a result, the BB1 method may need many iterations to reduce the corresponding components of the gradients $g_{k}^{(9)}$ or $g_{k}^{(10)}$.

In huang2019asymptotic , we showed that a family of gradient methods including SD and MG will asymptotically reduce their searches in a two-dimensional subspace and could be accelerated by exploiting certain orthogonal properties in this two-dimensional subspace. In a similar spirit, we could also accelerate the convergence of the class of gradient methods (8) in a lower dimensional subspace if certain orthogonal properties hold.

Suppose that, for a given $k>0$, there exists a $q_{k}$ satisfying

Since this $q_{k}$ is also invariant under translations and rotations, for later analysis we may still assume $A$ in (12) is diagonal as in (9). The following lemma shows a generalized orthogonal property for $q_{k}$ and $g_{k+1}$, which is a key property for deriving our new stepsize in the next subsection.

In this subsection, we derive a new stepsize based on the iterations of gradient method (8). We show that combining the new stepsize with gradient method (8), we can achieve finite termination for minimizing two-dimensional strongly convex quadratic functions.

By Lemma 1, we have that $g_{k}^{T}\Psi(A)q_{k-1}=0$ for $k>0$. Now, suppose both $\Psi^{r}(A)q_{k-1}$ and $\Psi^{1-r}(A)g_{k}$ are nonzero vectors, where $r\in\mathbb{R}$. Let us consider to minimize the function $f$ in a two-dimensional subspace spanned by $\frac{\Psi^{r}(A)q_{k-1}}{\|\Psi^{r}(A)q_{k-1}\|}$ and $\frac{\Psi^{1-r}(A)g_{k}}{\|\Psi^{1-r}(A)g_{k}\|}$, and let

where

and

Denote the components of $H_{k}$ by $H_{k}^{(ij)}$, $i,j=1,2$ and notice that $B_{k}B_{k}^{T}=I$ by $g_{k}^{T}\Psi(A)q_{k-1}=0$. Then, we have the following finite termination theorem.

Notice that by setting $k_{0}=2$ in the above Theorem 2.1, the new gradient method in Theorem 2.1 will find the exact minimizer in at most $5$ iterations when minimizing a two-dimensional strongly convex quadratic function. In fact, since $\Delta=\lambda_{1}\lambda_{2}$ and $H_{k}^{(11)}+H_{k}^{(22)}=\lambda_{1}+\lambda_{2}$, the equation (19) has two positive roots $1/\lambda_{1}$ and $1/\lambda_{2}$. This observation allows us to use the stepsize $\tilde{\alpha}_{k_{0}}$ with some retards as stated in the following corollary, which would lead us a more convenient way for choosing stepsizes when the objective function is not quadratic.

By setting $\Psi(A)=I$, $\Psi(A)=A$ and $r=1/2$ in (16), and setting $k_{0}=k$ in (17), we can derive the following two stepsizes:

and

respectively, where

By (23) and (24), we have

Hence, both $\tilde{\alpha}_{k}^{BB1}$ and $\tilde{\alpha}_{k}^{BB2}$ are short monotone steps for reducing the value and gradient norm of the objective function, respectively. And it follows from Theorem 2.1 that by inserting the monotone steps $\tilde{\alpha}_{k}^{BB1}$ and $\tilde{\alpha}_{k}^{BB2}$ into the BB1 and BB2 methods, respectively, the gradient method will have finite termination for minimizing two-dimensional strongly convex quadratic functions.

To numerically verify this finite termination property, we apply the method (8) with $\Psi(A)=I$ (i.e., the BB1 method) and $\tilde{\alpha}_{2}^{BB1}$ given by (23) to minimize a two-dimensional quadratic function (4) with

We run the algorithm for five iterations using ten random starting points. The averaged values of $\|g_{5}\|$ and $f(x_{5})$ are presented in Table 1. Moreover, we also run BB1 method for a comparison purpose. We can observe that for different values of $\lambda$, the values of $\|g_{5}\|$ and $f(x_{5})$ obtained by BB1 method with $\tilde{\alpha}_{2}^{BB1}$ given by (23) are numerically very close to zero. However, even for the case $\lambda=10$, $\|g_{5}\|$ and $f(x_{5})$ obtained by pure BB1 method are far away from zero. These numerical results coincide with our analysis and show that the nonmonotone method (8) can be significantly accelerated by incorporating proper monotone steps.

As mentioned in Section 2.1, the stepsizes $\alpha_{k}^{BB1}$ generated by the BB1 method may be far away from the reciprocals of the largest eigenvalues of the Hessian matrix $A$ of the quadratic function (4). In other words, the stepsize $\alpha_{k}^{BB1}$ may be too large to effectively decrease the components of gradient $g_{k}$ corresponding to the first several largest eigenvalues, which, by (11), can be greatly reduced when small stepsizes are employed. In addition, it has been observed by many works in the recent literature that gradient methods using long and short stepsizes adaptively generally perform much better than those using only one type of stepsizes, for example see dai2003altermin ; de2014efficient ; di2018steplength ; gonzaga2016steepest ; huang2019gradient ; huang2019asymptotic . So, we would like to develop gradient methods that combines the two nonmonotone BB stepsizes with the short monotone stepsize given by (17).

We first extend the orthogonal property developed in Lemma 1 and the finite termination result given in Theorem 2.1.

Based on Lemma 2 and using the same arguments as those in the proof of Theorem 2.1, we can obtain the following finite termination result even different function $\Psi$’s are used in (8) to obtain the stepsizes.

Theorem 3.1 allows us to incorporate the nonmonotone BB stepsizes $\alpha_{k}^{BB1}$ and $\alpha_{k}^{BB2}$, and the short monotone stepsize $\tilde{\alpha}_{k}^{BB2}$ in one gradient method. Alternate or adaptive scheme has been employed for choosing long and short stepsizes in BB-type methods dai2005projected ; zhou2006gradient . And recent studies show that adaptive strategies are more preferred than the alternate scheme dhl2018 ; zhou2006gradient . Hence, we would like develop adaptive strategies to choose proper stepsizes for our new gradient methods. In particular, our adaptive nonmonotone gradient method (ANGM) takes the long BB stepsize $\alpha_{k}^{BB1}$ when $\alpha_{k}^{BB2}/\alpha_{k}^{BB1}\geq\tau_{1}$ for some $\tau_{1}\in(0,1)$. Otherwise, a short stepsize $\alpha_{k}^{BB2}$ or $\tilde{\alpha}_{k}^{BB2}$ will be taken depending on the ratio $\|g_{k-1}\|/\|g_{k}\|$. Notice that $\alpha_{k}^{BB2}$ minimizes the gradient in the sense that

So, when $\|g_{k-1}\|/\|g_{k}\|>\tau_{2}$ for some $\tau_{2}>1$, i.e. the gradient norm decreases, the previous stepsize $\alpha_{k-1}$ is often a reasonable approximation of $\alpha_{k}^{BB2}$. By our numerical experiments, when BB method is applied the searches are often dominated in some two-dimensional subspaces. And the new gradient method in Theorem 3.1 would have finite convergence for minimizing two-dimensional quadratic function when the new stepsize $\tilde{\alpha}_{k}^{BB2}$ is applied after some BB2 steps. Hence, our ANGM would employ the new monotone stepsize $\tilde{\alpha}_{k}^{BB2}$ when $\|g_{k-1}\|\geq\tau_{2}\|g_{k}\|$; otherwise, certain BB2 steps should be taken. In practice, we find that when $\|g_{k-1}\|\leq\tau_{2}\|g_{k}\|$, ANGM often has good performance by taking the stepsize $\min\{\alpha_{k}^{BB2},\alpha_{k-1}^{BB2}\}$. To summarize, our ANGM applies the following adaptive strategies for choosing stepsizes:

Notice that the calculation of $\tilde{\alpha}_{k}^{BB2}$ needs to compute $\alpha_{k}^{MG}$ which is not easy to obtain when the objective function is not quadratic. In stead, the calculation of $\tilde{\alpha}_{k-1}^{BB2}$ will just require $\alpha_{k}^{BB2}$, which is readily available even for general objective function. Moreover, it is found in recent research that gradient methods using retard stepsizes can often lead better performances friedlander1998gradient . Hence, in the first variant of ANGM, we simply replace $\tilde{\alpha}_{k}^{BB2}$ in (30) by $\tilde{\alpha}_{k-1}^{BB2}$, i.e. the stepsizes are chosen as

We call the gradient method using stepsize (31) ANGR1. On the other hand, since the calculation of $\tilde{\alpha}_{k-1}^{BB2}$ also needs $\hat{\alpha}_{k-2}$ and $\Gamma_{k-1}$ and by (26),

to simplify ANGR1, we may further replace $\tilde{\alpha}_{k-1}^{BB2}$ in (31) by its upper bound in (32). As a result, we have the second variant of ANGM, which chooses stepsizes as

We call the gradient method using stepsize (33) ANGR2.

In terms of global convergence for minimizing quadratic function (4), by (26), we can easily show the $R$-linear global convergence of ANGM since it satisfies the property in dai2003alternate . Similarly, $R$-linear convergence of ANGR1 and ANGR2 can be also established. See the proof of Theorem 3 in dhl2018 for example.