The long time behavior and the rate of convergence of symplectic convex algorithms obtained via splitting discretizations of inertial damping systems

6 years ago

(preprint), Optimization (numerical), paper

Abstract

In this paper we propose new numerical algorithms in the setting of unconstrained optimization problems and we give proof for the rate of convergence in the iterates of the objective function.

Furthermore, our algorithms are based upon splitting and symplectic methods and they preserve the energy properties of the inherent continuous dynamical system that contains a Hessian perturbation.

At the same time, we show that Nesterov gradient method is equivalent to a Lie-Trotter splitting applied to a Hessian driven damping system.

Finally, some numerical experiments are presented in order to validate the theoretical results.

Authors

Cristian Daniel Alecsa
Tiberiu Popoviciu Institute of Numerical Analysis Romanian Academy Cluj-Napoca, Romania
Romanian Institute of Science and Technology, Cluj-Napoca, Romania

Keywords

unconstrained optimization; rate of convergence; inertial algorithms; Nesterov; convex function; Hessian; dynamical system; splitting operator; vanishing damping; Lie-Trotter.

Paper coordinates

Cristian-Daniel Alecsa, The long time behavior and the rate of convergence of symplectic convex algorithms obtained via splitting discretizations of inertial damping systems, arXiv:2001.10831

PDF

https://arxiv.org/pdf/2001.10831

About this paper

Journal

Publisher Name

DOI

Print ISSN

Online ISSN

google scholar link

[1] C.D. Alecsa, S.C. Laszlo, T. Pinta, An extension of the second order dynamical system that models Nesterov’s convex gradient method, 2019,arxiv:1908.02574

[2] C.D. Alecsa, I. Boros, T. Pinta, New optimization algorithms for neural network training using operator splitting techniques, 2019, arxiv:1904.12952

[3] F. Alvarez, H. Attouch, J. Bolte, P. Redont, A second-order gradient like dissipative dynamical system with Hessian-driven damping. Application to optimization and mechanics, J. Math. Pures Appl. 81 (2002), no. 8, pp. 747-779.

[4] H. Attouch, Z. Chbani, J. Fadili, H. Riahi, First-order optimization algorithms via inertial systems with Hessian driven damping, 2019, arxiv:1907.10536

[5] H. Attouch, A. Cabot, Asymptotic stabilization of inertial gradient dynamics with time-dependent viscosity, J. Differential Equations 263 (2017), pp. 5412-5458, doi:10.1016/j.jde.2017.06.024

[6] H. Attouch, A. Cabot, Convergence rates of inertial forward-backward algorithms, SIAM J. Optim (28) 2018, no.1, pp. 849-874, https://doi.org/10.1137/17M1114739

[7] H. Attouch, X. Goudou, P. Redont, The heavy ball with friction method. I. The continuous dynamical system: global exploration of the local minima of a real-valued function by asymptotic analysis of a dissipative system, Commun. Contemp. Math. 2 (2000), pp. 1-34.

[8] H. Attouch, J. Peypouquet, P. Redont, Fast convex minimization via inertial dynamics with Hessian driven damping, J. Differential Equations 261 (2016), pp. 5734-5783.

[9] N. Bansal, A. Gupta, Potential-function proofs for gradient methods, Theory of Computing (2019), pp. 1-32.

[10] A. Batkai, P. Csomos, B. Farkas, G. Nickel, Operator splitting for non-autonomous evolution equations, J. Functional Analysis 260 (2011), pp. 2163-2190.

[11] A. Beck, Introduction to Nonlinear Optimization: Theory, Algorithms, and Applications with MATLAB, vol. 19, SIAM, 2014.

[12] S. Blanes, F. Casas, A. Murua, Splitting and composition methods in the numerical integration of differential equations, Bol. Soc. Esp. Mat. Apl. 45 (2008), pp. 89-145.

[13] R.I. Bot, R. Csetnek, S.C. Laszlo, A second order dynamical approach with variable damping to nonconvex smooth minimization, Applic. Analysis, 2018, doi:10.1080/00036811.2018.1495330

[14] A. Cabot, H. Engler, S. Gadat, On the long time behavior of second-order differential equations with asymptotically small dissipation, Trans. Amer. Math. Soc. 361 (2009), pp. 5983-6017.

[15] A. Cabot, H. Engler, S. Gadat, Second order differential equations with asymptotically small dissipation and piecewise flat potentials, Electron J. Differential Equations 17 (2009), pp. 33-38.

[16] C. Castera, J. Bolte, C. Fevotte, E. Pauwels, An inertial Newton algorithm for deep learning, 2019, arXiv:1905.12278.

[17] A. Chambolle, C. Dossal, On the convergence of the iterates of FISTA, HAL Id : hal-01060130, 2014, https://hal.inria.fr/hal-01060130v3https://hal.inria.fr/hal-01060130v3.

[18] L. Chen, H. Luo, First order optimization methods based on Hessian-driven Nesterov accelerated gradient flow, 2019, arxiv:1912.09276.

[19] A. Defazio, On the curved geometry of accelerated optimization, Adv. Neural Inf. Processing Syst. 33 (NIPS 2019), 2019.

[20] K. Feng, On difference schemes and symplectic geometry, Proceedings of the 5-th Intern. Symposium on differential geometry & differential equations, August 1984, Beijing (1985), pp. 42-58.

[21] G. Franca, D.P. Robinson, R. Vidal, Gradient flows and accelerated proximal splitting methods, 2019, arXiv:1908.00865.

[22] G. Franca, J. Sulam, D.P. Robinson, R. Vidal, Conformal symplectic and relativistic optimization, 2019, arXiv:1903.04100.

[23] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration. Structure-Preserving Algorithms for ODE’s, second edition, Springer-Verlag Berlin Heidelberg, 2006, doi:10.1007/3-540-30666-8.

[24] E. Hairer, G. Wanner, Euler methods, explicit, implicit and symplectic. In : Encyclopedia of Applied and Computational Mathematics, B. Engquist ed., Springer, Berlin Heidelberg, 2015, pp. 451-455.

[25] E. Hansen, F. Kramer, A. Ostermann, A second-order positivity preserving scheme for semilinear parabolic problems, Appl. Num. Math. 62 (2012), pp. 1428-1435.

[26] E. Hansen, A. Ostermann, Exponential splitting for unbounded operators, Math. of Comput. 78 (2009), no. 267, pp. 1485-1496.

[27] S.C. Laszlo, Convergence rates for an inertial algorithm of gradient type associated to a smooth nonconvex minimization, 2018, arXiv:1807.00387.

[28] G. Marchuk, Some applications of splitting-up methods to the solution of mathematical physics problems, Aplikace Matematiky 13 (1968), pp. 103-132.

[29] M. Muehlebach, M.I. Jordan, A dynamical systems perspective on Nesterov acceleration, 2019, arXiv:1905.07436.

[30] Y. Nesterov, A method for solving a convex programming problem with convergence rate O(1/k 2 ) (in Russian), Soviet Math. Dokl. 27 (1983), no. 2, pp. 543-547.

[31] Y. Nesterov, Introductory lectures on convex optimization : A basic course, vol. 87 of Applied Optimization, Kluwer Academic Publishers, Boston, MA, 2004.

[32] N.C. Nguyen, P. Fernandez, R.M. Freund, J. Peraire, Accelerated residual methods for iterative solution of systems of equations, SIAM J. Sci. Computing 40 (2018), pp. A3157-A3179.

[33] M. Laborde, A.M. Oberman, A Lyapunov analysis for accelerated gradient methods : from deterministic to stochastic case, 2019, arXiv:1908.07861.

[34] B.T. Polyak, Some methods of speeding up the convergence of iteration methods, USSR Comput. Math. Math. Phys. (4) 1964, pp. 1-17.

[35] J.M. Sanz-Serna, Symplectic methods. In : Encyclopedia of Applied and Computational Mathematics, B. Engquist ed., Springer, Berlin Heidelberg, 2015, pp. 1451-1458.

[36] B. Shi, S.S. Du, W.J. Su, M. I. Jordan, Acceleration via symplectic discretization of high-resolution differential equations, 2019, arXiv:1902.03694.

[37] B. Shi, S.S. Iyengar, A conservation law method based on optimization. In : Mathematical Theories of Machine Learning – Theory and Applications, Springer, Cham, 2019, pp. 63-85 https://doi.org/10.1007/978-3-030-17076-9_8.

[38] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968), no. 3, pp. 506-517.

[39] W.J. Su, S. Boyd, E.J. Candes, A differential equation for modeling Nesterov’s accelerated gradient method : Theory and insights, J. Machine Learning Res. (17) 2016, pp. 1-43.

[40] T. Sun, P. Yin, D. Li, C. Huang, L. Guan, H. Jiang, Non-ergodic convergence analysis of heavy-ball algorithms, AAAI, 2019.

[41] H.F. Trotter, On the product of semi-groups of operators, Proc. Am. Math. Soc. 10 (1959), no. 4, pp. 545-551.

[42] P. Tseng, On accelerated proximal-gradient methods for convex-concave optimization, Manuscript, May 21 2018, mit.edu/~dimitrib/PTseng/papers/apgm.pdf.

[43] R. de Vogelaere, Methods of integration which preserve the contact transformation property of the Hamiltonian equations, Report no. 4, Dept. Math., Univ. of Notre Dame, Notre Dame, Ind., 1956.

Paper (preprint) in HTML form

The long time behavior and the rate of convergence of symplectic convex algorithms obtained via splitting discretizations of inertial damping systems.

Cristian Daniel Alecsa
Tiberiu Popoviciu Institute of Numerical Analysis
Romanian Academy
Cluj-Napoca, Romania
cristian.alecsa@ictp.acad.ro

Romanian Institute of Science and Technology
Cluj-Napoca, Romania
alecsa@rist.ro Corresponding author. This work was supported by a grant of Ministry of Research and Innovation, CNCS - UEFISCDI, project number PN-III-P1-1.1-TE-2016-0266, within PNCDI III.

Abstract. In this paper we propose new numerical algorithms in the setting of unconstrained optimization problems and we study the rate of convergence in the iterates of the objective function. Furthermore, our algorithms are based upon splitting and symplectic methods and they preserve the energy properties of the inherent continuous dynamical system that contains a Hessian perturbation. In particular, we show that Nesterov gradient method is equivalent to a Lie-Trotter splitting applied to a Hessian driven damping system. Finally, some numerical experiments are presented in order to validate the theoretical results.

Key Words. unconstrained optimization problems, rate of convergence, inertial algorithms, Nesterov convex gradient method, convex function, Hessian driven damping, dynamical system, Lie-Trotter splitting

1 Preliminaries

The theory of accelerated optimization algorithms initiated from the work of Polyak (see [34]) and it peaked with the convex gradient method developed by Nesterov in [30]. The theory of dynamical systems associated to Nesterov accelerated method was introduced first in [39] by Su, Boyd and Candès, where they pointed out that the method they have developed was informal. Also for the case of convex objective functions, the first true model for Nesterov discretization was considered very recently by Muehlebach and Jordan in [29] and in the generalization given in [1] by Alecsa, László and Pinţa, where under the gradient was considered a perturbation of the exact solution of the dynamical system. On the other hand, also recently, in [19] Defazio showed that the accelerated gradient method can be seen as a proximal point method in a Riemannian manifold, but the theoretical analysis was done under the assumption of strong convexity. The first contribution of the present paper is to give a rigorous alternative to these approaches and show that Nesterov method is in fact a Lie-Trotter discretization of a Hessian-driven damping dynamical system (see Theorem 2 and Subsection 2.1). Our second important contribution is to show that many well-known optimization algorithms are in fact a combination between splitting and symplectic methods and we further show that the algorithms based upon Hessian-driven damping discretizations (similar to those in [4]) have a continuous dynamical system that is also constructed with respect to the aforementioned methods (see Subsection 2.2 and Theorem 4). On the other hand, we introduce some generalized optimization algorithms associated to Hessian-driven damping dynamical systems that have many advantages from a practical point of view(see LT-S-IGAHD and Section 4), and this represents our third and last contribution of the present article.

We consider the following unconstrained optimization problem

\displaystyle\min\limits_{x\in\mathbb{R}^{N}}f(x)

(1)

where $f:\mathbb{R}^{N}\to\mathbb{R}$ is a convex and Fréchet differentiable function, for which its gradient $\nabla f:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is $L$ -Lipschitz, i.e. there exists $L>0$ such that for each $x,y\in\mathbb{R}^{N}$

\displaystyle\|\nabla f(x)-\nabla f(y)\|\leq L\|x-y\|\ .

(2)

For simplicity, we consider that all of the second order equations that appear in this paper are coupled with some initial conditions $x(t_{0})=x_{0}\in\mathbb{R}^{N}$ and $\dot{x}(t_{0})=v_{0}\in\mathbb{R}^{N}$ , respectively (where $t_{0}>0$ ). This assumption is imposed only for the clarity of the paper, since (as in [4] and other similar papers), we are interested in the asymptotic behavior of a trajectory of the underlying dynamical system. On the other hand, this is also true for the algorithmic counterparts, where, given $x_{0}\in\mathbb{R}^{N}$ and taking $x_{1}\in\mathbb{R}^{N}$ arbitrary, the accelerated inertial methods start from $n=1$ in order to compute $x_{2}$ .
Now, we give a brief presentation regarding the discrete and the continuous analysis of optimization problems. In [34], Polyak has introduced the heavy ball with friction dynamical system, defined as

\displaystyle\ddot{x}(t)+\gamma\dot{x}(t)=-\nabla f(x(t))\ .

(HBF_γ)

Here $\gamma>0$ is a fixed viscous damping parameter that represents the coefficient of friction of the dynamical system (HBF_γ). Furthermore, the combination between the second derivative of the solution and $\gamma$ is used in order to accelerate the gradient methods using inertial coefficients. For the (HBF_γ) the objective function $f$ does not decrease along the trajectories. On the other hand, the global energy (kinetic + potential) $\mathcal{E}(t)=\frac{1}{2}\|\dot{x}(t)\|^{2}+f(x(t))$ is monotone decreasing, by being a continuous Lyapunov-type functional. Furthermore, we remind that Attouch. et. al. (in [7]) have studied the asymptotic properties of this dynamical system.
The natural discretization of (HBF_γ) is Polyak’s inertial method (see [34]), i.e.

\displaystyle\begin{cases}y_{n}=x_{n}+(1-h\gamma)(x_{n}-x_{n-1})\\ x_{n+1}=y_{n}-h^{2}\nabla f(x_{n})\end{cases}

(PIM_γ)

where, in contrast with descent methods, it has an additional inertial term $x_{n}-x_{n-1}$ . If we divide this difference between two consecutive iteration values by the stepsize $h$ , then we obtain a basic approximation to the first order derivative of the solution $x(t)$ of (HBF_γ). As a side note, we mention that in [40], the authors made the first non-ergodic analysis of the rate of convergence for a variant of Polyak’s algorithm (PIM_γ), that contains a non-constant stepsize, related to the inertial term, under the assumption that the objective function is coercive. Now, we shift our attention to the most popular inertial algorithm, developed by Y. Nesterov in [30]. The following update rule is often called Nesterov’s second accelerated method and it takes the following form :

\displaystyle\begin{cases}y_{n}=x_{n}+\alpha_{n}(x_{n}-x_{n-1})\\ x_{n+1}=y_{n}-h^{2}\nabla f(y_{n})\end{cases}

(AGM2_α)

where, in practical applications, $\alpha_{n}=\dfrac{n}{n+\alpha}$ (see the work of Chambolle and Dossal [17]). A different set of update rules is given by the so-called Nesterov’s first accelerated method (briefly AGM1), which gives the best convergence rate in the objective function, and it is defined by a sequence $(\lambda_{n})_{n\in\mathbb{N}}$ that has the general term defined recursively as $\lambda_{n}=\left(1+\sqrt{1+4\lambda^{2}_{n-1}}\right)/2$ . The equivalence between (AGM1) and (AGM2_α) for $\alpha=2$ has been used in [9] regarding a Lyapunov analysis of Nesterov’s algorithm. Also, in [42] some unified generalizations were made using Bregman distances. In contrast with Polyak’s method (PIM_γ), in [30] it was shown that Nesterov’s algorithm has an improved rate of convergence, namely $O\left(1/n^{2}\right)$ . Further, we point out that, in a nonconvex setting, a two-step inertial method similar to (AGM2_α) has been recently studied in [27]. In [39], Su et. al. introduced the continuous counterpart of Nesterov’s dynamical system of the form

\displaystyle\ddot{x}(t)+\dfrac{\alpha}{t}\dot{x}(t)=-\nabla f(x(t))

(AVD_α)

The dynamical system with asymptotic vanishing damping was introduced in the framework of convex optimization problems. For $\alpha\geq 3$ , it was shown that each trajectory of (AVD_α) satisfies the rate of convergence : $f(x(t))-\min f=O\left(1/t^{2}\right)\text{ as }t\to\infty$ . The relationship between Nesterov’s algorithm (AGM2_α) and the dynamical system (AVD_α) is based upon vanishing stepsize arguments, and so the Nesterov’s optimization scheme can not be considered a proper discretization of the associated differential equation (even the authors of [39] pointed out that their derivation of the dynamical system is informal). A general version of AVD_α can be found in the seminal papers [14], [15] and in their extensions given in [5], where $\dfrac{\alpha}{t}$ is replaced by a general function $a=a(t)$ .
In [13], a perturbed variant of (AVD_α) was analyzed in the setting of nonconvex minimization problems, using the Kurdyka-Lojasiewicz property imposed on a regularization of the objective function. Very recently, in [29], it was shown that Nesterov convex gradient method is the natural discretization of a modified second order dynamical system. Further, a generalized version of the one from [29] was studied in [1], in which the authors proved an $O\left(1/t^{2}\right)$ rate of convergence as $t\to\infty$ for the trajectories of the dynamical system. that contains a perturbation under the gradient (see [1]):

\displaystyle\ddot{x}(t)+\dfrac{\alpha}{t}\dot{x}(t)=-\nabla f\left(x\left(t\right)+\left(\gamma-\dfrac{\alpha\gamma^{2}}{t}\right)\dot{x}(t)\right)

(Extended-AVD)

On the other hand, for this continuous version of Nesterov’s algorithm some theoretical results were derived related to the regularity of the solution of the associated initial value problem.
Now, we shift our attention to generalizations of the aforementioned dynamical systems, in which the derivative of the solution is perturbed with the value of the Hessian along the trajectories. In [4], Attouch et. al. introduced the Hessian-driven damping dynamical system, i.e.

\displaystyle\ddot{x}(t)+\dfrac{\alpha}{t}\dot{x}(t)+\beta(t)\nabla^{2}f(x(t))\dot{x}(t)=-b(t)\nabla f(x(t))\ .

(DIN-AVD

{}_{\tfrac{\alpha}{t},\beta(t),b(t)}

)

Furthermore, the inertial system with constant damping that represents a generalization of the heavy ball with friction was first studied in [3]. Also see [8] for the case of $\left(\text{DIN-AVD}_{\tfrac{\alpha}{t},\beta,1}\right)$ . The idea behind the Hessian dynamics is the fact that the oscillations that are present both in the continuous and in the discrete case are neutralized in the presence of convex objective functions. Now, we present some particular cases that play a key role in the present paper. The Hessian system $\left(\text{DIN-AVD}_{\tfrac{3}{t},0,1}\right)$ was considered in [33], where M. Laborde and A.M. Oberman pointed out that Nesterov’s discretization is in fact a forward Euler discretization for the Hessian dynamical system. In fact, they showed the method is a forward Euler discretization with a predictor value under the gradient, but not a true discretization. So, in this paper, we claim that Nesterov accelerated gradient method is indeed a proper discretization scheme of $\left(\text{DIN-AVD}_{\tfrac{\alpha}{t},\beta,1+\tfrac{\beta}{t}}\right)$ , namely a full Lie-Trotter splitting discretization composed of forward Euler methods. Attouch et. al. proposed a Nesterov-like algorithm but with two consecutive gradients that represent a discretization of the Hessian. This is linked with the term $\nabla^{2}f(x(t))\dot{x}(t)$ that acts like an additional dissipation term. We shall show a link between the true discretization of a modified DIN-AVD ${}_{\tfrac{\alpha}{t},\beta(t),b(t)}$ type system that is a Polyak-type method and the algorithm introduced by Attouch et. al. On the other hand, another particular case is $\left(\text{DIN-AVD}_{\alpha,\beta,1}\right)$ , where in [16], the forward Euler method was constructed as an optimizer to deep learning neural networks using this dynamical system. At the same time, this discretization scheme is a forward Euler applied to a first order dynamical system that is equivalent to the second order differential equation, taking account the regularity of $f$ , i.e.

\displaystyle\begin{cases}\dot{x}(t)+\beta\nabla f(x(t))+\left(\alpha-\dfrac{1}{\beta}\right)x(t)+\dfrac{1}{\beta}v(t)=0\\ \dot{v}(t)+\left(\alpha-\dfrac{1}{\beta}\right)x(t)+\dfrac{1}{\beta}v(t)=0\end{cases}

(3)

We mention that the particular case of the vanishing damping Hessian system $\left(\text{DIN-AVD}_{\tfrac{3}{t},\beta,1+\tfrac{\beta}{t}}\right)$ can also be written as a first order dynamical system. The following form is much simpler than the former first order system. This will be used frequently in our section of main results.

\displaystyle\begin{cases}\dot{x}(t)=\dfrac{2}{t}\left(v\left(t\right)-x\left(t\right)\right)-\beta\nabla f(x(t))\\ \dot{v}(t)=-\dfrac{t}{2}\nabla f(x(t))\end{cases}

(4)

Finally, we end this section by pointing out that, for $\beta=0$ , the dynamical system (4) reduces to system that is equivalent to the second order evolution equation (AVD_α) with $\alpha=3$ . For other details, we refer to [33].

1.1 Some notes on splitting operator techniques

In this subsection we present the most used splitting techniques for a general autonomous evolution equation of the form

\displaystyle\dfrac{d}{dt}\varphi(t)=F(\varphi(t))

(5)

where $\varphi:\mathbb{R}^{N}\to\mathbb{R}$ , $\varphi=\varphi(t)$ is the flow of the dynamical system. If the vector field from the right hand side can be decomposed into two vector fields, i.e. $F=F^{[1]}+F^{[2]}$ , then we attach the associated sub-problems of (5), namely

\displaystyle(P_{1}):\dfrac{d}{dt}\varphi(t)=F^{[1]}(\varphi(t))\quad\text{ and }\quad(P_{2}):\dfrac{d}{dt}\varphi(t)=F^{[2]}(\varphi(t))

(6)

Let $\varphi^{[1]}=\varphi^{[1]}(t)$ be the solution of $(P_{1})$ and $\varphi^{[2]}=\varphi^{[2]}(t)$ be the solution of $(P_{2})$ . Then, the Lie-Trotter splitting or sequential splitting (see [41]) after one stepsize $h$ is :

\displaystyle\psi_{h}=\varphi^{[1]}_{h}\circ\varphi^{[2]}_{h}

(LTS₁)

Another Lie-Trotter splitting operator is to first apply the exact solution of $(P_{1})$ as an initial condition to $(P_{2})$ , i.e.

\displaystyle\psi^{\ast}_{h}=\varphi^{[2]}_{h}\circ\varphi^{[1]}_{h}

(LTS₂)

Here we have employed the usual notation: $\varphi^{[i]}_{t}:=\varphi^{[i]}(t)$ for $i\in\{1,2\}$ . For a discrete time $t_{n}=nh$ , (LTS₁), becomes $\psi(t_{n})=(\varphi^{[1]}\circ\varphi^{[2]})(t_{n})=\left(\varphi^{[1]}_{h}\circ\varphi^{[2]}_{h}\right)^{n}(h)$ and (LTS₂) is equivalent with $\psi^{\ast}(t_{n})=\left(\varphi^{[2]}_{h}\circ\varphi^{[1]}_{h}\right)^{n}(h)$ . A more advanced type of splitting technique is the so-called Strang or Marchuk splitting (see [28] and [38]). From a numerical point of view, this has an advandage over the Lie-Trotter splitting by the fact that is second-order accurate. Again, giving a stepsize $h$ , the Strang splitting at $t=h$ , is defined as

\displaystyle\psi_{h}=\varphi^{[1]}_{h/2}\circ\varphi^{[2]}_{h}\circ\varphi^{[1]}_{h/2}

(SS₁)

Again, the Strang splitting operator can appear into a different symmetric form, namely

\displaystyle\psi^{\ast}_{h}=\varphi^{[2]}_{h/2}\circ\varphi^{[1]}_{h}\circ\varphi^{[2]}_{h/2}

(SS₂)

Applying recursively, we obtain the splitting operator formula for the discrete time $t_{n}=nh$ . Quite interestingly, if, for example, the exact flows $\varphi^{[1]}(t)$ and $\varphi^{[2]}(t)$ are not available, then one can use consistent approximations to these solutions. In this manner, from now on, $\tilde{\varphi}^{[1]}_{h}$ and $\tilde{\varphi}^{[2]}_{h}$ represents numerical discretization of the flows on the two sub-problems $(P_{1})$ and $(P_{2})$ at time $t_{n}=nh$ , respectively. Then, for example the Lie-Trotter splitting (LTS₁) can be defined as

\displaystyle\psi_{h}=\tilde{\varphi}^{[1]}_{h}\circ\tilde{\varphi}^{[2]}_{h}

(7)

This is valid also for the Strang splitting operator formula (see [35]). On the other hand, for the case of multiple vector fields, when $F=F^{[1]}+\ldots+F^{[m]}$ , we consider $m$ sub-problems, i.e.

\displaystyle(P_{1}):\dfrac{d}{dt}\varphi(t)=F^{[1]}(\varphi(t))\quad\ldots\ldots\quad(P_{m}):\dfrac{d}{dt}\varphi(t)=F^{[m]}(\varphi(t))

(8)

In this case, the Lie-Trotter splitting becomes (see [23])

\displaystyle\psi_{h}=\varphi^{[1]}_{h}\circ\ldots\circ\varphi^{[m]}_{h}

(9)

For example, since in the next sections we will use the case of $m=3$ , we observe that the formula (9) is in fact

\displaystyle\psi_{h}=\left(\varphi^{[1]}_{h}\circ\varphi^{[2]}_{h}\right)\circ\varphi^{[3]}_{h}

(10)

This is equivalent to the fact that we divide the evolution equation into two-subproblems and we use a numerical approximation to the solution of the second sub-problem. Furthermore, the first sub-problem will also be divided into another two-subproblems that can also be decomposed according to the Lie-Trotter splitting.
The case of unbounded linear operator has been studied in [26], when $F(\varphi(t))=A\varphi(t)+B\varphi(t)$ . Furthermore, in [25], Hansen et. al. considered the splitting for the case of semilinear evolution equations, when one of the sub-problems is not linear, i.e. $F(\varphi(t))=A\varphi(t)+F^{[2]}(\varphi(t))$ . For both this cases, Lie-Trotter splitting and Strang-Marchuk splitting ca be defined with the above formulas, using the exact flows or some approximations to the exact solution on the sub-problems.
Now, we turn our attention to the case of non-autonomous evolution equations of the form

\displaystyle\dfrac{d}{dt}\varphi(t)=F(t,\varphi(t))

(11)

There are some approaches for this kind of problem. In [10], Batkai et. al. have used a similar approach as Hansen, Kramer and Ostermann. They have considered the case when $F=F^{[1]}+F^{[2]}$ . In their paper, they applied the Lie-Trotter splitting, when $F^{[1]}(\varphi(t))=A(t)\varphi(t)$ and $F^{[2]}(t)=B(t)\varphi(t)$ . For brevity, we recall the sequential splitting (or Lie-Trotter splitting) on the the interval of discretization $(t_{n-1},t_{n}]$ :

