New optimization algorithms for neural network training using operator splitting techniques

TL;DR

This paper introduces novel optimization algorithms for neural network training based on operator splitting techniques, demonstrating their effectiveness through empirical convergence analysis and validation on standard datasets.

Contribution

It proposes a new class of optimization algorithms using operator splitting for neural networks, with empirical validation of convergence and performance.

Findings

Algorithms converge towards local minima

Validated on MNIST, MNIST-Fashion, CIFAR-10 datasets

Show improved or comparable accuracy and loss reduction

Abstract

In the following paper we present a new type of optimization algorithms adapted for neural network training. These algorithms are based upon sequential operator splitting technique for some associated dynamical systems. Furthermore, we investigate through numerical simulations the empirical rate of convergence of these iterative schemes toward a local minimum of the loss function, with some suitable choices of the underlying hyper-parameters. We validate the convergence of these optimizers using the results of the accuracy and of the loss function on the MNIST, MNIST-Fashion and CIFAR 10 classification datasets.

Tables

Table 1. Table 1: Scores of the accuracy metric versus step-size and epochs for different optimizers on MNIST dataset (left) and on MNIST-Fashion dataset (right), using N ​ L ​ L ​ L ​ o ​ s ​ s 𝑁 𝐿 𝐿 𝐿 𝑜 𝑠 𝑠 NLLLoss .

	MNIST			MNIST-Fashion
	20	50	100	20	50	100	$h$
SGD							${10}^{-3}$
							${10}^{-2}$
							${10}^{-1}$
NaG (0.5)							${10}^{-3}$
							${10}^{-2}$
							${10}^{-1}$
SSA 1							${10}^{-3}$
							${10}^{-2}$
							${10}^{-1}$
SSA 2							${10}^{-3}$
							${10}^{-2}$
							${10}^{-1}$
RMSProp							${10}^{-3}$
							${10}^{-2}$
							${10}^{-1}$
Adam							${10}^{-3}$
							${10}^{-2}$
							${10}^{-1}$
Adadelta							$1$
SSA 1 - Ada							$1$

Table 2. Table 2: Scores of the accuracy metric versus step-size and epochs for different optimizers on CIFAR 10 dataset, using Cross Entropy Loss .

	20	50	100	150	$h$
SGD					${10}^{-3}$
NaG (0.9)					${10}^{-3}$
SSA 1					${10}^{-3}$
SSA 2					${10}^{-3}$
RMSProp					${10}^{-3}$
Adam					${10}^{-3}$
Adagrad					${10}^{-3}$
Adadelta					$1$
SSA 1 - Ada					$1$

Table 3. Table 3: MNIST execution time.

	mean	std	min	25%	50%	75%	max	sum
Adadelta	$10.25$	$1.51$	$7.63$	$8.97$	$10.45$	$11.44$	$12.96$	$1025.33$
Adam	$10.22$	$1.70$	$7.39$	$8.80$	$10.06$	$11.55$	$14.48$	$1022.28$
RMSProp	$9.14$	$1.58$	$7.55$	$7.73$	$8.51$	$10.33$	$12.82$	$914.26$
NaG	$10.24$	$1.82$	$7.44$	$8.81$	$10.08$	$11.80$	$14.40$	$1024.89$
SGD	$8.39$	$1.12$	$7.52$	$7.73$	$7.89$	$8.69$	$12.63$	$839.54$
SSA1	$10.01$	$1.53$	$7.61$	$8.58$	$10.11$	$11.25$	$13.29$	$1001.22$
SSA2	$9.57$	$1.50$	$7.48$	$8.20$	$9.39$	$10.64$	$13.99$	$957.61$
SSA1-Ada	$11.61$	$1.42$	$8.83$	$10.62$	$11.54$	$12.58$	$15.40$	$1161.33$

Table 4. Table 4: MNIST-Fashion execution time.

	mean	std	min	25%	50%	75%	max	sum
Adadelta	$9.84$	$1.62$	$7.59$	$8.29$	$9.85$	$11.03$	$13.96$	$984.44$
Adam	$10.25$	$1.77$	$7.48$	$8.62$	$10.14$	$11.67$	$14.18$	$1025.30$
RMSProp	$10.03$	$1.93$	$7.49$	$8.29$	$9.53$	$11.44$	$14.49$	$1003.04$
NaG	$9.83$	$1.73$	$7.60$	$8.40$	$9.58$	$10.84$	$14.26$	$983.52$
SGD	$8.15$	$0.87$	$7.46$	$7.65$	$7.82$	$8.43$	$13.37$	$815.31$
SSA1	$9.68$	$1.52$	$7.65$	$8.58$	$9.23$	$10.47$	$13.75$	$968.59$
SSA2	$9.62$	$1.61$	$7.47$	$8.12$	$9.65$	$10.69$	$13.22$	$962.08$
SSA1-Ada	$11.70$	$1.25$	$8.85$	$10.93$	$11.77$	$12.62$	$14.52$	$1170.42$

Table 5. Table 5: CIFAR execution time.

	mean	std	min	25%	50%	75%	max	sum
Adadelta	$104.74$	$0.99$	$102.65$	$104.27$	$104.60$	$105.18$	$110.52$	$10474.36$
Adagrad	$91.82$	$0.42$	$91.39$	$91.47$	$91.72$	$92.07$	$93.11$	$9182.64$
Adam	$97.30$	$0.58$	$96.14$	$96.87$	$97.17$	$97.56$	$98.80$	$9730.21$
RMSProp	$94.52$	$0.50$	$93.11$	$94.17$	$94.40$	$94.83$	$96.21$	$9452.95$
NaG	$91.00$	$0.46$	$90.43$	$90.61$	$90.92$	$91.16$	$92.73$	$9100.56$
SGD	$85.92$	$0.29$	$85.55$	$85.68$	$85.84$	$86.12$	$87.18$	$8592.01$
SSA1	$112.23$	$2.25$	$107.88$	$110.18$	$112.33$	$114.33$	$115.79$	$11223.46$
SSA2	$106.81$	$2.55$	$102.76$	$104.63$	$105.65$	$109.61$	$110.96$	$10681.89$
SSA1-Ada	$179.89$	$4.54$	$168.86$	$176.53$	$179.98$	$183.60$	$186.96$	$17989.26$

Equations

$$u_{n+1}=u_n-h\nabla f(u_n).$$

$$\dot{u}(t)=-\nabla f(u(t)).$$

