Direct Nonlinear Acceleration

TL;DR

This paper introduces direct nonlinear acceleration (DNA), a new method that improves optimization speed by directly minimizing the function value at extrapolated points, outperforming previous RNA methods in various datasets.

Contribution

The paper proposes a novel direct nonlinear acceleration technique that minimizes the function value at extrapolated points, offering better performance than existing RNA methods.

Findings

DNA significantly outperforms RNA on synthetic datasets

DNA accelerates neural network training effectively

Computational cost of DNA is comparable to RNA

Abstract

Optimization acceleration techniques such as momentum play a key role in state-of-the-art machine learning algorithms. Recently, generic vector sequence extrapolation techniques, such as regularized nonlinear acceleration (RNA) of Scieur et al., were proposed and shown to accelerate fixed point iterations. In contrast to RNA which computes extrapolation coefficients by (approximately) setting the gradient of the objective function to zero at the extrapolated point, we propose a more direct approach, which we call direct nonlinear acceleration (DNA). In DNA, we aim to minimize (an approximation of) the function value at the extrapolated point instead. We adopt a regularized approach with regularizers designed to prevent the model from entering a region in which the functional approximation is less precise. While the computational cost of DNA is comparable to that of RNA, our direct approach significantly outperforms RNA on both synthetic and real-world datasets. While the focus of this paper is on convex problems, we obtain very encouraging results in accelerating the training of neural networks.

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }f(x),$$

$$x_{k+1}=x_k-{\alpha }_k\nabla f(x_k),$$

$$c^{\star }\approx \arg \underset{c}{\min }\left\{\Vert \nabla f\left(\sum _{k=0}^K{c}_k{x}_k\right)\Vert :c\in {\mathbb{R}}^{K+1},\sum _{k=0}^K{c}_k=1\right\},$$

$$\boxed{\underset{c\in {\mathbb{R}}^{K+1}}{\min }f\left(\sum _{k=0}^K{c}_k{x}_k\right)+\lambda g(c),}$$

$$c^{\star }=\arg \underset{c}{\min }\left\{\Vert \sum _{k=0}^K{c}_k\nabla f\left(x_k\right)\Vert =\Vert \sum _{k=0}^K{c}_k{\tilde{R}}_k\Vert :c\in {\mathbb{R}}^{K+1},\sum _{k=0}^K{c}_k=1\right\},$$

$$\boxed{\underset{c\in {\mathbb{R}}^{K+1}}{\min }f\left(Xc\right),}$$

$$\nabla f(Xc)$$

$$\boxed{\underset{c\in {\mathbb{R}}^{K+1}}{\min }f(Xc)+\lambda 1_S(c)=\underset{c}{\min }\left\{f\left(Xc\right):c\in {\mathbb{R}}^{K+1},\sum _{k=0}^K{c}_k=1\right\},}$$

$$\underset{\begin{array}{c}c\in {\mathbb{R}}^{K+1}\end{array}}{\min }f\left(Xc\right)+\frac{\lambda }{2}\Vert Xc-y{\Vert }^2,$$

$$(X^{\top }R+\lambda X^{\top }X)c=\lambda X^{\top }y-X^{\top }\nabla f(0).$$

$$\underset{\begin{array}{c}c\in {\mathbb{R}}^{K+1}\end{array}}{\min }f\left(Xc\right)+\frac{\lambda }{2}\Vert c-e{\Vert }^2,$$

$$\underset{x}{\min }f(x):=\frac{1}{2}\Vert Ax-y{\Vert }^2$$

$$\underset{x}{\min }f(x):=\frac{1}{2}\Vert Ax-y{\Vert }^2+\frac{1}{2n}\Vert x{\Vert }^2,$$