\displaystyle\begin{cases}\dfrac{d}{dt}\varphi^{[1]}(t)=B(t)\varphi^{[1]}(t)\quad,\quad t\in(t_{n-1},t_{n}]\quad\text{ with }\varphi^{[1]}(t_{n-1})=\varphi^{LT}(t_{n-1})\\ \dfrac{d}{dt}\varphi^{[2]}(t)=A(t)\varphi^{[2]}(t)\quad,\quad t\in(t_{n-1},t_{n}]\quad\text{ with }\varphi^{[2]}(t_{n-1})=\varphi^{[1]}(t_{n-1})\end{cases}

(12)

where $\varphi^{LT}(t_{n-1})$ represent the previous discretization solution that was applied to $(t_{n-2},t_{n-1}]$ . Also, from the above numerical scheme we get that $\varphi^{LT}(t_{n})=\varphi^{[2]}(t_{n})$ . On the other hand, taking $t=t_{n-1}$ , we obtain the full discretization of Lie-Trotter splitting, given by the formula (LTS₂). In the following sections, we will employ this method, even though, Batkai et. al. have provided another formula when the evolution equation $\dfrac{d}{dt}\varphi(t)=A(r)\varphi(t)+B(r)\varphi(t)$ can be solved for each $r\in\mathbb{R}$ . The only difference is that in (12), we can evaluate $B$ exactly in the current discretization time element $t_{n}$ . Since we consider the updated value of $t_{n}$ much more appropriately for the discretization schemes, we will use the previous formula. Also, we provide a reference for some splitting techniques that can be applied to optimization problems. For this, see the work of Blanes et. al. [12]. We emphasize the fact that, recently, in [2] the authors have used theses splitting techniques in order to devise some optimizers used in neural network training. The analysis is in contrast with the present article by the fact that in [2], the optimization algorithms were developed in an heuristic manner based upon algebraic manipulations. Nevertheless, this shows the applicability of the splitting-type methods to machine learning. Finally, we point out that splitting techniques have been used recently also in [21], where the authors have employed balanced and rebalanced splitting methods in connection with the equilibrium states of some ODE’s associated to proximal algorithms. This differs from our analysis by the fact that we use only the Lie-Trotter discretization scheme.

1.2 A brief review of symplectic methods

We consider the setting of autonomous evolution equation given by the formula (5). Let $\varphi=(x,v):\mathbb{R}^{N}\to\mathbb{R}^{N}$ and consider $H:\mathbb{R}^{N}\to\mathbb{R}$ to be the Hamiltonian. We write the first order differential equation as a Hamiltonian dynamical system, using $F$ from (5) with $F(\varphi(t))=(\nabla_{v}H(x(t),v(t)),-\nabla_{x}H(x(t),v(t)))^{T}$ in the following form (see [22]):

\displaystyle\begin{pmatrix}\dfrac{d}{dt}x(t)\\ \\ \dfrac{d}{dt}v(t)\end{pmatrix}=\begin{pmatrix}\nabla_{v}H(x(t),v(t))\\ \\ -\nabla_{x}H(x(t),v(t))\end{pmatrix}

(HS)

For the pioneering works on differential geometry and basic Hamiltonian systems, we refer to [20] and [43]. Furthermore, for brevity, we recall (see [23] and [24]) the symplectic Euler methods related to (HS)

\displaystyle\begin{cases}x_{n+1}=x_{n}+h\nabla_{v}H(x_{n+1},v_{n})\\ v_{n+1}=v_{n}-h\nabla_{x}H(x_{n+1},v_{n})\end{cases}

(SE₁)

respectively

\displaystyle\begin{cases}x_{n+1}=x_{n}+h\nabla_{v}H(x_{n},v_{n+1})\\ v_{n+1}=v_{n}-h\nabla_{x}H(x_{n},v_{n+1})\end{cases}

(SE₂)

Other symplectic-type algorithms can be constructed with combinations of the form $(x_{n},v_{n+1})$ and $(x_{n+1},v_{n})$ . Also from the work of Hairer, Lubich and Wanner [23], we recall a more advanced type of implicit-explicit symplectic integrators, namely the Stormer-Verlet schemes :

\displaystyle\begin{cases}x_{n+1/2}=x_{n}+\dfrac{h}{2}\nabla_{v}H(x_{n+1/2},v_{n})\\ v_{n+1}=v_{n}-\dfrac{h}{2}\left(\nabla_{x}H(x_{n+1/2},v_{n})+\nabla_{x}H(x_{n+1/2},v_{n+1})\right)\\ x_{n+1}=x_{n+1/2}+\dfrac{h}{2}\nabla_{v}H(x_{n+1/2},v_{n+1})\\ \end{cases}

(SV₁)

and a similar symmetric method

\displaystyle\begin{cases}v_{n+1/2}=v_{n}-\dfrac{h}{2}\nabla_{x}H(x_{n},v_{n+1/2})\\ x_{n+1}=x_{n}+\dfrac{h}{2}\left(\nabla_{v}H(x_{n},v_{n+1/2})+\nabla_{v}H(x_{n+1},v_{n+1/2})\right)\\ v_{n+1}=v_{n+1/2}-\dfrac{h}{2}\nabla_{x}H(x_{n+1},v_{n+1/2})\\ \end{cases}

(SV₂)

Fro brevity, we point out that $x_{n+1/2}$ from (SV₁) and $v_{n+1/2}$ from (SV₂) represent just a notation for intermediate steps that are used in the construction of the algorithms. Interestingly, under general assumptions, for a symplectic integrator of convergence order $\mu$ , one has that $H(x_{n},v_{n})=H(x_{0},v_{0})+O(h^{\mu})$ , where the constant does not depend on the iteration $n$ over exponentially long time intervals, namely $nh\leq e^{h_{0}/2h}$ . This leads to a very good energy behavior of symplectic methods and so these methods represent a good choice for the discretization of a conservative Hamiltonian system.
Also, it is known that if the Hamiltonian is separable, i.e. $H(x,v)=T(v)+U(x)$ , then (HS) becomes

\displaystyle\begin{pmatrix}\dfrac{d}{dt}x(t)\\ \\ \dfrac{d}{dt}v(t)\end{pmatrix}=\begin{pmatrix}\nabla_{v}T(v(t))\\ \\ -\nabla_{x}U(x(t))\end{pmatrix}

(13)

Furthermore, dividing (13) into two sub-problems

\displaystyle(P_{1}):\begin{cases}\dfrac{d}{dt}x(t)=0\\ \\ \dfrac{d}{dt}v(t)=-\nabla_{x}U(x(t))\end{cases}\text{ and }\quad(P_{2}):\begin{cases}\dfrac{d}{dt}x(t)=\nabla_{v}T(v(t))\\ \\ \dfrac{d}{dt}v(t)=0\end{cases}

(14)

and applying the explicit forward Euler discretizations to both the sub-problems, then the symplectic Euler methods (SE₁) and (SE₂) lead to the Lie-Trotter full discretizations (LTS₁) and (LTS₂). The same analysis can be depicted for the Stormer Verlet methods (SV₁) and (SV₂), where these methods are equivalent to the Strang-Marchuk discretization schemes.
Recently, in [22], França et. al. devised some optimization schemes based on some Hamiltonian systems perturbed by a dissipation term $\gamma>0$ , of the form

\displaystyle\begin{pmatrix}\dfrac{d}{dt}x(t)\\ \\ \dfrac{d}{dt}v(t)\end{pmatrix}=\begin{pmatrix}\nabla_{v}H(x(t),v(t))\\ \\ -\nabla_{x}H(x(t),v(t))-\gamma v(t)\end{pmatrix}

(HSD_γ)

where the separable Hamiltonian $H:\mathbb{R}^{N}\to\mathbb{R}$ is defined such as $H(x,v)=\dfrac{\|v\|^{2}}{2m}+f(x)$ , with $m>0$ related to the inertial positive coefficient $\gamma$ . Here, the first term $T(v)=\dfrac{\|v\|^{2}}{2m}$ is the kinetic energy and $U(x)=f(x)$ is the potential energy of the system. Furthermore, the heavy ball with friction system (HBF_γ) is a particular case of (HSD_γ) and can be written as

\displaystyle\begin{pmatrix}\dfrac{d}{dt}x(t)\\ \\ \dfrac{d}{dt}v(t)\end{pmatrix}=\begin{pmatrix}v(t)\\ \\ -\nabla f(x(t))-\gamma v(t)\end{pmatrix}

(15)

and the Hamiltonian of the problem is a Lyapunov function and satisfies $\dfrac{d}{dt}H(x(t),v(t))=-\gamma\|v(t)\|^{2}\leq 0$ along the trajectories. Furthermore, it was shown that Polyak’s method (PIM_γ) is a conformal symplectic Euler method for (15), but Nesterov’s method is not. Also, Polyak’s algorithm can be shown to be equivalent to a Lie-Trotter splitting (in connection with the decomposition of HSD_γ), where for the conservative subproblem one applies a symplectic Euler method and for the dissipative subproblem an explicit forward Euler for the exact solution is used (for details see [22]). Finally, we remind that for heavy ball system (HBF_γ) without friction (when $\gamma=0$ ), the symplectic Euler and the symplectic Stormer-Verlet algorithms can be considered as the natural discretizations of the underlying ODE (see [23] and [37]).

1.3 The organization of the paper

The outline of the present paper is the following. In Subsection 2.1, we show that Nesterov’s algorithm can be viewed as a Lie-Trotter or sequential splitting for a Hessian driven damping system. Furthermore, we also present an alternative version of writing the momentum based methods with asymptotic vanishing damping that resembles the continuous dynamical system. Also, we present some well known results for the heavy ball system and the extended dynamical system with asymptotic vanishing damping in connection to symplectic Euler integrator. In Subsection 2.3 and Subsection 2.4, we consider some results regarding various algorithms linked to symplectic integrators for the extended AVD system. In Section 3, we prove the rate of convergence of the Hessian-based algorithms in the framework of convex , differentiable objective functions, such that the gradient satisfies the Lipschitz continuity property. Finally, in Section 4, we validate our theoretical results with some numerical experiments involving various choices concerning the parameters of our newly introduced algorithms in the framework of classical convex test functions. For conciseness, we point out a brief analysis that will be done in the following sequel.

In Subsection 2.2, we show some connections between the continuous version and the discrete version of some dynamical systems. Then dynamical system (Extended-AVD) can be seen as a generalization of (AVD_α) with an extra perturbation under the gradient that is similar to the updates in Nesterov’s algorithm iterations. On the other hand, the (AVD_α) dynamical system can be generalized through the Hessian driven damping evolution equation (DIN-AVD ${}_{\tfrac{\alpha}{t},\beta(t),b(t)}$ ), by using a perturbation term consisting on the product of the Hessian matrix and the gradient of the underlying solution. For conciseness, we remind that the difference between Polyak’s method (PIM_γ) and (AGM2_α) lies on the fact that the momentum method is a splitting method for a system similar to (AVD_α) but with constant damping, while Nesterov’s method is also a splitting based-discretization combined with a symplectic Euler method, but for the extended continuous system. Quite suprinsingly, our novelty was given in the second section concerning the fact that Nesterov’s method can be seen also a discretization (by using a triple splitting technique) but for the Hessian driven damping system (DIN-AVD ${}_{\tfrac{\alpha}{t},\beta(t),b(t)}$ ) with $\beta(t)=\beta$ and $b(t)=1+\beta/t$ . At the same time, regarding the Hessian-based dynamical system, as in the Nesterov’s case, IGAHD algorithm (see [4]), is also a splitting-based method where on the potential forces we apply the symplectic Euler method. The key role in our analysis is the fact that the symplectic methods are playing a special role in the discretizations of the continuous dynamical systems. We emphasize the fact that we have four types of symplectic methods, namely (SE₁), (SE₂), (SV₁) and (SV₂). By a trial and error approach, these methods, in combination with splitting discretizations, applied on the Hessian driven damping system (DIN-AVD ${}_{\tfrac{\alpha}{t},\beta(t),b(t)}$ ) lead to unfeasible algorithms due to the fact that the continuous evolution equation is of rather complicated nature. Also, using these techniques, one obtains algorithms containing at least with three gradient evaluations or with a product on the form $\nabla^{2}f(x(t))\cdot\nabla f(x(t))$ that is a computational burden (even after a forward-type discretization). So, our aim Subsection 2.3 and Subsection 2.4, is to focus only on the (Extended-AVD) dynamical system. We also have observed that the alternative variant of the symplectic Euler in contrast with the one used in [29] leads to the same empirical convergence as Nesterov’s algorithm but with an additional gradient evaluation. On the other hand, Stormer-Verlet algorithms leads to three gradient evaluations, which is unfeasible. In this way, we will analyze in an intuitive manner only the Stormer-Verlet integrator (SV₂) and the symplectic Euler integrator (SE₁) that induces two gradient evaluations. This is similar to the two gradient scheme IGAHD. Moreover, we will introduce another dynamical system that for a Lie-Trotter discretization we will obtain the perturbed Nesterov-type algorithm introduced in [32].

2 The discretizations of the dynamical systems

2.1 Nesterov’s accelerated method is a Lie-Trotter splitting method applied to a Hessian-based dynamical system

In [33], it was shown that the system $\left(\text{DIN-AVD}_{\tfrac{3}{t},\beta,1+\tfrac{\beta}{t}}\right)$ with the vanishing damping $3/t$ can also be written as a first order dynamical system, namely (4). For brevity, in our first result, we will show that the general vanishing damping system with Hessian perturbation $\left(\text{DIN-AVD}_{\tfrac{3}{t},\beta,1+\tfrac{\beta}{t}}\right)$ can be written as a first order system. This result will be used also in this section when we will prove that Nesterov’s method is actually a Lie-Trotter splitting applied to three sub-problems of this dynamical system.

Proposition 1.

Let $f\in C^{2}(\mathbb{R}^{n},\mathbb{R})$ and let $\beta\in\mathbb{R}$ and $\alpha>1$ . Further, consider the $\left(\text{DIN-AVD}_{\tfrac{\alpha}{t},\beta,1+\tfrac{\beta}{t}}\right)$ second order evolution equation, i.e.

\displaystyle\ddot{x}(t)+\dfrac{\alpha}{t}\dot{x}(t)+\beta\nabla^{2}f(x(t))\dot{x}(t)=-\left(1+\dfrac{\beta}{t}\right)\nabla f(x(t))

(16)

On the other hand, for $t>0$ consider the non-autonomous first order dynamical system :

\displaystyle\begin{cases}\dot{x}(t)=\dfrac{\alpha-1}{t}(v(t)-x(t))-\beta\nabla f(x(t))\\ \dot{v}(t)=-\dfrac{t}{\alpha-1}\nabla f(x(t))\end{cases}

(17)

Then, the dynamical systems (16) and (17) are equivalent, in the sense that each solution of the second order dynamical system is a solution of the first order system and viceversa.

Proof .

Let $(x(t),v(t))\in\mathbb{R}^{n}\times\mathbb{R}^{N}$ be a solution for (17). We will show that $x=x(t)$ also satisfies the Hessian driven damping system (16). Firstly, using the fact that

v(t)=x(t)+\dfrac{t}{\alpha-1}\left(\dot{x}(t)+\beta\nabla f(x(t))\right),

by differentiation we obtain that

\dot{v}(t)=\dot{x}(t)+\dfrac{1}{\alpha-1}\left(\dot{x}(t)+\beta\nabla f(x(t))\right)+\dfrac{t}{\alpha-1}\left(\ddot{x}(t)+\beta\nabla^{2}f(x(t))\dot{x}(t)\right).

Since $\dot{v}(t)=-\dfrac{t}{\alpha-1}\nabla f(x(t))$ , one obtains that

-\dfrac{t}{\alpha-1}\nabla f(x(t))=\dot{x}(t)+\dfrac{1}{\alpha-1}\dot{x}(t)+\dfrac{\beta}{\alpha-1}\nabla f(x(t))+\dfrac{t}{\alpha-1}\ddot{x}(t)+\dfrac{\beta t}{\alpha-1}\nabla^{2}f(x(t))\dot{x}(t)\ ,

which is equivalent to (16).
On the other hand, let $x=x(t)$ be a solution to the Hessian damping system (16). We will show that, if we define $v=v(t)$ as in the first equation of (17), then the derivative of the function $v$ satisfies the second equation of (17). In this sense, by the definition of $v=v(t)$ , we have that

v(t)=x(t)+\dfrac{t}{\alpha-1}\left(\dot{x}(t)+\beta\nabla f(x(t))\right)\ .

By differentiation, as before, we obtain that

\dot{v}(t)=\dot{x}(t)+\dfrac{1}{\alpha-1}\left(\dot{x}(t)+\beta\nabla f(x(t))\right)+\dfrac{t}{\alpha-1}\left(\ddot{x}(t)+\beta\nabla^{2}f(x(t))\dot{x}(t)\right).

Using the fact that $x=x(t)$ satisfies (16), it follows that

\dot{v}(t)=\dot{x}(t)+\dfrac{1}{\alpha-1}\dot{x}(t)+\dfrac{\beta}{\alpha-1}\nabla f(x(t))+\dfrac{t}{\alpha-1}\left(-\dfrac{\alpha}{t}\dot{x}(t)-\left(1+\dfrac{\beta}{t}\right)\nabla f(x(t))\right)\ .

Finally, this equivalent exactly with the second line of (17). $\blacksquare$

Now, the following theorem is the main result of this subsection. In [33], the authors showed the equivalence between (AGM2_α) and (NaG) for the case when $\alpha=3$ . Moreover, we emphasize that their method closely resembles the forward Euler method applied to the dynamical system (17). This intuitive idea was pointed out in the authors work. On the other hand, since under the gradient we use a predictor $y_{n}$ and not the true value $x_{n}$ at the discretization time $t_{n}$ , this method can not be interpreted as a proper discretization to the first order dynamical system. Following this idea, we will show that Nesterov’s method (NaG) can be viewed as a discrete counterpart of (17) using the idea of Lie-Trotter splitting. In fact, we show that Nesterov’s accelerated gradient method is a Lie-Trotter splitting applied to a decomposition consisting of three subproblems at time $t_{n}=nh$ .

Theorem 2.

Nesterov’s accelerated gradient method (AGM2_α) is the numerical discretization of Lie-Trotter splitting applied to two non-autonomous subproblems and one autonomous subproblem at the discretization time $t_{n}=nh$ . Furthermore, for each continuous subproblem, an explicit forward Euler method is applied to each of them.

Proof .

(I.) The first step of the proof is inspired by the work of A.M. Oberman and M. Laborde [33]. We will present an alternative form of the Nesterov’s algorithm (AGM2_α) with general vanishing damping with coefficient $\alpha>1$ , that, we will observe, it is very similar to the continuous counterpart (17). For this, let $\alpha>1$ and consider the Nesterov’s discretization scheme (AGM2_α). At the same time, consider the following alternate form of the Nesterov’s accelerated gradient method :

\displaystyle\begin{cases}x_{n+1}=x_{n}+\dfrac{(\alpha-1)h}{t_{n}}(v_{n}-x_{n})-\tilde{\beta}\nabla f(y_{n})\\ v_{n+1}=v_{n}-\dfrac{ht_{n}}{\alpha-1}\nabla f(y_{n})\end{cases}

(NaG)

Moreover, let $\tilde{\beta}=h^{2}$ and $y_{n}=x_{n}+\dfrac{h(\alpha-1)}{t_{n}}(v_{n}-x_{n})$ . Then, we show that the algorithms (AGM2_α) and (NaG) are equivalent. Furthermore, $t_{n}=h(n+\alpha)$ is equivalent to $\alpha_{n}=\dfrac{n}{n+\alpha}$ and $t_{n}=hn$ is equivalent to the choice $\alpha_{n}=\dfrac{n-\alpha}{n}$ . We begin by rewritting the discretization scheme (NaG) as follows

\displaystyle\begin{cases}x_{n+1}-x_{n}-\dfrac{(\alpha-1)h}{t_{n}}v_{n}+\dfrac{(\alpha-1)h}{t_{n}}x_{n}=-h^{2}\nabla f(y_{n})\\ \\ v_{n+1}-v_{n}=-\dfrac{ht_{n}}{\alpha-1}\nabla f(y_{n})\end{cases}

Multiplying the second equation by $\dfrac{-h(\alpha-1)}{t_{n}}$ and adding it to the first one, we obtain that :

v_{n+1}=x_{n}+\dfrac{t_{n}}{h(\alpha-1)}(x_{n+1}-x_{n})\ .

Now, using the fact that $y_{n}=x_{n}+\dfrac{h(\alpha-1)}{t_{n}}v_{n}-\dfrac{h(\alpha-1)}{t_{n}}x_{n}$ and that $v_{n}=x_{n-1}+\dfrac{t_{n-1}}{h(\alpha-1)}(x_{n}-x_{n-1})$ , it follows that

y_{n}=x_{n}+(x_{n}-x_{n-1})\left(\dfrac{t_{n-1}}{t_{n}}-\dfrac{h(\alpha-1)}{t_{n}}\right)\ .

In this way, we obtain Nesterov’s algorithm (AGM2_α), with the inertial coefficient $\alpha_{n}=\dfrac{1}{t_{n}}\left(t_{n-1}-h(\alpha-1)\right)$ . Finally, if we take the discretization time to be $t_{n}=nh$ , where $h$ is the stepsize, then $\alpha_{n}=\dfrac{n-\alpha}{n}$ . On the other hand, the choice $t_{n}=h(n+\alpha)$ leads to $\alpha_{n}=\dfrac{n}{n+\alpha}$ . These two represent, in general, the classical choices for the inertial sequence $(\alpha_{n})_{n\in\mathbb{N}}$ .
(II.) Now, the second step is to consider splitting the first order dynamical system (17) in the following subproblems:

\displaystyle(P_{1}):\begin{cases}\dot{x}(t)=\dfrac{\alpha-1}{t}(v(t)-x(t))\\ \dot{v}(t)=-\dfrac{t}{\alpha-1}\nabla f(x(t))\end{cases}\quad(P_{2}):\begin{cases}\dot{x}(t)=-\beta\nabla f(x(t))\\ \dot{v(t)}=0\end{cases}

Using the notation from the first section, the splitting discrete solution is $\psi(t_{n})=\tilde{\phi}^{[2]}(t_{n})\circ\tilde{\phi}^{[1]}(t_{n})$ . We consider $(\tilde{x}_{n+1},\tilde{v}_{n+1})^{T}$ to be the numerical solution for the first subproblem $(P_{1})$ at time $t_{n}=nh$ . Then, applying a forward Euler discretization on $(P_{2})$ and taking $\beta=h$ , one obtains that

\displaystyle\begin{cases}x_{n+1}=\tilde{x}_{n+1}-h^{2}\nabla f(\tilde{x}_{n+1})\\ v_{n+1}=\tilde{v}_{n+1}\end{cases}

where $(x_{n+1},v_{n+1})^{T}$ is the current splitting solution at time $t_{n}=nh$ . Now, also using forward Euler discretization we find the numerical solution $(\tilde{x}_{n+1},v_{n+1})^{T}$ at time $t_{n}=nh$ for the continuous subproblem $(P_{1})$ . As before, we employ splitting the subproblem $(P_{1})$ into two non-autonomous subproblems, the first one containing constant velocity $v=v(t_{n})$ and the second one preserving the value of $x=x(t_{n})$ . So, let

\displaystyle(P_{1})_{1}:\begin{cases}\dot{x}(t)=\dfrac{\alpha-1}{t}(v(t)-x(t))\\ \dot{v}(t)=0\end{cases}\quad(P_{1})_{2}:\begin{cases}\dot{x}(t)=0\\ \dot{v}(t)=\dfrac{-t}{\alpha-1}\nabla f(x(t))\end{cases}

Then, the numerical Lie-Trotter splitting on $(P_{1})$ at the discrete time $t_{n}=nh$ is defined as $\tilde{\phi}^{[1]}(t_{n})=(\tilde{\phi}^{[12]}\circ\tilde{\phi}^{[11]})(t_{n})$ , where $\tilde{\phi}^{[1i]}$ is the forward Euler operator that approximates the flow of the subproblem $(P_{1})_{i}$ . If we apply the forward Euler numerical scheme on both the new subproblems, we obtain the following :

\displaystyle\begin{cases}\tilde{x}_{n+1/2}=x_{n}+\dfrac{h(\alpha-1)}{t_{n}}(v_{n}-x_{n})\\ \tilde{v}_{n+1/2}=v_{n}\end{cases}\text{ and }\begin{cases}\tilde{x}_{n+1}=\tilde{x}_{n+1/2}\\ \tilde{v}_{n+1}=\tilde{v}_{n+1/2}-\dfrac{ht_{n}}{\alpha-1}\nabla f(\tilde{x}_{n+1/2})\end{cases}

Then, it follows that the Lie-Trotter approximation to $(P_{1})$ can be written as follows

\displaystyle\begin{cases}\tilde{x}_{n+1}=x_{n}+\dfrac{h(\alpha-1)}{t_{n}}(v_{n}-x_{n})\\ \tilde{v}_{n+1}=v_{n}-\dfrac{ht_{n}}{\alpha-1}\nabla f\left(x_{n}+\dfrac{h(\alpha-1)}{t_{n}}\left(v_{n}-x_{n}\right)\right)\end{cases}

where $(x_{n},v_{n})^{T}$ is the approximation of the solution of the continuous problem at time $t_{n-1}=(n-1)h$ and since we are employing a sequential-type splitting, it is the numerical solution at time $t_{n-1}=(n-1)h$ of $(P_{2})$ . In this way, putting all together, we obtain that

\displaystyle\begin{cases}x_{n+1}=x_{n}+\dfrac{h(\alpha-1)}{t_{n}}(v_{n}-x_{n})-h^{2}\nabla f\left(x_{n}+\dfrac{h(\alpha-1)}{t_{n}}\left(v_{n}-x_{n}\right)\right)\\ v_{n+1}=v_{n}-\dfrac{ht_{n}}{\alpha-1}\nabla f\left(x_{n}+\dfrac{h(\alpha-1)}{t_{n}}\left(v_{n}-x_{n}\right)\right)\end{cases}

By taking $y_{n}=x_{n}+\dfrac{h(\alpha-1)}{t_{n}}(v_{n}-x_{n})$ and using the first step of the proof, we observe that the conclusion follows easily. $\blacksquare$

Finally, we give some observations regarding the Lie-Trotter splitting that was employed in the proof presented before.

Remark 3.

The numerical solution at time $t_{n}=nh$ can be written as $\psi(t_{n})=\left(\tilde{\phi}^{[2]}\circ\left(\tilde{\phi}^{[12]}\circ\tilde{\phi}^{[11]}\right)\right)(t_{n})$ , where the three operators are the approximation to the flows of the three subproblems considered above. Using the ideas that we have presented, we observe that all of the approximations are made at the discrete time $t_{n}=nh$ . Consequently, in contrast to the Strang-Marchuk splitting, we have used the stepsize $h$ and not $h/2$ in the numerical solutions of the subproblems. Finally, we observe that the splitting technique has the advantage that, at a fixed discrete time $t_{n}$ , the same stepsize and the same value of $t_{n}$ are used in all of the discretized subproblems.
Finally, we point out that, in the present proof, we have considered the following technique: we have taken the dynamical system (17) and then we applied the Lie-Trotter splitting with the choice $\beta=h$ . Then, we have showed the equivalence with NaG. But, we remind that in NaG we have chosen $\tilde{\beta}=h^{2}$ and this is valid since this choice does not correspond to a dynamical system, by the fact that it is a parameter from the discretization scheme.

2.2 The role of symplectic methods in optimization algorithms