$$\{\begin{array}{l}{u}_{n+1}=y_n-h\nabla f(u_n)\\ y_n=u_n+\gamma (u_n-u_{n-1}).\end{array}$$

$$\{\begin{array}{l}{u}_{n+1}=y_n-{\beta }_n\nabla f(u_n)\\ y_n=u_n+{\alpha }_n(u_n-u_{n-1}),\end{array}$$

$$\{\begin{array}{l}{u}_{n+1}=y_n-s\nabla f(y_n)\\ y_n=u_n+{\beta }_n(u_n-u_{n-1}).\end{array}$$

$$\ddot{u}(t)+\frac{3}{t}\dot{u}(t)=-\nabla f(u(t)).$$

$$\frac{u_{n+1}-2u_n+u_{n-1}}{h^2}+\frac{1-{\beta }_n}{h}\cdot \frac{u_n-u_{n-1}}{h}=-\nabla f(y_n).$$

$$\frac{1-{\beta }_n}{h}=\frac{3}{t+2h}.$$

$$\underset{h\to 0}{\lim }\frac{1-{\beta }_n}{h}=\frac{3}{t}.$$

$$\{\begin{array}{l}\dot{u}(t)=v(t)\\ \dot{v}(t)=-\frac{3}{t}v(t)-\nabla f(u(t)).\end{array}$$

$$\{\begin{array}{l}\frac{u_{n+1}-u_n}{h}=v_{n+1}\\ \frac{v_{n+1}-v_n}{h}=-\frac{v_n(1-{\beta }_n)}{h}-\nabla f(y_n),\end{array}$$

$$\{\begin{array}{l}{u}_{n+1}=u_n+h{v}_{n+1}\\ v_{n+1}={\beta }_n{v}_n-h\nabla f(y_n).\end{array}$$

$$\ddot{u}(t)+\left[\delta (t)-2\frac{\dot{\delta }(t)}{\delta (t)}\right]\dot{u}(t)=-{\delta }^2(t)\nabla f(u(t)).$$

$$\{\begin{array}{l}\dot{u}(t)=0\\ \dot{v}(t)=-\delta (t)v(t)\end{array}\text{ and }\{\begin{array}{l}\dot{u}(t)={\delta }^2(t)v(t)\\ \dot{v}(t)=-\nabla f(u(t)).\end{array}$$

$$\ddot{u}(t)+\delta (t)\dot{u}(t)=-\nabla f(u(t)).$$

$$\{\begin{array}{l}\dot{u}(t)=0\\ \dot{v}(t)=-\delta (t)v(t)\end{array}\text{ and }\{\begin{array}{l}\dot{u}(t)=v(t)\\ \dot{v}(t)=-\nabla f(u(t)).\end{array}$$

$$\dot{X}(t)=(A+B)X(t),$$

$$\dot{X}(t)=AX(t)\text{ and }\dot{X}(t)=BX(t).$$

$$e^{A+B}=e^A{e}^B,$$

$$\Vert e^{(A+B)h}-e^{Ah}{e}^{Bh}\Vert \le \mathcal{O}(h),$$

$$\ddot{u}(t)+\dot{u}(t)=-\nabla f(u(t)).$$

$$\{\begin{array}{l}\dot{u}(t)=0\\ \dot{v}(t)=-v(t).\end{array}$$

$$\{\begin{array}{l}\dot{u}(t)=v(t)\\ \dot{v}(t)=-\nabla f(u(t)).\end{array}$$

$$\{\begin{array}{l}{u}_{n+1/2}=u_n\\ \frac{v_{n+1/2}-v_n}{h}=-{\beta }_n{v}_n,\end{array}$$

$$\{\begin{array}{l}\frac{u_{n+1}-u_{n+1/2}}{h}={\beta }_n^2{v}_{n+1}\\ \frac{v_{n+1}-v_{n+1/2}}{h}=-\nabla f(y_n).\end{array}$$

$${\hat{v}}_{n+1}=v_{n+1}-h\cdot \left(\frac{1}{{\beta }_n^2}-1\right)\nabla f(y_n).$$

$$\{\begin{array}{l}\frac{u_{n+1}-u_{n+1/2}}{h}={\beta }_n^2{\hat{v}}_{n+1}\\ \frac{v_{n+1}-v_{n+1/2}}{h}=-\nabla f(y_n).\end{array}$$

$$\{\begin{array}{l}{u}_{n+1}=u_n+h{\beta }_n^2(1-h{\beta }_n)v_n-h^2\nabla f(y_n)\\ v_{n+1}=(1-h{\beta }_n)v_n-h\nabla f(y_n),\end{array}$$

$$y_n=u_n+h{\beta }_n{v}_n.$$

$$\{\begin{array}{l}{y}_n=u_n+h{\beta }_n{v}_n\\ v_{n+1}={\beta }_n^k\cdot ((1-h{\beta }_n)v_n-h\nabla f(y_n))\\ u_{n+1}=u_n+{\beta }_n(1-h{\beta }_n)(y_n-u_n)-h^2\nabla f(y_n),\end{array}$$

Full text

1 Introduction

The purpose of the present section is to revisit the basic algorithms that are used for unconstrained optimization problems. We focus our attention on the empirical order of convergence in the setting of deep and convolution neural networks. The rigorous mathematical proofs for the convex and non-convex cases will be given in a follow-up article. From now on, the objective function will be denoted by $f$, the underlying Lipschitz constant by $L$ and the step-size of the discretization scheme by $h$. Firstly, we recall that the best known numerical algorithm used in the setting of unconstrained optimization problems is the gradient descent scheme (see [ Nocedal and Wright (2006) ]):

[TABLE]

This simple algorithm is in fact the explicit forward Euler method and it represents a proper numerical discretization for the following continuous dynamical system:

[TABLE]

On the other hand, Polyak (see [ Polyak (1964) ]) introduced a new type of accelerated gradient-type method based upon the idea of an inertial term. This is a two-step method that is heavily used in optimization applications. Polyak’s momentum algorithm with constant coefficients is the following:

[TABLE]

From [ Ruder (2016) ] we remind that the gradient descent and the stochastic variant of gradient descent, namely SGD (see [ Bottou (2012) ] on Page 2) have trouble navigating landscapes that differ from one dimension to another. The momentum method algorithm 3 is used in order to damp the oscillations of the one-step methods like gradient descent and also in order to prevent the divergence of these type of algorithms. In the optimization literature, the case of the Polyak momentum algorithm with non-constant coefficients is thoroughly studied in the framework of both convex and non-convex minimization problems. This is given as follows:

[TABLE]

where $\alpha_{n}\in[0,1)$ and $\beta_{n}>0$. Quite interestingly, in [ Sun et. al. (2018) ], the authors proved an $\mathcal{O}\left(\tfrac{1}{n}\right)$ theoretical convergence rate, under the assumption that the objective function is coercive, $(\alpha_{n})_{n\in\mathbb{N}}$ is a decreasing sequence and the step-size parameter satisfies $\beta_{n}=\tfrac{2(1-\alpha_{n})c}{L}$, for a fixed element $c$ from $(0,1)$.

Also, the most important optimization algorithm is the inertial algorithm that was developed by Y. Nesterov in [ Nesterov (1983) ]. This accelerated inertial scheme has the following form:

[TABLE]

Moreover, Su, Boyd and Candés [ Su et. al. (2016) ] showed that the underlying dynamical system of 5 is the following second-order differential system, with non-constant gradient damping:

[TABLE]

It is known (see [ Ruder (2016) ]) that 5 is used in order to update the iteration values $u_{n}$ using the inertial coefficients. This will give us a practical way in order to approximate the next values of the parameters, using future updated values. Also, in machine learning applications, the inertial term $\beta_{n}$ is taken to be a constant with an approximate value between $0.5$ and $0.9$.

Further, it is known that 5 exhibits an $\mathcal{O}\left(\tfrac{1}{n^{2}}\right)$ convergence rate, under the strong assumption of convexity of the objective function $f$, when $s\leq\tfrac{1}{L}$. Last but not least, we recall that another competitive forward-backward type algorithm was proposed in [ Boţ et. al. (2016) ] by Boţ, Csetnek and László. It is worth noticing that this algorithm is an extension of the Polyak’s method if the nonsmooth term vanishes, in the setting of the non-convexity of the objective function. Finally, to wrap this introductory section up, we remind that in order to obtain the convergence rates of optimization schemes in the framework of non-convex minimization problems, the regularization function (which is, in fact, a discrete Lyapunov function) must obey the KL property (e.g. see [ Boţ et. al. (2016) ]).

1.1 Some notes on Nesterov’s accelerated method

In this subsection we consider some explanations regarding the asymptotic behavior of Nesterov’s accelerated gradient method. We emphasize the importance of the inertial damping sequence $(\beta_{n})_{n\in\mathbb{N}}$ that converges to $1$ as $n$ goes to $\infty$ and as the underlying step-size is fixed. Moreover, the additional iteration sequence $(y_{n})_{n\in\mathbb{N}}$ that, for every $n\in\mathbb{N}$, is similar to a predictor step in the direction of the critical points, can be written with respect to the discretized velocity. This key idea will be further used in the algorithms that we shall consider.

Remark 1

In 5, the momentum sequence $(\beta_{n})_{n\in\mathbb{N}}$ has the general term defined as $\beta_{n}=\tfrac{n}{n+3}$ or $\beta_{n}=\tfrac{n-1}{n+2}$, for $n\geq 1$ (for simplicity, we consider the last one). The constant term $s$ is taken to be $h^{2}$, where $h$ is the step-size of the numerical discretization. Furthermore, it is known that the 5 algorithm can be considered as a numerical discretization based upon central and forward finite differences of the continuous counterpart 6. In this case, one can write 5 in the following form:

[TABLE]

From a numerical point of view, taking $t=nh$, i.e. $t=n\sqrt{s}$, one obtains that

[TABLE]

Taking the limit as the step-size $h$ goes to [math] for the fixed value of $t$, we obtain that

[TABLE]

We are to highlight that the above limit will play a key role in the development of the splitting optimization algorithms that we will introduce after this section. Moreover, the inertial parameter $\beta_{n}$ converges to $1$ as $n$ goes to $+\infty$, since, with the above consideration, $\beta_{n}=\tfrac{t-h}{t+2h}$, as $n\geq 1$.

We can rewrite the 5 algorithm in a more intuitive form, that is based on the concept of velocity of the underlying second order differential system. One can see that the dynamical system 6 can be written as

[TABLE]

Relating to 7, the Euler-type discretization leads to the following numerical algorithm:

[TABLE]

where the sequence $(y_{n})_{n\in\mathbb{N}}$ will depend on the velocity term $v_{n}$, for $n\in\mathbb{N}$.

After some algebratic manipulations, the 5 algorithm can be written in the following manner:

[TABLE]

In the case of 5, the auxiliary iteration $y_{n}=u_{n}+\beta_{n}(u_{n}-u_{n-1})$ will lead to $y_{n}=u_{n}+h\beta_{n}v_{n}$, and consequently the optimization algorithms 5 and 8 are mathematically equivalent. This means that the momentum sequence $(y_{n})_{n\in\mathbb{N}}$ can be defined as a linear combination of the primary iteration $u_{n}$ and the underpinning velocity $v_{n}$.

1.2 The organization of the paper and the key ingredients of the splitting techniques

In this subsection we briefly present the main ideas concerning the concept of operator splitting. Furthermore, we will explain in detail the role of each section and subsection, respectively. At the same time, we shall infer the role of the splitting methods in optimization problems, where we emphasize the most important mathematical details that are used in the derivation of the formulas concerning the sequential splitting discretizations. The paper is organized as follows: in the present section, namely Section 1, we have presented a brief introduction into the most used optimization algorithms along with their dynamical systems. Furthermore, in Subsection 1.1, we have considered a crucial remark concerned with the asymptotical behavior of Nesterov method 5. This observation is related to the fact that the additional iteration $y_{n}$ can be written as $y_{n}=u_{n}+h\beta_{n}v_{n}$, where $v_{n}$ is the usual velocity at the discrete time $t=nh$. At the same discretization time $t$, the inertial sequence $(\beta_{n})_{n\in\mathbb{N}}$ is given by $\beta_{n}=(t-h)/(t+2h)$, so it converges assymptotically to $1$ as time grows to infinity. This observations will be the main ingredients in constructing sequential splitting-based optimization algorithms. On the other hand, in Section 2, we consider the methods and techniques that represent the backbone of our research. These are based on the idea that an abstract evolution equation of the form $\dot{X}(t)=(A+B)X(t)$ can be decomposed into two subproblems, i.e. $\dot{X}(t)=AX(t)$ and $\dot{X}(t)=BX(t)$, respectively. This leads to the natural idea that, instead of using a numerical method for the full dynamical system, we employ two discretizations, one for each subproblem. So, if $X(t)=(u(t),v(t))$, then at the discrete time $t=nh$, the first discretization operator maps $(u_{n},v_{n})$ into $(u_{n+1/2},v_{n+1/2})$, and the second discretization operator that is in fact related to the second subproblem maps $(u_{n+1/2},v_{n+1/2})$ into $(u_{n+1},v_{n+1})$. It is known that this simple approach of splitting a vector field into two subvector fields is keen on the idea of separating the parts of the dynamical systems that have different physical interpretations. Now, in Subsection 2.1 and in Subsection 2.2 we introduce two completely new optimization algorithms. The key idea is that each of them is a discretization, as time is large enough, to an assymptotical evolution equation in the finite dimensional space $\mathbb{R}^{N}$, namely 11. We shall explain briefly the structure of each of these numerical methods. Regarding Subsection 2.1, for the first splitting algorithm, we consider the following dynamical system:

[TABLE]

Then, using the notation $\dot{u}(t)=\delta^{2}(t)v(t)$, we consider the following attached subproblems:

[TABLE]

Then, using on the first subproblem 14 and on the second one 15, we obtain our first splitting algorithm, but, as we shall see, with some modifications on the velocity. The same technique is applied to a different dynamical system, namely:

[TABLE]

Now, for the full dynamical system we use the notation $\dot{u}(t)=v(t)$. As before, we have two dynamical systems that are represented by the sum of two vector fields, such that this sum is equal to the original vector field of 10. This is equivalent to

[TABLE]

Now, using the numerical discretization 14 on the first subproblem and 20 on the second one, we obtain our second splitting algorithm. Now, we point out that, in our construction, we have used the fact that at $t=nh$, $\beta_{n}=(t-h)/(t+2h)$, where $h>0$ is the chosen step-size. This fact motivates the idea that we can take $\delta(t)=(t-\delta)/(t+2\delta)$, such that it is the continuous counterpart of $\beta_{n}$. Furthermore, choosing $\delta>0$ equal to $h$, we obtain that $\delta(nh)=\beta_{n}$, for each $n\in\mathbb{N}$. On the other hand, as $t\to\infty$, as in the case of the inertial sequence $(\beta_{n})_{n\in\mathbb{N}}$, $\delta(t)\to 1$ as $t\to\infty$ and so, both the dynamical systems 9 and 10 are converging to 11 asymptotically. This is the first motivation why, in Subsection 2.1 and in Subsection 2.2 we have taken 11 to be a common ground to both our optimizers. The second motivation is that 11 is a linear autonomous evolution equation, so in this sense, it is much natural to apply a sequential splitting to this dynamical system.

Now, in Section 3, we have three subsections. In Subsection 3.1, we present the basic concepts of machine learning algorithms. We deal with Adam, Adadelta, Adagrad, SGD, NaG and RMSProp in terms of mini-batches, epochs and iterations. The algorithmic terminology is deeply inspired by Ruder’s work [ Ruder (2016) ]. In Subsection 3.2, we present the structure of the neural networks that are used in the training process on the following datasets: MNIST, MNIST-Fashion and CIFAR-10. In Subsection 3.3, we present our results regarding the comparison of our optimization algorithms with the classical ones belonging to the machine learning literature. Furthermore, we present the simulations concerning the loss, the accuracy and the training time of these optimizers. Finally, in Subsection 4 we point out the novelty, the conclusions and the open problems related to both the theoretical and the computational aspects of our splitting-based numerical methods.

2 Methods

In this section, we present our proposed algorithms that drawn upon the operator splitting technique. The idea of symmetric Strang splitting has its roots back to [ Strang (1968) ] and [ Marchuk (1968) ]. The numerical methods that are based upon the idea of Strang splitting or sequential splitting (see [ Farago (2007) ] on Page 444) are known to be second order and first order, respectively. They are often used in solving ordinary and partial differential equations (see [ Holden et. al. (2010) ]). Briefly, the operator splitting technique is a numerical method for the semi-discretization of a linear system of ODE’s of the following abstract form

[TABLE]

where $A$ and $B$ are two given matrices.

The sequential splitting method is given by the approximation of two linear sub-problems, namely

[TABLE]

The full discretizations of these continuous sub-problems are combined, such that the iteration values for the first discretization is the previous value in the second discretization scheme. Normally, one has that

[TABLE]

if and only if $A$ and $B$ commute. At the same time, in the field of numerical analysis of evolution equations, the major interest lies in finding the theoretical order of convergence of the semi-discretization splitting method, that is buttressed by the evaluation of the absolute error at the first discretization point between the exact and the numerical splitting solution. In the case of Lie (sequential) splitting, this satisfies

[TABLE]

One the other hand, in the framework of optimization problems, one needs to focus the attention on the asymptotic convergence, on infinite intervals, since one is interested in the limit of the solution of the underlying dynamical system. Furthermore, in [ Hansen and Ostermann (2009) ], Hansen et. al. have proved the convergence of the Strang semi-discretization method in the setting of linear evolution equations. The case of semilinear evolution equations can be found in [ Hansen et. al., (2012) ].

We observe that our algorithms start from a dynamical system that is similar to that of Nesterov’s, i.e. 6. But, our arguments are purely formal, in the sense that we adopt an intuitive way to infer a set of splitting-based numerical algorithms. We use some algebraic manipulations that will lead to some optimization schemes, using the idea of velocity given a second-order dynamical system. Firstly, we shall consider some remarks about the Nesterov’s dynamical system 6, that will be crucial for the constant-damping dynamical systems of our optimization schemes.

2.1 Sequential splitting algorithm I

We consider the following dynamical system with constant damping:

[TABLE]

In the vein of Nesterov’s algorithm 5, one can easily apply a forward Euler or a Crank-Nicolson type method to 11, such that it preserves the properties of the continuous dynamical system. Our approach is rather empirical: we split the continuous dynamical system into two sub-problems and discretize each of them with some basic numerical methods. After that, we merge these two discretization into one algorithm. This numerical procedure will lead to different velocity updates for the current iteration value.

Now, the first subproblem will contain the linear part of the velocity updates and is the most basic one:

[TABLE]

Further, the second subproblem will contain the continuous gradient of the objective function $f$, as follows:

[TABLE]

On the other hand, for the subproblem 12, we employ the explicit forward Euler discretization, as follows

[TABLE]

where $u_{n+1/2}$ is the iteration term for the discretization of the continuous differential system 12. Further, we know that $\beta_{n}$ converges to $1$ as $n$ goes to $+\infty$, so 14 is a proper numerical discretization, as time grows to infinity.

At the same time, for the second continuous subproblem 13, we consider the symplectic Euler numerical algorithm for the iteration value $u_{n}$ and for the velocity updates $v_{n}$, in the following way:

[TABLE]

Now, we introduce $\hat{v}_{n+1}$ which is a perturbed velocity under the value of the gradient. That is

[TABLE]

So, we modify our algorithm on the second sub-problem 13 in the following manner:

[TABLE]

If we combine 16 with 14, after some algebraic manipulations, we will eventually obtain the following algorithm:

[TABLE]

where the additional iteration $y_{n}$ will depend on $u_{n}$ and $v_{n}$. Now, the most intuitive approach for our discretization is to take $\hat{y}_{n}=u_{n}+h\beta_{n}^{2}(1-h\beta_{n})v_{n}$, in order to have $u_{n+1}=\hat{y}_{n}-h^{2}\nabla f(y_{n})$, which is similar to the case of Nesterov’s algorithm 5. From the preliminary results of the numerical simulations we infer that a better approach would be to take the same value for the momentum iteration updates $y_{n}$ as in the case of 5, i.e.

[TABLE]

Finally, we consider the modified form of the algorithm 17, using the addition of a hyper-parameter $k$ that will boost the velocity values in order to reach the minimum of the objective function at least as fast as 5:

[TABLE]

where the hyper-parameter $k$ is chosen such as $k\geq 0$.

2.2 Sequential splitting algorithm II

In the case of the second algorithm we propose, we also employ the sequential splitting. We consider the dynamical system 11, which we divide into the same two sub-problems as before, that is 12 and 13. For the first sub-problem, we employ the same discretization, namely the forward Euler method 14. On the other hand, for the second sub-problem we consider the following discrete counterpart, namely also the explicit Euler method:

[TABLE]

In this case, we observe that for the principal iterate values $u_{n+1}$, we have used the previous value updates of the velocity, namely $v_{n+1/2}$. In this sense, our sequential splitting algorithm is of backward type. Merging 14 and 20, we obtain the following optimization algorithm that is based upon the sequential splitting technique:

[TABLE]

where, as before, $y_{n}$ will depend on $u_{n}$ and on the velocity values $v_{n}$. One can notice that the natural choice for $y_{n}$ is $u_{n+1}$. As in 17, we shall opt for $y_{n}=u_{n}+h\beta_{n}v_{n}$. Subsequently, one obtains that the value of $hv_{n}$ is equal to $(y_{n}-u_{n})\cdot\dfrac{1}{\beta_{n}}$.

Taking these into account, we consider the following inertial algorithm that will be appropriate for unconstrained optimization problems:

[TABLE]

Also, in order to improve the empirical rate of convergence of our algorithm, we will add one hyper-parameter $k\geq 0$. From our numerical computations we set the default value of $k$ equal to $2.0$. So our proposed algorithm is the following:

[TABLE]

Finally, we observe that our new algorithm 22 differs from the first splitting algorithm 19 in the iterates values $u_{n+1}$, since the latter algorithm does not contain the gradient of the objective function $f$. In addition to this, the velocity values $v_{n}$ are updated with different parameters. In the framework of neural network training they depend on the current iteration, where the gradient is computed only on a mini-batch dataset.

Remark 2

If we compare our optimization schemes 19 and 22, we can see that the latter one is a proper discretization scheme of the associated asymptotic dynamical system. Further, the first splitting based algorithm has a different asymptotic dynamical system, by the fact that the velocity used in the second sub-problem 13 is perturbed with the value of the gradient in the point $u_{n}+h\beta_{n}v_{n}$.

3 Results

3.1 A gentle introduction to machine learning algorithms

In this section, we shall review some optimization algorithms that are often used in the training of deep and convolution neural networks. In our neural network training, we use only the stochastic variant of these algorithms, including epoch-training with mini-batches. Furthermore, we shall present algorithmically our optimization schemes 19 and 22, that will be also used with their stochastic counterpart. From now on, the first splitting-based algorithm 19 will be called SSA1 .

In the neural network training, the objective (also called loss function) must be minimized with respect to some parameters. In general, these parameters contain the weights and biases in the neural network. From now on, we consider $\theta$ to be the collection of weights and biases and $(x^{(i)},y^{(i)})$ to be the given training input values and target values, respectively.

In the first section, we have introduced the most basic optimization algorithm 1, namely gradient descent scheme. More often it is used the stochastic variant of this algorithm, namely stochastic gradient descent (see [ Bottou (2012) ] on Page 2 ). In the full version of this algorithm, namely mini-batch stochastic gradient descent, one updates the parameter iteration values for each shuffled mini-batch $m$ of training examples is given as follows

[TABLE]

This algorithm leads to a more stable convergence than plain 1, but updates the parameters with high variance (see [ Ruder (2016) ]) and will eventually lead to oscillations in the decrease of the objective function. The advantage is that at each epoch we dispose of redundant data input values, since at each iteration we compute the gradient on a mini-batch of a given size $m$. Despite the clear disadvantage of using the same learning rate (step-size) $h$ for every component of the vector $\theta_{n+1}$, the algorithm has lower chances than plain 1 to be stuck at a saddle point of the objective function.

Other iteration schemes that are deployed in neural network training are the adaptive algorithms. The most basic one is Adagrad (see [ Duchi et. al. (2011) ] on Page 2122). This is used for large-scale neural networks, but has the disadvantage that it has a very aggressive learning rate update. So, a modification of this algorithm is Adadelta (see [ Zeiler (2012) ] on Page 3), in which one stores previous values of the squared gradient of the objective function. Further, it is known that decreasing the effect of the past gradients lead to less sensitivity in choosing the hyper-parameters of the neural model. In Adadelta, the so-called running average is defined as

[TABLE]

where $\gamma$ is taken, in general, 0.9. Also, for each epoch and for a chosen mini-batch $(x^{(i:i+m)},y^{(i:i+m)})$ of a given suitable size $m$, another exponential decay average is defined, such that

[TABLE]

with the remark that the update values are given by the negative ratio of two adaptive rates. The coefficient from the nominator contains the square root of the exponential energy decay and the denominator is the square root of the past gradients. For more details, we let the reader follow [ Zeiler (2012) ].

Another popular algorithm in neural network training is RMSprop which was introduced in the Coursera lecture class by G. Hinton [ Tieleman and Hinton (2012) ] on Slide 26. Following [ Ruder (2016) ], RMSprop algorithm is in fact the same as Adadelta’s vector update, using the adaptive learning rate of the form

[TABLE]

In general, the default value for the starting learning rate is of order of magnitude $10^{-3}$. One can observe that in the case of RMSprop algorithm, the effective value for the adaptive learning rate depends on the magnitude of the squared energy of the past gradient values, for each vector component of the underlying parameter.

Last but not least, we recall the optimization algorithm Adam (see [ Kingma and Ba (2014) ] on Page 2), that can be considered as a modified gradient descent using two biased corrected moments, i.e.

[TABLE]

where, as before, $m$ is the size of the mini-batch used in the current iteration $n$. So, Adam is a modified adaptive gradient descent type-scheme using two moments that are corrected in each iteration. Also, in the machine learning community, Adam seems to be the most popular optimization algorithm and one of the most easy optimizers that require minimal tuning of the parameters (see [ Karpathy (2017) ]). From [ Dozat (2016) ] Page 3, we recall the Nadam algorithm, which basically represents the combination between Adam and the momentum method. In [ Sutskever et. al. (2013) ], it is mentioned that this algorithm is comparable to first order adaptive methods. The inertial parameter in the momentum method leads to higher quality updates of the weights and biases values, where in one updates the parameters after the momentum correction of the past gradients. Last but not least, we reassert that instead of using adaptive learning rate algorithms, one can use 23 using learning rate updates that depend on the current epoch training(for more details, follow [ Smith (2015) ]). For some references regarding the choice of the optimal learning rate, the optimizers for neural networks and for an empirical analysis of the training and classification problems, we refer to [ Bengio (2012) ], [ Bottou et. al. (2018) ], [ Zhang et. al. (2016) ] and [ Bishop (1995) ].

In classification-type tasks, stochastic gradient descent 23 is better the generalization errors obtained on the validation and test datasets than adaptive learning rate algorithms. For a full discussion of this problem, we refer to [ Wilson et. al. (2017) ]. Also, for some general notions concerning deep learning principles, classification and optimization problems we refer to [ Nielsen (2015) ], [ Mehta et. al. (2018) ], [ Higham and Higham (2018) ], [ LeCun et. al. (1998) ] and [ Goodfellow et. al. (2016) ]. Now, at the end of this section we consider the stochastic version of the SSA1 optimizer, that is 19. This is presented below, as Algorithm 1. Moreover, in our algorithms the current iteration is denoted by $n$ and the learning rate by $h$.

At the end of this sub-section, we present the mini-batch stochastic algorithms that are similar to the adaptive Adadelta optimizer, namely Algorithm 3. In the numerical computations, we shall briefly call it SSA1-Ada .

Remark 3

From our numerical results we have observed that the adaptive counterpart of SSA 1 algorithm, namely SSA 1-Ada is competitive with all other adaptive-type schemes. A similar adaptive algorithm can be given for SSA 2 . Further, also from our computations we have observed that this algorithm reaches almost $10\%$ accuracy on both MNIST and MNIST-Fashion. So, our results from the last sub-section does not include the adaptive mini-batch variant of our second proposed optimizer.

3.2 Description of the neural network

For our numerical simulations we have used some adaptive and non-adaptive stochastic algorithms with mini-batches in order to minimize the loss function of a convolution neural network. We have used both well known Hand Written digit Recognition Dataset (MNIST) from http://yann.lecun.com/exdb/mnist/ and the harder drop-in replacement MNIST-Fashion presented in https://github.com/zalandoresearch/fashion-mnist. For the neural network we have used the MNIST example provided by Pytorch at https://github.com/pytorch/examples/tree/master/mnist. The CNN has a convolution layer with a kernel size of 5 pixels which converts the single channel input to a 20 channel output. After RELU is applied to the output of this layer, it is max pooled with a kernel of 2 pixels and is passed as the input to another convolution layer wit the same kernel size, but whose number of output channels is 50. RELU and the same max pooling technique are applied again. This output is flattened and then is connected to a layer of 500 hidden neurons. The output of this hidden layer is passed through another RELU and then is connected to a 10 neuron output layer. The final output is obtained via the logsoftmax activation function of the output layer activation values, given by the formula

[TABLE]

where $\bf{x}^{[i]}$ represent the $i^{th}$ component of the entry vector $\bf{x}$.

Also, this model does not use any dropout and for the loss function we have used the Negative Log Likelihood given by the formula

[TABLE]

Furthermore, we have adapted the code provided by Pytorch’s example in order to facilitate the choice of different optimizing algorithms and different loss functions. The new algorithms we introduced were implemented as classes in the Pytorch optim package, using a coding style as similar as possible to the style used by the Pytorch developers. For the experiments we have employed the standard split for both MNIST and MNIST-Fashion, so the 60000 images were divided in two sets, one with 50000 images for training and one with 10000 images for testing. We did not use validation data because we desired to mimic the results obtained by the example provided by Pytorch. The data sets used were loaded via the Pytorch loaders and were normalized with mean $0.1307$ and standard deviation $0.3081$. We divided both the training and test data in mini-batches and then we did our training in a stochastic manner. For the shuffling of the mini-batches we used the random seed 1. We have also used the seed as Pytorch’s seed for the random initialization of the weights and biases. We saved the model that we started the experiments using Pytorch’s built-in save feature.

On the other hand, CIFAR 10 (https://www.cs.toronto.edu/~kriz/cifar.html) is a dataset of images from 10 different categories. There are 60000 color square images measuring 32 by 32 pixels. The dataset was proposed by the Canadian Institute for Advanced Research. In our computations, we have used the GoogLeNet model proposed in [ Szegedy et. al. (2015) ] by Szegedy et. al. The model contains 22 inception layers and has achieved an error rate of 6.67% in the top 5 categories. We chose this model because it is circumscribed by the 1.5 billion add-multiply operations budget, which in turn means a faster, shorter training time for our optimizer. The model’s implementation is available on GitHub at https://github.com/kuangliu/pytorch-cifar.git. In order to load the data, we have used the methods provided by Pytorch datasets package. On top of this, we emphasize the fact that our results on CIFAR 10 lie upon the use of the Cross Entropy Loss, which is a combination between Negative Log Likelihood loss and softmax activation function.

Moreover, for the simulation that entails CIFAR 10, we employed the following learning rates: for the adaptive algorithms Adadelta, Algorithm 3, we took the default learning rate 1.0. Also, for the other optimizers, i.e. RMSProp, Adagrad , Adam , 23, 5, Algorithm 1 and Algorithm 2, we took the learning rate 0.001, since this represents a good value for 200 epochs of training.

Last but not least, we emphasize the fact that our numerical computations were made on a GPU:NVIDIA Tesla V100-SXM2, which is a GPU data center that has 5120 NVIDIA Cuda cores, GPU memory 16 GB HBM2 and its double precision performance is around 7.5 TFlops. In a nutshell, we point out that in Keras and Tensorflow, the stochastic algorithms with mini-batches present heavy oscillations in the decrease of the loss function, since the graphical representations of the loss function are made with respect to each iteration concerning mini-batches. In our case, the plots are given in terms of accuracy and loss values at the end of each epoch, and not in terms of each particular iteration at different mini-batch datasets.

3.3 The neural network experiments

In Figure 1(a), we plotted the overall accuracy on the MNIST test dataset. The learning rate of 23, 5, Algorithm 1 and Algorithm 2 has been set to a high value, i.e. $0.1$ and the graph represents the accuracy values on 50 epochs. Our algorithms are comparable to the stochastic version of the stochastic gradient descent and with the Nesterov’s algorithm. On this graph we observe that our optimizers have a better accuracy on the MNIST dataset, but present higher oscillations due to their velocity updates that are based on the constant parameters $k$ and $q$.

Furthermore, since the decrease in the loss function presents the same behavior as in the increase of the accuracy from Figure 1(b), we have presented the decrease in the objective function. We have stuck to the same 50 epochs an d 0.1 learning rate of the chosen algorithms. We can observe that the loss function of the 23 showcases high oscillations over the first 20 epochs. Due to the inherent inertial moment, 5 alleviates these oscillations and stabilizes the decrease of the objective function. Conversely, the loss value of these algorithms increases after 20 epochs and this can be prevented with some early-stopping techniques. Last but not least, SSA 1 and SSA 2 present a similar value with respect to the loss function as 5, on the test dataset, and this is in correlation with the increased accuracy observed in the Figure 1(a).

In Figure 2, we considered the graph of the accuracy values of the same algorithms as before, but for the MNIST-Fashion test dataset. It is known that on the MNIST-Fashion dataset, the accuracy is not as high as in the classical MNIST dataset. Quite surprisingly, our newly introduced SSA 1 and SSA 2 seem to behave better than what we have seen in Figure 1(a). Due to their velocity updates in the gradient values, especially in the sequential splitting Algorithm 1, these optimizers do not oscillate with respect to the increase in epochs. While this stands true for our splitting-based schemes, 23 and 5 present powerful degradation in the accuracy on the same test dataset. Also, a similar analysis can be made for the decrease in the loss function, as in Figure 1(b).

In Figures 1(a), 1(b) and 2, we have the behavior of our stochastic non-adaptive algorithms in the presence of the training of the convolution neural network. Further, SSA 1 and SSA 2 optimizers have some additional hyper-parameter $k$. This is clearly an advantage in the velocity updates of these discretizations. On the other hand, for more challenging classification problems, one needs to choose the optimal values for these parameters. Even though it is of good practice to employ in the numerical simulations a Bayesian optimization technique, it is easily verifiable that it is time consuming for these two parameters, in addition with the learning rate values. Furthermore, on more complex classification problems we have observed, through our numerical simulations, that moderate values of both $k$ , i.e. $k=2.0$ can be considered as default values for this optimization schemes, yielding to a good convergence behavior of our proposed algorithms.

So far we have put an emphasis on our numerical results on some basic experiments. All of our computations are presented in Table 1 for both MNIST and MNIST-Fashion datasets. Also, since the phenomenon of overfitting is crucial in machine learning, we have illustrated the accuracy on the training set (lower diagonal) and the accuracy on the test set (upper diagonal). For our non-adaptive optimizers we have set the step-size to vary from $0.001$ to $0.1$ and for our adaptive counterparts we have set the learning rate to the default value of $1$, since SSA 1-Ada is a combination between Adadelta, Algorithm 1. We focus our attention on a few remarks regarding our results for the accuracy values on the MNIST dataset. Both stochastic gradient descent and Nesterov’s algorithm with momentum 0.5 have almost 0.99 accuracy on the test set, when the learning rate is equal to 0.1. Furthermore, SSA 1 and SSA 2 display almost the same behavior when the step-size is large enough, i.e. they have almost 0.992 and 0.993 accuracy value on the test set. On the other hand, a lower learning rate of approximate value of 0.001 leads to a difference in the train and test accuracy for the splitting algorithms, yet for 23 and 5 it does not lead to the overfitting problem when lower step-size values are used. In addition, it seems that 5 achieves better accuracy than both SGD and SSA 1 after 20 epochs until this value is stabilized. For SSA 1 and SSA 2 there seems to be a negligible difference in training versus testing accuracy after the first 50 epochs. Now, regarding our adaptive algorithms, for the MNIST dataset, the Adadelta algorithm achieves a very good accuracy before 20 epochs of learning. The same stands true for SSA 1-Ada , but these algorithms degrade in their accuracy and the need for early-stopping is a serious issue. Only 20 epochs are needed in order to train our adaptive optimizers. In light of stabilizing the decrease in the loss function, one can decrease after a number of epochs the value of the velocity parameter $k$ down to 0.1 in order to prevent high values in the objective function.

Now, we turn our attention to the results on the MNIST-Fashion dataset (the results are in the right part of the Table 1). We can observe that 5 with inertial parameter 0.5 achieves almost 0.875 accuracy value, taking into account that the step-size was chosen as 0.1. Further, SSA 1 and SSA 2 both achieve circa 0.92 accuracy on this test set. For lower values of the learning rate 5 is better, but it seems that the hyper-parameters $k$ and $q$ compensate well into the velocity updates in order to achieve almost the same values of accuracy. The advantage of 23 and 5 is that on the 20 epochs they achieve a lower value in the loss function, but after 50 the difference is negligible. Quite interestingly, in the case of the non-adaptive optimizers, Algorithm 3 does not degrade after 20 epochs, even though MNIST-Fashion is harder to classify in comparison with the classical MNIST dataset. On the other hand, our adaptive algorithms seem to provide lower accuracy values at the end of the $100^{th}$ epoch and, in order to compensate, one can tune parameter $k$ to a greater value that will boost the value of the iterations. Finally, for RMSProp and Adam, we infer that they are not really suited for a step-size greater than $0.01$. Their advantage is that they do not overfit, which can be noticed from Table 1, where the difference between the accuracy on training and test datasets appears imperceptible.

Concerning Figure 3(a) and Figure 3(b), we considered the comparison on CIFAR 10 test dataset for different optimization algorithms, using 80 epochs and 150 epochs, respectively. It can be derived that Algorithm 1 does not present oscillations but they need more epochs in order to converge to a local minimum of the Cross Entropy Loss function. Much better results are obtained by using adaptive variants, due to the tuning of the learning rate at each iteration on the mini-batches. It is worth noticing that the other adaptive algorithms RMSProp and Adam present natural oscillations determined by their aggressive step-size tuning. On the other hand, we recall that these results depend heavily on the randomization of weights and biases, since the local minimum of the loss function at which they converge depends on each and every optimizer. The full results are given in Table 2. We ran our algorithms to 150 epochs and we have chosen the optimal step-size for each algorithm. SSA 1 and SSA 2 are both comparable to 5 (with momentum 0.9), but in the first 20 epochs they present lower values in accuracy. On the other hand, RMSProp is comparable with Adadelta and Adagrad, but the latter one reaches lower values of accuracy after 20 epochs, due to the aggressive nature of the learning-rate tuning. Last but not least, we mention that in the last two lines of Table 2, SSA 1-Ada yield better results compared to all the algorithms we have presented so far. Needless to say, a rigorous analysis of these optimization schemes entails numerous simulations, with different random values for the weights and biases of the neural network, in order to reach an average value for the accuracy on the test set. As far our splitting-based schemes go, the hyper-parameter $k$ play a major part, as in they must tuned by way of the validation set. This aims at achieving the lowest possible values in the chosen loss function.

Remark 4

It is interesting to notice that in Table 1 and in Table 2 we have the values for the accuracy on the train and test datasets for different optimizers. For our simulations we have used a fixed random seed. For practical purposes, one must generate different simulations on these datasets, since the convergence of these algorithms depend heavily on the initial random values of the weights and biases in the neural network.

Remark 5

In the 5 implementation from Pytorch, the momentum coefficient $\beta_{n}$ is set to a fixed value between $0.5$ and $0.9$. Thus, for a smooth implementation of our optimization schemes, one can also use a fixed inertial coefficient, instead of the non-constant sequence, defined as $\beta_{n}=\tfrac{n}{n+3}$ or $\beta_{n}=\tfrac{n-1}{n+2}$ (see Remark 1), where $n$ represents the iteration employed in the computation of the gradient on the current mini-batch.

Remark 6

The algorithms SSA 1 , SSA 2 and SSA 1-Ada are based upon the hyper-parameter $k$, with the default value of 2.0. The increase of the loss function values and the decrease of the accuracy on the test dataset, call for the tuning $k$ with respect to the current iteration, where the gradient of the mini-batches is computed. For example, if the accuracy does not change from epoch to epoch, then one can boost the value of $k$ from 2.0 to 10.0. On the flip side, if the chances of overfitting are high, then one can tune $k$ from 2.0 down to 0.1. The parameter $k$ can be non-constant, i.e. $k=k(n)$. If the algorithm indeed overfits, one can set it to decay exponentially, such as $k(n)=e^{k(-1/n)}$, where $k$ stands for the default value.

At last, we present a crucial remark concerning the adaptive counterpart of the first splitting algorithm.

Remark 7

A method for improving SSA 1-Ada , given by Algorithm 3 is to use the current value of the gradient $\nabla f_{\theta}(z_{n+1};x^{(i:i+m)},y^{(i:i+m)})$ and the squared gradient $\nabla^{2}f_{\theta}(z_{n+1};x^{(i:i+m)},y^{(i:i+m)})$ at the intermediate point $z_{n+1}=\theta_{n}+h\beta_{n}v_{n}$. A simple way for implementing this is to compute the additional iteration $z_{n+1}$ before the accumulation of gradient and before the computation updates. Simply put, this means moving line 18 after line 9 in the algorithm and then use the gradient and the squared gradient in this newly updated iteration.

Now, before our next remark consisting of some explanations behind the training time of the neural network, we emphasize the fact that SSA1-Ada can be implemented as in Algorithm 3 or as specified in Remark 7. The main idea is that different constructions of the adaptive splitting-based algorithms represent only a guideline for deep neural network implementations of SSA1-Ada .

Remark 8

Now, we present some facts about the training time regarding the machine learning algorithms. As above, we have considered the three datasets, namely MNIST, MNIST-Fashion and CIFAR-10. On each of them we have compared the time to train the neural network with the same optimizers as before: Adadelta, Adam, RMSProp, 5, SGD, SSA1 , SSA2 and SSA 1-Ada . The results are presented in Table 3, Table 4 and Table 5, respectively. The procedure is in a sense statistical: we have stored the training time on 100 epochs and then we have constructed an array consisting of the 100 simulations. As an implementation we have used the classical Pandas toolbox from Python, in which we have stored the mean, the standard deviation, the minimum value, the maximum value, the quantiles and the sum of all these values, for the array containing the results of the training time corresponding to 100 epochs. Now, we will briefly present the comparative results for each of the datasets. In Table 3 we observe that SSA1 and SSA2 present an mean value close to the other algorithms. The standard deviation is, on average, as the same as Adadelta and RMSProp. In this case only SGD has a much lower value in the standar deviation, suggesting that is more stable than the other optimizers. The maximum values suggests also that our algorithms stand close to Adadelta, Adam and 5. On the other hand, we point out that adaptive SSA1-Ada has values that are close enough to Adadelta. This is linked to the fact that both Adadelta and our algorithm are adaptive counterparts of inertial optimizers. Now, we turn our attention to Table 4. In the case of the MNIST-Fashion dataset, the mean of the algorithms is lower in comparison with the case of the classical MNIST. Furthermore, the behavior of SSA1 and SSA2 is similar to the previous experiment in the case of the mean and of the standard deviation. This is connected to the fact that the same network architecture was used in the case of both datasets. Finally, we will present the computational results regarding Table 5. In the case of complex architectures like that used in CIFAR-10 dataset, our splitting algorithms have a high standard deviation, which is in resonance with the instability regarding the training time over a long period of epochs. Moreover, the mean of SSA1 and SSA2 are close to the values of Adadelta and Adagrad. Also, it seems that 5 and SGD are faster than our splitting algorithms. Now, the adaptive SSA1-Ada requires much more training time and has a higher variance in comparison with the other stochastic algorithms. Even though this algorithm requires a greater training time, the advantage is that it compensates in accuracy and it does not lead to overfitting (see the previous tables). Finally, we point out that, in order to make SSA1-Ada much faster, one can employ the Remark 7, where we have discussed a different implementation. This will lead to a faster algorithm that encompasses both the inertial effects and also the stability in the loss function.

Remark 9

Nesterov algorithm 5 developed in [ Nesterov (1983) ] has been analyzed theoretically in different areas of the optimization community. The algorithm is based upon the inertial coefficient given as $\beta_{n}=\tfrac{n-1}{n+2}$ or $\beta_{n}=\tfrac{n}{n+3}$, as in Remark 1. From a practical point of view, in PyTorch 111 https://pytorch.org/docs/stable/optim.html and in Tensorflow 222 https://www.tensorflow.org/api_docs/python/tf/compat/v1/train/MomentumOptimizer , the momentum coefficient $\beta_{n}$ is considered as a constant, with values between $0.5$ and $0.9$. The implementation of the momentum (see [ Sutskever et. al. (2013) ]) as a constant term $\beta_{n}=\beta$ gives a boost in the empirical rate of convergence, due to the fact that $\lim\limits_{n\to\infty}\beta_{n}=1$. Hence, we have the natural theoretical interpretation of 5 with $\beta_{n}=\tfrac{n-1}{n+2}$ or $\beta_{n}=\tfrac{n}{n+3}$ and the practical version, when $\beta_{n}=\beta$. We gave a theoretical motivation for both our splitting-based algorithms with a nonconstant inertial sequence $(\beta_{n})_{n\in\mathbb{N}}$ and we have compared with adaptive or non-adaptive algorithms (e.g. 5 and 23). Now, taking into account the above explanation, we have considered in Figure 4(a) and in Figure 4(b), some numerical experiments on the CIFAR 10 test dataset, over 150 epochs. We have considered the following algorithms: Adadelta (with constant learning $1.0$), Adam and RMSprop algorithms (with learning rate $0.003$), SGD and Nesterov (with learning rate $0.5$) and both SSA1 and SSA2 , with $0.001$. Furthermore, the constant momentum term of Nesterov was set to $0.5$ as in the case of the constant variants of our algorithms (also for these algorithms we have set the learning rate to $0.5$). We observe, that Adam and RMSprop are unstable with respect to the learning rate due to their inherent adaptativity. On the other hand, we have chosen a larger learning rate for the constant momentum variants of both SSA1 and SSA2 , in order to show the potential smoothing effect available also for high values of the learning rate, with respect to the optimizers oscillations (due to the stochasticity).

We show that, similar to the computational implementation of 5, the constant momentum algorithms SSA 1 and SSA2 admit a better rate of convergence when the momentum coefficient is $\beta_{n}=\beta\in[0.5,0.9]$. Furthermore, we can observe that our algorithms have a faster rate of convergence compared with all the algorithms we have considered so far (these algorithms have the name SSA1 CONST and SSA2 CONST , respectively), in the sense that we have obtained a better accuracy. The comparison reveals two computational aspects:

As in the case of 5, we have constructed a non-adaptive optimizer with a momentum term $\beta_{n}=\tfrac{n-1}{n+2}$ or $\beta_{n}=\tfrac{n}{n+3}$. From a theoretical point of view, our algorithms resemble 23 and 5 where $\beta_{n}$ depends on the iteration values, but from a practical point of view, we have considered that a fair comparison needs to be made when $\beta_{n}=\beta$ with values between $0.5$ and $0.9$ for both SSA1 and SSA2 , and also for SGD and Nesterov accelerated gradient method (see also the implementation of momentum method in PyTorch as reminded above). In this case, obviously the underlying dynamical systems will simplify and will contain a constant damping term, but our purpose of the present remark is to give a fundamental approach that gives us a good and practical behavior related to speeding up the empirical convergence rate of the optimizers.

Last but not least, we mention that our algorithms have a similar property to non-adaptive optimizers like SGD and 5, in the sense that they generalize better than adaptive algorithms (see for example [ Wilson et. al. (2017) ]). This can be seen in Figure 4(a) and in Figure 4(b), respectively. We mention that although the property that we have mentioned above is present for both inertial algorithms and also for SSA1 and SSA2 , our algorithms are not inertial optimizers, but they are multi-step type methods since they contain both $y_{n}$ and $y_{n-1}$. For example, 19 can be written as follows:

[TABLE]

Similarly, 22 can be reduced to the following alternate form:

[TABLE]

Finally, we emphasize the fact that one can reach the above algorithms (equivalent with SSA1 and SSA2 ), using the value of discrete velocity, i.e. this consists in expressing $v_{n+1}$ in terms of $u_{n+1}$, $u_{n}$, $v_{n}$ and $y_{n}$.

4 Discussions

In this paper we have introduced new types of optimization algorithms that are competitive in the neural network training with the momentum-based algorithms and with the numerical schemes. These are based upon the idea of adaptive learning rates. The approach of using the operator splitting technique is new in the field of machine learning and one can see that it is an efficient way of developing inertial optimizers. Concerning inertial algorithms, we follow [Wilson et. al., 2017], in which one infers that SGD and 5 generalize much better than the adaptive algorithms like Adam and Adagrad. Also, even when these adaptive algorithms achieve the same or better training accuracy, the non-adaptive ones give, in the long run, much better accuracy values on the test data. We have observed empirically that different types of inertial algorithms, as in the case of SSA1 and SSA2 generalize better than the adaptive optimization methods. At the same time, in [ Keskar and Socher (2017) ], the authors have developed a method in order to make Adam generalize as good as SGD, that is based on a triggering condition and this represents a key idea in order to make adaptive algorithms to share the same behavior as the inertial ones. This is the first motivation to the fact that we have emphasized the key role of inertial algorithms that are closely related, but ultimately different, than 5. On the other hand, we point out that, in the optimization community, the sequential splitting method has drawn attention due to the recent article of M. Muehlebach and M.I. Jordan, [ Muehlebach and Jordan (2019) ]. The authors have showed that Nesterov’s method 5 is a combination between the sequential splitting method and a symplectic Euler method. Recently, the authors of [ França et. al. (2019) ] have also used splitting techniques of ODE’s, but to proximal-type algorithms. This differs from other papers from the fact that they have used balanced and rebalanced splitting methods that are in deep connection with the equilibrium states of the ODE’s subproblems. This represents our second motivation, in the sense that it is a strong argument that operator splitting techniques are not only an elegant approach, but also it gives some insights concerning the evolution equations associated to the discretization methods. Elseways, the comparison of our splitting-based algorithms were done with some classical optimizers, like in the paper [ Dozat (2016) ], [ Kingma and Ba (2014) ], [ Wilson et. al. (2017) ], [ Zeiler (2012) ] and [ Sutskever et. al. (2013) ]. Our full implementation is based on Python package PyTorch. In a future work, we aim to add our algorithms to Keras and Tensorflow. We consider adding these algorithms to Autograd, so as to be used as optimizers for particular objective functions, which, in their turn, can be deployed for the approximation of solutions of multi-dimensional initial value problems and boundary value problems. We also intend to present mathematically rigorous proof for the convergence rate of these optimization algorithms for convex and non-convex objective functions. For the convex functions, our aim is to employ the technique concerning suitable Lyapunov functions. Such work is underpinned by [ Da Silva and Gazeau (2018) ], where Lyapunov functions were used for the associated dynamical systems of some adaptive optimizers. To wrap things up, non-convex optimization functions call for the use of the KL property of the underlying regularization of the objective function, as in [ Boţ et. al. (2016) ].

Acknowledgments

We would like to thank DirectMailers for allowing us to use their NVIDIA DGX-1 server with eight Tesla V100-SXM2 GPUs, which enabled us to do our simulations for the training of the neural network. Furthermore, we give our appreciation to Luigi Malagò from Romanian Institute of Science and Technology (RIST) from Cluj-Napoca for some valuable suggestions regarding the effect of the learning rate with respect to the splitting-type optimization algorithms. Finally, we give our appreciation to anonymous reviewers who have helped us organize the present research article, since their comments improved greatly the experimental part of our research.