$$f_i(x)=\log (1+\exp (-y_i⟨A(:,i),x⟩)+\frac{1}{2m}\Vert x{\Vert }^2.$$

$$x_{k+1}=\Phi (x_k).$$

$$\boxed{{\hat{c}}^{(k)}=\arg \underset{c}{\min }\Vert F^{k}c\Vert \mathrm{subject}\mathrm{to}\sum _i{c}_i=1,}$$

$$\nabla h(c)$$

$$\nabla f(Xc)$$

$$L(c,\lambda )=h(c)+\lambda \left(\sum _{k=0}^K{c}_k-1\right),$$

$$\nabla L_x(x,\lambda )$$

$$\nabla L_{\lambda }(x,\lambda )$$

$$\nabla f\left(Xc\right)=\sum _{k=0}^K{c}_k\nabla f\left(x_k\right).$$

$$\nabla f\left(x_k\right)=\frac{x_k-x_{k+1}}{{\alpha }_k}.$$

$$\nabla f\left(Xc\right)=\sum _{k=0}^K{c}_k\frac{x_k-x_{k+1}}{{\alpha }_k}=\tilde{R}c.$$

$$c=-\lambda {\left(X^{\top }\tilde{R}\right)}^{-1}\text{mathbbm}1.$$

$$\lambda =\frac{-1}{\text{mathbbm}{1}^{\top }{\left(X^{\top }\tilde{R}\right)}^{-1}\text{mathbbm}1},$$

$$c=\frac{{\left(X^{\top }\tilde{R}\right)}^{-1}\text{mathbbm}1}{\text{mathbbm}{1}^{\top }{\left(X^{\top }\tilde{R}\right)}^{-1}\text{mathbbm}1}.$$

$$f(Xc)=\frac{1}{2}{c}^{\top }{X}^{\top }AXc,\text{ for DNA}{c}_{\mathrm{D}}=0,\text{and for DNA-1}{c}_{D1}=\frac{z}{\text{mathbbm}{1}^{\top }z}\text{ where }z=(X^{\top }\tilde{R}{)}^{-1}\text{mathbbm}1.$$

$$f_D=0,\text{ and }{f}_{D1}=\frac{\text{mathbbm}{1}^{\top }({\tilde{R}}^{\top }X{)}^{-1}{X}^{\top }AX(X^{\top }\tilde{R}{)}^{-1}\text{mathbbm}1}{2(\text{mathbbm}{1}^{\top }z{)}^2},$$

$$f_{D1}=\frac{1}{2\text{mathbbm}{1}^{\top }(X^{\top }AX{)}^{-1}\text{mathbbm}1}.$$

$$c_{\mathrm{R}}=\frac{({\tilde{R}}^{\top }\tilde{R}{)}^{-1}\text{mathbbm}1}{\text{mathbbm}{1}^{\top }({\tilde{R}}^{\top }\tilde{R}{)}^{-1}\text{mathbbm}1}.$$

Full text

1 Introduction

In this paper we consider the generic unconstrained minimization problem

[TABLE]

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a smooth objective function and bounded from below. One of the most fundamental methods for solving (1) is gradient descent (GD), on which many state-of-the-art methods are based. Given current iterate $x_{k}\in\mathbb{R}^{n}$, the update rule of GD is

[TABLE]

where $\alpha_{k}>0$ is a stepsize. The efficiency of GD depends on further properties of $f$. Assuming $f$ is $L$–smooth and $\mu$–strongly convex, for instance, the iteration complexity of GD is ${\cal O}(\kappa{\rm log}(1/\epsilon))$, where $\kappa=L/\mu$ and $\epsilon$ is the target error tolerance. However, it is known that GD is not the “optimal” gradient type method: it can be accelerated.

The idea of accelerating converging optimization algorithms can track its history back to 1964 when Polyak proposed his “heavy ball” method (Polyak, 1964). In 1983, Nesterov proposed his accelerated version for general convex optimization problems. Comparing with Polyak’s method, Nesterov’s method gives acceleration for general convex and smooth problems and the iteration complexity improves to ${{\cal O}(1/\sqrt{\epsilon})}$ (Nesterov, 1983). In 2009, Beck and Teboulle proposed fast iterative shrinkage thresholding algorithm (FISTA) (Beck and Teboulle, 2009) that uses Nesterov’s momentum coefficient and accelerates proximal type algorithms to solve a more complex class of objective functions that combine a smooth, convex loss function (not necessarily differentiable) and a strongly convex, smooth penalty function (also see (Nesterov, 2007, 2013)). To develop further insights into Nesterov’s method, Su et al. (Su et al., 2014) examined a continuous time 2nd-order ODE which at its limit reduces to Nesterov’s accelerated gradient method. In addition, Lin et al. (Lin et al., 2015) introduced a generic approach known as catalyst that minimizes a convex objective function via an accelerated proximal point algorithm and gains acceleration in Nesterov’s sense. (Bubeck et al., 2015) proposed a geometric alternative to gradient descent that is inspired by ellipsoid method and produces acceleration with complexity ${{\cal O}(1/\sqrt{\epsilon})}$. Recently, (Zhu and Orecchia, 2017) used a linear coupling of gradient descent and mirror descent and claimed to attend acceleration in Nesterov’s sense as well. In contrast, the sequence acceleration techniques accelerate a sequence independently from the iterative method that produces this sequence. In other words, these techniques take a sequence $\{x_{k}\}$ and produce an accelerated sequence based on the linear combination of $x_{k}{\rm s}$ such that the new accelerated sequence converges faster than the original. In the same spirit, recently, Scieur et al. (Scieur et al., 2016, 2018) proposed an acceleration technique called regularized nonlinear acceleration (RNA). Scieur et al.’s idea is based on Aitken’s $\Delta^{2}$-algorithm (Aitken, 1927) and Wynn’s $\epsilon$-algorithm (Wynn, 1956) (or recursive formulation of generalized Shanks transform (Shanks, 1955; Wynn, 1956; Brezinski et al., 2018)). To achieve acceleration, Scieur et al. considered a technique known as minimum polynomial approximation and they assumed a linear model for the iterates near the optimum. They also proposed a regularized variant of their method to stabilize it numerically. The intuition behind the regularized nonlinear acceleration of Scieur et al. is very natural. To minimize $f$ as in (1), they considered the sequence of iterates $\{x_{k}\}_{k\geq 0}$ is generated by a fixed-point map. If $x^{\star}$ is a minimizer of $f$, $\nabla f(x^{\star})=0$, and hence through extrapolation one can find:

[TABLE]

such that the next (accelerated) point can be generated as a linear combination of $K+1$ previous iterates: $x=\sum_{i=0}^{K}c_{i}^{\star}x_{i}.$ We review RNA in detail in Section 2.

Notation. We denote the $\ell_{2}$-norm of a vector $x$ by $\|x\|$ and define $\|x\|_{M}$ by $\|x\|_{M}:=\sqrt{x^{\top}Mx}.$

1.1 Contributions

We highlight our main contributions in this paper as follows:

Direct nonlinear acceleration (DNA). Inspired by Anderson’s acceleration technique (Anderson, 1965) (see Appendix for a brief description of Anderson’s acceleration) and the work of Scieur et al. (Scieur et al., 2016), we propose an extrapolation technique that accelerates a converging iterative algorithm. However, in contrast to (Scieur et al., 2016), we find the extrapolation coefficients $c^{\star}$ by directly minimizing the function at the linear combination of $K+1$ iterates $\{x_{k}\}_{k\geq 0}^{K}$ with respect to $c\in\mathbb{R}^{K+1}$. In particular, for a given sequence of iterates $\{x_{k}\}_{k\geq 0}^{K}$ we propose to approximately solve:

[TABLE]

where $\lambda>0$ is a balancing parameter and $g$ is a penalty function. As our approach tries to minimize the functional value directly, we call it as direct nonlinear acceleration (DNA). We also note that our formulation shares some similarities with (Riseth, 2019; Zhang et al., 2018). However, unlike (Riseth, 2019), we do not require line search and check a decrease condition at each step of our algorithm. On the other hand, Zhang et al. (Zhang et al., 2018) do not consider a direct acceleration scheme as they deal with a fixed-point problem.

Regularization. We propose several versions of DNA by varying the penalty function $g(c)$. This helps us to deal with the numerical instability in solving a linear system as well as to control errors in gradient approximation. In our first version, we let $g(c)=1_{S}(c)$, where $S:=\{c\;:\;\sum_{i}c_{i}=1\}$ and $1_{S}(c)=0$ if $c\in S$, while $1_{S}(c)=+\infty$ otherwise. Later, we propose two regularized constraint-free versions to find a better minimum of the function $f$ by expanding the search space of extrapolating coefficients to $\mathbb{R}^{K+1}$ rather than restricting them over the space $S$. To this end, the first constraint-free version adds a quadratic regularization ${g(c)=\left\|\sum_{i=0}^{K}c_{i}x_{i}-y\right\|^{2}}$ to the objective function, where $y$ is a reference point and $g(c)$ controls how far we want the linear combination $\sum_{i}c_{i}x_{i}$ to deviate from $y$. In the second constraint-free version, we add the regularization directly on $c$. We add a quadratic term of the form $g(c)={\|c-e\|^{2}}$ to the objective function, where $e$ is a reference point to $c$ and $g(c)$ controls how far we want $c$ to deviate from $e$. In contrast, the regularized version of RNA only considers a ridge regularization $\|c\|^{2}$ for numerical stability. Trivially, we note that by setting $e=0$, we recover the regularization proposed in RNA. We argue that by using a different penalty function $g(c)$ as regularizer our DNA is more robust than RNA.

Quantification between RNA and DNA in minimizing quadratic functions by using GD iterates. If $g(c)=0$ or $g(c)=1_{\sum_{i}c_{i}=1}$, in terms of the functional value, we always obtain a better accelerated point than RNA. Moreover, the acceleration obtained by DNA can be theoretically directly implied from the existing results of Scieur et al.(Scieur et al., 2016). If $g(c)=0$, we show by a simple example on quadratic functions that DNA outperforms RNA by an arbitrary large margin. If $g(c)=1_{\sum_{i}c_{i}=1}$, we also quantify the functional values obtained from both RNA and DNA for quadratic functions and provide a bound on how DNA outperforms RNA in this setup.

Numerical results. Our empirical results show that for smooth and strongly convex functions, minimizing the functional value converges faster than RNA. In practice, our acceleration techniques are robust and outperform that of Scieur et al. (Scieur et al., 2016) by large margins in almost all experiments on both synthetic and real datasets. To further push the robustness of our methods, we test them on nonconvex problems as well. As a proof of concept, we trained a simple neural network classifier on MNIST dataset (LeCun et al., 2010) via GD and accelerate the GD iterates via the online scheme in (Scieur et al., 2016) for both RNA and DNA. Next, we train ResNet18 network (He et al., 2016) on CIFAR10 dataset (Krizhevsky and Hinton, 2009) by SGD and accelerate the SGD iterates via the online scheme in (Scieur et al., 2016) for both RNA and DNA. In both cases, DNA outperform RNA in lowering the generalization errors of the networks.

2 Regularized Nonlinear Acceleration

In RNA, one solves (3) by assuming that the gradient can be approximated by linearizing it in the neighborhood of $\{x_{k}\}_{k=0}^{K}.$ Thus, by assuming $\sum_{k=0}^{K}c_{k}=1$, the relation ${\left\|\nabla f\left(\sum_{k=0}^{K}c_{k}x_{k}\right)\right\|\approx\left\|\sum_{k=0}^{K}c_{k}\nabla f\left(x_{k}\right)\right\|}$ holds. Hence, one can approximately solve (3) via:

[TABLE]

where $\tilde{R}_{k}$ is the $k^{\rm th}$ column of the matrix $\tilde{R}$, which holds $\nabla f\left(x_{k}\right).$ Moreover (5) does not need an explicit access to the gradient and it can be seen as an approximated minimal polynomial extrapolation (AMPE) as in (Cabay and Jackson, 1976; Scieur et al., 2016, 2018). In this context, we should note that Scieur et al. indicated that the summability condition $c^{\top}\mathbbm{1}=1$ is not restrictive, where $\mathbbm{1}$ is a vector of all 1s. If the sequence $\{x_{k}\}$ is generated via GD (as in (2)), then $\tilde{R}=\left[\nicefrac{{(x_{0}-x_{1})}}{{\alpha_{0}}},\ldots,\nicefrac{{(x_{K}-x_{K+1})}}{{\alpha_{K}}}\right]$. Also, if $\tilde{R}^{\top}\tilde{R}$ is nonsingular, then the minimizer of (5) is explicitly given as: ${c^{\star}=\frac{(\tilde{R}^{\top}\tilde{R})^{-1}\mathbbm{1}}{\mathbbm{1}^{\top}(\tilde{R}^{\top}\tilde{R})^{-1}\mathbbm{1}}}$. If $\tilde{R}^{\top}\tilde{R}$ is singular then $c$ is not necessarily unique. Any $c$ of the form $\frac{z}{z^{\top}\mathbbm{1}}$, where $z$ is a solution of $\tilde{R}^{\top}\tilde{R}z=\mathbbm{1}$, is a solution of (5). To deal with the numerical instabilities and the case when the matrix $\tilde{R}^{\top}\tilde{R}$ is singular, Scieur et al. proposed to add a regularizer of the form $\lambda\|c\|^{2}$ to their problem, where $\lambda>0$. As a result, $c^{\star}$ is unique and given as ${c^{\star}=\frac{(\tilde{R}^{\top}\tilde{R}+\lambda I)^{-1}\mathbbm{1}}{\mathbbm{1}^{\top}(\tilde{R}^{\top}\tilde{R}+\lambda I)^{-1}\mathbbm{1}}}.$ The numerical procedure of RNA is given in Alg 1. For further details about RNA we refer the readers to (Scieur et al., 2016, 2018). Scieur et al. also explained several acceleration schemes to use with Algorithm 1.

3 Direct Nonlinear Acceleration

Instead of minimizing the norm of the gradient, we propose to minimize the objective function $f$ directly to obtain the coefficients $\{c_{k}\}$. We set $g(c)=0$ in (4) and we propose to solve the unconstrained minimization problem

[TABLE]

where $X=[x_{0},\ldots,x_{K}]$. We call problem (6) as direct nonlinear acceleration (DNA) without any constraint. If $f$ is quadratic, then we have the following lemma:

Lemma 1**.**

Let the objective function $f$ be quadratic and let $\{x_{k}\}$ be the iterates produced by (2) to minimize $f$. Then $c$ is a solution of the linear system $X^{\top}Rz=-X^{\top}\nabla f(0),$ where $R\in\mathbb{R}^{n\times(K+1)}$ is a matrix such that its $i^{\rm th}$ column is $R_{i}=\frac{x_{i}-x_{i+1}}{\alpha_{i}}-\nabla f(0)$ and $X=[x_{0},\ldots,x_{K}]$.

If $f$ is non-quadratic then we can approximately solve problem (6) by approximating its gradient by a linear model. In fact, we use the following approximation $\nabla f(x)\approx A(x-y_{x})+\nabla f(y_{x}),$ where we assume that $x$ is close to $y_{x}$ and $A$ is an approximation of the Hessian. Therefore, by setting $x=Xc$ and $y_{x}=y$ in the above, we have $\nabla f(Xc)\approx A(Xc-y)+\nabla f(y)=\sum_{i}c_{i}Ax_{i}-Ay+\nabla f(y),$ where $y$ is a referent point that is assumed to be in the neighborhood of $Xc$. For instance, one may choose $y$ to be $x_{K}$. Let $x_{i-1}$ be a referent point for $x_{i}$, that is, assume that $\nabla f(x_{i})\approx A(x_{i}-x_{i-1})+\nabla f(x_{i-1})$. Then one can show that $Ax_{i}=\nabla f(x_{i})-\nabla f(0)$. As a result, we have

[TABLE]

Therefore, from the first optimality condition and by using (7), we conclude that the solutions of (6) can be approximated by the solutions of the linear system $X^{\top}Rz=-X^{\top}\nabla f(0)$. We describe the numerical procedure in Alg 1 in the Appendix.

Comments on the convergence of DNA.

Let $H$ be the Hessian of $f$, where we assume that $f$ is quadratic. Also let $\lambda_{\rm max}(H)$ be the maximum eigenvalue of $H$

Lemma 2**.**

Let $c_{D}$ and $c_{R}$ be the extrapolation coefficients produced by DNA and RNA, respectively. Then ${\|Xc_{D}-x^{\star}\|_{H}^{2}\leq\|Xc_{R}-x^{\star}\|_{H}^{2}\leq\lambda_{\rm max}(H)\|Xc_{R}-x^{\star}\|^{2}.}$

Theorem 1**.**

Denote $\xi=\nicefrac{{(\sqrt{L}-\sqrt{\mu})}}{{(\sqrt{L}+\sqrt{{\mu}})}}.$ It is given in (Scieur et al., 2016) that for the coefficients $c_{R}$ produced by Alg 1: ${\|Xc_{R}-x^{\star}\|^{2}\leq\kappa(H)\tfrac{2\xi^{k}}{1+\xi^{2k}}\|x_{0}-x^{\star}\|^{2}.}$ Further, for Alg 3 we have ${\|Xc_{D}-x^{\star}\|_{H}^{2}\leq\lambda_{\rm max}(H)\kappa(H)\tfrac{2\xi^{k}}{1+\xi^{2k}}\|x_{0}-x^{\star}\|^{2}.}$

Remark 1*.*

The convergence rate for DNA for quadratic functions in Theorem 1 is the same as that for Krylov subspace methods (for example, conjugate gradient algorithm) up to a multiplicative scalar.

However, numerically DNA is unstable like RNA without regularization. In fact, the matrix $X^{\top}R$ can be very ill-conditioned and can lead to large errors in computing $c^{\star}$. Moreover, we accumulate errors in approximating the gradient via linearization as our approximation of the gradient is valid only in the neighborhood of the iterates $x_{0},\ldots,x_{K}$. To solve these problems, we propose three regularized versions of DNA by using three different regularizers in the form of $g(c)$ and show that they work well in practice. But one can explore different forms of $g(c)$ as regularizer. We explain them in the following sections.

3.1 DNA-1

This regularized version of DNA is directly influenced by Scieur et al. (Scieur et al., 2016). Here, we generate the extrapolated point $x$ as a linear combination of the set of $K+1$ iterates such that, $x=\sum_{k}c_{k}x_{k}$. Additionally, as in (Scieur et al., 2016, 2018), we assume the sum of the coefficients $c_{k}$ to be equal to 1. Therefore, for $c\in\mathbb{R}^{K+1}$ with sum of its elements equal to 1, we set $g(c)=1_{\sum_{i}c_{i}=1}$ in (4) and consider the following constrained problem:

[TABLE]

where $X=[x_{0},\ldots,x_{K}]$. We call this version of DNA as DNA-1.

Lemma 3**.**

If the objective function $f$ is quadratic and $X^{\top}\tilde{R}$ is nonsingular then $c=(X^{\top}\tilde{R})^{-1}\mathbbm{1}/\delta,$ where $\delta=\mathbbm{1}^{\top}(X^{\top}\tilde{R})^{-1}\mathbbm{1},$ and $\mathbbm{1}$ is the vector of dimension $K+1$ with all the components equal to 1 and ${\tilde{R}=\left[\nicefrac{{(x_{0}-x_{1})}}{{\alpha_{0}}},\ldots,\nicefrac{{(x_{K}-x_{K+1})}}{{\alpha_{K}}}\right]}$.

Similar to RNA, if $X^{\top}\tilde{R}$ is singular then $c$ is not necessarily unique. Any $c$ of the form $\frac{z}{z^{\top}\mathbbm{1}}$, where $z$ is a solution of $X^{\top}\tilde{R}z=\mathbbm{1}$, is a solution of (8). DNA-1 is described in Alg 3.

Comparison with RNA on simple quadratic functions.

Denote the functional value obtained by DNA, DNA-1 (Alg 3) and RNA (Alg 1) at an extrapolated point as $f_{D}$, $f_{D1}$ and $f_{R}$, respectively.

Proposition 1**.**

Let $A\in\mathbb{R}^{n\times n}$ be symmetric and positive definite and $f(x)=\frac{1}{2}x^{\top}Ax$ be a quadratic objective function. Let $X=[x_{0}\;\;x_{1}\;\cdots x_{k}]$ be a matrix generated by stacking $k$ iterates of GD to minimize $f$. Then the functional value of DNA, DNA-1 and RNA at the accelerated point are: $f_{D}=0$, ${f_{\rm D1}=\frac{1}{2\mathbbm{1}^{\top}(X^{\top}AX)^{-1}\mathbbm{1}}},$ and ${f_{\rm R}=\frac{\mathbbm{1}^{\top}(X^{\top}A^{2}X)^{-1}X^{\top}AX(X^{\top}A^{2}X)^{-1}\mathbbm{1}}{2(\mathbbm{1}^{\top}(X^{\top}A^{2}X)^{-1}\mathbbm{1})^{2}}},$ respectively.

We conclude that for this simple objective function, DNA reaches the optimal solution after the first acceleration. Moreover, one can choose the matrix $A$ such that $f_{R}$ is arbitrary large, and this example shows that DNA may outperform RNA by a large margin. The comparison between DNA-1 and RNA on the previous example is given in the following lemma and theorem.

Lemma 4**.**

We assume that the matrix $\tilde{R}$ has full column rank. With the notations used in Proposition 1, we have ${\nicefrac{{f_{{\rm R}}}}{{f_{{\rm D1}}}}=\|z\|_{A^{-1}}^{2}\|y\|^{2}\|z\|^{-4},}$ where ${z:=(\tilde{R}^{\dagger})^{\top}\mathbbm{1}=((AX)^{\dagger})^{\top}\mathbbm{1}}$ and ${y:=((A^{1/2}X)^{\dagger})^{\top}\mathbbm{1}}$. We have $y^{\top}A^{-1/2}z=z^{\top}z$; then, by using Cauchy-Schwarz inequality, we conclude that ${{\|z\|_{A^{-1}}\|y\|_{2}}\geq y^{\top}A^{-1/2}z={\|z\|_{2}^{2}}}$ whence $\nicefrac{{f_{\rm R}}}{{f_{\rm D1}}}\geq 1$.

Note that the ratio ${\frac{f_{\rm R}}{f_{\rm D1}}\geq 1}$ can be directly concluded from the definition of ${f_{\rm R}}$ and ${f_{\rm D1}}$. The main goal of the previous lemma is to exactly quantify the ratio between these two quantities. The following theorem gives more insight.

Theorem 2**.**

We have ${\frac{f_{{\rm R}}}{f_{{\rm D1}}}\leq U_{R}:=\|z\|_{A^{-1}}^{2}\|z\|_{A}^{2}\|z\|^{-4},}$ and $U_{R}\in\left[\nicefrac{{1}}{{2}}+\nicefrac{{\kappa(A)}}{{2}},\kappa(A)\right],$ where $\kappa(A)$ is the condition number of $A$.

The above theorem tells us, for a simple quadratic function, the ratio of the objective function values of DNA-1 and RNA may attain an order of $\kappa(A)$, but it never exceeds it. The theoretical quantification of the acceleration obtained by DNA and its different versions compared to RNA in more general problems is left for future work. Although DNA-1 can be seen as a regularized version of DNA, we still need to remedy the fact that the linearization of the gradient is not a good approximation in the entire space, and that the matrix $X^{\top}\tilde{R}$ may be singular. To this end, we impose some regularization such that the new extrapolated point stays near to some reference point. We propose two different ways in the following two sections.

3.2 DNA-2

We set $g(c)=\|Xc-y\|^{2}$ in (4) and consider a regularized version of problem (6):

[TABLE]

where $\lambda>0$ is a balancing parameter and $y$ is a reference point (a point supposed to be in the neighborhood of $Xc$). By taking the derivative of the objective in (9) with respect to $c$ and setting it to 0, we find ${X^{\top}\nabla f(Xc)+\lambda X^{\top}(Xc-y)=0},$ which after using the approximation (7) becomes ${X^{\top}(Rc+\nabla f(0))+\lambda X^{\top}(Xc-y)=0.}$ Finally, $c^{\star}$ is given as a solution to the linear system

[TABLE]

In general, $X^{\top}R$ is not necessarily symmetric. To justify the regularization further, one might symmetrize $X^{\top}R$ by its transpose. In our experiments, we obtained good performance without this.

We call this method DNA-2 (see Alg 4). Note that $X^{\top}R+\lambda X^{\top}X$ can be singular, especially near the optimal solution. To remedy this, we propose either to add another regularization to the problem (9), or to consider a direct regularization on $c$ instead of $Xc$. We explain this next.

3.3 DNA-3

We set $g(c)=\|c-e\|^{2}$ in (4) and consider a regularized version of (6) as

[TABLE]

where $\lambda>0$ and $e$ is a reference point for $c$. By taking the derivative with respect to $c$ and setting it to 0, we find ${X^{\top}\nabla f(Xc)+\lambda(c-e)=0},$ which after using the approximation (7) becomes ${X^{\top}(Rc+\nabla f(0))+\lambda X^{\top}(c-e)=0.}$ Therefore, $c^{\star}$ is given as a solution to the linear system: ${(X^{\top}R+\lambda I)c=\lambda e-X^{\top}\nabla f(0).}$ We call this method DNA-3, and describe it in Alg 5.

4 Numerical illustration

We evaluate our techniques and compare against RNA and GD by using both synthetic data as well as real-world datasets. Overall, we find that DNA outperforms RNA in most settings by large margins.

Experimental setup. Our experimental setup comprises of 3 typical problems, least squares, ridge regression, and logistic regression, for which the optimal solution $x^{\star}$ is either known or can be evaluated using a numerical solver. We apply the online acceleration scheme in (Scieur et al., 2016) and compare 3 versions of DNA against RNA and GD. Our results show the difference between the functional values at the extrapolated point and at the optimal solution on a logarithmic scale (the lower the better), as the iterations progress. The primary objective of our simulations is to show the effectiveness of DNA and its different versions to accelerate a converging, deterministic optimization algorithm. Therefore, we do not report any computation time of the algorithms and we do not claim these implementations are optimized. Note that the computation bottleneck of all algorithms (including RNA) is solving the linear system to calculate $c$, and because the dimensionality of the linear systems is the same in RNA and DNA, the extra cost is the same in both approaches. In our experiments, we consider a fixed stepsize $\alpha_{k}=1/L$ for GD, where $L$ is the Lipschitz constant of $\nabla f$. We note that for DNA-1 and 2 we need to use the stepsize explicitly to construct $R$ as defined in Lem 1.

Least Squares. We consider a least squares regression problem of the form

[TABLE]

where $A\in\mathbb{R}^{m\times n}$ with $m>n$ is the data matrix, $y\in\mathbb{R}^{m}$ is the response vector. For $m>n$ and ${\rm rank}(A)=n$, the objective function $f$ in (12) is strongly convex. The optimal solution $x^{\star}$ to (12) is given by ${x^{\star}=\operatorname{argmin}_{x}f(x)=(A^{\top}A)^{-1}A^{\top}y.}$ For least squares we only consider the overdetermined systems, that is, $m>n$.

Ridge Regression. The classic ridge regression problem is of the form:

[TABLE]

where $A\in\mathbb{R}^{m\times n}$ is the data matrix, $y\in\mathbb{R}^{n}$ is the response vector. The optimal solution $x^{\star}$ to (13) is given by $x^{\star}=\operatorname{argmin}_{x}f(x)=(A^{\top}A+\tfrac{1}{2n}I)^{-1}A^{\top}y.$

Logistic Regression. In logistic regression with $\ell_{2}$ regularization, the objective function $f(x)$ is the summation of $n$ loss function of the form:

[TABLE]

We use the MATLAB function fminunc to numerically obtain the minimizer of $f$ in this case.

Synthetic Data. To compare the performance of different methods under different acceleration schemes, we are interested in the case where matrix $A$ has a known singular value distribution and we consider the cases where $A$ has varying condition numbers. We note that the condition number of $A$ is defined as $\kappa(A):=\nicefrac{{\lambda_{\max}(A)}}{{\lambda_{\min}(A)}}$, where $\lambda$ is the eigenvalue of $A$. We first generate a random matrix and let $U\Sigma V^{\top}$ be its SVD. Next we create a vector $S\in\mathbb{R}^{\min\{m,n\}}$ with entries $s_{i}\in\mathbb{R}^{+}$ arranged in an nonincreasing order such that $s_{1}$ is maximum and $s_{\min\{m,n\}}$ is minimum. Finally, we form the test matrix $A$ as $A=U{\rm diag}(s_{1}\;\;s_{2}\;\cdots s_{\min\{m,n\}})V^{\top}$ such that $A$ will have a higher condition number if $\nicefrac{{s_{1}}}{{s_{\min\{m,n\}}}}$ is large and smaller condition number if $\nicefrac{{s_{1}}}{{s_{\min\{m,n\}}}}$ is small. We create the vector $y$ as a random vector.

Real Data. We use 15 different real-world datasets from the LIBSVM repository (Chang and Lin, 2011). We set apart the datasets with $y$-labels as $\{-1,1\}$ for Logistic regression and used the remaining 12 multi-label datasets for least squares and ridge regression problems. We use the matrix $A$ in its crude form, that is, without any scaling/normalizing or centralizing its rows or columns.

Acceleration Results. We use GD as our baseline algorithm and the online acceleration scheme as explained in (Scieur et al., 2016) for both RNA and DNAs to accelerate the GD iterates. For synthetic data (see Fig 1), we see that for smaller condition numbers and for quadratic objective functions, DNA-1 and RNA has almost similar performance and DNA-2 and 3 show faster decrease of $f(x_{k})-f(x^{\star})$, but all of them are very competitive. As the condition number of the problems becomes huge, DNA-2 and 3 outperform RNA by large margins. However, for logistic regression problems we see performance gains for all versions of DNA compared to RNA. Though for huge condition numbers, for logistic regression problems, the performance of DNA-2 depends on the hyperparameter $\lambda$. We argue with experimental evidence as in Fig 1 that for huge condition numbers the sensitivity of the performance of DNA-2 is problem specific. We use an additional regularizer $\epsilon\|c\|^{2}$, where $\epsilon\approx 10^{-14}$ to find a stable solution to (10) of DNA-2. Next, on real-world datasets in Figs  2 and 3, we see that all versions of DNA outperform RNA, except in a few cases, where RNA and DNA-1 have almost similar performance. We indicate the oscillating nature of DNA-2 in some plots is due to its problem-specific sensitivity to the regularizer. In Fig 4, we find for logistic regression problems on real datasets, DNA outperfoms RNA. We owe the success of DNA on non-quadratic problems to its adaptive gradient approximation. We note that the performance of all algorithms on the offline scheme of (Scieur et al., 2016) are similar to online scheme of (Scieur et al., 2016). However, on the second online scheme used in (Scieur et al., 2018), all the algorithms perform extremely poorly. Therefore, we do not report the results in this paper.

Application to the non-convex world: Accelerating neural network training. Modern deep learning requires optimization algorithms to work in a nonconvex setup. Although this is not the main goal of this paper, nevertheless, we implement our acceleration techniques for training neural networks and obtain surprisingly promising results. We only use DNA-1 for experiments in this section. Tuning the hyperparameter $\lambda$ for the other versions of DNAs requires more time, and we leave this for future research. The Pytorch implementation of RNA is based on (Scieur et al., 2018). Finally, see Fig 7 in Appendix for more results.

MNIST Classification. First, we trained a simple two-layer neural network classifier on MNIST dataset (LeCun et al., 2010) via GD and accelerate the GD iterates via the online scheme in (Scieur et al., 2016) for both RNA and DNA-1. The two-layer neural network is wildely adopted in most tutorials that use MNIST dataset 111https://github.com/pytorch/examples/blob/master/mnist/main.py. In Fig 5 (a), DNA-1 gains acceleration by using GD iterates with a window size $k=5$. However, RNA fails to accelerate the GD iterates. This motivated us to train the same network on MNIST dataset classification (LeCun et al., 2010) via SGD as baseline algorithm and accelerate the SGD iterates via the online scheme in (Scieur et al., 2016) for both RNA and DNA-1 (as in Fig 5 (b)). Again, with window size $k=5$, DNA-1 achieves better acceleration than RNA.

ResNet18 on CIFAR10. Finally, we train the ResNet18 network (He et al., 2016) on CIFAR10 dataset (Krizhevsky and Hinton, 2009) by SGD. Each epoch of SGD consists of multiple iterations and each iteration applies to $128$ training samples. The size of the training set is $5\times 10^{4}$ and the size of validation set is $10^{4}$. Each sample is a $32\times 32$ resolution color image and they are categorized into 10 classes. We accelerate the SGD iterates via the online scheme in (Scieur et al., 2016) for both RNA and DNA-1. Again DNA-1 outperforms RNA in lowering the generalization error of the network (see Fig 5 (c)).

Appendix

Appendix A Anderson’s Acceleration (Anderson, 1965)

There are several acceleration techniques that have been proposed in the literature and they pose a lot of similarities. We quote the authors from (Brezinski et al., 2018) – “Methods for accelerating the convergence of various processes have been developed by researchers in a wide range of disciplines, often without being aware of similar efforts undertaken elsewhere.” In 1965 Anderson’s acceleration was designed to accelerate Picard iteration for electronic structure computations. Because it is relevant in our current work, we give a brief description of it for completeness.

For a given sequence of iterate $\{x_{k}\}$ with $x_{k}\in\mathbb{R}^{n}$ and a mapping $\Phi(\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, the fixed-point algorithm generates a recursive update of the iterates as:

[TABLE]

Let there be $m_{k}+1$ evaluations of the fixed point map $\phi$. Anderson’s acceleration technique computes a new iteration as a linear combination of the previous $m_{k}+1$ evaluations. We explain it formally in Alg 6. In Alg 6, $m$ is considered as a hyperparameter that sets the quantity $m_{k}$ as $\min\{m,k\}$, where $k$ is the iteration counter and $m$ is known as the depth. This is used to determine the window size to compute $\hat{c}$–the coefficients for linear combination of the fixed point evaluations. In other words, in each iteration, by solving the optimization problem:

[TABLE]

one can obtain the extrapolation coefficients $\hat{c}^{(k)}$ that help to determine the accelerated point $x_{k+1}.$ Toth and Kelley pointed out that, in principle, any norm can be used in the minimization step (Toth and Kelley, 2015).

The summability of the coefficients $c_{i}$ or the normalization condition was not explicitly mentioned in the original work of Anderson. Because $c_{i}$’s can be determined up to a multiplicative scalar, one can impose the normalization condition. However, it does not restrict generality. We refer the readers to (Kelley, 2018; Toth and Kelley, 2015; Walker and Ni, 2011; Anderson, 1965; Brezinski et al., 2018) for a comprehensive idea of Anderson’s acceleration technique.

Appendix B Acceleration Schema

In this section, we explain different acceleration schema used by Scieur et al. in (Scieur et al., 2016, 2018) for completeness.

Remark 2*.*

For GD, the updates of the iterates are done via the simple update rule (2), which is explained in Fig 6(a).

B.0.1 Online Scheme 1

We explain the acceleration scheme proposed in (Scieur et al., 2016) herein. First, we run $k$ iterations of GD to produce the sequence of iterates $\{x_{i}\}_{i=1}^{k}$ and then use extrapolation to generate a new point $x_{k}^{\prime}$.We use $x_{k}^{\prime}$ as the initial point of GD and produce a set of next $k$ iterates via GD. At this end, we further use the extrapolation scheme to produce a second offline update $x_{k+1}^{\prime}$ which is used as next the initial point of GD, and this process continues. See Fig 6(b).

B.0.2 Online Scheme 2

The acceleration scheme proposed in (Scieur et al., 2018), is more involved than the one propsed in (Scieur et al., 2016). First, we run $k$ iterations of GD to produce a sequence of iterates $\{x_{i}\}_{i=1}^{k}$ and then use extrapolation to generate $x_{k}^{\prime}$. Next, we use $x_{k}^{\prime}$ as the starting point of GD to produce $x_{k+1}$. Now, we start from the second iterate $x_{2}$ and consider a set of $k$ iterates $\{x_{2},x_{3},\cdots,x_{k},x_{k+1}\}$ to produce the second offline update $x_{k+1}^{\prime}$ via extrapolation which is to be used as the next starting point of GD, and this process continues. See Fig 6(c).

B.0.3 Offline scheme

Lastly, we describe an offline update scheme, as illustrated in Fig 6(d). First, we run the GD to produce the sequence of iterates $\{x_{k}\}$ and then use the acceleration on the set of first $k$ iterates to produce the first offline update $x_{k}^{\prime}$ and concatenate it with the previous $(k-1)$ GD updates. Next, we start from the second iterate $x_{2}$ and consider a set of $k$ iterates to produce the second offline update $x_{k+1}^{\prime}$ via acceleration and this process continues. As a result, the offline accelerated updates are generated as $\{x_{1},x_{2}\cdots,x_{k}^{\prime},x_{k+1}^{\prime},\cdots\}$.

Proof.

of Lemma 1. Let $h(c)=f(Xc)$, from the first order optimality condition we have

[TABLE]

For quadratic objective function the gradient is affine, i.e

[TABLE]

By using the relation between the iterates of GD method we find $\nabla f\left(x_{k}\right)=\frac{x_{k}-x_{k+1}}{\alpha_{k}}.$ Hence $\nabla f(Xc)=Rc+\nabla f(0).$ By injecting this in the first order optimality condition we get the result. ∎

Proof.

of Lemma 2.

Since $f$ is quadratic then $f(x)=f(x^{\star})+\|x-x^{\star}\|_{H}^{2}$. Therefore, from the definition of $c_{R}$ and $c_{D}$ and using Proposition 2.2 in (Scieur et al., 2016) we conclude the result. ∎

Proof.

of Lemma 3.

The Lagrangian of the problem (5) is

[TABLE]

where $\lambda>0$ is the Lagrange multiplier. The first order optimality conditions are

[TABLE]

For quadratic objective functions, the gradient is affine and because $\sum_{k=0}^{K}c_{k}=1$ we have

[TABLE]

By using the relation between the iterates of GD method we find

[TABLE]

By using the above expression in equation (18) we further get

[TABLE]

Substituting (19) in the first optimality condition and solving for $c$ we get

[TABLE]

Next we use it in the second optimality condition and solve it for $\lambda$ to find

[TABLE]

and therefore the final expression for $c$ is

[TABLE]

∎

Appendix C Example with Quadratic function.

Let $f(x)=\frac{1}{2}x^{\top}Ax$, where $A$ is symmetric and positive definite. We know $\nabla f(x)=Ax.$ By using the extrapolation we find the coefficients $c_{i}$s such that $x=\sum_{i}c_{i}x_{i}=Xc$, where $X=[x_{0}\;\;x_{1}\;\cdots x_{k}]$ is a matrix generated by stacking $k$ iterates as its column and $c\in\mathbb{R}^{k}$ is a vector of coefficients. We know

[TABLE]

Therefore, we find

[TABLE]

which for $\tilde{R}=[\nabla f(x_{0})\;\;\nabla f(x_{1})\;\cdots\nabla f(x_{k})]=[Ax_{0}\;\;Ax_{1}\;\cdots Ax_{k}]=AX$ further reduces to

[TABLE]

Similarly, we find for RNA

[TABLE]

Therefore,

[TABLE]

which further reduces to

[TABLE]

In order to prove Lemma 4 we need the following Lemma.

Lemma 5**.**

*If the sequence of iterates $\{x_{k}\}$ are linearly independent then we have:

(i) the matrices $AX$ and $A^{\frac{1}{2}}X$ have full column ranks.

(ii) $(A^{1/2}X)^{\dagger}A^{-1/2}((AX)^{\dagger})^{\top}=(X^{\top}A^{2}X)^{-1}.$*

Proof.

(i) Since $A$ is symmetric and positive definite, ${\rm rank}(A)={\rm rank}(A^{1/2})=n$. As the iterates $\{x_{k}\}$ are linearly independent, $X=[x_{0},\ldots,x_{K}]\in\mathbb{R}^{n\times(K+1)}$ has full column rank. Therefore, the matrices $AX$ and $A^{\frac{1}{2}}X$ have full column ranks.

(ii) We know if a matrix $B$ is of full column rank then $B^{\dagger}=(B^{\top}B)^{-1}B^{\top}.$ By using the above and (i) and we find

[TABLE]

Hence the result. ∎

Proof.

of Lemma 2. We have $\tilde{R}=AX$. Since $A$ is symmetric and positive definite, it is invertible and $X=A^{-1}R.$ Set $y=((A^{1/2}X)^{\dagger})^{\top}\mathbbm{1}$. and let $\tilde{R}^{\dagger}$ be the pseudo-inverse of $R$. Therefore, $\tilde{R}^{\dagger}$ can be computed as

[TABLE]

and $(\tilde{R}^{\dagger})^{\top}$ is

[TABLE]

We also note that $\tilde{R}^{\dagger}\tilde{R}=I_{k}$ and $\tilde{R}^{\top}(\tilde{R}^{\dagger})^{\top}=I_{k}$, where $I_{k}$ is an identity matrix of size $k$. Therefore, we have

[TABLE]

Since, $(X^{\top}A^{1/2})^{\top}=A^{1/2}X$, by using the property of pseudo-inverse, we can write

[TABLE]

and the above expression becomes

[TABLE]

Similarly we find, $(X^{\top}A)^{\top}=AX$, and again by using the property of pseudo-inverse, we can write

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

Combining (20) and (21) we obtain the ratio between $f_{R}$ and $f_{D}$.

From lemma 5 we have $y^{\top}A^{-1/2}z=z^{\top}z$ then by uisng Cauchy Swartz inequality we conclude that ${\|z\|_{A^{-1}}\|y\|_{2}}\geq y^{\top}A^{-1/2}z={\|z\|_{2}^{2}}$ whence $\frac{f_{\rm R}}{f_{\rm D1}}\geq 1$. ∎

Proof.

of Theorem 2. Recall that $\frac{f_{{\rm R}}}{f_{{\rm D1}}}=\frac{\|z\|_{A^{-1}}^{2}\|y\|_{2}^{2}}{\|z\|_{2}^{4}}.$ Also recall that $z:=(\tilde{R}^{\dagger})^{\top}\mathbbm{1}=((AX)^{\dagger})^{\top}\mathbbm{1}$ and $y:=((A^{1/2}X)^{\dagger})^{\top}\mathbbm{1}$. Therefore, $A^{1/2}z$ is the minimum norm solution to the linear system: $X^{\top}A^{1/2}A^{1/2}z=\mathbbm{1}$ and similarly, $y$ is the minimum norm solution to the linear system: $X^{\top}A^{1/2}y=\mathbbm{1}$. By using the above fact, we find $\|y\|_{2}\leq\|A^{1/2}z\|_{2}=\|z\|_{A}$ and we can rewrite the ratio as:

[TABLE]

From (22) the quantity $\displaystyle{\max_{z\neq 0}\frac{\|z\|_{A^{-1}}^{2}\|z\|_{A}^{2}}{\|z\|_{2}^{4}}}$ is equivalent to $\displaystyle{\max_{\|z\|_{2}=1}{\|z\|_{A^{-1}}^{2}\|z\|_{A}^{2}}}\leq\displaystyle{\max_{\|z\|_{2}=1}\|z\|_{A^{-1}}^{2}\max_{\|z\|_{2}=1}\|z\|_{A}^{2}}:=U_{R}.$

Note that $U_{R}\leq{\lambda_{\max}(A^{-1})}\lambda_{\max}(A)=\frac{\lambda_{\max}(A)}{\lambda_{\min}(A)}=\kappa(A).$

Let $U\Sigma U^{\top}$ be an eigenvalue decomposition of $A$ then $U_{R}=\displaystyle{\max_{\|z\|_{2}=1}\|z\|_{\Sigma^{-1}}^{2}\max_{\|z\|_{2}=1}\|z\|_{\Sigma}^{2}}$. By considering a vector with $1/\sqrt{2}$ at the first and last position and zero everywhere else we conclude that

[TABLE]

∎

Appendix D Reproducible research

See the LIBSVM dataset from the repository online: https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/. See the source code of RNA from: https://github.com/windows7lover/RegularizedNonlinearAcceleration. For MATLAB and Pytorch code that is used to produce all the results for our DNA, please email the authors.