In this subsection we will make a brief review concerning the role of the sympletic integrators in the setting of unconstrained optimization problems. Furthermore, we point out that many classical inertial algorithms that are based upon dynamical systems with vanishing or with constant damping can be interpreted as a combination between the Lie-Trotter splitting and choices of symplectic methods. Firstly, let’s consider the heavy ball with friction model (HBF_γ). Also, consider its natural discretization, namely the Polyak’s algorithm given by (PIM_γ). In [22], it was pointed out that the momentum method can be interpreted as a splitting discretization between conservative and dissipative vector fields. Indeed, consider the subproblems associated to (HBF_γ) as follows

\displaystyle(P_{1}):\begin{cases}\dot{x}(t)=0\\ \dot{v}(t)=-\gamma v\end{cases}\text{ and }(P_{2}):\begin{cases}\dot{x}(t)=v\\ \dot{v}(t)=-\nabla f(x(t))\end{cases}

where the Hamiltonian $H(x,v)=\dfrac{\|v\|^{2}}{2}+f(x)$ is the energy of the entire dynamical system. The first subproblem $(P_{1})$ represents the dissipative system and the second one, i.e. $(P_{2})$ is the conservative system. Then, if one applies the explicit forward Euler method on $(P_{1})$ and the symplectic Euler method (SE₂) on $(P_{2})$ , it follows that the Lie-Trotter discretization scheme at time $t_{n}=nh$ , i.e. $\psi(t_{n})=\tilde{\phi}^{[2]}(t_{n})\circ\tilde{\phi}^{[1]}(t_{n})$ becomes

\displaystyle\begin{cases}y_{n}=x_{n}+h(1-h\gamma)v_{n}\\ x_{n+1}=y_{n}-h^{2}\nabla f(x_{n})\end{cases}

which is an alternative form of writting the momentum method (PIM_γ), using the idea that, from the symplectic discretization of the second subproblem, $v_{n+1}=\dfrac{x_{n+1}-x_{n}}{h}$ is the discretization of the velocity of the underlying system.
Now, we turn our attention to the dynamical system that models Nesterov’s accelerated gradient method (AGM2_α). This dynamical system is the true continuous counterpart of Nesterov’s algorithm, since under the gradient it contains a continuous version of the inertial term $y_{n}$ , i.e. this inertial term is a discretization of a combination between the exact flow and its first order derivative. For brevity, we recall the idea behind [29], where the authors considered Nesterov’s method as a combination of a Lie-Trotter splitting and a symplectic Euler method. For this, we consider a more generalized vanishing damping system Extended-AVD.
It is obvious that taking $\gamma=h$ , one obtains the Nesterov’s algorithm (AGM2_α). Furthermore, (Extended-AVD) can be written as a first order evolution equation, using the idea of velocity, as we have pointed out in the case of the heavy ball system :

\displaystyle\begin{cases}\dot{x}(t)=v(t)\\ \dot{v}(t)=-\dfrac{\alpha}{t}v(t)-\nabla f\left(x\left(t\right)+\left(\gamma-\dfrac{\alpha\gamma^{2}}{t}\right)v(t)\right)\end{cases}

(18)

By the fact that the above system has no inherent energetic structure, in [29], the term $\nabla f(x(t))$ was added and substracted, such that one can use, as in the case of the heavy ball with friction system, the dissipative and conservative parts of the system. Thus, the following subproblems were considered :

\displaystyle(P_{1}):\begin{cases}\dot{x}(t)=0\\ \dot{v}(t)=-\dfrac{\alpha}{t}v+\nabla f(x(t))-\nabla f\left(x\left(t\right)+\left(\gamma-\dfrac{\alpha\gamma^{2}}{t}\right)v(t)\right)\end{cases}\text{ and }(P_{2}):\begin{cases}\dot{x}(t)=v\\ \dot{v}(t)=-\nabla f(x(t))\end{cases}

Here, the first subproblem $(P_{1})$ represents the non-potential forces and the second subproblem $(P_{2})$ stands for the potential forces. Moreover, as in the case of the heavy ball with friction system, applying a forward Euler method to the first subproblem and the symplectic Euler integrator (SE₂) to the second one leads to the Nesterov’s method (AGM2_α).
Now, we shift our focus to the Hessian-driven damping dynamical system $\left(\text{DIN-AVD}_{\tfrac{\alpha}{t},\beta,1+\tfrac{\beta}{t}}\right)$ . First of all, we will present a natural discretization of this dynamical system that resembles a Polyak-type counterpart to the IGAHD algorithm from [4]. Second of all, we show that a discretization resembling the IGAHD algorithm can be obtained, as in the case of Nesterov’s method, by a combination of a Lie-Trotter splitting method with a symplectic Euler integrator, but for a Hessian-based system endowed with a perturbation under the gradient as in (Extended-AVD). We consider the $\left(\text{DIN-AVD}_{\tfrac{\alpha}{t},\beta,1+\tfrac{\beta}{t}}\right)$ . The natural discretization this system is a Polyak-type method, where at time $t_{n}=nh$ under the gradient appears the old iteration $x_{n}$ . We approximate $\nabla^{2}f(x(t))\dot{x}(t)$ by the difference of consecutive gradients as in the work of Attouch et. al. [4] and we consider the discretization of the continuous system, as follows :

	$\displaystyle\dfrac{x_{n+1}-2x_{n}+x_{n-1}}{h^{2}}+\dfrac{\alpha}{t_{n}}\dfrac{x_{n}-x_{n-1}}{h}+\dfrac{\beta}{h}\left[\nabla f(x_{n})-\nabla f(x_{n-1})\right]=-\nabla f(x_{n})-\dfrac{\beta}{t_{n}}\nabla f(x_{n})$
	$\displaystyle\dfrac{x_{n+1}-2x_{n}+x_{n-1}}{h^{2}}+\dfrac{\alpha}{nh}\dfrac{x_{n}-x_{n-1}}{h}+\dfrac{\beta}{h}\left[\nabla f(x_{n})-\nabla f(x_{n-1})\right]=-\nabla f(x_{n})-\dfrac{\beta}{nh}\nabla f(x_{n})$

By some easy computations, it follows that

x_{n+1}=x_{n}+\left(1-\dfrac{\alpha}{n}\right)(x_{n}-x_{n-1})-\beta h\left[\nabla f(x_{n})-\nabla f(x_{n-1})\right]-\dfrac{\beta h}{n}\nabla f(x_{n})-h^{2}\nabla f(x_{n})

Define $\alpha:=1-\dfrac{\alpha}{n}$ . Then, we obtain that

\displaystyle\begin{cases}&y_{n}=x_{n}+\alpha_{n}(x_{n}-x_{n-1})-\beta h\left[\nabla f(x_{n})-\nabla f(x_{n-1})\right]-\dfrac{\beta h}{n}\nabla f(x_{n})\\ &x_{n+1}=y_{n}-h^{2}\nabla f(x_{n})\end{cases}

(Polyak-IGAHD)

Now, our only result of this subsection consists in showing that a similar algorithm to the IGAHD (see [4]) with a predictor $y_{n}$ under the gradient can be obtained as in the case of (Extended-AVD) by applying a sequential splitting method with a symplectic Euler integrator on the second subproblem. Furthermore, the only difference between our algorithm and IGAHD is that the latter one contains $\nabla f(x_{n-1})$ instead $\nabla f(x_{n})$ . The choice $\nabla f(x_{n-1})$ seems unnatural, by the fact that at the discrete time $t_{n}=nh$ , the gradient of the flow is approximated by $\nabla f(x_{n-1})$ and not with the current gradient value. Nevertheless, by using the same technique of approximation, we can easily obtain a similar discretization scheme with the IGAHD algorithm.

Theorem 4.

Let $f\in C^{2}(\mathbb{R}^{N},\mathbb{R})$ and consider the following Hessian driven damping second order system :

\displaystyle\hskip-14.22636pt\ddot{x}(t)+\dfrac{\alpha}{t}\dot{x}(t)+\beta\nabla^{2}f(x(t))\dot{x}(t)=-\dfrac{\beta}{t}\nabla f(x(t))-\nabla f\left(x(t)+\gamma\left(1-\dfrac{\alpha\gamma}{t}\right)\dot{x}(t)-\beta\gamma^{2}\nabla^{2}f(x(t))\dot{x}(t)-\dfrac{\beta\gamma^{2}}{t}\nabla f(x(t))\right)

Then, a Lie-Trotter splitting discretization, with a forward Euler method on the first subproblem and a symplectic integrator on the second one leads to an IGAHD-type algorithm.

Proof .

As in the case of (Extended-AVD) we consider two subproblems, one containing the so-called non-potential forces and the other one containing the potential forces :

\displaystyle(P_{1}):\begin{cases}\dot{x}=0\\ \dot{v}=-\dfrac{\alpha}{t}v(t)-\beta\nabla^{2}f(x(t))v(t)-\dfrac{\beta}{t}\nabla f(x(t))-\nabla f\left(\mathcal{A}(t,x(t),v(t))\right)+\nabla f(x(t))\end{cases}

where

\mathcal{A}(t,x(t),v(t))=x(t)+\gamma\left(1-\dfrac{\alpha\gamma}{t}\right)v(t)-\beta\gamma^{2}\nabla f(x(t))v(t)-\dfrac{\beta\gamma^{2}}{t}\nabla f(x(t))

and the second subproblem as

\displaystyle(P_{2}):\begin{cases}\dot{x}=v\\ \dot{v}=-\nabla f(x(t))\end{cases}

Taking $\gamma=h$ and defining $\alpha_{n}:=1-\dfrac{\alpha h}{t_{n}}$ we obtain that

\displaystyle\begin{cases}x_{n+1/2}=x_{n}\\ v_{n+1/2}=\alpha_{n}v_{n}-\beta h\nabla^{2}f(x_{n})v_{n}-\dfrac{\beta h}{t_{n}}\nabla f(x_{n})+h\nabla f(x_{n})-h\nabla f\left(x_{n}+h\alpha_{n}v_{n}-\beta h^{2}\nabla^{2}f(x_{n})v_{n}-\dfrac{\beta h^{2}}{t_{n}}\nabla f(x_{n})\right)\end{cases}

Furthermore, for the second subproblem $(P_{2})$ , we apply the symplectic Euler integrator (SE₂) and obtain that

\displaystyle\begin{cases}x_{n+1}=x_{n+1/2}+hv_{n+1}\\ v_{n+1}=v_{n+1/2}-h\nabla f(x_{n+1/2})\end{cases}

Using the fact that for $t_{n}=nh$ , the splitting solution is the discrete solution of $(P_{2})$ , namely $(x_{n+1},v_{n+1})^{T}\in\mathbb{R}^{N}\times\mathbb{R}^{N}$ , we obtain from the first equation given above, that $v_{n}=\dfrac{x_{n}-x_{n-1}}{h}$ . By using this and the approximation involving the Hessian, i.e. $\nabla^{2}f(x_{n})v_{n}$ at time $t_{n}=nh$ by two consecutive gradients, it follows that

\displaystyle\begin{cases}&y_{n}=x_{n}+\alpha_{n}(x_{n}-x_{n-1})-\beta h\left[\nabla f(x_{n})-\nabla f(x_{n-1})\right]-\dfrac{\beta h}{n}\nabla f(x_{n})\\ &x_{n+1}=y_{n}-h^{2}\nabla f(y_{n})\end{cases}

(IGAHD-type)

$\blacksquare$

We end this section by emphasizing the role of choosing the appropriate first order evolution equation for the Hessian driven damping system in the splitting methods.

Remark 5.

1.) At the second order dynamical system from Theorem 4, under the gradient at $\mathcal{A}(t,x(t),\dot{x}(t))$ , we can approximate $\dfrac{\beta\gamma^{2}}{t}\nabla f(x(t))$ as $\dfrac{\beta h^{2}}{t_{n}}\nabla f(x_{n-1})$ . In this way, we obtain the classical IGAHD algorithm from [4].
2.) The idea behind Theorem 4 is the fact that we separate the perturbed gradient $-\dfrac{\beta}{t}\nabla f(x(t))$ from the unperturbed gradient value. Moreover, the latter one is similar to the (Extended-AVD) since, under the gradient, the value contains a combination of both the exact flow $x=x(t)$ and its first order derivative $\dot{x}=\dot{x}(t)$ .
3.) In [36], Shi et. al. have studied a symplectic Euler method that was applied to a full Hessian driven damping system similar to DIN-AVD ${}_{\tfrac{\alpha}{t},\beta(t),b(t)}$ . The key difference between our analysis and theirs is that our algorithms are also based upon a Lie-Trotter discretization in combination with symplectic methods.

2.3 Discussion on some splitting techniques for the extended AVD dynamical system

Until now, we have used the symplectic Euler method (SE₂) but, in the present subsection, we consider a Lie-Trotter splitting discretization scheme by combining a forward Euler method on the first subproblem and the symplectic Euler integrator (SE₁) on the second one, on the dynamical system (Extended-AVD). In this way we consider the following two subproblems:

\displaystyle(P_{1}):\begin{cases}\dot{x}(t)=0\\ \dot{v}(t)=-\dfrac{\alpha}{t}v+\nabla f(x(t))-\nabla f\left(x\left(t\right)+\left(\gamma-\dfrac{\alpha\gamma^{2}}{t}\right)v(t)\right)\end{cases}\text{ and }(P_{2}):\begin{cases}\dot{x}(t)=v\\ \dot{v}(t)=-\nabla f(x(t))\end{cases}

By using a forward Euler discretization on $(P_{1})$ and setting $\gamma=h$ , one obtains that

\displaystyle\begin{cases}&x_{n+1/2}=x_{n}\\ &v_{n+1/2}=\left(1-\dfrac{\alpha h}{t_{n}}\right)v_{n}-h\nabla f\left(x_{n}+h\left(1-\dfrac{\alpha h}{t_{n}}\right)v_{n}\right)+h\nabla f(x_{n})\end{cases}

Also, applying the symplectic Euler method (SE₁) on the subproblem $(P_{2})$ , it follows that

\displaystyle\begin{cases}&x_{n+1}=x_{n+1/2}+hv_{n+1/2}\\ &v_{n+1}=v_{n+1/2}-h\nabla f(x_{n+1})\end{cases}

Now, as before, using the discretization time $t_{n}=nh$ and denoting $\alpha_{n}:=1-\dfrac{\alpha}{n}$ , we introduce the auxiliary term $y_{n}=x_{n}+h\alpha_{n}v_{n}$ . So, we get that

\displaystyle\begin{cases}&x_{n+1}=x_{n}+h\alpha_{n}v_{n}-h^{2}\nabla f(y_{n})+h^{2}\nabla f(x_{n})\\ &v_{n+1}=\alpha_{n}v_{n}-h\nabla f(y_{n})+h\nabla f(x_{n})-h\nabla f(x_{n+1})\end{cases}

By multiplying the second equation by $h$ and substracting them, we obtain the perturbed velocity which differs by the one obtained in the case of (SE₂) by an extra gradient, namely:

v_{n}=\dfrac{x_{n}-x_{n-1}}{h}-h\nabla f(x_{n})\,.

Finally, we obtain our first algorithm, i.e.

\displaystyle\begin{cases}&y_{n}=x_{n}+\alpha_{n}(x_{n}-x_{n-1})-h^{2}\alpha_{n}\nabla f(x_{n})\\ &x_{n+1}=y_{n}-h^{2}\nabla f(y_{n})+h^{2}\nabla f(x_{n})\end{cases}

(LT-SE1)

Now, we consider a similar approach, namely on subproblem $(P_{1})$ we apply the forward Euler method, and on $(P_{2})$ we consider the Stormer-Verlet algorithm (SV₂). As before, for the first subproblem we have that

\displaystyle\begin{cases}&x_{n+1/2}=x_{n}\\ &v_{n+1/2}=\left(1-\dfrac{\alpha h}{t_{n}}\right)v_{n}-h\nabla f\left(x_{n}+h\left(1-\dfrac{\alpha h}{t_{n}}\right)v_{n}\right)+h\nabla f(x_{n})\end{cases}

Furthermore, using (SV₂) on the second sub-problem, it follows that

\displaystyle\begin{cases}&\tilde{v}_{n+1/2}=v_{n+1/2}-\dfrac{h}{2}\nabla f(x_{n+1/2})\\ &x_{n+1}=x_{n+1/2}+h\tilde{v}_{n+1/2}\\ &v_{n+1}=\tilde{v}_{n+1/2}-\dfrac{h}{2}\nabla f(x_{n+1})\end{cases}

As before, using $y_{n}=x_{n}+h\alpha_{n}v_{n}$ , where $t_{n}$ is the discretization time and $\alpha_{n}$ represents the intertial momentum coefficient, it follows that

\displaystyle\begin{cases}&x_{n+1}=x_{n}+h\alpha_{n}v_{n}-h^{2}\nabla f(y_{n})+h^{2}\nabla f(x_{n})-\dfrac{h^{2}}{2}\nabla f(x_{n})\\ &v_{n+1}=\alpha_{n}v_{n}-h\nabla f(y_{n})+h\nabla f(x_{n})-\dfrac{h}{2}\nabla f(x_{n})-\dfrac{h}{2}\nabla f(x_{n+1})\end{cases}

In a similar spirit with the determination of the algorithm (LT-SE1), we obtain our second algorithm, i.e.

\displaystyle\begin{cases}&y_{n}=x_{n}+\alpha_{n}(x_{n}-x_{n-1})-\dfrac{h^{2}}{2}\alpha_{n}\nabla f(x_{n})\\ &x_{n+1}=y_{n}-h^{2}\nabla f(y_{n})+\dfrac{h^{2}}{2}\nabla f(x_{n})\end{cases}

(LT-SV2)

We emphasize that our newly introduced algorithms (LT-SE1) and (LT-SV2) differ by a $\dfrac{h}{2}$ term at the gradient $\nabla f(x_{n})$ , in the case of the Stormer-Verlet integrator the perturbed velocity contains also a $\dfrac{h}{2}$ term, which in an intuitive sense, it annihilates much more properly large gradient values.
Now, we shift our attention to an algorithm that was developed by N.C. Nguyen et.al. in [32], that was used as a residual method for solving systems of equations. Using a variant related to the constant stepsize $h$ , the residual method is the following:

\displaystyle\begin{cases}&y_{n}=x_{n}+\alpha_{n}(x_{n}-x_{n-1})-h^{2}(1+\alpha_{n})\nabla f(x_{n})\\ &x_{n+1}=y_{n}-h^{2}\nabla f(y_{n})\end{cases}

(ARDM)

where the method reads as follows: accelerated residual descent method. Also, N.C. Nguyen et.al. have considered a restarting strategy that it is out of the scope of the present paper. Furthermore, their arguments into developing (ARDM) are based upon the stability theory related to finite difference approximations and are also based upon an approach similar to that of Su et.al. (see [39]). Our next purpose is to introduce a new dynamical system, such that the optimization algorithm (ARDM) can be considered as a Lie-Trotter splitting for a combination between the forward Euler and the symplectic Euler method (SE₂), as in the case of (Extended-AVD) dynamical system. For this, we consider the following second order evolution equation:

\displaystyle\ddot{x}(t)+\dfrac{\alpha}{t}\dot{x}(t)=-\nabla f\left(x\left(t\right)+\left(\gamma-\dfrac{\alpha\gamma^{2}}{t}\right)\dot{x}(t)-\gamma^{2}\left(2-\dfrac{\alpha\gamma}{t}\right)\nabla f(x(t))\right)-\left(2-\dfrac{\alpha\gamma}{t}\right)\nabla f(x(t))

(DynSys-ARDM)

As before, we divide this evolution equation into dissipative and conservative parts, i.e. into non-potential and potential forces:

\displaystyle\hskip-14.22636pt(P_{1}):\begin{cases}\dot{x}(t)=0\\ \dot{v}(t)=-\dfrac{\alpha}{t}v(t)+\nabla f(x(t))-\nabla f\left(x\left(t\right)+\left(\gamma-\dfrac{\alpha\gamma^{2}}{t}\right)v(t)-\gamma^{2}\left(2-\dfrac{\alpha\gamma}{t}\right)\nabla f(x(t))\right)-\left(2-\dfrac{\alpha\gamma}{t}\right)\nabla f(x(t))\end{cases}

and

\displaystyle(P_{2}):\begin{cases}\dot{x}(t)=v\\ \dot{v}(t)=-\nabla f(x(t))\end{cases}

Using a forward Euler explict method on $(P_{1})$ and taking $\gamma=h$ , we obtain that

\displaystyle\begin{cases}&x_{n+1/2}=x_{n}\\ &v_{n+1/2}=\left(1-\dfrac{\alpha h}{t_{n}}\right)v_{n}-\left(2-\dfrac{\alpha h}{t_{n}}\right)\nabla f(x_{n})-h\nabla f\left(x_{n}+h\left(1-\dfrac{\alpha h}{t_{n}}\right)v_{n}-h^{2}\left(2-\dfrac{\alpha h}{t_{n}}\right)\nabla f(x_{n})\right)+h\nabla f(x_{n})\end{cases}

On the other hand, using the symplectic Euler method (SE₂) on $(P_{2})$ , we arrive at the following equalities:

\displaystyle\begin{cases}&x_{n+1}=x_{n+1/2}+hv_{n+1}\\ &v_{n+1}=v_{n+1/2}-h\nabla f(x_{n+1/2})\end{cases}

Now, in contrast with the previous algorithm, at the discrete time $t_{n}=nh$ , we consider the additional iteration $y_{n}=x_{n}+h\alpha_{n}v_{n}-h^{2}(1+\alpha_{n})\nabla f(x_{n})$ . Using this notation, we get the following:

\displaystyle\begin{cases}&x_{n+1}=x_{n}+h\alpha_{n}v_{n}-h^{2}\nabla f(y_{n})-h^{2}(1+\alpha_{n})\nabla f(x_{n})\\ &v_{n+1}=\alpha_{n}v_{n}-h\nabla f(y_{n})-h(1+\alpha_{n})\nabla f(x_{n})\end{cases}

As in our splitting-based algorithms that were presented above, we multiply the second equation by $h$ and by substracting them, we get the algorithm (ARDM). Finally, we end this subsection by considered a remark that will represent a key role in our theoretical analysis.

Remark 6.

So far, we have considered three algorithms, namely (ARDM), (LT-SE1) and (LT-SV2). The difference between the first one and the symplectic methods applied to the (Extended-AVD) dynamical system is that the accelerated residual method contains a perturbation $h^{2}(1+\alpha_{n})$ only in the additional inertial term $y_{n}$ . The sequential Euler and Stormer-Verlet algorithms for the extended dynamical system with asymptotically vanishing damping contains a perturbation at $y_{n}$ and also at the current iteration value $x_{n+1}$ . In our terms, if an element appears at $x_{n+1}$ and it contains the gradient $\nabla f(x_{n})$ , we say that the algorithm has a bias correction term. This is the case for the (LT-SE1) and (LT-SV2) algorithms. On the other hand, we can say that (ARDM) is bias corrected free. This difference can be seen in the definition of the inertial momentum term, where at (ARDM), $y_{n}$ contains $h^{2}(1+\alpha_{n})\nabla f(x_{n})$ . So, now we consider a general form of inertial type algorithms that contain two gradient evaluations and that encompasses (ARDM), (LT-SE1) and (LT-SV2):

\displaystyle\begin{cases}&y_{n}=x_{n}+\alpha_{n}(x_{n}-x_{n-1})-\omega_{n}\nabla f(x_{n})\\ &x_{n+1}=y_{n}-h^{2}\nabla f(y_{n})+\gamma_{n}\nabla f(x_{n})\end{cases}

(LT-Symplectic)

where $\alpha_{n}=1-\dfrac{\alpha}{n}$ and $t_{n}=nh$ . Also, we point out that $\omega_{n}$ is the underlying weight that appears in the Lie-Trotter discretization and $\gamma_{n}$ is the bias correction term. Last, we observe that for (LT-SE1) we have that $\omega_{n}=h^{2}\alpha_{n}$ and $\gamma_{n}=h^{2}$ . Furthermore, for (LT-SV2), these values are halved. On the other hand, (ARDM) can be obtained by setting $\omega_{n}=h^{2}(1+\alpha_{n})$ and $\gamma_{n}=0$ , for each $n\in\mathbb{N}$ .

2.4 Fixed points, spurious solutions and exploding gradients : Heuristic interpretation

In this section we present some intuitive explanations regarding the general algorithm of type LT-Symplectic. We consider the asymptotic behavior, the critical points and the equilibrium states related to our optimization methods. Also, regarding the phenomenons that we observe, we shall present also some remedies.
First of all, we shall start with Nesterov algorithm AGM2_α. Let’s suppose that there exists $x^{\ast}\in\mathbb{R}^{N}$ such that the sequence $(x_{n})_{n\in\mathbb{N}}$ is convergent to $x^{\ast}$ . Then, let $\lim\limits y_{n}=y^{\ast}$ . We observe that $y^{\ast}=x^{\ast}$ and so $\nabla f(x^{\ast})=0$ . This is equivalent to the fact that $x^{\ast}\in\operatorname*{argmin}f$ . At the same time, define the operator $T_{h}:\mathbb{R}^{N}\to\mathbb{R}^{N}$ , as $T_{h}(y)=y-h^{2}\nabla f(y)$ . From all of this it follows that if $\lim\limits_{n\to\infty}x_{n}=x^{\ast}$ , then we have that $x^{\ast}\in\operatorname*{argmin}f$ , i.e. $x^{\ast}=T_{h}(x^{\ast})$ . So, the behavior of the Nesterov scheme is related to the fact that if the principal sequence converges, then the limit is also a minimizer of the convex function $f$ and at the same time it is a fixed point of the associated operator $T_{h}$ .
Now, we turn our attention to the IGAHD-type algorithm. Keeping in mind our notation from LT-Symplectic, let $\omega_{n}=\dfrac{\beta h}{n}$ . In this case $\lim\limits_{n\to\infty}\omega_{n}=0$ . The fact that the limit of the sequence $(\omega_{n})_{n\in\mathbb{N}}$ is $0$ will be crucial in our analysis. We shall use the fact that two consecutive iterates of the gradient represent a natural approximation of the Hessian, and this is related to $\lim\limits_{n\to\infty}\left[\nabla f(x_{n})-\nabla f(x_{n-1})\right]=0$ , since $\lim\limits_{n\to\mathbb{N}}x_{n}=x^{\ast}$ . As before, we have considered that $x^{\ast}$ the limit of $(x_{n})_{n\in\mathbb{N}}$ . Again, denoting by $y^{\ast}$ the limit of $(y_{n})_{n\in\mathbb{N}}$ , we find that, as in the case of AGM2_α, $x^{\ast}\in\operatorname*{argmin}f$ (here an important property relies on the fact that the coefficient of the difference between consecutive gradient $\beta h$ is bounded). So, the same conclusion is valid also for the IGAHD-type algorithm, regarding the fixed points of the operator $T_{h}$ and the minimizers of $f$ . So, if the major iteration of the algorithm converges, then the limit is also a minimizer of the convex objective function $f$ .
We have said that the property that $\lim\limits_{n\to\infty}\omega_{n}=0$ is crucial. We will call this phenomenon by the name annihilating the exploding gradient, since if the limit was strictly positive then the gradient at $x_{n}$ would not be disappearing at all, and so we are left with a residual that, as time grows, will tend to the gradient of $x^{\ast}$ , which could be large enough if $x^{\ast}$ was not a minimizer. For instance, let $\omega=\lim\limits_{n\in\mathbb{N}}\omega_{n}$ , such that $\omega>0$ . As before, let $\lim\limits_{n\in\mathbb{N}}x_{n}=x^{\ast}$ and $\lim\limits_{n\in\mathbb{N}}y_{n}=y^{\ast}$ . It is trivial to show that $y^{\ast}=x^{\ast}-\omega\nabla f(x^{\ast})$ , and so $x^{\ast}$ verifies the following equation:

\displaystyle\nabla f(x^{\ast})=-\dfrac{h^{2}}{\omega}\nabla f(x^{\ast}-\omega\nabla f(x^{\ast}))

(19)

It is easy to see that a minimizer satisfies 19. On the other hand, if $x^{\ast}$ is not a minimizer then we have the so called spurious roots. This means that not all the solutions of 19 are minimizers of the objective function $f$ . As in the previous cases, by the fact that $x_{n+1}=y_{n}-h^{2}\nabla f(y_{n})$ , let $T_{h}$ defined as before. A simple analysis shows that the limit $x^{\ast}$ of $(x_{n})_{n\in\mathbb{N}}$ satisfies the relation 19. Also, we have that $y^{\ast}=T_{h}(y^{\ast})$ if and only if $y^{\ast}=x^{\ast}-\omega\nabla f(x^{\ast})\in\operatorname*{argmin}f$ . This is the reason we also call these elements spurious fixed points, by the fact that they are fixed points of $T_{h}$ , but they can be represented as a perturbation of the limit point $x^{\ast}$ .
Now, let’s consider the ARDM algorithm. Using, the same notation as before, we have that $\omega=2h^{2}$ , which is strictly positive. By making the same analysis we reach the following equation between $x^{\ast}$ and $y^{\ast}$ , i.e

\displaystyle\nabla f(x^{\ast})=-\dfrac{1}{2}\nabla f(x^{\ast}-2h^{2}\nabla f(x^{\ast}))

(20)

This comes from 19 by taking $\omega=2h^{2}$ . Also, as for the case of the previous discussions regarding the modification of $\omega$ in IGAHD-type, we obtain the same conclusion regarding the spurious roots of ARDM.
The algorithm LT-SE1 meets the requirement that $\omega=h^{2}>0$ . Also, letting $\lim\limits_{n\to\infty}x^{\ast}$ and $\lim\limits_{n\to\infty y_{n}}=y^{\ast}$ , one can show that we obtain the equation

\displaystyle\nabla f(x^{\ast})=-\nabla f(x^{\ast}-h^{2}\nabla f(x^{\ast}))

(21)

So, in all of the three cases from above, namely 19, 20 and 21 we have that the limit point $x^{\ast}$ is not necessarily a minimizer of the convex function $f$ . We want also to point out that, for LT-SE1, we have that $x_{n+1}=y_{n}-h^{2}\nabla f(y_{n})+h^{2}\nabla f(x_{n})$ .In this case, let $T_{h}:\mathbb{R}^{N}\times\mathbb{R}^{N}\to\mathbb{R}^{N}$ be defined as $T_{h}(y,x)=y-h^{2}\nabla f(y)+h^{2}\nabla f(x)$ . In this case, by simple computations, we obtain that, for the case when $x^{\ast}$ is the limit of $(x_{n})_{n\in\mathbb{N}}$ and $y^{\ast}$ is the limit of $(y_{n})_{n\in\mathbb{N}}$ , we have that $T_{h}(y^{\ast},x^{\ast})=T_{h}(x^{\ast}-h^{2}\nabla f(x^{\ast}),x^{\ast})=x^{\ast}-h^{2}\nabla f(y^{\ast})$ . So, we have that $x^{\ast}=T_{h}(y^{\ast},x^{\ast})$ is equivalent with the fact that $y^{\ast}=x^{\ast}-h^{2}\nabla f(x^{\ast})\in\operatorname*{argmin}f$ . Nevertheless, let us notice that $x^{\ast}=T_{h}(x^{\ast},x^{\ast})$ means that $x^{\ast}$ is a minimizer. So, it is worth pointing out that the limit point $x^{\ast}$ minimizes $f$ if and only if is a fixed point of $T_{h}$ on the cartesian product $\mathbb{R}^{N}\times\mathbb{R}^{N}$ , i.e. a fixed point in the classical sense (i.e. on the diagonal set).
Bearing all these in mind, we present a remedy for the algorithms that belong to the class LT-Symplectic, such that they preserve fixed points, do not have the property of gradient explosion and do not introduce spurious roots. For LT-Symplectic, suppose that $\lim\limits_{n\to\infty}\omega_{n}=0$ , $\lim\limits_{n\to\infty}\gamma_{n}=\gamma\geq 0$ and that $\gamma\neq h^{2}$ . As before, let $x^{\ast}$ be the limit of $(x_{n})_{n\in\mathbb{N}}$ and $y^{\ast}$ be the limit of $(y_{n})_{n\in\mathbb{N}}$ . In this case, by the fact that $\omega_{n}$ converges to $0$ such that the gradient $\nabla f(x^{\ast})$ is annihilated, we obtain that $y^{\ast}=x^{\ast}$ . This leads to the fact that $(\gamma-h^{2})\nabla f(x^{\ast})=0$ , so $x^{\ast}\in\operatorname*{argmin}f$ . In this case, $T_{h}(y,x)=y-h^{2}\nabla f(y)+\gamma\nabla f(x)$ . In this case, $T_{h}(y^{\ast},x^{\ast})=x^{\ast}$ is equivalent with $x^{\ast}\in\operatorname*{argmin}f$ . So, with the above conditions on the coefficients, we have a natural extension of the asymptotic behavior of Nesterov algorithm AGM2_α, where the operator $T_{h}(y)$ is extended to the cartesian product by the mean of $T_{h}(y,x)$ .
We end this section, by considered another remedy but a little bit more particular, namely we modify explicitly LT-SE1, such that it has the properties presented above. In the case of LT-SE1, by our symplectic approach, we have the perturbed velocity $v_{n}=\dfrac{x_{n}-x_{n-1}}{h}-h\nabla f(x_{n})$ . This is different from the symplectic approach of Nesterov’s algorithm AGM2_α by the fact that, in our proposed scheme, $\lim\limits_{n\to\infty}v_{n}=-h\nabla f(x^{\ast})$ . Regarding the Extended-AVD dynamical system, the velocity $v$ satisfies $\lim\limits_{t\to\infty}v(t)=0$ , i.e. the velocity of the continuous system approaches $0$ as time grows. We want to reflect this property also to our discretization method. So, we want to have that $v_{n}=\dfrac{x_{n}-x_{n-1}}{h}-h\theta_{n-1}\nabla f(x_{n})$ , with $\lim\limits_{n\to\infty}\theta_{n}=0$ . We recall that in the construction of LT-SE1, we have used the technique of adding and subtracting $\nabla f(x(t))$ in order to split the potential and non-potential forces. We shall make this trick, but for $\theta(t)\nabla f(x(t))$ , where $\theta$ vanishes at $\infty$ . Furthermore, we consider the discretized version, i.e. $\theta(t_{n})\approx\theta_{n}$ .
As for LT-SE1, we consider a Lie-Trotter splitting discretization scheme by combining a forward Euler method on the first subproblem and the symplectic Euler integrator (SE₁) on the second one, on the dynamical system (Extended-AVD). Let the following two subproblems:

\displaystyle(P_{1}):\begin{cases}\dot{x}(t)=0\\ \dot{v}(t)=-\dfrac{\alpha}{t}v+\theta(t)\nabla f(x(t))-\nabla f\left(x\left(t\right)+\left(\gamma-\dfrac{\alpha\gamma^{2}}{t}\right)v(t)\right)\end{cases}\text{ and }(P_{2}):\begin{cases}\dot{x}(t)=v\\ \dot{v}(t)=-\theta(t)\nabla f(x(t))\end{cases}

By using a forward Euler discretization on $(P_{1})$ , it follows that

\displaystyle\begin{cases}&x_{n+1/2}=x_{n}\\ &v_{n+1/2}=\left(1-\dfrac{\alpha h}{t_{n}}\right)v_{n}-h\nabla f\left(x_{n}+h\left(1-\dfrac{\alpha h}{t_{n}}\right)v_{n}\right)+h\theta_{n}\nabla f(x_{n})\end{cases}

Also, applying the symplectic Euler method (SE₁) on the subproblem $(P_{2})$ , it follows that

\displaystyle\begin{cases}&x_{n+1}=x_{n+1/2}+hv_{n+1/2}\\ &v_{n+1}=v_{n+1/2}-h\theta_{n}\nabla f(x_{n+1})\end{cases}

Using the discrete time $t_{n}=nh$ and keeping in mind that $\alpha_{n}:=1-\dfrac{\alpha}{n}$ and that $y_{n}=x_{n}+h\alpha_{n}v_{n}$ , we get:

\displaystyle\begin{cases}&x_{n+1}=x_{n}+h\alpha_{n}v_{n}-h^{2}\nabla f(y_{n})+h^{2}\theta_{n}\nabla f(x_{n})\\ &v_{n+1}=\alpha_{n}v_{n}-h\nabla f(y_{n})+h\theta_{n}\nabla f(x_{n})-h\theta_{n}\nabla f(x_{n+1})\end{cases}

By multiplying the second equation by $h$ and substracting them, we obtain the modified velocity, i.e.

v_{n}=\dfrac{x_{n}-x_{n-1}}{h}-h\theta_{n-1}\nabla f(x_{n})\,.

In this way, we obtain the third symplectic-type method:

\displaystyle\begin{cases}&y_{n}=x_{n}+\alpha_{n}(x_{n}-x_{n-1})-h^{2}\alpha_{n}\theta_{n-1}\nabla f(x_{n})\\ &x_{n+1}=y_{n}-h^{2}\nabla f(y_{n})+h^{2}\theta_{n}\nabla f(x_{n})\end{cases}

(LT-SE3)

In this case, we have that $\lim\limits_{n\to\infty}\omega_{n}=\lim\limits_{n\to\infty}h^{2}\alpha_{n}\theta_{n-1}=0$ and also $\lim\limits_{n\to\infty}\gamma_{n}=\lim\limits_{n\to\infty}h^{2}\theta_{n}=0$ . Finally, we observe that we are in the case of our previous discussion regarding the remedy concerning the exploding gradients and the spurious roots.

3 Proposed algorithms and their rate of convergence

Even though we have introduced new algorithms, namely LT-Symplectic via an intuitive geometric approach by the means of symplectic methods, we will focus our attention on a generalized Hessian driven damping system and not on the extended AVD system. The discussion about our choice is postponed to Remark 11. Hence, our aim is to consider a theoretical and an empirical convergence analysis regarding Nesterov’s algorithm (AGM2_α), IGAHD algorithm from [4] and our approach via symplectic integrators, i.e. algorithms of the form LT-S-IGAHD.

3.1 Modified Hessian-driven damping continuous and discrete systems

In the previous sections, we have presented the key idea behind the splitting technique applied to gradient-type dynamical systems. We have considered LT-Symplectic and LT-SE3 algorithms that have a geometric interpretation that is related to symplectic integrators. We observe that all these optimization algorithms have two coefficients in front of the gradient value at $x_{n}$ , namely $\gamma_{n}$ and $\omega_{n}$ , for each $n\in\mathbb{N}$ . On the other hand, there are the Hessian-driven damping discretizations like IGAHD-type, which contain the difference of two consecutive gradients, i.e $\nabla f(x_{n})-\nabla f(x_{n-1})$ along with $\nabla f(x_{n})$ , under the gradient of the inertial element $y_{n}$ . By the fact that the effect of the Hessian is to dampen the oscillations, it is more natural to consider a Hessian-driven damping algorithm. So, in the present subsection we shall modify some of the previous dynamical systems in order to make consistent modification into our proposed algorithms. At the same time, our algorithms that we propose shadow the influence of the spurious roots and exploding gradients (explained in the Subsection 2.4), preserve the Hessian structure and are derived from a symplectic perspective as the other methods derived above. In order to do this, for $n\geq 1$ we consider the following generalization of IGAHD-type, namely the Lie-Trotter-Symplectic-IGAHD, i.e.

\displaystyle\begin{cases}&y_{n}=x_{n}+\alpha_{n}(x_{n}-x_{n-1})-\lambda_{n}\left[\nabla f(x_{n})-\nabla f(x_{n-1})\right]-\omega_{n}\nabla f(x_{n})\\ &x_{n+1}=y_{n}-h^{2}\nabla f(y_{n})+\gamma_{n}\nabla f(x_{n})\end{cases}

(LT-S-IGAHD)

where $\lambda_{n}$ and the coefficients $\omega_{n}$ and $\gamma_{n}$ are consistent with the notations from LT-Symplectic. Furthermore, for $\gamma_{n}=0$ , $\lambda_{n}=\beta h$ and $\omega_{n}=\dfrac{\beta h}{n}$ , we obtain as a particular case IGAHD-type. On the other hand, taking into account Subsection 2.4, it is naturally to consider $\omega_{n}$ such as $\lim\limits_{n\to\infty}\omega_{n}=0$ . Also, for the case of IGAHD-type this property is easily verifiable, with $\omega_{n}=O\left(\dfrac{1}{n}\right)$ as $n\to\infty$ . Now, we present our result of this subsection, in which we introduce a new dynamical system such that LT-S-IGAHD is a Lie-Trotter splitting discretization.

Theorem 7.

Let $f\in C^{2}(\mathbb{R}^{N},\mathbb{R})$ and consider the following Hessian driven damping second order system :

\displaystyle\hskip-14.22636pt\ddot{x}(t)+\dfrac{\alpha}{t}\dot{x}(t)+\dfrac{\lambda(t)}{\gamma}\nabla f^{2}(x(t))\dot{x}(t)=-\dfrac{1}{\gamma^{2}}\left(\omega(t)-\tilde{\gamma}(t)\right)\nabla f(x(t))-\nabla f\left(\mathcal{A}(t,x(t),\dot{x}(t))\right)

where $\tilde{\gamma}=\tilde{\gamma}(t)$ and $\omega=\omega(t)$ are given functions. Also, for a pair of the form $(t,x(t),y(t))$ we have used the notation

\mathcal{A}(t,x(t),y(t)):=x(t)+\gamma\left(1-\dfrac{\alpha\gamma}{t}\right)y(t)-\gamma\lambda(t)\nabla^{2}f(x(t))y(t)-\omega(t)\nabla f(x(t))

Then, a Lie-Trotter splitting discretization, with a forward Euler method on the first subproblem and a symplectic integrator on the second one leads to LT-S-IGAHD algorithm.

Proof .

As in the case of (IGAHD-type) we consider two subproblems, one containing the so-called non-potential forces and the other one containing the potential forces :

\displaystyle\hskip-14.22636pt(P_{1}):\begin{cases}\dot{x}(t)=0\\ \dot{v}(t)=-\dfrac{\alpha}{t}v(t)-\dfrac{\lambda(t)}{\gamma}\nabla f^{2}(x(t))v(t)-\dfrac{\omega(t)}{\gamma^{2}}\nabla f(x(t))+\dfrac{\tilde{\gamma}(t)}{\gamma^{2}}\nabla f(x(t))-\nabla f\left(\mathcal{A}(t,x(t),v(t))\right)+\nabla f(x(t))\end{cases}

and the second subproblem as

\displaystyle(P_{2}):\begin{cases}\dot{x}(t)=v(t)\\ \dot{v}(t)=-\nabla f(x(t))\end{cases}

Taking $\gamma=h$ , approximating $\tilde{\gamma}(t_{n})\approx\gamma_{n}$ , $\lambda(t_{n})\approx\lambda_{n}$ and $\omega(t_{n})\approx\omega_{n}$ , and defining $\alpha_{n}:=1-\dfrac{\alpha h}{t_{n}}$ we obtain that

\displaystyle\begin{cases}x_{n+1/2}=x_{n}\\ v_{n+1/2}=\alpha_{n}v_{n}-\lambda_{n}\nabla^{2}f(x_{n})v_{n}-\dfrac{\omega_{n}}{h}\nabla f(x_{n})+\dfrac{\gamma_{n}}{h}\nabla f(x_{n})+h\nabla f(x_{n})-h\nabla f\left(y_{n}\right)\end{cases}

where we have the notation for the auxiliary iteration element $y_{n}:=\mathcal{A}(t_{n},x_{n},v_{n})$ , i.e. $y_{n}=x_{n}+h\alpha_{n}v_{n}-h\lambda_{n}\nabla^{2}f(x_{n})v_{n}-\omega_{n}\nabla f(x_{n})$ . Furthermore, for the second subproblem $(P_{2})$ , we apply the symplectic Euler integrator (SE₂) and obtain that

\displaystyle\begin{cases}x_{n+1}=x_{n+1/2}+hv_{n+1}\\ v_{n+1}=v_{n+1/2}-h\nabla f(x_{n+1/2})\end{cases}

Using the fact that for $t_{n}=nh$ , the splitting solution is the discrete solution of $(P_{2})$ , namely $(x_{n+1},v_{n+1})^{T}\in\mathbb{R}^{N}\times\mathbb{R}^{N}$ , we obtain, from the definition of $x_{n+1}$ , that $v_{n}=\dfrac{x_{n}-x_{n-1}}{h}$ . Using the approximation technique for $\nabla^{2}f(x_{n})v_{n}$ as in Theorem 4, it follows that our algorithm is identical with LT-S-IGAHD. $\blacksquare$

Now, we will briefly give a chain of observations that will give insights into the construction of LT-S-IGAHD type algorithms.

Remark 8.

The same principle as derived in Theorem 7 can be applied to the DynSys-ARDM. It is easy to see that if we add a perturbation of the form $\dfrac{\tilde{\gamma}(t)}{\gamma^{2}}\nabla f(x(t))$ , the we will obtain an additional element $\gamma_{n}\nabla f(x_{n})$ in the discrete version of the exact solution, i.e. $x_{n+1}$ . Furthermore, we will explain in detail the technical reason why the new added iteration is significant in the process of creating new splitting-type algorithms. For example, we refer to the particular case of the dynamical system from Theorem 7, for the case when $\omega(t)=\dfrac{\beta\gamma^{2}}{t}$ . This is the system corresponding to the IGAHD dynamical system from Theorem 4, but with an additional gradient with coefficient $\dfrac{\tilde{\gamma}(t)}{\gamma^{2}}$ . Even though we have two gradients in the right hand side, namely $\dfrac{\beta}{t}\nabla f(x(t))$ and $\dfrac{\tilde{\gamma}(t)}{\gamma^{2}}\nabla f(x(t))$ , the first one is proportional to the one under the gradient of $\mathcal{A}(t,x(t),\dot{x}(t))$ , i.e. $\dfrac{\beta\gamma^{2}}{t}\nabla f(x(t))$ . This basic observation (also valid for the term involving the Hessian, which also has $\gamma^{2}$ as a coefficient of proportionality) leads to the crucial idea in which we can construct algorithms of the form LT-S-IGAHD. This is equivalent to the fact that if a term proportional to the gradient at the current value does not appear under the gradient of $\mathcal{A}(t,x(t),x,\dot{x}(t))$ , then this term will not appear in the inertial term $y_{n}$ but will appear in $x_{n+1}$ , i.e. it is added to $y_{n}$ and $\nabla f(y_{n})$ .

Remark 9.

In the Subsection 2.4 we have pointed out that in order to make the algorithms have the natural behavior in the sense that they do not possess exploding gradients or spurious fixed points, a natural technique is to make the division into potential and non-potential forces and add and substract $\theta(t)\nabla f(x(t))$ as in the case of LT-SE3. This method was used for the LT-SE3 algorithm by the fact that we have employed the SE₁ symplectic integrator. In the case of LT-S-IGAHD (as in the case of Nesterov’s method AGM2_α), the symplectic numerical scheme SE₂ does not require such a process.

Remark 10.

One might argue that our Hessian-type algorithms LT-S-IGAHD have a major drawback, namely we require the computation of $\nabla f(x_{n})$ at $x_{n+1}$ , along with $\nabla f(y_{n})$ . This is not true since $\nabla f(x_{n})$ was already computed in order to find the value of the inertial term $y_{n}$ so, the algorithms require the same computational time.

Remark 11.

One might wonder why we do not analyzed theoretically LT-Symplectic type algorithms and why we have introduced an additional gradient value $\nabla f(x_{n-1})$ . The motivation is two-fold. First of all, the product $\nabla^{2}f(x(t))\dot{x}(t)$ is equal to $\dfrac{d}{dt}\left(\nabla f(x(t))\right)$ and so it is natural to consider the effects of the Hessian perturbation. Also, the Hessian-type systems lead to a much clearer interpretation of numerical methods in optimization (for this see Subsection 2.1). Second of all, the previous value under the gradient $\nabla f(x_{n-1})$ leads to an intuitive discrete Lyapunov function as we will see in the main result of the next subsection. Without this term, through tedious computations, one can see that the Lyapunov function can not be naturally chosen.

Remark 12.

Recently, in [18], L. Chen and H. Luo have presented some Hessian-driven damping dynamical systems. They have incorporated the effects of a dynamically changing damping coefficient. Also, they have considered alternative proofs concerning the continuous and the discrete energy functions using the so-called strong Lyapunov property. On the other hand, their optimization methods are in connection to the first order system associated to the second order evolution equation containing the Hessian. This idea is related to the technique of A.M. Oberman and M. Laborde [33], in which Nesterov type accelerated gradient methods can be interpreted as explicit discretizations of Hessian systems, but with an extra gradient descent step. We point out that the algorithms developed in [18] come from an interpretation different to that developed by us in Subsection 2.1. Our algorithms have a beautiful geometric interpretation related to the idea of symplectic integrators that almost conserve the energy and splitting methods that separate the potential and non-potential forces. Also, they have presented that Nesterov’s method is an explicit method for a Hessian type system but with an addition gradient step. This can be seen in contrast with the present work where we actually show in a rigorous manner how AGM2_α can be constructed through the lens of the Lie-Trotter splitting.

3.2 The rate of convergence of the optimization methods

In the present section, we shall consider some necessary tools for the convergence analysis of the Hessian-driven damping discretization LT-S-IGAHD. For this, we recall that an important assumption concerning the theoretical rate of convergence of optimization algorithms that the obhective function has a Lipschitz continuous gradient. An interesting consequence is the so-called descent lemma (see inequality 2.1.6 from page 56 of [31]), which is presented below.

Lemma 13 (Descent lemma).

Let $f:\mathbb{R}^{N}\to\mathbb{R}$ be Fŕechet differentiable, such that $\nabla f:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is $L$ - Lipschitz continuous. Then, one has that

\displaystyle f(y)\leq f(x)+\langle\nabla f(x),y-x\rangle+\dfrac{L}{2}\|y-x\|^{2}\quad,\quad\forall x,y\in\mathbb{R}^{N}\quad.

(DL)

It is worth noticing that if the objective function is convex, then we have actually an equivalence between DL and the $L$ - Lipschitz continuity of $\nabla f$ . Furthermore, Lemma DL presents a quadratic upper bound of the underlying objective function. Now, we shall briefly recall from [4] extended variant of Lemma 13, where a modified variant of DL contains also the norm of the difference between two values of $\nabla f$ (as a side note, see also inequality 2.17 from page 57 of [31]).

Lemma 14 (Extended descent lemma).

Let $f:\mathbb{R}^{N}\to\mathbb{R}$ be convex and Fréchet differentiable, such that $\nabla f:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is $L$ - Lipschitz continuous. Consider a parameter $s$ , such that $0<s\leq\dfrac{1}{L}$ . Then one has that

\displaystyle f(y-s\nabla f(y))\leq f(x)+\langle\nabla f(y),y-x\rangle-\dfrac{s}{2}\|\nabla f(y)\|^{2}-\dfrac{s}{2}\|\nabla f(x)-\nabla f(y)\|^{2}\quad,\quad\forall x,y\in\mathbb{R}^{N}\quad.

(EDL)

In the following sequel, we present a technical lemma that will play a crucial role in our main result regarding the asymptotic convergence of the Hessian-type algorithms. So, it is the time to ’extend’ the ’Extended descent lemma (Lemma 14). This represents our first result of this section.

Lemma 15 (Three points extended lemma).

Let $f:\mathbb{R}^{N}\to\mathbb{R}$ be a convex and Fréchet differentiable function, such that $\nabla f:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is $L$ - Lipschitz continuous. Let $s\in\left(0,\dfrac{1}{L}\right]$ . For simplicity, consider the following notations: $\mathcal{A}_{1}(s):=s$ , $\mathcal{A}_{2}(s):=-\mathcal{A}_{1}(s)$ , $\mathcal{A}_{3}(s;\gamma,L):=-\gamma(1-Ls)$ , $\mathcal{A}_{4}(s):=\mathcal{A}_{1}(s)-\dfrac{s}{2}$ and $\mathcal{A}_{5}(s;\gamma):=-\gamma^{2}/(2s)$ , where $\gamma\geq 0$ . Then $\forall x,y,z\in\mathbb{R}^{N}$ , one has that $f$ satisfies the following inequality:

$\displaystyle f(y-s\nabla f(y)+\gamma\nabla f(z))$	$\displaystyle\leq f(x)+\langle\nabla f(y),y-x\rangle-\mathcal{A}_{1}(s)\\|\nabla f(y)\\|^{2}$	(E-EDL)
	$\displaystyle-\left[\mathcal{A}_{2}(s)\langle\nabla f(y),\nabla f(x)\rangle+\mathcal{A}_{3}(s;\gamma,L)\langle\nabla f(y),\nabla f(z)\rangle\right]$
	$\displaystyle-\left[\mathcal{A}_{4}(s)\\|\nabla f(x)\\|^{2}+\mathcal{A}_{5}(s;\gamma)\\|\nabla f(z)\\|^{2}\right]\quad.$

Proof .

Let $x,y,z\in\mathbb{R}^{N}$ . For simplicity, consider $y^{+}:=y-s\nabla f(y)+\gamma\nabla f(z)$ . Applying DL for $y^{+}$ and $y$ , it follows that

\displaystyle f(y^{+})\leq f(y)+\langle\nabla f(y),y^{+}-y\rangle+\dfrac{L}{2}\|y^{+}-y\|^{2}\quad.

This leads to

	$\displaystyle f(y^{+})\leq f(y)+\langle\nabla f(y),-s\nabla f(y)+\gamma\nabla f(z)\rangle+\dfrac{L}{2}\\|-s\nabla f(y)+\gamma\nabla f(z)\\|^{2}$
	$\displaystyle f(y^{+})\leq f(y)-s\\|\nabla f(y)\\|^{2}+\gamma\langle\nabla f(y),\nabla f(z)\rangle+\dfrac{L}{2}\langle-s\nabla f(y)+\gamma\nabla f(z),-s\nabla f(y)+\gamma\nabla f(z)\rangle$
	$\displaystyle f(y^{+})\leq f(y)-s\\|\nabla f(y)\\|^{2}+\gamma\langle\nabla f(y),\nabla f(z)\rangle+\dfrac{L}{2}s^{2}\\|\nabla f(y)\\|^{2}+\dfrac{L}{2}\gamma^{2}\\|\nabla f(z)\\|^{2}-Ls\gamma\langle\nabla f(y),\nabla f(z)\rangle$

So, it follows that

\displaystyle f(y^{+})\leq f(y)-\dfrac{s}{2}\left(2-Ls\right)\|\nabla f(y)\|^{2}+\dfrac{L}{2}\gamma^{2}\|\nabla f(z)\|^{2}+\gamma(1-Ls)\langle\nabla f(y),\nabla f(z)\rangle

(22)

Taking into account that $\nabla f$ is $L$ - Lipschitz continuous is equivalent to the fact that the conjugate function $f^{\ast}$ is $\dfrac{1}{L}$ - strongly convex and applying this property for the pair $(\nabla f(y),\nabla f(x))$ , we get that

\displaystyle f^{\ast}(\nabla f(y))\geq f^{\ast}(\nabla f(x))+\langle\nabla f^{\ast}(\nabla f(x)),\nabla f(y)-\nabla f(x)\rangle+\dfrac{1}{2L}\|\nabla f(y)-\nabla f(x)\|^{2}

Using the fact that $(\nabla f)^{-1}=\nabla f^{\ast}$ , it simplifies to

\displaystyle f^{\ast}(\nabla f(y))\geq f^{\ast}(\nabla f(x))+\langle x,\nabla f(y)-\nabla f(x)\rangle+\dfrac{1}{2L}\|\nabla f(y)-\nabla f(x)\|^{2}

Using Fenchel identity for $y$ , namely $f(y)=\langle\nabla f(y),y\rangle-f^{\ast}(\nabla f(y))$ , we get

\displaystyle\langle\nabla f(y),y\rangle-f(y)\geq f^{\ast}(\nabla f(x))+\langle x,\nabla f(y)-\nabla f(x)\rangle+\dfrac{1}{2L}\|\nabla f(y)-\nabla f(x)\|^{2}

So, this is equivalent to

	$\displaystyle f(y)\leq\langle\nabla f(y),y\rangle-f^{\ast}(\nabla f(x))-\langle x,\nabla f(y)-\nabla f(x)\rangle-\dfrac{1}{2L}\\|\nabla f(y)-\nabla f(x)\\|^{2}$
	$\displaystyle f(y)\leq\left[\langle x,\nabla f(x)\rangle-f^{\ast}(\nabla f(x))\right]+\langle\nabla f(y),y-x\rangle-\dfrac{1}{2L}\\|\nabla f(y)-\nabla f(x)\\|^{2}$

Again using Fenchel identity but applied on $x$ , namely $f(x)=\langle x,\nabla f(x)\rangle-f^{\ast}(\nabla f(x))$ , it follows that

\displaystyle f(y)\leq f(x)+\langle\nabla f(y),y-x\rangle-\dfrac{1}{2L}\|\nabla f(y)-\nabla f(x)\|^{2}\quad.

(23)

Taking into account inequalities 22 and 23, thus we obtain

	$\displaystyle f(y^{+})$	$\displaystyle\leq f(x)+\langle\nabla f(y),y-x\rangle-\dfrac{1}{2L}\\|\nabla f(y)-\nabla f(x)\\|^{2}-\dfrac{s}{2}\left(2-Ls\right)\\|\nabla f(y)\\|^{2}$		(24)
		$\displaystyle+\dfrac{L}{2}\gamma^{2}\\|\nabla f(z)\\|^{2}+\gamma(1-Ls)\langle\nabla f(z),\nabla f(y)\rangle\quad.$

From the theorem’s assumptions, we know that $\dfrac{s}{2}\leq s-\dfrac{Ls^{2}}{2}$ , $\dfrac{1}{L}\geq s$ and $-\dfrac{1}{2L}\leq-\dfrac{s}{2}$ . From the inequality 24 and from the condition on $L$ , one obtains that

	$\displaystyle f(y^{+})$	$\displaystyle\leq f(x)+\langle\nabla f(y),y-x\rangle-\dfrac{s}{2}\\|\nabla f(y)-\nabla f(x)\\|^{2}-\dfrac{s}{2}\\|\nabla f(y)\\|^{2}$
		$\displaystyle+\dfrac{L\gamma^{2}}{2}\\|\nabla f(z)\\|^{2}+\gamma(1-Ls)\langle\nabla f(z),\nabla f(y)\rangle$

Through some simple computations, it follows

	$\displaystyle f(y^{+})$	$\displaystyle\leq f(x)+\langle\nabla f(y),y-x\rangle-\dfrac{s}{2}\\|\nabla f(y)\\|^{2}-\dfrac{s}{2}\\|\nabla f(x)\\|^{2}$
		$\displaystyle+s\langle\nabla f(y),\nabla f(x)\rangle-\dfrac{s}{2}\\|\nabla f(y)\\|^{2}+\dfrac{L\gamma^{2}}{2}\\|\nabla f(z)\\|^{2}+\gamma(1-Ls)\langle\nabla f(z),\nabla f(y)\rangle$

Hence,

	$\displaystyle f(y^{+})$	$\displaystyle\leq f(x)+\langle\nabla f(y),y-x\rangle-\left[\dfrac{s}{2}+\dfrac{s}{2}\right]\\|\nabla f(y)\\|^{2}-\left[\dfrac{s}{2}\\|\nabla f(x)\\|^{2}-\dfrac{L\gamma^{2}}{2}\\|\nabla f(z)\\|^{2}\right]$
		$\displaystyle-\left[-s\langle\nabla f(y),\nabla f(x)\rangle-\gamma(1-Ls)\langle\nabla f(y),\nabla f(z)\rangle)\right]$

Using the fact $\dfrac{L}{2}\leq\dfrac{1}{2s}$ , the last inequality leads to

	$\displaystyle f(y^{+})$	$\displaystyle\leq f(x)+\langle\nabla f(y),y-x\rangle-s\\|\nabla f(y)\\|^{2}-\left[\dfrac{s}{2}\\|\nabla f(x)\\|^{2}-\dfrac{\gamma^{2}}{2s}\\|\nabla f(z)\\|^{2}\right]$
		$\displaystyle-\left[-s\langle\nabla f(y),\nabla f(x)\rangle-\gamma(1-Ls)\langle\nabla f(y),\nabla f(z)\rangle\right]$

Using the notations for $\mathcal{A}_{1}(s)$ , $\mathcal{A}_{2}(s)$ , $\mathcal{A}_{3}(s;\gamma,L)$ , $\mathcal{A}_{4}(s)$ and $\mathcal{A}_{5}(s;\gamma)$ , the inequality E-EDL follows easily. $\blacksquare$

Corollary 16.

Consider $x,y,z\in\mathbb{R}^{N}$ . By the fact that the stepsize $s$ lies in $\left(0,\dfrac{1}{L}\right]$ , it follows that the major bound E-EDL for $f(y^{+})$ reduces to EDL, with $\gamma=0$ .

Before embarking into our main result of the present paper, we consider some technical lemmas that will be used in the forthcoming sequel. These results involve inequalities that are related to quadratic inequalities with nonconstant coefficients. Even though they are elementary results, for brevity we give them along with their proofs.

Lemma 17.

Let $a,b,c:\mathbb{R}\to\mathbb{R}$ be given functions. If there exists $x\in\mathbb{R}$ , such that the following inequalities are satisfied: $a(x)>0$ and $b^{2}(x)-4a(x)c(x)\leq 0$ , then $x$ satisfies $a(x)x^{2}+b(x)x+c(x)\geq 0$ .

Proof .

First of all, we consider the following computations:

	$\displaystyle a(x)\left(x+\dfrac{b(x)}{2a(x)}\right)^{2}-\dfrac{b^{2}(x)-4a(x)c(x)}{4a(x)}$	$\displaystyle=a(x)\left[x^{2}+\dfrac{b^{2}(x)}{4a^{2}(x)}+x\dfrac{b(x)}{a(x)}\right]-\dfrac{b^{2}(x)}{4a(x)}+c(x)$
		$\displaystyle=a(x)x^{2}+\dfrac{b^{2}(x)}{4a^{2}(x)}a(x)+x\dfrac{b(x)}{a(x)}a(x)-\dfrac{b^{2}(x)}{4a(x)}+c(x)\,.$

This is equivalent to the fact that

a(x)x^{2}+b(x)x+c(x)=a(x)\left(x+\dfrac{b(x)}{2a(x)}\right)^{2}-\dfrac{b^{2}(x)-4a(x)c(x)}{4a(x)}\,.

Now, $a(x)x^{2}+b(x)x+c(x)\geq 0\Longleftrightarrow a(x)\left(x+\dfrac{b(x)}{2a(x)}\right)^{2}\geq\dfrac{b^{2}(x)-4a(x)c(x)}{4a(x)}$ . Taking into account that $b^{2}(x)-4a(x)c(x)\geq 0$ and dividing by $a(x)>0$ , it follows that $\left(x+\dfrac{b(x)}{2a(x)}\right)^{2}\geq\dfrac{b^{2}(x)-4a(x)c(x)}{4a^{2}(x)}$ . Likewise, let $\left(x+\dfrac{b(x)}{2a(x)}\right)^{2}\geq\dfrac{b^{2}(x)-4a(x)c(x)}{4a^{2}(x)}$ . As before, we observe that these terms are positive and multiplying be $a(x)>0$ , we obtain that $a(x)x^{2}+b(x)x+c(x)\geq 0$ . So, wrapping things up, we need to show that $\left(x+\dfrac{b(x)}{2a(x)}\right)^{2}\geq\dfrac{b^{2}(x)-4a(x)c(x)}{4a^{2}(x)}$ , since this is equivalent to the quadratic inequality with nonconstant coefficients. Because $b^{2}(x)-a(x)c(x)\leq 0$ , then it follows that $\dfrac{b^{2}(x)-4a(x)c(x)}{4a^{2}(x)}\leq 0$ . Combining this with the fact that the left hand side of the inequality is positive, the conclusion stands up. $\blacksquare$

Lemma 18.

Let $a,b,c:\mathbb{R}\to\mathbb{R}$ be given functions. If there exists $x\in\mathbb{R}$ , such that the following inequalities are satisfied: $a(x)>0$ , $b^{2}(x)-4a(x)c(x)\geq 0$ and $x\in I_{1}\cup I_{2}$ , where: $I_{1}:=\left(-\infty,\dfrac{-b(x)-\sqrt{b^{2}(x)-4a(x)c(x)}}{2a(x)}\right]$ and $I_{2}:=\left[\dfrac{-b(x)+\sqrt{b^{2}(x)-4a(x)c(x)}}{2a(x)},+\infty\right)$ , then $x$ satisfies $a(x)x^{2}+b(x)x+c(x)\geq 0$ .

Proof .

As before, we must show that $\left(x+\dfrac{b(x)}{2a(x)}\right)^{2}\geq\dfrac{b^{2}(x)-4a(x)c(x)}{4a^{2}(x)}$ , so we must analyze two cases:

1)

$x\in I_{1}$ , that is equivalent to $x\leq-\dfrac{b(x)}{2a(x)}-\dfrac{\sqrt{b^{2}(x)-4a(x)c(x)}}{2a(x)}$ . Thus, we have that $\dfrac{\sqrt{b^{2}(x)-4a(x)c(x)}}{2a(x)}\leq-x-\dfrac{b(x)}{2a(x)}=-\left(x+\dfrac{b(x)}{2a(x)}\right)$ . But, it is easy to see that $x\leq-\dfrac{b(x)}{2a(x)}$ , i.e. $-\left(x+\dfrac{b(x)}{2a(x)}\right)\geq 0$ . So, it follows that $\dfrac{\sqrt{b^{2}(x)-4a(x)c(x)}}{2a(x)}\leq-\left(x+\dfrac{b(x)}{2a(x)}\right)=\Big{|}x+\dfrac{b(x)}{2a(x)}\Big{|}$ . Squaring the positive terms, it follows that $\dfrac{b^{2}(x)-4a(x)c(x)}{4a^{2}(x)}\leq\left(x+\dfrac{b(x)}{2a(x)}\right)^{2}$ .
2)

$x\in I_{2}$ , i.e. $x\geq-\dfrac{b(x)}{2a(x)}+\dfrac{\sqrt{b^{2}(x)-4a(x)c(x)}}{2a(x)}$ , which is equivalent to $\dfrac{\sqrt{b^{2}(x)-4a(x)c(x)}}{2a(x)}\leq x+\dfrac{b(x)}{2a(x)}$ . But, since $x+\dfrac{b(x)}{2a(x)}\geq\dfrac{\sqrt{b^{2}(x)-4a(x)c(x)}}{2a(x)}$ , it follows that $x+\dfrac{b(x)}{2a(x)}\geq 0$ , so $x+\dfrac{b(x)}{2a(x)}=\Big{|}x+\dfrac{b(x)}{2a(x)}\Big{|}$ . Hence, using that $\dfrac{\sqrt{b^{2}(x)-4a(x)c(x)}}{2a(x)}\leq\Big{|}x+\dfrac{b(x)}{2a(x)}\Big{|}$ and squaring the positive terms, we get the conclusion that $\dfrac{b^{2}(x)-4a(x)c(x)}{4a^{2}(x)}\leq\left(x+\dfrac{b(x)}{2a(x)}\right)^{2}$ .

$\blacksquare$

Until our main result, we consider some observations regarding the functions that appear in Lemma 17 and Lemma 18.

Remark 19.

Related to Lemma 17 and Lemma 18, we have the following:

•

For a given pair of functions $(a(x),b(x),c(x))$ , if the element $x\in\mathbb{R}$ exists, then the quadratic inequality takes place. But this does not mean that such a point exists for all pair of functions $(a(x),b(x),c(x))$ .
•

It is easy to observe that, in Lemma 17, the conditions $a(x)>0$ and $b^{2}(x)-4a(x)c(x)\leq 0$ implies that $c(x)\geq 0$ . Furthermore, we do not have a positivity-type inequality on the function $b=b(x)$ .
•

In Lemma 18, we observe that $a(x)>0$ and if $c(x)\leq 0$ , then it implies that $4a(x)c(x)\leq b^{2}(x)$ is trivially satisfied. As before, we do not need the positivity of $b=b(x)$ .

In what follows we have our main result of this section. Given a convex objective function endowed with Lipschitz continuous gradient, we show that the rate of convergence of the algorithms given by the general discretization LT-S-IGAHD is $O\left(\dfrac{1}{n^{2}}\right)$ as $n\to\infty$ , under some suitable conditions over the parameters. Given $x^{\ast}\in\operatorname*{argmin}f$ (under the assumption that $\operatorname*{argmin}f\neq\emptyset$ ), the Lyapunov analysis that we employ is based upon the discrete energy defined as:

		$\displaystyle\mathcal{E}_{n}:=t_{n}^{2}\left(f(x_{n})-f(x^{\ast})\right)+\dfrac{1}{2s}\\|z_{n}\\|^{2}$		(Discrete-Energy)
		$\displaystyle z_{n}:=(x_{n-1}-x^{\ast})+t_{n}(x_{n}-x_{n-1})+\lambda_{n}t_{n+1}\nabla f(x_{n-1})\,,$

where $t_{n+1}=\dfrac{n}{\alpha-1}$ , with $\alpha\geq 3$ . We shall also use $\alpha_{n}=\dfrac{n-\alpha}{n}$ , so it follows that $t_{n}-1=\alpha_{n}t_{n+1}$ . Moreover, our technique for the proof is inspired by the ones given in [4] and [6].

Theorem 20.

Let $f:\mathbb{R}^{N}\to\mathbb{R}$ be a convex and Fréchet differentiable function whose gradient $\nabla f:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is $L$ - Lipschitz continuous. Also, suppose that $\operatorname*{argmin}f\neq\emptyset$ . Let $x_{n}$ be generated by the algorithm LT-S-IGAHD for every $n\geq 1$ and, for simplicity, denote $h^{2}$ as $s$ , which satisfies the condition $s\in\left(0,\dfrac{1}{L}\right]$ . Suppose that the sequences $(\gamma_{n})_{n\in\mathbb{N}}$ and $(\lambda_{n})_{n\in\mathbb{N}}$ are convergent. Likewise, assume that the following assumptions are satisfied:

i)

there exists $N^{\prime}\in\mathbb{N}$ such that $\left[s(1-\gamma_{n}L)-\dfrac{n+1}{n}\lambda_{n+1}\right]^{2}<s^{2}-\gamma_{n}^{2}\,,\forall n>max\{\alpha-1,N^{\prime}\}$ ;
ii)

$\gamma_{n}=(\lambda_{n}+\omega_{n})-\dfrac{n+1}{n}\lambda_{n+1}\,,\forall n\geq\alpha-1$ .

For $x^{\ast}\in\operatorname*{argmin}f$ , the sequence $(\mathcal{E})_{n\in\mathbb{N}}$ defined in Discrete-Energy is non-increasing and the following converge rate is satisfied:

f(x_{n})-f(x^{\ast})=O\left(\dfrac{1}{n^{2}}\right)\text{ as }n\to\infty

Proof .

( Step I ) The first step of the proof is to apply the so-called Three points extended lemma, namely Lemma 15 and then construct the difference $\mathcal{E}_{n+1}-\mathcal{E}_{n}$ of the discrete energy Discrete-Energy.
We apply the inequality E-EDL with the pair $x=z=x_{n}$ , $y=y_{n}$ and $\gamma=\gamma_{n}$ . Then, we obtain that:

	$\displaystyle f(y_{n}-s\nabla f(y_{n})+\gamma_{n}\nabla f(x_{n}))$	$\displaystyle\leq f(x_{n})+\langle\nabla f(y_{n}),y_{n}-x_{n}\rangle-\mathcal{A}_{1}(s)\\|\nabla f(y_{n})\\|^{2}$
		$\displaystyle-\left[\mathcal{A}_{2}(s)\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle+\mathcal{A}_{3}(s;\gamma_{n},L)\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle\right]$
		$\displaystyle-\left[\mathcal{A}_{4}(s)\\|\nabla f(x_{n})\\|^{2}+\mathcal{A}_{5}(s;\gamma_{n})\\|\nabla f(x_{n})\\|^{2}\right]\,.$

This means that

$\displaystyle f(x_{n+1})$	$\displaystyle\leq f(x_{n})+\langle\nabla f(y_{n}),y_{n}-x_{n}\rangle-\mathcal{A}_{1}(s)\\|\nabla f(y_{n})\\|^{2}$	(25)
	$\displaystyle-\left[\mathcal{A}_{2}(s)+\mathcal{A}_{3}(s;\gamma_{n},L)\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$
	$\displaystyle-\left[\mathcal{A}_{4}(s)+\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}\,.$

As before, we will apply E-EDL, but with the choices $x=x^{\ast}$ , $y=y_{n}$ , $z=x_{n}$ and $\gamma=\gamma_{n}$ . It follows that:

	$\displaystyle f(y_{n}-s\nabla f(y_{n})+\gamma_{n}\nabla f(x_{n}))$	$\displaystyle\leq f(x^{\ast})+\langle\nabla f(y_{n}),y_{n}-x^{\ast}\rangle-\mathcal{A}_{1}(s)\\|\nabla f(y_{n})\\|^{2}$
		$\displaystyle-\left[\mathcal{A}_{2}(s)\langle\nabla f(y_{n}),\nabla f(x^{\ast})\rangle+\mathcal{A}_{3}(s;\gamma_{n},L)\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle\right]$
		$\displaystyle-\left[\mathcal{A}_{4}(s)\\|\nabla f(x^{\ast})\\|^{2}+\mathcal{A}_{5}(s;\gamma_{n})\\|\nabla f(x_{n})\\|^{2}\right]\,.$

Using the fact that $x^{\ast}\in\operatorname*{argmin}f$ , so $\nabla f(x^{\ast})=0$ , we obtain

	$\displaystyle f(x_{n+1})$	$\displaystyle\leq f(x^{\ast})+\langle\nabla f(y_{n}),y_{n}-x^{\ast}\rangle-\mathcal{A}_{1}(s)\\|\nabla f(y_{n})\\|^{2}$		(26)
		$\displaystyle-\mathcal{A}_{3}(s;\gamma_{n},L)\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle-\mathcal{A}_{5}(s;\gamma_{n})\\|\nabla f(x_{n})\\|^{2}\,.$

Using the definition of $t_{n+1}$ , we consider

\displaystyle n\geq\alpha-1\,.

(27)

hence $t_{n+1}-1\geq 0$ . Now, we multiply 25 by $(t_{n+1}-1)$ and then add it to 26. One obtains that

	$\displaystyle(t_{n+1}-1)f(x_{n+1})+f(x_{n+1})$	$\displaystyle\leq(t_{n+1}-1)f(x_{n})+f(x^{\ast})+(t_{n+1}-1)\langle\nabla f(y_{n}),y_{n}-x_{n}\rangle+\langle\nabla f(y_{n}),y_{n}-x^{\ast}\rangle$
		$\displaystyle-(t_{n+1}-1)\left[\mathcal{A}_{2}(s)+\mathcal{A}_{3}(s;\gamma_{n},L)\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle-\mathcal{A}_{1}(s)(t_{n+1}-1)\\|\nabla f(y_{n})\\|^{2}$
		$\displaystyle-(t_{n+1}-1)\left[\mathcal{A}_{4}(s)+\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}-\mathcal{A}_{1}(s)\\|\nabla f(y_{n})\\|^{2}$
		$\displaystyle-\mathcal{A}_{3}(s;\gamma_{n},L)\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle-\mathcal{A}_{5}(s;\gamma_{n})\\|\nabla f(x_{n})\\|^{2}\,.$

By grouping up the terms, it follows that

	$\displaystyle t_{n+1}f(x_{n+1})-(t_{n+1}-1)f(x_{n})-f(x^{\ast})$	$\displaystyle\leq(t_{n+1}-1)\langle\nabla f(y_{n}),y_{n}-x_{n}\rangle+\langle\nabla f(y_{n}),y_{n}-x^{\ast}\rangle$
		$\displaystyle-\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$
		$\displaystyle-\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}-\mathcal{A}_{1}(s)t_{n+1}\\|\nabla f(y_{n})\\|^{2}\,.$

We rewrite the left hand side of the inequality from above as follows:

	$\displaystyle t_{n+1}f(x_{n+1})-(t_{n+1}-1)f(x_{n})-f(x^{\ast})$	$\displaystyle=t_{n+1}f(x_{n+1})-(t_{n+1}-1)f(x_{n})-(t_{n+1}-t_{n+1}+1)f(x^{\ast})$
		$\displaystyle=t_{n+1}f(x_{n+1})-(t_{n+1}-1)f(x_{n})-t_{n+1}f(x^{\ast})-(1-t_{n+1})f(x^{\ast})$
		$\displaystyle=t_{n+1}\left[f(x_{n+1})-f(x^{\ast})\right]-(t_{n+1}-1)\left[f(x_{n})-f(x^{\ast})\right]$

Then, the above inequality takes the following form:

$\displaystyle t_{n+1}\left[f(x_{n+1})-f(x^{\ast})\right]$	$\displaystyle\leq(t_{n+1}-1)\left[f(x_{n})-f(x^{\ast})\right]+\langle\nabla f(y_{n}),(t_{n+1}-1)(y_{n}-x_{n})+(y_{n}-x^{\ast})\rangle$	(28)
	$\displaystyle-\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$
	$\displaystyle-\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}-\mathcal{A}_{1}(s)t_{n+1}\\|\nabla f(y_{n})\\|^{2}\,.$

Using the fact that $t_{n+1}\geq 0$ , for every $n\geq 0$ , hence for every $n$ satisfying 27, we multiply the inequality 28 by $t_{n+1}$ and it follows that:

	$\displaystyle t_{n+1}^{2}\left[f(x_{n+1})-f(x^{\ast})\right]$	$\displaystyle\leq(t_{n+1}^{2}-t_{n+1})\left[f(x_{n})-f(x^{\ast})\right]+t_{n+1}\langle\nabla f(y_{n}),(t_{n+1}-1)(y_{n}-x_{n})+(y_{n}-x^{\ast})\rangle$
		$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$
		$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}-\mathcal{A}_{1}(s)t_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}\,.$

Since $(t_{n+1}^{2}-t_{n+1})\left[f(x_{n})-f(x^{\ast})\right]=(t_{n+1}^{2}-t_{n+1}-t_{n}^{2})\left[f(x_{n})-f(x^{\ast})\right]+t_{n}^{2}\left[f(x_{n})-f(x^{\ast})\right]$ and using the fact that $t_{n+1}^{2}-t_{n+1}-t_{n}^{2}\leq 0$ and that $f(x_{n})-f(x^{\ast})\geq 0$ (since $x^{\ast}\in\operatorname*{argmin}f$ ) for every $n\geq 0$ and thus for each $n$ satisfying 27, it follows

	$\displaystyle t_{n+1}^{2}\left[f(x_{n+1})-f(x^{\ast})\right]$	$\displaystyle\leq t_{n}^{2}\left[f(x_{n})-f(x^{\ast})\right]+t_{n+1}\langle\nabla f(y_{n}),(t_{n+1}-1)(y_{n}-x_{n})+(y_{n}-x^{\ast})\rangle$
		$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$
		$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}-\mathcal{A}_{1}(s)t_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}\,.$

Using the definition of $\mathcal{E}_{n}$ given in Discrete-Energy, one can easily see that

$\displaystyle\mathcal{E}_{n+1}-\mathcal{E}_{n}$	$\displaystyle\leq t_{n+1}\langle\nabla f(y_{n}),(t_{n+1}-1)(y_{n}-x_{n})+(y_{n}-x^{\ast})\rangle$	(29)
	$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$
	$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}-\mathcal{A}_{1}(s)t_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}$
	$\displaystyle+\dfrac{1}{2s}\\|z_{n+1}\\|^{2}-\dfrac{1}{2s}\\|z_{n}\\|^{2}\,.$

In computing the last term involving $z_{n}$ and $z_{n+1}$ , we shall use the following well-known equality:

\displaystyle\dfrac{1}{2}\|z_{n+1}\|^{2}-\dfrac{1}{2}\|z_{n}\|^{2}=\langle z_{n+1}-z_{n},z_{n+1}\rangle-\dfrac{1}{2}\|z_{n+1}-z_{n}\|^{2}\,.

(30)

Armed with the definition of $z_{n}$ given in Discrete-Energy, we have the following estimations:

	$\displaystyle z_{n+1}-z_{n}$	$\displaystyle=x_{n}-x^{\ast}+t_{n+1}(x_{n+1}-x_{n})+\lambda_{n+1}t_{n+2}\nabla f(x_{n})-x_{n-1}+x^{\ast}-t_{n}(x_{n}-x_{n-1})-\lambda_{n}t_{n+1}\nabla f(x_{n-1})$
		$\displaystyle=(x_{n}-x_{n-1})+t_{n+1}(x_{n+1}-x_{n})-t_{n}(x_{n}-x_{n-1})+\left[\lambda_{n+1}t_{n+2}\nabla f(x_{n})-\lambda_{n}t_{n+1}\nabla f(x_{n-1})\right]$
		$\displaystyle=t_{n+1}(x_{n+1}-x_{n})-(t_{n}-1)(x_{n}-x_{n-1})+\left[\lambda_{n+1}t_{n+2}\nabla f(x_{n})-\lambda_{n}t_{n+1}\nabla f(x_{n-1})\right]\,.$

But, we know that $t_{n}-1=\alpha_{n}t_{n+1}$ for each $n\geq 0$ , consequently for every $n$ satisfying 27, we get that

	$\displaystyle z_{n+1}-z_{n}$	$\displaystyle=t_{n+1}(x_{n+1}-x_{n})-\alpha_{n}t_{n+1}(x_{n}-x_{n-1})+\left[\lambda_{n+1}t_{n+2}\nabla f(x_{n})-\lambda_{n}t_{n+1}\nabla f(x_{n-1})\right]$
		$\displaystyle=t_{n+1}\left[x_{n+1}-(x_{n}+\alpha_{n}(x_{n}-x_{n-1}))\right]+\left[\lambda_{n+1}t_{n+2}\nabla f(x_{n})-\lambda_{n}t_{n+1}\nabla f(x_{n-1})\right]\,.$

From the definition of LT-S-IGAHD, it follows that $x_{n}+\alpha_{n}(x_{n}-x_{n-1})=y_{n}+\lambda_{n}\left[\nabla f(x_{n})-\nabla f(x_{n-1})\right]+\omega_{n}\nabla f(x_{n})\,.$ Hence, it follows that

	$\displaystyle z_{n+1}-z_{n}$	$\displaystyle=t_{n+1}\left[x_{n+1}-y_{n}-\lambda_{n}\left(\nabla f(x_{n})-\nabla f(x_{n-1})\right)-\omega_{n}\nabla f(x_{n})\right]+\left[\lambda_{n+1}t_{n+2}\nabla f(x_{n})-\lambda_{n}t_{n+1}\nabla f(x_{n-1})\right]$
	$\displaystyle z_{n+1}-z_{n}$	$\displaystyle=t_{n+1}(x_{n+1}-y_{n})+\left[-\lambda_{n}t_{n+1}-\omega_{n}t_{n+1}+\lambda_{n+1}t_{n+2}\right]\nabla f(x_{n})\,.$

But, from LT-S-IGAHD, we have that $x_{n+1}-y_{n}=-s\nabla f(y_{n})+\gamma_{n}\nabla f(x_{n})$ . This leads to

\displaystyle z_{n+1}-z_{n}=-st_{n+1}\nabla f(y_{n})+\left[\gamma_{n}t_{n+1}+\lambda_{n+1}t_{n+2}-\lambda_{n}t_{n+1}-\omega_{n}t_{n+1}\right]\nabla f(x_{n})\,.

(31)

It is easy to remark that assumption (ii) is equivalent to $\gamma_{n}=(\lambda_{n}+\omega_{n})-\dfrac{t_{n+2}}{t_{n+1}}\lambda_{n+1}$ for every $n\geq\alpha-1$ , by taking into account that $t_{n+1}=\dfrac{n}{\alpha-1}$ and so $t_{n+2}=\dfrac{n+1}{\alpha-1}$ . Now, multiplying by $t_{n+1}$ , it follows that $\gamma_{n}t_{n+1}=(\lambda_{n}+\omega_{n})t_{n+1}-\lambda_{n+1}t_{n+2}$ . Using this relationship, 27 and 31, for $n\geq\alpha-1$ , we have that

\displaystyle z_{n+1}-z_{n}=-st_{n+1}\nabla f(y_{n})

(32)

At the same time, by using 30 and 32, we get

\displaystyle\dfrac{1}{2s}\|z_{n+1}\|^{2}-\dfrac{1}{2s}\|z_{n}\|^{2}=-\langle t_{n+1}\nabla f(y_{n}),(x_{n}-x^{\ast})+t_{n+1}(x_{n+1}-x_{n})+\lambda_{n+1}t_{n+2}\nabla f(x_{n})\rangle-\dfrac{s}{2}t_{n+1}^{2}\|\nabla f(y_{n})\|^{2}

This means that

	$\displaystyle\dfrac{1}{2s}\\|z_{n+1}\\|^{2}-\dfrac{1}{2s}\\|z_{n}\\|^{2}$	$\displaystyle=-\lambda_{n+1}t_{n+1}t_{n+2}\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle-\dfrac{s}{2}t_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}$		(33)
		$\displaystyle-t_{n+1}\langle\nabla f(y_{n}),(x_{n}-x^{\ast})+t_{n+1}(x_{n+1}-x_{n})\rangle\,.$

Combining 29 with 33, we obtain

	$\displaystyle\mathcal{E}_{n+1}-\mathcal{E}_{n}$	$\displaystyle\leq t_{n+1}\langle\nabla f(y_{n}),(t_{n+1}-1)(y_{n}-x_{n})+(y_{n}-x^{\ast})\rangle$
		$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$
		$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}-\mathcal{A}_{1}(s)t_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}$
		$\displaystyle-\lambda_{n+1}t_{n+1}t_{n+2}\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle-\dfrac{s}{2}t_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}-t_{n+1}\langle\nabla f(y_{n}),(x_{n}-x^{\ast})+t_{n+1}(x_{n+1}-x_{n})\rangle\,.$

On the other hand, we make the following estimations:

	$\displaystyle\langle\nabla f(y_{n}),(t_{n+1}-1)(y_{n}-x_{n})+(y_{n}-x^{\ast})\rangle-\langle\nabla f(y_{n}),(x_{n}-x^{\ast})+t_{n+1}(x_{n+1}-x_{n})\rangle=$
	$\displaystyle\langle\nabla f(y_{n}),(t_{n+1}-1)(y_{n}-x_{n})+(y_{n}-x^{\ast})-(x_{n}-x^{\ast})-t_{n+1}(x_{n+1}-x_{n})\rangle=$
	$\displaystyle\langle\nabla f(y_{n}),t_{n+1}(y_{n}-x_{n+1})\rangle\,.$

From LT-S-IGAHD, we know that $y_{n}-x_{n+1}=s\nabla f(y_{n})-\gamma_{n}\nabla f(x_{n})$ and this leads to

	$\displaystyle\langle\nabla f(y_{n}),(t_{n+1}-1)(y_{n}-x_{n})+(y_{n}-x^{\ast})\rangle-\langle\nabla f(y_{n}),(x_{n}-x^{\ast})+t_{n+1}(x_{n+1}-x_{n})\rangle=$
	$\displaystyle\langle\nabla f(y_{n}),st_{n+1}\nabla f(y_{n})-\gamma_{n}t_{n+1}\nabla f(x_{n})\rangle=st_{n+1}\\|\nabla f(y_{n})\\|^{2}-\gamma_{n}t_{n+1}\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle\,.$

So, it follows that

	$\displaystyle\mathcal{E}_{n+1}-\mathcal{E}_{n}$	$\displaystyle\leq-\mathcal{A}_{1}(s)t_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}+st_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}-\dfrac{s}{2}t_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}$
		$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$
		$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}$
		$\displaystyle-\lambda_{n+1}t_{n+1}t_{n+2}\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle-\gamma_{n}t_{n+1}^{2}\langle f(y_{n}),f(x_{n})\rangle\,.$

Using the fact that $\mathcal{A}_{1}(s)=s$ , one obtains

	$\displaystyle\mathcal{E}_{n+1}-\mathcal{E}_{n}$	$\displaystyle\leq-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)+\lambda_{n+1}t_{n+2}+\gamma_{n}t_{n+1}\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$
		$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}-\dfrac{s}{2}t_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}\,.$

Recalling the fact that we have already pointed out that assumption (ii) is equivalent to $\lambda_{n+1}t_{n+2}+\gamma_{n}t_{n+1}=(\lambda_{n}+\omega_{n})t_{n+1}$ for every $n\geq\alpha-1$ , we find that

	$\displaystyle\mathcal{E}_{n+1}-\mathcal{E}_{n}$	$\displaystyle\leq-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)+(\lambda_{n}+\omega_{n})t_{n+1}\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$
		$\displaystyle-t_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}-\dfrac{s}{2}t_{n+1}^{2}\\|\nabla f(y_{n})\\|^{2}\,.$

So, the difference in the discrete energy functional can be written as

\displaystyle\mathcal{E}_{n+1}-\mathcal{E}_{n}\leq-t_{n+1}\tilde{F}_{n}\,,

(34)

where

	$\displaystyle\tilde{F}_{n}$	$\displaystyle:=\dfrac{s}{2}t_{n+1}\\|\nabla f(y_{n})\\|^{2}+\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)+(\lambda_{n}+\omega_{n})t_{n+1}\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle$		(35)
		$\displaystyle+\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}\,.$

( Step II ) The second step of our proof consists in applying Lemma 17 and Lemma 18 to two quadratic inequalities with nonconstant coefficients and show the conditions on the index $n$ such that they are satisfied.
For simplicity, let’s denote $Y:=\nabla f(y_{n})$ and $X:=\nabla f(x_{n})$ . Hence, we obtain that

	$\displaystyle\tilde{F}_{n}=\dfrac{s}{2}t_{n+1}\\|Y\\|^{2}$	$\displaystyle+\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)+(\lambda_{n}+\omega_{n})t_{n+1}\right]\langle Y,X\rangle$
		$\displaystyle+\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|X\\|^{2}\,.$

Using Cauchy-Schwarz inequality, it follows that

	$\displaystyle\tilde{F}_{n}\geq\dfrac{s}{2}t_{n+1}\\|Y\\|^{2}$	$\displaystyle-\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)+(\lambda_{n}+\omega_{n})t_{n+1}\right]\\|Y\\|\\|X\\|$		(36)
		$\displaystyle+\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|X\\|^{2}:=F_{n}\,.$

The next step of the present proof is to apply Lemma 17 and to show that for $n$ large enough, $F_{n}\geq 0$ . Using Lemma 17, we observe that the coefficient of $\|Y\|^{2}$ is positive, namely $\dfrac{s}{2}t_{n+1}>0$ for $n\geq\alpha-1$ . Now, the only thing left for us is to show that

\displaystyle\left[(1-t_{n+1})\mathcal{A}_{2}(s)-t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)-(\lambda_{n}+\omega_{n})t_{n+1}\right]^{2}\leq 2st_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\,.

(37)

This means that

	$\displaystyle(1-t_{n+1})^{2}\mathcal{A}_{2}^{2}(s)+t_{n+1}^{2}\mathcal{A}_{3}^{2}(s;\gamma_{n},L)-2t_{n+1}(1-t_{n+1})\mathcal{A}_{2}(s)\mathcal{A}_{3}(s;\gamma_{n},L)+(\lambda_{n}+\omega_{n})^{2}t_{n+1}^{2}$
	$\displaystyle-2(\lambda_{n}+\omega_{n})t_{n+1}\left[(1-t_{n+1})\mathcal{A}_{2}(s)-t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)\right]\leq 2st_{n+1}\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right],.$

Expanding the terms, it follows that

	$\displaystyle t_{n+1}^{2}\mathcal{A}_{2}^{2}(s)+\mathcal{A}_{2}^{2}(s)-2t_{n+1}\mathcal{A}_{2}^{2}(s)+t_{n+1}^{2}\mathcal{A}_{3}^{2}(s;\gamma_{n},L)-2t_{n+1}\mathcal{A}_{2}(s)\mathcal{A}_{3}(s;\gamma_{n},L)+2t_{n+1}^{2}\mathcal{A}_{2}(s)\mathcal{A}_{3}(s;\gamma_{n},L)$
	$\displaystyle+(\lambda_{n}+\omega_{n})^{2}t_{n+1}^{2}-2\mathcal{A}_{2}(s)(\lambda_{n}+\omega_{n})t_{n+1}+2\mathcal{A}_{2}(s)(\lambda_{n}+\omega_{n})t_{n+1}^{2}+2(\lambda_{n}+\omega_{n})t_{n+1}^{2}\mathcal{A}_{3}(s;\gamma_{n},L)$
	$\displaystyle\leq 2s\mathcal{A}_{4}(s)t_{n+1}^{2}-2s\mathcal{A}_{4}(s)t_{n+1}+2st_{n+1}^{2}\mathcal{A}_{5}(s;\gamma_{n})\,.$

Grouping by $t_{n+1}$ and $t_{n+1}^{2}$ implies that

	$\displaystyle t_{n+1}^{2}[\mathcal{A}_{2}^{2}(s)+\mathcal{A}_{3}^{2}(s;\gamma_{n},L)+2\mathcal{A}_{2}(s)\mathcal{A}_{3}(s;\gamma_{n},L)+(\lambda_{n}+\omega_{n})^{2}+2\mathcal{A}_{2}(s)(\lambda_{n}+\omega_{n})+2(\lambda_{n}+\omega_{n})\mathcal{A}_{3}(s;\gamma_{n},L)$
	$\displaystyle-2s\mathcal{A}_{4}(s)-2s\mathcal{A}_{5}(s;\gamma_{n})]+t_{n+1}[-2\mathcal{A}_{2}^{2}(s)-2\mathcal{A}_{2}(s)\mathcal{A}_{3}(s;\gamma_{n},L)-2\mathcal{A}_{2}(s)(\lambda_{n}+\omega_{n})+2s\mathcal{A}_{4}(s)]+\mathcal{A}_{2}^{2}(s)\leq 0\,.$

Finally, one has that

		$\displaystyle t_{n+1}^{2}\left[(\mathcal{A}_{2}(s)+\mathcal{A}_{3}(s;\gamma_{n},L))^{2}-2s(\mathcal{A}_{4}(s)+\mathcal{A}_{5}(s;\gamma_{n}))+(\lambda_{n}+\omega_{n})^{2}+2(\lambda_{n}+\omega_{n})(\mathcal{A}_{2}(s)+\mathcal{A}_{3}(s;\gamma_{n},L))\right]+$		(38)
		$\displaystyle t_{n+1}\left[-2\mathcal{A}_{2}^{2}(s)-2\mathcal{A}_{2}(s)\mathcal{A}_{3}(s;\gamma_{n},L)-2(\lambda_{n}+\omega_{n})\mathcal{A}_{2}(s)+2s\mathcal{A}_{4}(s)\right]+\mathcal{A}_{2}^{2}(s)\leq 0\,.$

We rewrite 38 as

\displaystyle G_{n}t_{n+1}^{2}-H_{n}t_{n+1}-I_{n}\geq 0\,,

(39)

where the coefficients are

\displaystyle\begin{cases}G_{n}&:=-(\mathcal{A}_{2}(s)+\mathcal{A}_{3}(s;\gamma_{n},L))^{2}+2s(\mathcal{A}_{4}(s)+\mathcal{A}_{5}(s;\gamma_{n}))-(\lambda_{n}+\omega_{n})^{2}-2(\lambda_{n}+\omega_{n})(\mathcal{A}_{2}(s)+\mathcal{A}_{3}(s;\gamma_{n},L))\\ H_{n}&:=-2\mathcal{A}_{2}^{2}(s)-2\mathcal{A}_{2}(s)\mathcal{A}_{3}(s;\gamma_{n},L)-2(\lambda_{n}+\omega_{n})\mathcal{A}_{2}(s)+2s\mathcal{A}_{4}(s)\\ I_{n}&:=\mathcal{A}_{2}^{2}(s)\end{cases}

Since $I_{n}=s^{2}>0$ , in order to use Lemma 18 and taking into account Remark 19, we only need to show that for $n$ large enough $G_{n}>0$ , hence $H_{n}^{2}+4G_{n}I_{n}\geq 0$ . Now, this means that 39 is equivalent to the fact that $(-G_{n})t_{n+1}^{2}+H_{n}t_{n+1}+I_{n}\leq 0$ . Based upon 38, we compute $-G_{n}$ as follows:

	$\displaystyle-G_{n}$	$\displaystyle=\left[s+\gamma_{n}(1-Ls)\right]^{2}-s^{2}+\gamma_{n}^{2}+(\lambda_{n}+\omega_{n})^{2}+2(\lambda_{n}+\omega_{n})(-s-\gamma_{n}(1-Ls))$
		$\displaystyle=\left[s+\gamma_{n}(1-Ls)\right]^{2}+(\gamma_{n}^{2}-s^{2})+(\lambda_{n}+\omega_{n})^{2}-2(\lambda_{n}+\omega_{n})(s+\gamma_{n}(1-Ls))$
		$\displaystyle=\left[s+\gamma_{n}(1-Ls)\right]^{2}+(\gamma_{n}^{2}-s^{2})+(\lambda_{n}+\omega_{n})^{2}-(\lambda_{n}+\omega_{n})(s+\gamma_{n}(1-Ls))-(\lambda_{n}+\omega_{n})(s+\gamma_{n}(1-Ls))$
		$\displaystyle=\left[s+\gamma_{n}(1-Ls)\right]\cdot\left[s+\gamma_{n}(1-Ls)-(\lambda_{n}+\omega_{n})\right]+(\gamma_{n}^{2}-s^{2})+(\lambda_{n}+\omega_{n})\left[(\lambda_{n}+\omega_{n})-(s+\gamma_{n}(1-Ls))\right]$
		$\displaystyle=(\gamma_{n}^{2}-s^{2})+\left[s+\gamma_{n}(1-Ls)-(\lambda_{n}+\omega_{n})\right]\cdot\left[s+\gamma_{n}(1-Ls)\right]-(\lambda_{n}+\omega_{n})\left[s+\gamma_{n}(1-Ls)-(\lambda_{n}+\omega_{n})\right]$
		$\displaystyle=(\gamma_{n}^{2}-s^{2})+\left[(s+\gamma_{n}(1-Ls))-(\lambda_{n}+\omega_{n})\right]\cdot\left[(s+\gamma_{n}(1-Ls))-(\lambda_{n}+\omega_{n})\right]$
		$\displaystyle=(\gamma_{n}^{2}-s^{2})+\left[(s+\gamma_{n}(1-Ls))-(\lambda_{n}+\omega_{n})\right]^{2}\,.$

Considering the above computations for $G_{n}$ and the form of $\lambda_{n}+\omega_{n}$ from assumption (ii), we obtain, after some basic evaluations, exactly assumption (i), where it is satisfied for each $n>max\{\alpha-1,N^{\prime}\}$ .
( Step III ) In the last step of the present proof we recollect the results from above and apply Lemma 17 and Lemma 18 and actually show the $O\left(\dfrac{1}{n^{2}}\right)$ asymptotic rate of converge. At the same time, we actually determine the index $N\in\mathbb{N}$ , such that for $n\geq N$ , the rate of convergence in the objective function is obtained.
From now on, we consider only when $n>max\{N_{1},N_{2},N^{\prime}\}$ , where $N_{1}:=\alpha-1$ , $N_{2}:=(\alpha-1)\dfrac{H_{n}+\sqrt{H_{n}^{2}+4G_{n}I_{n}}}{2G_{n}}$ and $N^{\prime}$ the index from the hypothesis of the theorem. Now, we know from the theorem’s hypothesis that $(\gamma_{n})_{n\in\mathbb{N}}$ and $(\lambda_{n})_{n\in\mathbb{N}}$ are convergent. By using assumption (ii), it follows that $(\omega_{n})_{n\in\mathbb{N}}$ is also convergent. This implies that also $(H_{n})_{n\in\mathbb{N}}$ and $(G_{n})_{n\in\mathbb{N}}$ are convergent. Since these two sequences are convergent, they are bounded. Hence, the sequence $\left((\alpha-1)\dfrac{H_{n}+\sqrt{H_{n}^{2}+4G_{n}I_{n}}}{2G_{n}}\right)_{n\in\mathbb{N}}$ is also bounded. So, let $M>0$ such that $(\alpha-1)\dfrac{H_{n}+\sqrt{H_{n}^{2}+4G_{n}I_{n}}}{2G_{n}}<M$ . Then, we can take $n\geq M>N_{2}$ (in this sense $n>N_{2}$ is well defined).
Consider 39, with the nonconstant coefficients $G_{n}$ , $H_{n}$ and $I_{n}$ . We apply Lemma 18 with $x=t_{n+1}$ , $a(x)=a(t_{n+1})=G_{n}$ , $b(x)=b(t_{n+1})=-H_{n}$ and $c(x)=c(t_{n+1})=-I_{n}$ . Since $n>N_{2}$ , then we consider find an infinity of points $x=t_{n+1}\geq\dfrac{H_{n}+\sqrt{H_{n}^{2}+4G_{n}I_{n}}}{2G_{n}}$ that obviously satisfy $c(t_{n+1})=-s^{2}<0$ , since $s>0$ . Also, the condition $b^{2}(t_{n+1})-4a(t_{n+1})c(t_{n+1})\geq 0$ is equivalent to $H_{n}^{2}+4G_{n}I_{n}\geq 0$ . Because $n>max\{N_{1},N_{2},N^{\prime}\}$ , assumption (i) is in fact $G_{n}>0$ , i.e. $a(t_{n+1})>0$ (this stands true for $n\geq\max\{\alpha-1,N^{\prime}\}$ ). This means that indeed we have $H_{n}^{2}+4G_{n}I_{n}\geq 0$ , for each $n>max\{N_{1},N_{2},N^{\prime}\}$ .
After these computations, now we can consider the quadratic inequality with nonconstant coefficients $F_{n}\geq 0$ , with $F_{n}$ given in 36 and where $n>\max\{N_{1},N_{2},N^{\prime}\}$ . We apply Lemma 17 with $x=\|\nabla f(y_{n})\|$ , $a(x)=a(\|\nabla f(y_{n})\|)=\dfrac{s}{2}t_{n+1}$ , $b(x)=b(\|\nabla f(y_{n})\|)=\left[(1-t_{n+1})\mathcal{A}_{2}(s)-t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)-(\lambda_{n}+\omega_{n})t_{n+1}\right]\|\nabla f(x_{n})\|$ and $c(x)=c(\|\nabla f(y_{n})\|)=\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\|\nabla f(x_{n})\|^{2}$ . So, in this case we consider $x=t_{n+1}$ with $n>\max\{N_{1},N_{2},N^{\prime}\}$ . The condition $a(\|\nabla f(y_{n})\|)>0$ is obviously satisfied for $n>\max\{N_{1},N_{2},N^{\prime}\}$ . Furthermore, the condition is equivalent to 37. Moreover, this is in fact 39, that we have shown that is valid for $n>max\{N_{1},N_{2},N^{\prime}\}$ . We conclude that $F_{n}\geq 0$ for every $\nabla f(y_{n})$ with $n>max\{N_{1},N_{2},N^{\prime}\}$ .
By 34, we find that $\mathcal{E}_{n+1}-\mathcal{E}_{n}\leq 0$ , so $\mathcal{E}_{n}\leq\mathcal{E}_{N}$ , where $N:=\max\{N_{1},N_{2},N^{\prime}\}$ . By the fact that $t_{n}=\dfrac{n-1}{\alpha-1}$ and utilizing the definition Discrete-Energy, we have that $t_{n}^{2}\left[f(x_{n})-f(x^{\ast})\right]\leq\mathcal{E}_{N}$ . Finally, we obtain that $f(x_{n})-f(x^{\ast})\leq\dfrac{\mathcal{E}_{N}(\alpha-1)^{2}}{(n-1)^{2}}$ , for $n>N$ and the proof is over. $\blacksquare$

Remark 21.

Take the case when the sequence $(\lambda_{n})_{n\in\mathbb{N}}$ is constant, with $\lambda_{n}=\beta\sqrt{s}$ . Then, the value of $z_{n}$ given in Discrete-Energy is $z_{n}=(x_{n-1}-x^{\ast})+t_{n}(x_{n}-x_{n-1})+\beta\sqrt{s}t_{n+1}\nabla f(x_{n-1})$ . This discrete value differs significantly from that given in [4], where $z_{n}=(x_{n-1}-x^{\ast})+t_{n}\left[(x_{n}-x_{n-1})+\beta\sqrt{s}\nabla f(x_{n-1})\right]$ , by the fact that we have used $\nabla f(x_{n})$ and not $\nabla f(x_{n-1})$ in LT-S-IGAHD (see Remark 5 and Subsection 2.2). Last of all, we point out that even though Discrete-Energy is similar to that used in [4], there are subtle differences consisting in the fact that our choice of the discrete functional is optimal, since in $z_{n+1}$ we got rid of $\nabla f(x_{n-1})$ . This means that every other types of discrete Lyapunov functionals that can be constructed, contain $\nabla f(x_{n-1})$ , that can’t be eliminated by the fact that they have a nonconstant term in front of them.

Remark 22.

We observe that assumption (ii) of Theorem 20 reduces to the fact that $\gamma^{2}_{n}<s^{2}$ , namely $\gamma_{n}\in(-s,s)$ for each $n>max\{\alpha-1,N^{\prime}\}$ . Furthermore, following the last step of the previous proof and considering Lemma 17 and Remark 19, the following sequence of equivalences stand true:

	$\displaystyle\gamma_{n}^{2}<s^{2}\left(1-\dfrac{\alpha-1}{n}\right)$	$\displaystyle\Longleftrightarrow s-\dfrac{\gamma_{n}^{2}}{s}>\dfrac{s(\alpha-1)}{n}$
		$\displaystyle\Longleftrightarrow\dfrac{n}{\alpha-1}\left(\dfrac{s}{2}-\dfrac{\gamma_{n}^{2}}{2s}\right)>\dfrac{s}{2}\Longleftrightarrow t_{n+1}\left[\mathcal{A}_{4}(s)+\mathcal{A}_{5}(s;\gamma_{n})\right]>\mathcal{A}_{4}(s)\,,$

So, in fact we have a stronger bound for $\gamma_{n}$ , namely $-s^{2}\left(1-\dfrac{\alpha-1}{n}\right)<\gamma_{n}<s^{2}\left(1-\dfrac{\alpha-1}{n}\right)$ , and this is valid for all $n>max\{\alpha-1,N^{\prime}\}$ .

Remark 23.

In Theorem 20 we have considered the strict inequality belonging to assumption (ii) to be valid for each $n>\max\{\alpha-1,N^{\prime}\}$ . But, if we look closer at the third step of the proof from the same theorem, we observe that we could have taken $n>\max\{\alpha-1,N_{2},N^{\prime}\}$ , where $N_{2}$ was defined through $H_{n}$ , $G_{n}$ and $I_{n}$ . Hence, in this way, we obtain a more relaxed bound on the restriction on the index $n$ .
On the other hand, we observe that the assumptions that the sequences $(\gamma_{n})_{n\in\mathbb{N}}$ and $(\lambda_{n})_{n\in\mathbb{N}}$ are bounded, so also $(\omega_{n})_{n\in\mathbb{N}}$ , can be linked to our reasoning from Subsection 2.4, where we have supposed the convergence of the sequences in our heuristic analysis.

We end this section by giving some examples of algorithms of tpye LT-S-IGAHD that satisfy the assumptions from Theorem 20 and we also show that indeed exists an index $N^{\prime}$ for each of them. In our first example, we consider a general class of Hessian driven optimization algorithms with asymptotic vanishing damping, that include three parameters $a$ , $\mu$ and $b$ that will play a crucial role in the numerical behavior, in the sense that, from a scientific computing point of view, we have a much larger freedom and flexibility into the hypertuning of the parameters.

Example 24.

We consider the choice of the coefficients as follows:
$\gamma_{n}=s\sqrt{\dfrac{\alpha-1}{n+a}}$ , $\lambda_{n+1}=s\dfrac{n}{n+1}+\mu\dfrac{n}{(n+1)(n+b)}$ and $\omega_{n}=s\sqrt{\dfrac{\alpha-1}{n+a}}+\dfrac{s}{n}+\mu\left[\dfrac{1}{n+b}-\dfrac{n-1}{n(n+b-1)}\right]$ , where $a\geq 0$ , $b\geq 0$ and $\mu\geq 0$ . Now, we validate our algorithm as follows:
With the choices of the sequences $(\gamma_{n})_{n\in\mathbb{N}}$ and $(\lambda_{n})_{n\in\mathbb{N}}$ , we obtain the following chain of equivalences:

	$\displaystyle\omega_{n}$	$\displaystyle=s\sqrt{\dfrac{\alpha-1}{n+a}}+\dfrac{s}{n}+\mu\left[\dfrac{1}{n+b}-\dfrac{n-1}{n(n+b-1)}\right]$
	$\displaystyle\omega_{n}$	$\displaystyle=s\sqrt{\dfrac{\alpha-1}{n+a}}-s\left(\dfrac{n-1}{n}-1\right)+\mu\left[\dfrac{1}{n+b}-\dfrac{n-1}{n(n+b-1)}\right]$
	$\displaystyle\omega_{n}$	$\displaystyle=s\sqrt{\dfrac{\alpha-1}{n+a}}-s\dfrac{n-1}{n}-\mu\dfrac{n-1}{n(n+b-1)}+s+\dfrac{\mu}{n+b}$
	$\displaystyle\omega_{n}$	$\displaystyle=\gamma_{n}-\lambda_{n}+\dfrac{n+1}{n}\lambda_{n+1}$

so assumption (ii) of Theorem 20 is satisfied for each $n\geq 1$ .
We must show assumption (i) of Theorem 20. For this, we have the following equivalences:

	$\displaystyle\left[s(1-\gamma_{n}L)-\dfrac{n+1}{n}\lambda_{n+1}\right]^{2}<s^{2}-\gamma_{n}^{2}$	$\displaystyle\Longleftrightarrow\left[s\gamma_{n}L+\dfrac{\mu}{n+b}\right]^{2}<s^{2}-\gamma_{n}^{2}$
		$\displaystyle\Longleftrightarrow s^{2}\gamma_{n}^{2}L^{2}+2s\mu L\gamma_{n}\dfrac{1}{n+b}+\dfrac{\mu^{2}}{(n+b)^{2}}<s^{2}-\gamma_{n}^{2}$
		$\displaystyle\Longleftrightarrow s^{4}\dfrac{\alpha-1}{n+a}L^{2}+2s^{2}\mu L\dfrac{\sqrt{\alpha-1}}{(n+b)\sqrt{n+a}}+\dfrac{\mu^{2}}{(n+b)^{2}}<s^{2}-s^{2}\dfrac{\alpha-1}{n+a}$

Multiplying by $1/s^{2}$ , we get that

\displaystyle s^{2}\dfrac{\alpha-1}{n+a}L^{2}+2\mu L\dfrac{\sqrt{\alpha-1}}{(n+b)\sqrt{n+a}}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+b)^{2}}+\dfrac{\alpha-1}{n+a}<1\,.

(40)

Let $\sqrt{n+a}<n+b$ . So $n+a-n^{2}-2bn-b^{2}<0$ , which means that $n^{2}+(2b-1)n+(b^{2}-a)>0$ , with $\Delta=(2b-1)^{2}-4(b^{2}-a)$ . We consider the following cases :

•

If $\Delta\geq 0$ , namely $b-a\leq\dfrac{1}{4}$ then we can take $n>\dfrac{1-2b+\sqrt{4(a-b)+1}}{2}$ , so the inequality is satisfied.
•

If $\Delta<0$ , i.e. $b-a>\dfrac{1}{4}$ , then we can take $n\geq 1$ .

Since

s^{2}\dfrac{\alpha-1}{n+a}L^{2}+2\mu L\dfrac{\sqrt{\alpha-1}}{(n+b)\sqrt{n+a}}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+b)^{2}}+\dfrac{\alpha-1}{n+a}<s^{2}\dfrac{\alpha-1}{n+a}L^{2}+2\mu L\dfrac{\sqrt{\alpha-1}}{n+a}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{n+a}+\dfrac{\alpha-1}{n+a}\,,

and taking into account 40, then it is enough to show that:

	$\displaystyle s^{2}\dfrac{\alpha-1}{n+a}L^{2}+2\mu L\dfrac{\sqrt{\alpha-1}}{n+a}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{n+a}+\dfrac{\alpha-1}{n+a}<1$	$\displaystyle\Longleftrightarrow\dfrac{1}{n+a}\left[(\alpha-1)(s^{2}L^{2}+1)+2\mu L\sqrt{\alpha-1}+\left(\dfrac{\mu}{s}\right)^{2}\right]<1$
		$\displaystyle\Longleftrightarrow n>(\alpha-1)(s^{2}L^{2}+1)+2\mu L\sqrt{\alpha-1}+\left(\dfrac{\mu}{s}\right)^{2}-a\,.$

Finally, if $b-a>\dfrac{1}{4}$ then we can consider

N^{\prime}:=(\alpha-1)(s^{2}L^{2}+1)+2\mu L\sqrt{\alpha-1}+\left(\dfrac{\mu}{s}\right)^{2}-a\,,

and if $b-a\leq\dfrac{1}{4}$ , then we can take

N^{\prime}:=\max\Big{\{}\dfrac{1-2b+\sqrt{4(a-b)+1}}{2},(\alpha-1)(s^{2}L^{2}+1)+2\mu L\sqrt{\alpha-1}+\left(\dfrac{\mu}{s}\right)^{2}-a\Big{\}}.

The second example is also based upon three coefficient $\beta$ , $b$ and $\mu$ . It is worth noticing that for the particular case when $\mu=0$ we obtain the IGAHD-type algorithm that differs from that given in [4] (see Remark 5 and Subsection 2.2).

Example 25.

We consider the following choice of the coefficients:
$\gamma_{n}=0$ , $\lambda_{n+1}=\beta\sqrt{s}+\dfrac{\mu}{n+b}$ and $\omega_{n}=\dfrac{\beta\sqrt{s}}{n}+\mu\left[\dfrac{n+1}{n(n+b)}-\dfrac{1}{n+b-1}\right]$ , where $0<\beta<2\sqrt{s}$ , $b>0$ and $\mu\geq 0$ . Now, we validate our algorithm as follows:
With the choices of the sequences $(\gamma_{n})_{n\in\mathbb{N}}$ and $(\lambda_{n})_{n\in\mathbb{N}}$ , we obtain the following chain of equivalences:

	$\displaystyle\omega_{n}$	$\displaystyle=\dfrac{\beta\sqrt{s}}{n}+\mu\left[\dfrac{n+1}{n(n+b)}-\dfrac{1}{n+b-1}\right]$
	$\displaystyle\omega_{n}$	$\displaystyle=\dfrac{n+1}{n}\beta\sqrt{s}+\mu\dfrac{n+1}{n(n+b)}-\beta\sqrt{s}-\dfrac{\mu}{n+b-1}$
	$\displaystyle\omega_{n}$	$\displaystyle=\dfrac{n+1}{n}\lambda_{n+1}-\lambda_{n}$
	$\displaystyle\omega_{n}$	$\displaystyle=\gamma_{n}-\lambda_{n}+\dfrac{n+1}{n}\lambda_{n+1}$

We must show assumption (i) of Theorem 20. For this, we have the following equivalences:

	$\displaystyle\left[s(1-\gamma_{n}L)-\dfrac{n+1}{n}\lambda_{n+1}\right]^{2}<s^{2}-\gamma_{n}^{2}$	$\displaystyle\Longleftrightarrow\left[s-\dfrac{n+1}{n}\lambda_{n+1}\right]^{2}<s^{2}$
		$\displaystyle\Longleftrightarrow\left[s-\dfrac{n+1}{n}\beta\sqrt{s}-\mu\dfrac{n+1}{n(n+b)}\right]^{2}<s^{2}$
		$\displaystyle\Longleftrightarrow-\left[\dfrac{n+1}{n}\beta\sqrt{s}+\mu\dfrac{n+1}{n(n+b)}\right]\left[2s-\dfrac{n+1}{n}\beta\sqrt{s}-\mu\dfrac{n+1}{n(n+b)}\right]<0$

By the fact that $\dfrac{n+1}{n}\beta\sqrt{s}+\mu\dfrac{n+1}{n(n+b)}$ is positive and since $\beta>0$ , this term can’t be $0$ . So, the inequality from above is satisfied if and only if $\dfrac{n+1}{n}\beta\sqrt{s}+\mu\dfrac{n+1}{n(n+b)}<2s$ . Multiplying by $\dfrac{n}{\sqrt{s}(n+1)}$ , it follows that $\beta+\dfrac{\mu}{\sqrt{s}(n+b)}<2\sqrt{s}\dfrac{n}{n+1}$ . Hence

\displaystyle\beta+\dfrac{\mu}{\sqrt{s}(n+b)}<2\sqrt{s}\dfrac{n}{n+1}\,.

(41)

Since $\beta<2\sqrt{s}$ , there exists $\tau>0$ such that $2\sqrt{s}-\beta>\tau$ . We consider $\epsilon=\dfrac{\tau}{\dfrac{\mu}{\sqrt{s}}+2\sqrt{s}}>0$ . From the Archimedean property, we know that for $\epsilon>0$ there exists $N\in\mathbb{N}$ such that $\epsilon>\dfrac{1}{N}$ , thus the strict inequality $2\sqrt{s}-\beta>\left(\dfrac{\mu}{\sqrt{s}}+2\sqrt{s}\right)\epsilon$ implies that $2\sqrt{s}-\beta>\left(\dfrac{\mu}{\sqrt{s}}+2\sqrt{s}\right)\dfrac{1}{N}$ . For $n+1\geq N$ and $n+b\geq N$ which are equivalent to $\dfrac{1}{N}\geq\dfrac{1}{n+1}$ and $\dfrac{1}{N}\geq\dfrac{1}{n+b}$ , we obtain that

\left(\dfrac{\mu}{\sqrt{s}}+2\sqrt{s}\right)\dfrac{1}{N}\geq\dfrac{\mu}{\sqrt{s}}\dfrac{1}{n+b}+2\sqrt{s}\dfrac{1}{n+1}\Rightarrow 2\sqrt{s}-\beta>\dfrac{\mu}{\sqrt{s}}\dfrac{1}{n+b}+2\sqrt{s}\dfrac{1}{n+1}.

Second to last, we can take $N^{\prime}:=max\{N-1,N-b\}$ and we obtain that 41 is satisfied for every $n\geq N^{\prime}$ . Finally, one can find in an explicit way the index $N^{\prime}$ . In order to do this, observe that 41 leads to a quadratic inequality, namely:

(2\sqrt{s}-\beta)n^{2}+\left[2\sqrt{s}b-\beta(b+1)-\dfrac{\mu}{\sqrt{s}}\right]n-\left(\beta b+\dfrac{\mu}{\sqrt{s}}\right)>0.

In this case, one can see that $N^{\prime}=\dfrac{\left[\beta(b+1)+\dfrac{\mu}{\sqrt{s}}-2\sqrt{s}b\right]+\sqrt{\Delta}}{2(2\sqrt{s}-\beta)}$ , where $\Delta:=\left(2b\sqrt{s}-\beta(b+1)-\dfrac{\mu}{\sqrt{s}}\right)^{2}+4(2\sqrt{s}-\beta)\left(\beta b+\dfrac{\mu}{\sqrt{s}}\right)\geq 0\,.$

Our last example consists in choosing $\gamma_{n}$ to take negative values, but $\lambda_{n}$ to have only positive values. Thus, $\omega_{n}$ is determined through assumption (ii) of Theorem 20.

Example 26.

We consider the choice of the coefficients as follows:
$\gamma_{n}=-\dfrac{s}{n+a}$ , $\lambda_{n+1}=s\dfrac{n}{n+1}+\mu\dfrac{n}{(n+1)(n+b)}$ and $\omega_{n}=-\dfrac{s}{n+a}+\dfrac{s}{n}+\mu\left[\dfrac{1}{n+b}-\dfrac{n-1}{n(n+b-1)}\right]$ , where $a\geq 0$ , $b\geq 0$ and $\mu\geq 0$ . Now, we validate our algorithm as follows:
With the choices of the sequences $(\gamma_{n})_{n\in\mathbb{N}}$ and $(\lambda_{n})_{n\in\mathbb{N}}$ , we obtain the following chain of equivalences:

	$\displaystyle\omega_{n}$	$\displaystyle=-\dfrac{s}{n+a}+\dfrac{s}{n}+\mu\left[\dfrac{1}{n+b}-\dfrac{n-1}{n(n+b-1)}\right]$
	$\displaystyle\omega_{n}$	$\displaystyle=-\dfrac{s}{n+a}-s\left(\dfrac{n-1}{n}-1\right)+\mu\left[\dfrac{1}{n+b}-\dfrac{n-1}{n(n+b-1)}\right]$
	$\displaystyle\omega_{n}$	$\displaystyle=-\dfrac{s}{n+a}-s\dfrac{n-1}{n}-\mu\dfrac{n-1}{n(n+b-1)}+s+\dfrac{\mu}{n+b}$
	$\displaystyle\omega_{n}$	$\displaystyle=\gamma_{n}-\lambda_{n}+\dfrac{n+1}{n}\lambda_{n+1}$

so assumption (ii) of Theorem 20 is satisfied for each $n\geq 1$ .
Now, regarding assumption (i), we consider the following equivalences:

	$\displaystyle\left[s(1-\gamma_{n}L)-\dfrac{n+1}{n}\lambda_{n+1}\right]^{2}<s^{2}-\gamma_{n}^{2}$	$\displaystyle\Longleftrightarrow\left[s\gamma_{n}L+\dfrac{\mu}{n+b}\right]^{2}<s^{2}-\gamma_{n}^{2}$
		$\displaystyle\Longleftrightarrow s^{2}\gamma_{n}^{2}L^{2}+2s\mu L\gamma_{n}\dfrac{1}{n+b}+\dfrac{\mu^{2}}{(n+b)^{2}}<s^{2}-\gamma_{n}^{2}$
		$\displaystyle\Longleftrightarrow s^{4}L^{2}\dfrac{1}{(n+a)^{2}}-2s^{2}\mu L\dfrac{1}{(n+a)(n+b)}+\dfrac{\mu^{2}}{(n+b)^{2}}<s^{2}-s^{2}\dfrac{1}{(n+a)^{2}}$

Multiplying by $1/s^{2}$ , we get that

\displaystyle\dfrac{(sL)^{2}+1}{(n+a)^{2}}-\dfrac{2\mu L}{(n+a)(n+b)}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+b)^{2}}<1\,.

(42)

Now, we analyze two cases:
1) Suppose that $a\leq b$ . Then, $\dfrac{1}{n+b}\leq\dfrac{1}{n+a}$ . We consider $n>N^{\prime}:=\sqrt{(sL)^{2}+1+\left(\dfrac{\mu}{s}\right)^{2}}-a$ . This leads to the fact that $(n+a)^{2}>(sL)^{2}+1+\left(\dfrac{\mu}{s}\right)^{2}$ . Thus, $\dfrac{(sL)^{2}+1}{(n+a)^{2}}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+a)^{2}}<1$ . At the same time, we have the following inequalities:

\displaystyle\dfrac{(sL)^{2}+1}{(n+a)^{2}}-\dfrac{2\mu L}{(n+a)(n+b)}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+b)^{2}}

\displaystyle\leq\dfrac{(sL)^{2}+1}{(n+a)^{2}}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+b)^{2}}\leq\dfrac{(sL)^{2}+1}{(n+a)^{2}}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+a)^{2}}<1\,.

2) Suppose that $b\leq a$ . Then, $\dfrac{1}{n+a}\leq\dfrac{1}{n+b}$ . We consider $n>N^{\prime}:=\sqrt{(sL)^{2}+1+\left(\dfrac{\mu}{s}\right)^{2}}-b$ . This leads to the fact that $(n+b)^{2}>(sL)^{2}+1+\left(\dfrac{\mu}{s}\right)^{2}$ . Thus, $\dfrac{(sL)^{2}+1}{(n+b)^{2}}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+b)^{2}}<1$ . At the same time, we have the following inequalities:

\displaystyle\dfrac{(sL)^{2}+1}{(n+a)^{2}}-\dfrac{2\mu L}{(n+a)(n+b)}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+b)^{2}}

\displaystyle\leq\dfrac{(sL)^{2}+1}{(n+a)^{2}}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+b)^{2}}\leq\dfrac{(sL)^{2}+1}{(n+b)^{2}}+\left(\dfrac{\mu}{s}\right)^{2}\dfrac{1}{(n+b)^{2}}<1\,.

Finally, we end this example emphasizing that we can take $N^{\prime}=\sqrt{3+\left(\dfrac{\mu}{s}\right)^{2}}-a$ if $a\leq b$ and $N^{\prime}=\sqrt{3+\left(\dfrac{\mu}{s}\right)^{2}}-b$ if $b\leq a$ , and in this way we obtain an index that does not depend on the Lipschitz constant $L$ .

4 Numerical experiments involving convex functions

The role of the present section is to validate our theoretical results through some computational simulations, involving some convex test functions. Furthermore, we consider the particular cases of LT-S-IGAHD given in Example 24, Example 25 and Example 26 and we study the role of the parameters for each of them. Also, our codes for the optimization algorithms are implemented in Matlab and the analysis of the convexity of the chosen functions is based upon the a fundamental theory that can be found in [11].

4.1 Algorithmic description and objective functions

We consider the algorithms of the form LT-S-IGAHD (with the stepsize $s=h^{2}$ ), starting from $n\geq 1$ . We initialize the starting point $x_{0}\in\mathbb{R}^{2}$ for functions $f:\mathbb{R}^{2}\to\mathbb{R}$ . Furthermore, even though in many simulations $x_{1}\in\mathbb{R}^{2}$ is considered arbitrary, for consistency we choose $x_{1}=x_{0}-s\nabla f(x_{0})$ , namely the first iteration at $n=1$ is given by a step of gradient descent. Throughout the present section we count the total number of iterations starting from $n=0$ in order to encompass the initial value $x_{0}$ (even though our algorithms start from $n=1$ ). Furthermore, since for $x_{2}$ we need only the values $x_{0}$ , $x_{1}$ and $y_{1}$ , for brevity we set $y_{0}=x_{0}$ . Using tolerance values of low order $\epsilon=10^{-10}$ , we employ the stopping criteria as $\|f(x_{n})-f(x_{n-1})\|_{2}\leq\epsilon$ or as $\|f(x_{n})-f(x^{\ast})\|_{2}\leq\epsilon$ in the case when the objective function admits a unique minimum point $x^{\ast}\in\mathbb{R}^{2}$ . For each of the functions that we consider we compute the Lipschitz constant in order to present the behavior of the discretizations with respect to the parameters. Furthermore, we consider the stepsize $s<\dfrac{1}{L}$ by the fact that in the optimal case when $s=\dfrac{1}{L}$ no interpretation can be given since the algorithms stop just after a few iterations. We consider the following convex functions:
1) The case of convex, but not strongly convex function (see [4]), namely $f:\mathbb{R}^{2}\to\mathbb{R}$ , defined as

\displaystyle f(x_{1},x_{2})=(x_{1}+x_{2})^{2}\,,\,\forall(x_{1},x_{2})\in\mathbb{R}^{2}\,,

(F1)

It is easy to observe that $\nabla f(x_{1},x_{2})=\begin{pmatrix}2(x_{1}+x_{2})\\ 2(x_{1}+x_{2})\end{pmatrix}$ . Furthermore, the Hessian of $f$ is $\nabla^{2}f(x_{1},x_{2})=\begin{pmatrix}2&2\\ 2&2\end{pmatrix}$ . Taking $z=\begin{pmatrix}z_{1}\\ z_{2}\end{pmatrix}$ , we infer that $z^{T}\nabla^{2}f(x_{1},x_{2})z=2(z_{1}+z_{2})^{2}\geq 0$ , so the Hessian matrix $\nabla^{2}f\succeq 0$ . Using the second order characterization of convexity (see Theorem 7.12 from [11]), we find that the function is convex. Furthermore, computing the solution of the equation $\nabla f(x_{1},x_{2})=\begin{pmatrix}0\\ 0\end{pmatrix}$ and using the fact from above that the objective function is convex, by using Theorem 7.9 of [11], we obtain that the set of stationary points is given by $\operatorname*{argmin}f=\{(x_{1},x_{2})\in\mathbb{R}^{2}\,/\,x_{1}=-x_{2}\}$ . By the fact that we use the $l_{2}$ norm on the gradient, we aim at evaluating the subordinate matrix norm. By the fact that $\|\nabla^{2}f(x_{1},x_{2})\|_{2}$ consists of the largest square root of the absolute value of the matrix eigenvalues, we find that $\|\nabla^{2}f(x_{1},x_{2})\|=4$ , so that the Lipschitz constant of the gradient $\nabla f$ is also $L=4$ . On the other hand, since the $\operatorname*{argmin}$ set contains an infinity of elements, we use the stopping criteria $\|f(x_{n})-f(x_{n-1})\|_{2}\leq\epsilon$ , for a given tolerance $\epsilon>0$ sufficiently small. Moreover, it is also natural to adapt the Discrete-Energy function for the case when we do not actually now to what minimizer our algorithm converges. For this, we consider

	$\displaystyle\mathcal{E}_{n}:=t_{n}^{2}\left(f(x_{n})-f(x_{M})\right)+\dfrac{1}{2s}\\|\tilde{z}_{n}\\|^{2}$		(43)
	$\displaystyle\tilde{z}_{n}:=(x_{n-1}-x_{M})+t_{n}(x_{n}-x_{n-1})+\lambda_{n}t_{n+1}\nabla f(x_{n-1})\,,$

where $M$ is determined by the stopping condition, i.e. when we reach the tolerance factor $\epsilon$ , we consider the last iteration value $x_{n}=x_{M}$ by the fact that for $\epsilon$ small enough, $x_{M}$ is close to $x^{\ast}\in\operatorname*{argmin}f$ .
2) The hybrid norm function $f:\mathbb{R}^{2}\to\mathbb{R}$ , defined as

\displaystyle f(x_{1},x_{2})=\sqrt{1+x_{1}^{2}}+\sqrt{1+x_{2}^{2}}\,\,,\,\forall(x_{1},x_{2})\in\mathbb{R}^{2}\,.

(F2)

The gradient of the objective function has the form $\nabla f(x_{1},x_{2})=\begin{pmatrix}\dfrac{x_{1}}{\sqrt{1+x_{1}^{2}}}\\ \dfrac{x_{2}}{\sqrt{1+x_{2}^{2}}}\end{pmatrix}$ . On the other hand, the Hessian matrix associated to $f$ is given by $\nabla^{2}f(x_{1},x_{2})=\begin{pmatrix}\dfrac{1}{(1+x_{1}^{2})^{3/2}}&0\\ 0&\dfrac{1}{(1+x_{2}^{2})^{3/2}}\end{pmatrix}$ . Taking $z=\begin{pmatrix}z_{1}\\ z_{2}\end{pmatrix}$ we obtain that $z^{T}\nabla^{2}f(x_{1},x_{2})z=z_{1}^{2}\dfrac{1}{(1+x_{1}^{2})^{3/2}}+z_{2}^{2}\dfrac{1}{(1+x_{2}^{2})^{3/2}}\geq 0$ , so the Hessian satisfies $\nabla^{2}f\succeq 0$ . As for the case of F1, we use Theorem 7.12 from [11] and we reach the convexity of F2. Furthermore, computing the solution of the equation $\nabla f(x_{1},x_{2})=\begin{pmatrix}0\\ 0\end{pmatrix}$ and using the fact from above that the objective function is convex, we obtain that the set of stationary points is given by $\operatorname*{argmin}f=\{(0,0)\}$ (this follows from Theorem 7.9 of [11], namely the necessity and sufficiency of stationary points). By the fact that we use the $l_{2}$ norm on the gradient, we aim at evaluating the subordinate matrix norm. It is easy to observe that $\|\nabla^{2}f(x_{1},x_{2})\|_{1}=\|\nabla^{2}f(x_{1},x_{2})\|_{\infty}=\max\Big{\{}\dfrac{1}{(1+x_{1}^{2})^{3/2}},\dfrac{1}{(1+x_{2}^{2})^{3/2}}\Big{\}}\leq 1$ . Since $\nabla^{2}f(x_{1},x_{2})\in\mathbb{R}^{2\times 2}$ , by elementary algebra computations, we have that $\|\nabla^{2}f(x_{1},x_{2})\|_{2}\leq\sqrt{2}\|\nabla^{2}f(x_{1},x_{2})\|_{\infty}$ , so the tightest Lipschitz constant of the gradient is $L=\sqrt{2}$ . At the same time, since the $\operatorname*{argmin}$ set contains just one global minimum point, we use the stopping criteria $\|f(x_{n})-f(x^{\ast})\|_{2}\leq\epsilon$ , where $x^{\ast}$ is the unique element of $\operatorname*{argmin}f$ . Also, we consider the Discrete-Energy for the Lyapunov functional.

Refer to caption — (a) The representation of F1 along with
the set of minimum points $\operatorname*{argmin}f$

Following the proof of Theorem 20, we consider $N_{1}=\alpha-1$ and $N_{2}$ given by $G_{n}$ , $H_{n}$ and $I_{n}$ as in the proof of the aforementioned theorem. Now we turn our attention to some choices of the parameters that appear in the three examples from the previous section, concerning the index $N^{\prime}$ . However, we specify that even though $G_{n}$ , $H_{n}$ and $I_{n}$ depend on the Lipschitz constant of the gradient $\nabla f$ (this is also valid for the case of Nesterov’s method AGM2_α), the empirical convergence is in fact determined by the index $N^{\prime}$ that contain some predetermined values of the inherent parameters.

4.2 Some discussions over the algorithms parameters

In this subsection we consider some situations over the parameters referring to Example 24, Example 25 and Example 26. For each of them, we analyze in an empiric manner the key role of the underlying parameters over the convergence of the optimization algorithms of type LT-S-IGAHD.
Likewise, for two parameters $A$ and $B$ , we use the heuristic notations ” $A\ll B$ ” for the situation when $A$ is much smaller than $B$ and ” $A\gg B$ ” for the case $A$ is much larger than $B$ , respectively. Therefore, for our analysis, we consider a heuristic analysis for our algorithms.

The case involving Example 24) We aim at considering two subcases:
A.) When $\dfrac{\mu}{s}\ll 1$ . This is, in fact, the practical case when the parameters have small value. Because, in general, $s$ is relatively small, this allow us to take $\mu\ll 1$ . For this, we consider the following situations:

A1.)

When $a+\dfrac{1}{4}<b$ . Now, if $a\gg 1$ ( $a$ has a big value), then we can take $a\gg\left(\dfrac{\mu}{s}\right)^{2}$ , so $N^{\prime}\ll 0$ . On the other hand, when $a\ll 1$ , then we can take $a\approx\left(\dfrac{\mu}{s}\right)^{2}$ , so $N^{\prime}>0$ , i.e. $N^{\prime}=O\left(\mu L\right)$ , but $\mu\ll 1$ annihilates the Lipschitz constant $L$ in the index value $N^{\prime}$ .
A2.)

When $b\leq a+\dfrac{1}{4}$ . If $b\gg 1$ , then $a\gg 1$ , so we may take $a\gg\left(\dfrac{\mu}{s}\right)^{2}$ . This is in fact the same situation as before when $a$ can take a big value, because the difference under the square root, namely $a-b$ , can be close to $0$ , in the sense that for $b\approx a$ , $N^{\prime}\approx\dfrac{1}{2}-b$ , i.e. $N^{\prime}\approx 0$ . On the other hand, for the case when $b\ll 1$ and $a\gg 1$ , we find that $N^{\prime}\approx\dfrac{\sqrt{\Delta}}{2}$ , hence $N^{\prime}=O(\sqrt{a})$ , if we can take $a$ large enough such that $a\gg\left(\dfrac{\mu}{s}\right)^{2}$ .

B.) When $\dfrac{\mu}{s}\gg 1$ . By the fact that in many situations $s\ll 1$ , we can take $\mu>1$ or $\mu\gg 1$ . For this, we consider the following situations:

B1.)

When $a+\dfrac{1}{4}<b$ . If $a\gg 1$ , then $b\gg 1$ . We can take $a\approx\left(\dfrac{\mu}{s}\right)^{2}$ , so $N^{\prime}=O(\mu L)$ , which can be eventually large. On the other hand, if $b\gg 1$ and $a\ll 1$ , we can take $a\ll\left(\dfrac{\mu}{s}\right)^{2}$ , so $N^{\prime}=O\left(\left(\dfrac{\mu}{s}\right)^{2}\right)$ , i.e. $N^{\prime}\gg 1$ .
B2.)

When $b\leq a+\dfrac{1}{4}$ . If $a\gg 1$ , thus $b\gg 1$ , then we are in a similar situations as presented above, since we can take $a$ close to $b$ such that the difference $a-b$ in the square root in the definition of $N^{\prime}$ is close to $0$ . On the other hand, if $a\gg 1$ and $b\ll 1$ , then we can take $a\approx\left(\dfrac{\mu}{s}\right)^{2}$ , thus $N^{\prime}=O\left(\sqrt{a}\right)$ , which is equivalent to $N^{\prime}\gg 1$ .

The case involving Example 25) For this example, we consider the only valid situation when $2\sqrt{s}>\beta$ :
Since the stepsize $s$ has small values in many practical applications, we must consider $\beta<1$ or $\beta\approx 0$ . For small values of $\beta$ , we observe that $N^{\prime}=O\left(\dfrac{\dfrac{\mu}{\sqrt{s}}-2\sqrt{s}b+\sqrt{\Delta}}{2(2\sqrt{s}-\beta)}\right)$ . Now, we consider the following subcases:
C1.) When $\dfrac{\mu}{\sqrt{s}}\gg 1$ and $b$ is small or has moderate values, we have that

N^{\prime}=O\left(\dfrac{\dfrac{\mu}{\sqrt{s}}+\sqrt{(2\sqrt{s}-\beta)\dfrac{\mu}{\sqrt{s}}+\left(\dfrac{\mu}{\sqrt{s}}\right)^{2}}}{2(2\sqrt{s}-\beta)}\right)\gg 1.

C2.) If we are in one of the cases when $\dfrac{\mu}{\sqrt{s}}<1$ or $\dfrac{\mu}{\sqrt{s}}\ll 1$ , then, by the analysis made above, we obtain that $N^{\prime}\approx 0$ , since $\mu<\sqrt{s}$ or $\mu\ll 1$ .
C3.) If $b\gg 0$ and dominates the other parameters, then it is easy to see that $N^{\prime}=O\left((2\sqrt{s}-\beta)\beta\sqrt{b}\right)$ , i.e. $N^{\prime}\gg 1$ . Furthermore, in order to have $N^{\prime}\approx 0$ , then we must take $\beta\ll b$ .

The case involving Example 26) This is similar with the first example involving the non-negative coefficients $a$ , $b$ and $\mu$ . For brevity, we consider the following cases for our heuristic interpretation:
D.) When $\dfrac{\mu}{s}\ll 1$ . Since many times $s\ll 1$ , we have that $\mu\ll 1$ . Here, we consider the following subcases:

D1.)

When $a\leq b$ . If $a\ll 1$ and $b\gg 1$ , i.e. $a$ small and $b$ has large values, we find that $N^{\prime}\approx\sqrt{3+\left(\dfrac{\mu}{s}\right)^{2}}=O\left(\left(\dfrac{\mu}{s}\right)^{2}\right)\approx 0$ . On the other hand, if $a\gg 1$ and $b\gg 1$ , we obtain that $N^{\prime}\approx 0$ .
D2.)

When $b\leq a$ . If $a\gg 1$ , then $b\gg 1$ and, as before, $N^{\prime}\approx 0$ . On the other hand, if $a\gg 1$ and $b\ll 1$ , it follows that $N^{\prime}\approx\sqrt{3+\left(\dfrac{\mu}{s}\right)^{2}}$ , i.e. the index $N^{\prime}$ is small.

E.) When $\dfrac{\mu}{s}\gg 1$ , then $\sqrt{(sL)^{2}+1+\left(\dfrac{\mu}{s}\right)^{2}}\gg 0$ . Likewise, in practical situations we are obliged to take $\mu\gg 1$ . We consider the following subcases:

E1.)

When $a\leq b$ . In order to obtain the optimal case when $N^{\prime}\leq 0$ , we can take $a\gg 0$ , more precisely $a\geq\sqrt{3+\left(\dfrac{\mu}{s}\right)^{2}}\geq\sqrt{(sL)^{2}+1+\left(\dfrac{\mu}{s}\right)^{2}}$ . It is trivial to see that since $a$ is large, then $b$ is also large.
E2.)

When $b\leq a$ . As in the previous case, we consider the optimal case when $N^{\prime}\leq 0$ , so we can take easily take $b\geq\sqrt{3+\left(\dfrac{\mu}{s}\right)^{2}}$ . In this situation $b\gg 1$ , so $a\gg 1$ .

4.3 A computational perspective over the optimization algorithms

In the present subsection, we consider some numerical experiments over the algorithms of type LT-S-IGAHD, taking into account the theoretical rate of convergence in the convex objective function obtained in the proof of Theorem 20. At the same time, we precisely consider the cases and the subcases presented in Subsection 4.2.
In Figure 2(a), we have considered the convex objective function F2 along with the initial condition $x_{0}=\begin{pmatrix}1\\ -2\end{pmatrix}$ Also, we have set $\alpha=3$ and the tolerance $\epsilon=10^{-10}$ . We have the behavior of the discretization schemes along the trajectories when the stepsize is taken as $s=0.001$ in order to simulate the underyling dynamical systems. In the Subfigure 2(a), we have considered the IGAHD algorithm $\beta=10^{-2}$ and $b=1$ (black color) and the case when $\beta=10^{-1}$ and $b=1$ (blue color). Furthermore, we have considered also the LT-S-IGAHD from Example 25, by the fact that when $\mu\approx 0$ , then the algorithms becomes IGAHD-type. We have chosen the following parameters : $\mu=10^{-1}$ , $\beta=10^{-2}$ , $b=1$ (red color), $\mu=10^{-1}$ , $\beta=10^{-1}$ , $b=1$ (saddle brown color, $[139,69,19]$ in RGB), $\mu=0.25$ , $\beta=10^{-2}$ , $b=1$ (yellow color) and $\mu=0.25$ , $\beta=10^{-1}$ , $b=1$ (rustic red color, $[75,0,13]$ in RGB). With this choices of coefficients, we observe that our algorithms have trajectories very close to the ones of IGAHD. This is valided by the theory where when $\mu$ is taken small, then the only difference between the algorithms gave in Example 25 and IGAHD is between $\nabla f(x_{n})$ and $\nabla f(x_{n-1})$ . On the other hand, in Subfigure 2(b), we have considered the algorithm from Example 26, with the following choices of the algorithm’s coefficients: $a=2$ , $\mu=10^{-15}$ , $b=0.1$ (red color), $a=2$ , $\mu=10^{-1}$ , $b=0.1$ (blue color) and finally, $a=25$ , $\mu=2$ , $b=3$ (black color). In this situation we observe the phenomenon of the high amplification of the coefficients, namely when the parameters have large values, the trajectories have a high acceleration toward the critical point. From this subfigure and also from our preliminary numerical simulations, we have observed that there are some optimal choices of the parameters that lead to a faster convergence, that is given by a low number of iterations for a given tolerance.

Now, we focus on the discrete velocity $v_{n}=(x_{n}-x_{n-1})/\sqrt{s}$ , where $\sqrt{s}=h>0$ . This is considered in Figure 3. We have taken $\epsilon=10^{-10}$ , $\alpha=3$ and $x_{0}=\begin{pmatrix}1\\ -2\end{pmatrix}$ For the algorithm LT-S-IGAHD from Example 25, we have set the values for the parameters as follows: $\mu=0.5$ , $\beta=10^{-5}$ and $b=2$ . In Subfigre 3(a) we have plotted the first component of velocity in $\mathbb{R}^{2}$ for $s=0.25$ , while in Subfigure 3(b) we have taken a smaller stepsize $s=0.025$ . We observe that, in the beginning, the velocity presents high oscillations and then it stabilizes through the value $0$ , and this is related to the fact that the velocity of the underlying evolution equation has the same behavior. The same stands true also for the second component of velocity and for other choices of the algorithm’s parameters. Moreover, we observe that the stabilization process acts much faster than the convergence of the numerical scheme, when the tolerance is of low value.
Second to last, we turn our attention to Table 1, Table 2, Table 3 and Table 4. Here, we have validated all of the cases from Subsection 4.2 involving both convex functions F1 and F2, respectively. First thing to note is that error represents the difference in two consecutive values in the objective function that is less than the given tolerance $\epsilon$ and along with this, it satisfies and additional stopping criteria, namely $n>N$ , where $N=\max\{N_{1},N_{2},N^{\prime}\}$ . We observe that our heuristic analysis stands very close to the simulations that we have made. For this, observe the values of $N^{\prime}$ given for each of the cases. An important difference is in the case $A_{2}$ where the value has an error by the fact that the coefficient $a$ is not much larger than $1$ . Furthermore, another difference is for the case $B_{2}$ when we do not have taken $a-b\approx 0$ and this shows the fact that the heuristic interpretation that we gave leads to instability in the coefficients, in the sense that the values must take extreme values (see also the situation $B_{1}$ ). This is linked to the fact that in Subsection 4.2 we have considered an asymptotic analysis when some coefficients have very small values while others have very large values. The cases $A_{2}$ , $B_{1}$ and $B_{2}$ justifies the difference between simulations and the heuristic analysis. On the other hand, for all of the other cases presented in the four tables, $N^{\prime}$ is very close to the actual theoretical value. Likewise, it is interesting to note that the empirical values $N_{2}$ and $N^{\prime}$ have an empiric relation, namely in the case of moderate or of high value coefficients, $N^{\prime}$ can be much greater than $N_{2}$ , so in this sense, in the majority of cases, the index $N^{\prime}$ has a large influence over the convergence of the algorithms of type LT-S-IGAHD. Nevertheless, there are some cases when $N_{2}$ has a larger influence over the rate of convergence, i.e. when the coefficients have extremely have numerical values, then $N_{2}$ increases (with respect to the values presented above in the other cases) and $N^{\prime}$ decreases. Even though $N^{\prime}$ is larger than $N_{2}$ , the difference is lower than in the previous cases. Finally, we point out that in our numerical simulations, the coefficient $G_{n}$ is strictly positive, when $n>\max\{N_{1},N^{\prime}\}$ . This remains valid for all of the cases and for all of the three examples that we gave and it is consistent with the arguments from the third step of the proof of Theorem 20.

Table 1: Analysis of the cases

A_{1}

A_{2}

B_{1}

and

B_{2}

of Subsection 4.2,
for the Example 24 and for the function F1

	$error$	$\varepsilon$	$\mu$	$a$	$b$	$N_{2}$	$N^{\prime}$	$N$
$A_{1}$	1.22e-11	1e-10	1e-02	4.0	10.0	3.91	-3.56	3.91
$A_{1}$	3.31e-11	1e-10	1e-02	1e-02	10.0	3.99	0.43	3.99
$A_{2}$	2.82e-11	1e-10	1e-02	10.0	4.0	3.80	-1.00	3.80
$A_{2}$	1.50e-12	1e-10	1e-05	3.0	1e-01	3.92	2.17	3.92
$B_{1}$	0.0	1e-10	1.0	100.0	102.0	3.74	11.63	11.63
$B_{1}$	0.0	1e-10	2.0	5e-01	6.0	3.48	422.45	422.45
$B_{2}$	0.0	1e-10	1.5	2.0	1.75	3.58	240.29	240.29
$B_{2}$	0.0	1e-10	1.5	225	1e-02	8.52	17.29	17.29

Table 2: Analysis of the cases

A_{1}

A_{2}

B_{1}

and

B_{2}

of Subsection 4.2,
for the Example 24 and for the function F2

	$error$	$\varepsilon$	$\mu$	$a$	$b$	$N_{2}$	$N^{\prime}$	$N$
$A_{1}$	7.03e-11	1e-10	1e-02	4.0	10.0	3.38	-3.91	3.38
$A_{1}$	1.92e-12	1e-10	1e-02	1e-02	10.0	3.38	0.08	3.38
$A_{2}$	3.47e-11	1e-10	1e-02	10.0	4.0	3.37	-1.0	3.37
$A_{2}$	7.53e-11	1e-10	1e-05	3.0	1e-01	3.38	2.17	3.38
$B_{1}$	0.0	1e-10	1	100	102	3.50	4.04	4.04
$B_{1}$	0.0	1e-10	2.0	5e-01	6.0	3.43	407.54	407.54
$B_{2}$	0.0	1e-10	1.5	2.0	1.75	3.5	229.04	229.04
$B_{2}$	0.0	1e-10	1.5	225	1e-02	4.46	15.50	15.50

Table 3: Analysis of the cases

D_{1}

D_{2}

E_{1}

and

E_{2}

of Subsection 4.2,
for the Example 26 and for the function F1

	$error$	$\varepsilon$	$\mu$	$a$	$b$	$N_{2}$	$N^{\prime}$	$N$
$D_{1}$	3.37e-11	1e-10	0.0	0.25	3.5	3.18	0.83	3.18
$D_{1}$	9.79e-11	1e-10	1e-3	1.25	5.5	3.18	-0.17	3.18
$D_{2}$	3.28e-11	1e-10	1e-3	5.5	1.25	3.19	-0.17	3.19
$D_{2}$	1.11e-11	1e-10	1e-3	3.5	0.25	3.19	0.83	3.19
$E_{1}$	0.0	1e-10	2.0	21.0	24.0	5.82	-0.97	5.82
$E_{2}$	0.0	1e-10	2.0	24.0	21.0	6.09	-0.97	6.09

Table 4: Analysis of the cases

D_{1}

D_{2}

E_{1}

and

E_{2}

of Subsection 4.2,
for the Example 26 and for the function F2

	$error$	$\varepsilon$	$\mu$	$a$	$b$	$N_{2}$	$N^{\prime}$	$N$
$D_{1}$	4.21e-12	1e-10	0.0	0.25	3.5	3.23	0.76	3.23
$D_{1}$	9.18e-11	1e-10	1e-03	1.25	5.5	3.23	-0.24	3.23
$D_{2}$	5.68e-11	1e-10	1e-03	5.5	1.25	3.23	-0.24	3.23
$D_{2}$	7.18e-11	1e-10	1e-03	3.5	0.25	3.23	0.76	3.23
$E_{1}$	0.0	1e-10	2.0	21.0	24.0	4.14	-0.97	4.14
$E_{2}$	0.0	1e-10	2.0	24.0	21.0	4.24	-0.97	4.24

Conclusions and Discussions

The novelty of the present paper consists in the following both theoretical and practical aspects regarding Hessian driven damping dynamical systems and also their numerical methods. We conclude the link between the discrete case and the continuous counterpart of the optimization algorithms, by pointing out our major contributions and some discussions that arise naturally from our research:

1)

In Subsection 2.1, we have shown that Nesterov’s algorithm AGM2_α is in fact a natural discretization of a Hessian-type dynamical system. Even though other authors (see [18] and [33]) have pointed out a link between continuous Hessian dynamical systems and Nesterov algorithm, we have introduced a rigorous perspective such that Nesterov method can be interpreted as a Lie-Trotter discretization of the underlying second evolution equation. This approach is totally new with respect to inertial algorithms, in the following sense: only recently (see [1] and [29]) it was shown that Nesterov’s method is a natural discretization of the Extended-AVD dynamical system. Before that, no other constructions were available such that Nesterov algorithm could be interpreted as an explicit numerical scheme of a continuous dynamical system without an extra-gradient step, so our method represents a totally new approach to optimization. Furthermore, the use of the splitting technique is somewhat new to optimization problems, so that this article strengthens the relationship between distinct mathematical fields.
2)

We have presented a geometrical-type perspective upon optimization methods in the setting of minimization problems that involve convex objective functions (see Subsection 2.3). Likewise, we gave intuitive explanations over the choice of our algorithms in Subsection 2.4. At the same time, in Subsection 3.2 we have constructed Lemma 15 that represents a generalization of Lemma 14 that was given in [4]. Also, our proof of Theorem 20 is based upon the proof from the same paper, namely [4], but in a much more general setting.
3)

We have introduced new algorithms that contain the discrete version of the Hessian system, namely LT-S-IGAHD. Also, the particular case given in Example 25 resembles the newly introduced algorithm IGAHD in [4]. Also, in Example 24 and Example 26 we have gave general algorithms that contain three arbitrary parameters (under some suitable conditions) that have numerical benefits in practical applications (see Subsections 4.2 and 4.3). On the other hand, we emphasize that our algorithms have been developed in an elegant manner based upon the combination of splitting methods with symplectic algorithms, similar to the simpler case of Nesterov inertial algorithm (see [29]).
4)

We have developed some new optimization algorithms based upon two gradient evaluations at each iteration. For each of these numerical schemes we have presented the associated continuous versions, such that these numerical schemes are constructed through various geometric principles like splitting and symplectic methods. The associated continuous dynamical systems DynSys-ARDM and the system from Theorem 7 are worth studying by the fact that their discretizations are obtained in a natural manner, as in the case with Nesterov’s method (see [1]and [29]).

Acknowledgments

I am gratefully indebted to Szilard Csaba László for his suggestions concerning the organization of this article and for the helpful comments regarding the dynamical systems and the optimization algorithms of the present research paper.

References

[1] C.D. Alecsa, S.C. László, T. Pinţa, An extension of the second order dynamical system that models Nesterov’s convex gradient method, 2019, arXiv:1908.02574.
[2] C.D. Alecsa, I. Boros, T. Pinţa, New optimization algorithms for neural network training using operator splitting techniques, 2019, arXiv:1904.12952.
[3] F. Álvarez, H. Attouch, J. Bolte, P. Redont, A second-order gradient like dissipative dynamical system with Hessian-driven damping. Application to optimization and mechanics, J. Math. Pures Appl. 81 (2002), no. 8, pp. 747-779.
[4] H. Attouch, Z. Chbani, J. Fadili, H. Riahi, First-order optimization algorithms via inertial systems with Hessian driven damping, 2019, arXiv:1907.10536.
[5] H. Attouch, A. Cabot, Asymptotic stabilization of inertial gradient dynamics with time-dependent viscosity, J. Differential Equations 263 (2017), pp. 5412-5458, doi:10.1016/j.jde.2017.06.024.
[6] H. Attouch, A. Cabot, Convergence rates of inertial forward-backward algorithms, SIAM J. Optim (28) 2018, no.1, pp. 849-874, https://doi.org/10.1137/17M1114739.
[7] H. Attouch, X. Goudou, P. Redont, The heavy ball with friction method. I. The continuous dynamical system: global exploration of the local minima of a real-valued function by asymptotic analysis of a dissipative system, Commun. Contemp. Math. 2 (2000), pp. 1-34.
[8] H. Attouch, J. Peypouquet, P. Redont, Fast convex minimization via inertial dynamics with Hessian driven damping, J. Differential Equations 261 (2016), pp. 5734-5783.
[9] N. Bansal, A. Gupta, Potential-function proofs for gradient methods, Theory of Computing (2019), pp. 1-32.
[10] A. Bátkai, P. Csomós, B. Farkas, G. Nickel, Operator splitting for non-autonomous evolution equations, J. Functional Analysis 260 (2011), pp. 2163-2190.
[11] A. Beck, Introduction to Nonlinear Optimization: Theory, Algorithms, and Applications with MATLAB, vol. 19, SIAM, 2014.
[12] S. Blanes, F. Casas, A. Murua, Splitting and composition methods in the numerical integration of differential equations, Bol. Soc. Esp. Mat. Apl. 45 (2008), pp. 89-145.
[13] R.I. Bot, R. Csetnek, S.C. László, A second order dynamical approach with variable damping to nonconvex smooth minimization, Applic. Analysis, 2018, doi:10.1080/00036811.2018.1495330.
[14] A. Cabot, H. Engler, S. Gadat, On the long time behavior of second-order differential equations with asymptotically small dissipation, Trans. Amer. Math. Soc. 361 (2009), pp. 5983-6017.
[15] A. Cabot, H. Engler, S. Gadat, Second order differential equations with asymptotically small dissipation and piecewise flat potentials, Electron J. Differential Equations 17 (2009), pp. 33-38.
[16] C. Castera, J. Bolte, C. Févotte, E. Pauwels, An inertial Newton algorithm for deep learning, 2019, arXiv:1905.12278.
[17] A. Chambolle, C. Dossal, On the convergence of the iterates of FISTA, HAL Id : hal-01060130, 2014, https://hal.inria.fr/hal-01060130v3.
[18] L. Chen, H. Luo, First order optimization methods based on Hessian-driven Nesterov accelerated gradient flow, 2019, arXiv:1912.09276.
[19] A. Defazio, On the curved geometry of accelerated optimization, Adv. Neural Inf. Processing Syst. 33 (NIPS 2019), 2019.
[20] K. Feng, On difference schemes and symplectic geometry, Proceedings of the 5-th Intern. Symposium on differential geometry $\&$ differential equations, August 1984, Beijing (1985), pp. 42-58.
[21] G. França, D.P. Robinson, R. Vidal, Gradient flows and accelerated proximal splitting methods, 2019, arXiv:1908.00865.
[22] G. França, J. Sulam, D.P. Robinson, R. Vidal, Conformal symplectic and relativistic optimization, 2019, arXiv:1903.04100.
[23] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration. Structure-Preserving Algorithms for ODE’s, second edition, Springer-Verlag Berlin Heidelberg, 2006, doi:10.1007/3-540-30666-8.
[24] E. Hairer, G. Wanner, Euler methods, explicit, implicit and symplectic. In : Encyclopedia of Applied and Computational Mathematics, B. Engquist ed., Springer, Berlin Heidelberg, 2015, pp. 451-455.
[25] E. Hansen, F. Kramer, A. Ostermann, A second-order positivity preserving scheme for semilinear parabolic problems, Appl. Num. Math. 62 (2012), pp. 1428-1435.
[26] E. Hansen, A. Ostermann, Exponential splitting for unbounded operators, Math. of Comput. 78 (2009), no. 267, pp. 1485-1496.
[27] S.C. László, Convergence rates for an inertial algorithm of gradient type associated to a smooth nonconvex minimization, 2018, arXiv:1807.00387.
[28] G. Marchuk, Some applications of splitting-up methods to the solution of mathematical physics problems, Aplikace Matematiky 13 (1968), pp. 103-132.
[29] M. Muehlebach, M.I. Jordan, A dynamical systems perspective on Nesterov acceleration, 2019, arXiv:1905.07436.
[30] Y. Nesterov, A method for solving a convex programming problem with convergence rate $O(1/k^{2})$ (in Russian), Soviet Math. Dokl. 27 (1983), no. 2, pp. 543-547.
[31] Y. Nesterov, Introductory lectures on convex optimization : A basic course, vol. 87 of Applied Optimization, Kluwer Academic Publishers, Boston, MA, 2004.
[32] N.C. Nguyen, P. Fernandez, R.M. Freund, J. Peraire, Accelerated residual methods for iterative solution of systems of equations, SIAM J. Sci. Computing 40 (2018), pp. A3157-A3179.
[33] M. Laborde, A.M. Oberman, A Lyapunov analysis for accelerated gradient methods : from deterministic to stochastic case, 2019, arXiv:1908.07861.
[34] B.T. Polyak, Some methods of speeding up the convergence of iteration methods, USSR Comput. Math. Math. Phys. (4) 1964, pp. 1-17.
[35] J.M. Sanz-Serna, Symplectic methods. In : Encyclopedia of Applied and Computational Mathematics, B. Engquist ed., Springer, Berlin Heidelberg, 2015, pp. 1451-1458.
[36] B. Shi, S.S. Du, W.J. Su, M. I. Jordan, Acceleration via symplectic discretization of high-resolution differential equations, 2019, arXiv:1902.03694.
[37] B. Shi, S.S. Iyengar, A conservation law method based on optimization. In : Mathematical Theories of Machine Learning - Theory and Applications, Springer, Cham, 2019, pp. 63-85 https://doi.org/10.1007/978-3-030-17076-9_8.
[38] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968), no. 3, pp. 506-517.
[39] W.J. Su, S. Boyd, E.J. Candès, A differential equation for modeling Nesterov’s accelerated gradient method : Theory and insights, J. Machine Learning Res. (17) 2016, pp. 1-43.
[40] T. Sun, P. Yin, D. Li, C. Huang, L. Guan, H. Jiang, Non-ergodic convergence analysis of heavy-ball algorithms, AAAI, 2019.
[41] H.F. Trotter, On the product of semi-groups of operators, Proc. Am. Math. Soc. 10 (1959), no. 4, pp. 545-551.
[42] P. Tseng, On accelerated proximal-gradient methods for convex-concave optimization, Manuscript, May 21 2018, mit.edu/~dimitrib/PTseng/papers/apgm.pdf.
[43] R. de Vogelaere, Methods of integration which preserve the contact transformation property of the Hamiltonian equations, Report no. 4, Dept. Math., Univ. of Notre Dame, Notre Dame, Ind., 1956.

2020

Positive periodic solutions for Lotka-Volterra systems with a general attack rate

April 29, 2022

On a class of Bernstein-type rational functions

March 13, 2023

Properties of discrete non-multiplicative operators

March 13, 2023

	$\displaystyle(t_{n+1}-1)f(x_{n+1})+f(x_{n+1})$	$\displaystyle\leq(t_{n+1}-1)f(x_{n})+f(x^{\ast})+(t_{n+1}-1)\langle\nabla f(y_{n}),y_{n}-x_{n}\rangle+\langle\nabla f(y_{n}),y_{n}-x^{\ast}\rangle$
		$\displaystyle-(t_{n+1}-1)\left[\mathcal{A}_{2}(s)+\mathcal{A}_{3}(s;\gamma_{n},L)\right]\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle-\mathcal{A}_{1}(s)(t_{n+1}-1)\\|\nabla f(y_{n})\\|^{2}$
		$\displaystyle-(t_{n+1}-1)\left[\mathcal{A}_{4}(s)+\mathcal{A}_{5}(s;\gamma_{n})\right]\\|\nabla f(x_{n})\\|^{2}-\mathcal{A}_{1}(s)\\|\nabla f(y_{n})\\|^{2}$
		$\displaystyle-\mathcal{A}_{3}(s;\gamma_{n},L)\langle\nabla f(y_{n}),\nabla f(x_{n})\rangle-\mathcal{A}_{5}(s;\gamma_{n})\\|\nabla f(x_{n})\\|^{2}\,.$

	$\displaystyle\tilde{F}_{n}\geq\dfrac{s}{2}t_{n+1}\\|Y\\|^{2}$	$\displaystyle-\left[(t_{n+1}-1)\mathcal{A}_{2}(s)+t_{n+1}\mathcal{A}_{3}(s;\gamma_{n},L)+(\lambda_{n}+\omega_{n})t_{n+1}\right]\\|Y\\|\\|X\\|$		(36)
		$\displaystyle+\left[(t_{n+1}-1)\mathcal{A}_{4}(s)+t_{n+1}\mathcal{A}_{5}(s;\gamma_{n})\right]\\|X\\|^{2}:=F_{n}\,.$

The long time behavior and the rate of convergence of symplectic convex algorithms obtained via splitting discretizations of inertial damping systems

Abstract

Authors

Keywords

Paper coordinates

PDF

About this paper

Journal

Publisher Name

DOI

Print ISSN

Online ISSN

References

Paper (preprint) in HTML form

The long time behavior and the rate of convergence of symplectic convex algorithms obtained via splitting discretizations of inertial damping systems.

1 Preliminaries

1.1 Some notes on splitting operator techniques

1.2 A brief review of symplectic methods

1.3 The organization of the paper

2 The discretizations of the dynamical systems

2.1 Nesterov’s accelerated method is a Lie-Trotter splitting method applied to a Hessian-based dynamical system

Proposition 1.

Proof .

Theorem 2.

Proof .

Remark 3.

2.2 The role of symplectic methods in optimization algorithms

Theorem 4.

Proof .

Remark 5.

2.3 Discussion on some splitting techniques for the extended AVD dynamical system

Remark 6.

2.4 Fixed points, spurious solutions and exploding gradients : Heuristic interpretation

3 Proposed algorithms and their rate of convergence

3.1 Modified Hessian-driven damping continuous and discrete systems

Theorem 7.

Proof .

Remark 8.

Remark 9.

Remark 10.

Remark 11.

Remark 12.

3.2 The rate of convergence of the optimization methods

Lemma 13 (Descent lemma).

Lemma 14 (Extended descent lemma).

Lemma 15 (Three points extended lemma).

Proof .

Corollary 16.

Lemma 17.

Proof .

Lemma 18.

Proof .

Remark 19.

Theorem 20.

Proof .

Remark 21.

Remark 22.

Remark 23.

Example 24.

Example 25.

Example 26.

4 Numerical experiments involving convex functions

4.1 Algorithmic description and objective functions

4.2 Some discussions over the algorithms parameters

4.3 A computational perspective over the optimization algorithms

Conclusions and Discussions

Acknowledgments

References

Related Posts