# Rethinking ResNets: Improved Stacking Strategies With High Order Schemes

Zhengbo Luo<EMAIL_ADDRESS>Zitang Sun<EMAIL_ADDRESS>Weilian Zhou<EMAIL_ADDRESS>Zizhang Wu
<EMAIL_ADDRESS>Sei-ichiro Kamata<EMAIL_ADDRESS>Graduate School of
IPS, Waseda University, Kitakyushu City, Fukuoka, Japan  Zongmu Technology,
Shanghai, China

###### Abstract

Various deep neural network architectures (DNNs) maintain massive vital
records in computer vision. While drawing attention worldwide, the design of
the overall structure lacks general guidance. Based on the relationship
between DNN design and numerical differential equations, we performed a fair
comparison of the residual design with higher-order perspectives. We show that
the widely used DNN design strategy, constantly stacking a small design
(usually 2-3 layers), could be easily improved, supported by solid theoretical
knowledge and with no extra parameters needed. We reorganise the residual
design in higher-order ways, which is inspired by the observation that many
effective networks can be interpreted as different numerical discretisations
of differential equations. The design of ResNet follows a relatively simple
scheme, which is Euler forward; however, the situation becomes complicated
rapidly while stacking. We suppose that stacked ResNet is somehow equalled to
a higher-order scheme; then, the current method of forwarding propagation
might be relatively weak compared with a typical high-order method such as
Runge-Kutta. We propose HO-ResNet to verify the hypothesis of widely used CV
benchmarks with sufficient experiments. Stable and noticeable increases in
performance are observed, and convergence and robustness are also improved.
Our stacking strategy improved ResNet-30 by 2.15 per cent and ResNet-58 by
2.35 per cent on CIFAR-10, with the same settings and parameters. The proposed
strategy is fundamental and theoretical and can therefore be applied to any
network as a general guideline.

###### keywords:

Deep neural networks , Neural Ordinary Differential Equations , Image
processing

††journal: Neurocomputing

##  1 Introduction

Deep neural networks (DNNs) have achieved many exemplary breakthroughs in
computer vision, image processing, and signal processing with powerful
learning representations from extremely deep structures and massive data.
Moreover, repeated simple functions are doing a surprisingly good job while
approximating complicated ones; even though this is not fully understood.

Many DNNs have been proposed for more specific tasks, which all perform well.
However, despite the tremendous success, we still lack a theoretical
understanding of DNNs. In this study, we examine DNNs from a continuous
perspective, regard DNNs as discrete dynamical systems, and then improve the
stacking strategy of DNNs with high-order numerical methods.

Figure 1: HO-ResNet, stacking the same layers in higher order: Given fixed
layers, with the same parameters, there are more advanced stacking strategies
supported by numerical methods. We implemented three of them to be HO-ResNet,
which achieved noticeable improvements in several aspects. Moreover, the
degradation problem occurs later with more advanced numerical schemes.

It is common for DNNs to have more than a hundred layers; however, one usually
designs a sub-network with around 2-3 layers and then repeats it multiple
times rather than designing a huge one directly. The fixed topology and
relationship between blocks lead to the development of dynamic systems.

Enormous possibilities and advantages are gained from this view. For example,
mathematical problems are easier to solve and provide direct links to physical
sciences and solid theoretical support from differential equations. This study
focuses on the depth of DNNs, specifically, the stacking strategy of blocks:
constructing DNNs with a given sub-net design in high-order ways.

We show that a certain DNN (ResNet as the example in this study) can be easily
improved by slight adjustments in the block stacking strategy, following
numerical methods, without any changes in width and depth. Moreover, the
improvements are directly proportional to the order of the numerical methods,
not only in terms of accuracy but also in terms of convergence and robustness.

###  1.1 Related Work

The relationship between DNNs and dynamical systems has been widely discussed
in recent years. In 2016, Liao et al. [1] stated the equivalence of ResNet [2]
and a specific recurrent network based on formulations from dynamical systems.
A systematic proposal was published by Weinan [3] in 2017, indicating that
DNNs could be thought of as a discretisation of continuous dynamical systems;
the proposal and neural ordinary differential equations motivated us to
rethink ResNet using advanced numerical methods.

####  1.1.1 Neural ODEs

In 2018, neural ordinary differential equations (NODEs [4]) took a step
forward, replacing the discrete DNNs with solvers of ordinary differential
equations (ODEs).

To highlight the equivalence, we take a residual network (ResNet) as an
example. Let model ResNet be
$\mathbf{RES}^{L}(\mathbf{x}):\mathbb{R}^{d_{\text{in
}}}\rightarrow\mathbb{R}^{d_{\text{out }}}$ as an $L$-layer residual neural
network. We assume that the hidden state $\mathcal{N}$ belongs to
$\mathbb{R}^{d_{\text{N}}}$, and the input layer is
$\mathcal{N}^{0}(\mathbf{x})=\mathbf{x}\in\mathbb{R}^{d_{\mathrm{in}}}$. The
residual design is described as follows:

$\mathcal{F}\left(\mathbf{x},\left\\{\bm{W}^{i}\right\\}\right)+\mathbf{x},\quad
s.t.\quad 1\leq i\leq L-1.$ (1)

$i$ is an integer only and
$\mathcal{F}\left(\mathbf{x},\left\\{\bm{W}^{i}\right\\}\right)$ represents
the residual mapping to be learned which we do not expand further.

Figure 2: Discrete and continuous approaches of feature mapping: Discrete
methods usually design a basic structure (2-3 layers), and repeat it N times
to have various depths. The size of model, computational cost and performance
changes with the depth. Meanwhile, the NODEs map features end-to-end with
fully continuous flow and O(1) para size.

A sequence of transformation to a hidden state $\mathcal{N}$ from depth $i$ to
$i+1$ in residual networks can be written as follows, where $i\in\\{1\ldots
L-1\\}$ and $\bm{W}$ refers to the weight matrix:

$\mathcal{N}^{i+1}=\mathcal{F}\left(\mathcal{N}^{i},\left\\{\bm{W}^{i}\right\\}\right)+\mathcal{N}^{i},\quad
s.t.\quad 1\leq i\leq L-1.$ (2)

The situation above always takes a step $\Delta i=1$, letting $\Delta
i\rightarrow 0$, we have

$\lim_{\Delta i\rightarrow 0}\frac{\mathcal{N}_{i+\Delta
i}-\mathcal{N}_{i}}{\Delta
i}=\frac{\mathrm{d}\mathcal{N}(i)}{\mathrm{d}i}=\mathcal{F}(\mathcal{N}(i),i),$
(3)

thereby, hidden states can be parameterized using an ODE. Data point $x$ can
be mapped into a set of features by solving the initial value problem (IVP):

$\frac{\mathrm{d}\mathcal{N}(i)}{\mathrm{d}i}=\mathcal{F}(\mathcal{N}(i),i),\quad\mathcal{N}(0)=\mathbf{x},$
(4)

at a certain position. $\mathcal{N}(i)$ corresponds to the features learned by
the model. NODEs map the input to output by solving an ODE starting from
$\mathcal{N}(0)$ and adjust the dynamics to fit the output of the system
closer to the label.

Networks | Corresponding ODE formula | ODE Scheme
---|---|---
ResNet [2] | $\mathcal{N}^{t+1}=\mathcal{N}^{t}+\mathcal{F}\left(\mathcal{N}^{t},\bm{W}^{t}\right)$ | Forward Euler
RevNet [5] |  $\mathcal{Y}^{t}=\mathcal{N}^{t}+\mathcal{F}\left(\mathcal{N}^{t+1},\bm{W}^{t}\right)$, $\mathcal{Y}^{t+1}=\mathcal{N}^{t+1}+\mathcal{G}\left(\mathcal{Y}^{t},\bm{W}^{t}\right)$ | Forward Euler
PolyNet [6] | $\mathcal{N}^{t+1}=\mathcal{N}^{t}+\mathcal{F}\left(\mathcal{N}^{t},\bm{W}^{t}\right)+\mathcal{F}\left(\mathcal{F}\left(\mathcal{N}^{t},\bm{W}^{t}\right)\right)$ | Backward Euler
LM-ResNet [7] | $\mathcal{N}^{t+1}=\left(1-\Delta t\right)\mathcal{N}^{t}+\Delta t\mathcal{F}\left(\mathcal{N}^{t-1},\bm{W}^{t-1}\right)+\mathcal{F}\left(\mathcal{N}^{t},\bm{W}^{t}\right)$ | Linear-MultiStep
Second-order CNNs [8] | $\mathcal{N}^{t+1}=2\mathcal{N}^{t}-\mathcal{N}^{t-1}+{\Delta t}^{2}\mathcal{F}\left(\mathcal{N}^{t},\bm{W}^{t}\right)$ | Second-order
Table 1: Connections between neural networks and ODE schemes: Some networks
have been proposed without awareness of their connection with ODE schemes such
as ResNet, RevNet, and PolyNet, while others have been proposed following the
guidelines of ODE schemes such as LM-ResNet and Second-Order CNNs. The form of
ResNet is widely discussed as the Euler forward scheme, which is the most
straightforward way to solve the initial value problem. RevNet is a reversible
network, which means that the dynamic can be simulated from the end time back
to the initial time. PolyNet includes polynomial compositions that can be
interpreted as approximations to one step of the implicit backward Euler
scheme. LM ResNet adopted the known linear multistep method in numerical ODEs
and achieved the same accuracy by half the parameter size. Moreover, the
second-order CNN considers higher-order states, which is another direction
that can be further explored.

In addition to ResNet, several outstanding networks are linked to different
numerical schemes (see Table.1). The NODEs are a family of DNNs that can be
interpreted as a continuous equivalent of a certain discrete DNN (see
Figure.2). Inspired by several results of NODEs, we believe that continuous
theoretical concepts could lead to equivalent improvements in discrete DNNs.

####  1.1.2 Common high-order numerical schemes for IVP

With respect to the design of DNNs, structures by NAS have been in a leading
position in recent years. However, recent work from the Google team revisited
ResNet to conclude that state-of-the-art NAS structures such as Efficient-Nets
[9] are not necessarily better than ResNets [10]. Another excellent piece of
work, RepVGG [11], also achieved state-of-the-art results with a simple VGG-
style sub-net. The comeback of ResNet and simple sub-net motivated us to
rethink the block-stacking strategy, i.e., how to stack these subnets in a
better way. We studied ResNet because complicated NAS designs did not seem
necessarily better than ResNet.

Euler’s scheme is a simple and elegant one-order method, but most networks
contain a few basic components. The true expression of the stacked network is
nested and in a considerably high order. Will advanced numerical methods make
things better? Our stacking designs are based on the corresponding high-order
methods for the differential equation. Adopting Euler forward as the baseline
ResNet [12], mid-point [13], Runge-Kutta 4th, and Verner’s RK-8(9)th Order
Embedded Runge-Kutta method [14, 15, 16] are the general guidelines for our
proposed stacking design.

####  1.1.3 Effective DNNs that related to ODE schemes

Table.1 shows several outstanding networks with their corresponding schemes.
Besides them, many powerful networks are partially related to the ODE scheme.
For example, dense connections are like shortcuts of high-order RK schemes
[17], a DenseNet block is nested high-order RK; the One-Shot Aggregation (OSA)
proposed by VoVNet [18] is very similar to how RK methods conclude the block
output from each state inside. Moreover, in NAS, which discovered several
complicated structures, its stacking links (regardless of sub-net differences)
may eventually converge to a specific high-order scheme.

###  1.2 Contribution

##### Motivation

Our strategy is motivated by multiple recent works mentioned above. The
comeback of the simple sub-net design made us focus more on the stacking of
blocks rather than on the design of a small sub-net. The improved performance
of NODEs with advanced solvers indicated to us that it would be of interest to
explore stacked discrete networks, following numerical schemes with various
orders.

##### Pros and Cons

Our method enables DNNs to achieve better performance, more stable loss
landscape (before and after training), faster convergence, and stronger
robustness, which we will show in Section.3. Moreover, one can hardly notice
the additional cost when adopting 2-4 order methods, such as mid-point and
RK-4, which we will analyse at the end of Section.2. However, the model
requires considerable memory to store middle states if fully following a very
high-order design, such as Verner’s RK-8(9)th scheme.

The main contributions of this study are listed as follows:

  * 1.

We propose the Higher-Order ResNet (HO-Net) to explore, with a given sub-net
design, how high-order methods help in terms of performance, convergence and
robustness.

  * 2.

We provide sufficient fair comparisons, stacking basic Euler scheme to the
same order to enable comparison with high-order methods, therefore no extra
parameters and steps.

  * 3.

Visualisations of loss landscapes are given for better understanding.

  * 4.

Corresponding theoretical supports, complexity analysis and sufficient
ablation studies.

##  2 High-Order Residual Networks

NODEs compute the gradients of a scalar-valued loss with respect to all inputs
of any ODE solver, without back-propagating through the operations of the
solver. Similar to different sub-net designs in discrete networks, the methods
of the solver make a difference under the same conditions. Thus, advanced
solvers make NODEs better, but incur extra computational costs. Figure.3 shows
a simple case for better understanding.

Figure 3: Various schemes with a simple case: From end-to-end, higher-order
methods with more steps lead to a closer fitting of the target function but
also additional steps and costs.

Will advanced numerical schemes make ResNet better in a discrete scheme with
the same cost and parameter size? In this study, baseline ResNet (Euler
forward scheme), mid-point scheme, RK-4, fixed Verner-8(9), and adaptive
Verner’s 8(9) have a basic sub-net with [2, 4, 8, 28, 28] layers,
respectively. When comparing ResNet with HO-ResNet, we maintained the same
depth and width.

For example, a ResNet-Euler-18 has eight (8 × 2 layers) Euler-style blocks,
one embedding layer, and one FC layer, where ResNet-Midpoint-18 and ResNet-
RK4-18 have four (4 × 4 layers) Midpoint-style blocks and two (2 × 8 layers)
RK4-style blocks. The higher-order design evidently impacts the network, even
without training (see Figure.4, in relation to visualisations of initial
conditions.

Figure 4: Loss landscape without training: When processing the same images,
ResNet-18 with Euler, mid-point, and RK-4 schemes have ranges of loss ([0,
43.16], [8.20, 18.79], [1.99, 3.08]), respectively. A high-order scheme
significantly improves the stability of models for unseen samples; the loss
will not shack too much, which leads to the robustness of the learning rate
and other properties, which will be discussed in a later section.

A few definitions are necessary before introducing HO-ResNet. The situations
discussed below are without the embedding layer and FC layers, which means
that the input and output are all features. We use $\mathcal{F}$ to describe
any sub-net/component, which could be any component such as a VGG-style two-
layer sub-net, three-layer bottle net, ViT block, self-attention layer, or
simple feed-forward network.

However, we used a VGG-style block in the experiments because the question
’what form shall $\mathcal{F}$ be?’ was not included in our discussion. Our
focus is on propagating features in various schemes with a given
$\mathcal{F}$. We refer to the input/output features as $\mathcal{N}^{in}$ and
$\mathcal{N}^{out}$, denote the middle states of the input as
$\mathcal{N}^{mid}$ if there are two layers, $\mathbf{k}_{i}$ is used to
indicate the hidden states of the input after a particular $i$th
$\mathcal{F}$.

###  2.1 Euler Forward Scheme

Let a network be ResNet-6 which means that two ResBlocks are stacked, and we
can write the transforms between the embedding layer and FC layer as follows:

$\mathcal{N}^{Mid}=\mathcal{F}\left(\mathcal{N}^{in},\left\\{\bm{W}^{i}\right\\}\right),$
(5)

$\mathcal{N}^{Out}=\mathcal{F}\left(\mathcal{N}^{Mid}+\mathcal{N}^{in},\left\\{\bm{W}^{i}\right\\}\right)+\mathcal{N}^{Mid}.$
(6)

$\mathcal{N}^{Mid}$ denotes the output of the first $\mathcal{F}$; one may
find that the behaviour of two stacked ResBlocks is very similar to a mid-
point ResBlock, as shown in Figure.5.

###  2.2 Mid-point Scheme

The Euler method updates it with the easiest one-step method, which is $\Delta
i\mathcal{F}\left(i,\mathcal{N}^{i}\right)$, while the midpoint method updates
it in a higher-order manner:

$\mathcal{N}^{i+\Delta i}=\mathcal{N}^{i}+\Delta
i\mathcal{F}\left(i+\frac{\Delta i}{2},\mathcal{N}^{i}+\frac{\Delta
i}{2}\mathcal{N}\left(i,\mathcal{N}^{i}\right)\right).$ (7)

If we compare ResBlock-Euler and ResBlock-Midpoint, in terms of network
design, we will have different stacking, as shown in Figure.5. Rather than a
design with two layers, we now have a block design consisting of four layers.
Nevertheless, we stack the Euler scheme to the same layers/orders for a fair
comparison.

There were two noteworthy differences: (1). When adding the shortcut to obtain
the mid-state, the mid-point design compressed the output of $\mathcal{F}$ by
half. (2). the second shortcut is from an earlier location, which is directly
from the input. If we write the mid-point ResBlock design as:

$\mathcal{N}^{Mid}=\mathcal{F}\left(\mathcal{N}^{in},\left\\{\bm{W}^{i}\right\\}\right),$
(8)

$\mathcal{N}^{Out}=\mathcal{F}\left(\frac{1}{2}\mathcal{N}^{Mid}+\mathcal{N}^{in},\left\\{\bm{W}^{i}\right\\}\right)+\mathcal{N}^{in}.$
(9)

One finds that it is very similar to the stacked ResBlock-Euler. Therefore, we
assume that stacked ResNet will be easily improved with higher-order methods
because numerical problems have solid theoretical backups, which could be more
accurate than the lower-order ones.

Figure 5: Mid-point ResBlock: Comparing two stacked ResBlock-Euler with a
single design that following the Mid-point method.

###  2.3 Fourth order Runge-Kutta Scheme

Figure 6: RK4 ResBlock: Comparing four stacked ResBlock-Euler with a single
design which following the 4th-order Runge-Kutta method.

Exploring this further, with a 4th order design, Zhu et al. [19] attempted an
RK style design. In this study, we use it as a guideline for stacking rather
than a certain RK-Net. The Runge-Kutta method (RK4) updates it in a 4-step
manner, and our implemented ResBlock-RK4 can be written as follows:

$\displaystyle\mathbf{k}_{1}$
$\displaystyle=\mathcal{F}\left(\mathcal{N}^{in}\right),$ (10)
$\displaystyle\mathbf{k}_{2}$
$\displaystyle=\mathcal{F}\left(\frac{1}{2}\mathbf{k}_{1}+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{3}$
$\displaystyle=\mathcal{F}\left(\frac{1}{2}\mathbf{k}_{2}+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{4}$
$\displaystyle=\mathcal{F}\left(\mathbf{k}_{3}+\mathcal{N}^{in}\right),$

And the network output is:

$\mathcal{N}^{Out}=\mathcal{N}^{in}+\frac{1}{6}\left(\mathbf{k}_{1}+2\mathbf{k}_{2}+2\mathbf{k}_{3}+\mathbf{k}_{4}\right).$
(11)

Implementing a single ResBlock-RK4 requires eight layers, and thus we compare
four stacked ResBlock-Euler with it. The stacked ResBlock-Euler can be written
in the same way for easy comparison:

$\displaystyle\mathbf{k}_{1}$
$\displaystyle=\mathcal{F}\left(\mathcal{N}^{in}\right),$ (12)
$\displaystyle\mathbf{k}_{2}$
$\displaystyle=\mathcal{F}\left(\mathbf{k}_{1}+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{3}$
$\displaystyle=\mathcal{F}\left(\mathbf{k}_{2}+\mathbf{k}_{1}\right),$
$\displaystyle\mathbf{k}_{4}$
$\displaystyle=\mathcal{F}\left(\mathbf{k}_{3}+\mathbf{k}_{2}\right),$

And the network output is:

$\mathcal{N}^{Out}=\mathbf{k}_{4}+\mathbf{k}_{3}.$ (13)

According to the summary above, we can easily understand the difference
between 4th fourth-order RK-ResNet and the four stacked ResBlocks. In
Figure.6, we find that RK4-ResBlock is somehow similar to the DenseNet design
[17]. However, DenseNet will shortcut every $\mathbf{k}$ to layers after it,
which the proposed RK4-ResNet does not. Moreover, HO-ResNet scales the
particular output of $\mathbf{k}$ with mathematical support.

Midpoint-ResBlock and RK4-ResBlock only shortcut $\mathbf{k}$ to the block
output; identity mappings inside the block are all from the input x of the
block. However, some more high-order methods will start to map $\mathbf{k}$
inside the block to other states. However, the coefficient is precisely
calculated, rather than simply added.

###  2.4 Verner’s 8(9)th order Runge-Kutta Scheme

Of course, one could keep seeking a more advanced scheme to construct ResBlock
in a higher-order manner, that contains more layers than just two or three.
There are many other versions to guide the network design. Many of these
methods are adaptive methods that use a scaling factor. This means that the
compression rate of $\mathbf{k}_{i}$ can be adaptive depending on the error
returned.

In this study, we implemented the mid-point and Rk4 methods directly, as they
are fixed methods. However, we did not fully implement the adaptive Verner’s
RK-8 [16], as shown in the caption of Figure.15 for reasons.

Figure 7: Verner’s RK-8(9)th order scheme: The full method has 16 steps and
extra processes for error estimation and scale factor adjustment. In our
implementation, Verner’s 8(9) th-order block has 28 layers instead of (8 × 2 ×
2) layers for two reasons: (1). $\mathbf{k}_{15}$ and $\mathbf{k}_{16}$ do not
impact the block output $\mathcal{N}^{Out}$ and are used to calculate the
adaptive factor $h$. (2). Moreover, one can update $h$ via deep learning
frameworks directly by setting it as a learnable tensor value. We compared
both fixed $h$ and adaptive $h$, where learnable $h$ is the only extra
parameter in HO-ResNet, usually no more than 50, making a rare difference
because most models have millions of parameters. $h$ equals one in a fixed
Verner block; for the adaptive version, one could use a shard $h$ for all
hidden states from $\mathbf{k}_{1}$ to $\mathbf{k}_{14}$ or every state could
have a different $h$.

$\displaystyle\mathbf{k}_{1}$
$\displaystyle=\mathcal{F}\left(\mathcal{N}^{in}\right),$ (14)
$\displaystyle\mathbf{k}_{2}$
$\displaystyle=\mathcal{F}\left(\frac{h}{12}\mathbf{k}_{1}+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{3}$
$\displaystyle=\mathcal{F}\left(\frac{h}{27}\left(\mathbf{k}_{1}+2\mathbf{k}_{2}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{4}$
$\displaystyle=\mathcal{F}\left(\frac{h}{24}\left(\mathbf{k}_{1}+3\mathbf{k}_{3}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{5}$
$\displaystyle=\mathcal{F}\left(\frac{h}{375}\left(\left(4+94\sqrt{6}\right)\mathbf{k}_{1}-\left(282+252\sqrt{6}\right)\mathbf{k}_{3}+\left(328+206\sqrt{6}\right)\mathbf{k}_{4}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{6}$
$\displaystyle=\mathcal{F}\left(h\left(\frac{9-\sqrt{6}}{150}\mathbf{k}_{1}+\frac{312+32\sqrt{6}}{1425}\mathbf{k}_{4}+\frac{69+29\sqrt{6}}{570}\mathbf{k}_{5}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{7}$
$\displaystyle=\mathcal{F}\left(h\left(\frac{927-347\sqrt{6}}{1250}\mathbf{k}_{1}+\frac{7328\sqrt{6}-16248}{9375}\mathbf{k}_{4}+\frac{179\sqrt{6}-489}{3750}\mathbf{k}_{5}+\frac{14268-5798\sqrt{6}}{9375}\mathbf{k}_{6}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{8}$
$\displaystyle=\mathcal{F}\left(\frac{h}{54}\left(4\mathbf{k}_{1}+\left(16-\sqrt{6}\right)\mathbf{k}_{6}+\left(16+\sqrt{6}\right)\mathbf{k}_{7}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{9}$
$\displaystyle=\mathcal{F}\left(\frac{h}{512}\left(38\mathbf{k}_{1}+\left(118-23\sqrt{6}\right)\mathbf{k}_{6}+\left(118+23\sqrt{6}\right)\mathbf{k}_{7}-18\mathbf{k}_{8}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{10}$
$\displaystyle=\mathcal{F}\left(h\left(\frac{11}{144}\mathbf{k}_{1}+\frac{266-\sqrt{6}}{864}\mathbf{k}_{6}+\frac{266+\sqrt{6}}{864}\mathbf{k}_{7}-\frac{1}{16}\mathbf{k}_{8}-\frac{8}{27}\mathbf{k}_{9}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{11}$
$\displaystyle=\mathcal{F}\left(h\left(\frac{5034-271\sqrt{6}}{61440}\mathbf{k}_{1}+\frac{7859-1626\sqrt{6}}{10240}\mathbf{k}_{7}+\frac{813\sqrt{6}-2232}{20480}\mathbf{k}_{8}+\frac{271\sqrt{6}-594}{960}\mathbf{k}_{9}+\frac{657-813\sqrt{6}}{5120}\mathbf{k}_{10}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{12}$
$\displaystyle=\mathcal{F}\left(h\left(-8.14164\mathbf{k}_{1}-574.436\mathbf{k}_{6}+847.88\mathbf{k}_{7}+113.719\mathbf{k}_{8}+626.94\mathbf{k}_{9}+605.73\mathbf{k}_{10}-328.69\mathbf{k}_{11}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{13}$
$\displaystyle=\mathcal{F}\left(h\left(0.0878\mathbf{k}_{1}+0.69337\mathbf{k}_{6}-1.9\mathbf{k}_{7}+0.23\mathbf{k}_{8}-0.69\mathbf{k}_{9}-0.077\mathbf{k}_{10}+2.49\mathbf{k}_{11}+0.0018\mathbf{k}_{12}\right)+\mathcal{N}^{in}\right),$
$\displaystyle\mathbf{k}_{14}$
$\displaystyle=\mathcal{F}\left(h\left(-0.1\mathbf{k}_{1}+5.575\mathbf{k}_{6}+7.486\mathbf{k}_{7}-6.23\mathbf{k}_{8}+2.27\mathbf{k}_{9}-4.89\mathbf{k}_{10}-4.86\mathbf{k}_{11}-0.0235\mathbf{k}_{12}+1.78\mathbf{k}_{13}\right)+\mathcal{N}^{in}\right),$

And the network output is:

$\mathcal{N}^{Out}=\mathcal{N}^{in}+h\left(0.06\mathbf{k}_{1}-0.19\mathbf{k}_{8}+0.72\mathbf{k}_{9}-0.72\mathbf{k}_{10}+0.75\mathbf{k}_{11}+0.0004\mathbf{k}_{12}+0.34\mathbf{k}_{13}+0.032\mathbf{k}_{14}\right)$
(15)

The full Verner’s RK-8(9) ResBlock design shall contain 16 $\mathcal{F}$ which
means 32 layers. We use $\mathbf{k}_{i}$ to represent the output of each
block, where $i$ is an integer in the range [1-16]. Unfortunately, no
flowchart is provided because of the complicated connections between the
hidden states.

The final output of the Verner-8(9) ResBlock depends only on the shortcut from
the start and $\mathbf{k}_{i}$, where $i$ in [1,8,9,10,11,12,13,14]. The scale
factor $h$ should be adjusted depending on $\mathbf{k}_{i}$, where $i$ in
[1,8,9,10,11,12,13,14,15,16]. Thereby, one can implement a fixed version with
an out-of-scale factor adjustment using 14 $\mathcal{F}$. However,
$\mathbf{k}_{15}$ and $\mathbf{k}_{16}$ are required to update $h$ if one
wants to fully follow Verner’s method, in which we update $h$ by directly
setting it to be learnable, in order to compare fixed and adaptive versions
with the same computational cost.

.

###  2.5 Complexity Analysis

We compare the situation in which the same size, number, and design of the
layers are used. The space complexity can be described by the

$\text{ Space }\sim O\left(\sum_{i=1}^{D}K_{i}^{2}\cdot C_{l-1}\cdot
C_{i}+\sum_{i=1}^{D}M^{2}\cdot C_{i}\right)$ (16)

where K is the kernel size, C refers to the number of channels, and D denotes
the depth of the network. Thus, the complexity of various ResNet designs is
the same, with no extra parameters, theoretically. However, in practice, extra
space is required to maintain the output inside the high-order block because
the block output depends on more states at particular layers rather than two
fixed positions, compared with the baseline ResNet.

Specifically, baseline ResNet only maintains a shortcut from the block input,
as does the mid-point block. For the Runge-Kutta 4 block, we need three more
shortcuts for the block output. However, this could be optimised to one more
shortcut.

For Verner’s RK-8 block, this situation becomes rather complicated. From the
state $\mathbf{k}_{5}$, every state starts to depend on other previous states,
besides the input from the last state; we shall maintain 16 more shortcuts for
the full version or 14 more shortcuts for our version. Furthermore, the space
requirement is extremely high while forwarding, which is a crucial shortage;
sufficient space should be provided when using such a design. Therefore, it is
not recommended to use Verner’s block; applying the adaptive step from
Verner’s block to a fourth-order scheme would be a better solution.

Our design does not affect the first term, $\sum_{i=1}^{D}K_{i}^{2}\cdot
C_{l-1}\cdot C_{i}$, but only the second term, $\sum_{i=1}^{D}M^{2}\cdot
C_{i}$. Baseline ResNet maintains only one feature map for shortcut as
$\sum_{i=1}^{D}M^{2}\cdot\left(C_{i}+1\right)$, similar to the midpoint block.
RK4-Block and Fixed-RK8-Block will take
$\sum_{i=1}^{D}M^{2}\cdot\left(C_{i}+2\right)$ and
$\sum_{i=1}^{D}M^{2}\cdot\left(C_{i}+15\right)$, respectively. We could say
that the extra space requirement could be ignored generally for second-fourth-
order schemes but is slightly too large for full Verner’s RK-8(9).

We may describe the time complexity by:

$\text{ Time }\sim O\left(\sum_{i=1}^{D}M_{i}^{2}\cdot K_{i}^{2}\cdot
C_{i-1}\cdot C_{i}\right)$ (17)

According to the explanation above, we know that HO-ResNet will frequently
scale the shortcuts, unavoidably creating extra multi/add processes. Compared
with the baseline ResNet, the midpoint block has one additional multiplication
operation for every four layers. Therefore, we shall add
$\frac{1}{4}\sum_{i=1}^{D}M^{2}\cdot C_{i}$.

The RK-4 block requires six more multiplication operations and three
additional operations every eight layers. Therefore, we shall add
$\frac{9}{8}\sum_{i=1}^{D}M^{2}\cdot C_{i}$. It is difficult to observe the
complexity changes among Euler, Midpoint, and RK4. However, Verner’s RK-8
requires a considerable number of storage spaces because of the highly
correlated hidden states and extensive additional multi/add scaling. The time
consumption is approximately 1.15x-1.2x while the previous three remain the
same.

Schemes | Euler | Midpoint | RK-4 | Verner-8(9)
---|---|---|---|---
Memory | 1x | 1x | 1.5x(1x) | 7x
Time | 1x | 1x | 1x | 1.15x-1.2x
Table 2: Multiple relationships of memory requirements and running time: The
second to fourth-order scheme has rare extra operations, and thus does not
reflect on the running time. However, more memory is required because RK-4
blends all hidden states for the output. Verner’s 8(9) requires huge amounts
of memory because it blends not only hidden states for output but also
previous states for hidden states in a later stage.

Table.2 shows the multiple relationships of memory and running time in
practice, taking the ResNet as the baseline. There is no room to optimize the
huge cost of Verner’s 8(9), but one could remove the blending of RK-4 and
output the $\mathcal{N}^{in}+\mathbf{k}_{4}$ only instead of
$\mathcal{N}^{in}+\frac{1}{6}\left(\mathbf{k}_{1}+2\mathbf{k}_{2}+2\mathbf{k}_{3}+\mathbf{k}_{4}\right)$.
RK-4 will have the same memory requirement as baseline Euler by doing so with
no evident performance drop.

Figure 8: Visualisation of untrained ResNet-58 with various schemes : The
range of loss exploded with a low-order scheme. Euler, Midpoint, RK-4, and
Verner’s 8(9) have loss ranges ([1e5, 1e7], [1028, 7097], [4.76, 13.04],
[1.47, 1.90]). The stability is significantly improved with higher order,
which means that gradient problems, such as vanishing or exploding, are less
likely to occur. Figure 9: ResNet-86 with various schemes against different
initial learning rates : Under the same condition, ResNet with higher-order
could be trained with a noticeable larger learning rate. Baseline Euler scheme
diverged with learning rate 0.13 and all lr larger than it, while Mid-point
could converge with lr close to 0.2, RK-4 could converge at lr 0.23 and
Verner’s-8(9) can converge at lr 0.4 (or even above). Figure 10: One midpoint
block vs. two Euler blocks: The second-order scheme shows noticeable
improvements in both accuracy and loss during training and testing. Figure 11:
One Rk4 block vs. two Midpoint blocks and four Euler Blocks: The fourth-order
RK-4 scheme outperforms the mid-point and baseline Euler as a basic sub-net.
Note that the advantages will become increasingly evident for deeper
situations. Figure 12: ResNet-58 with various schemes: Training accuracy and
loss on top and testing accuracy and loss at the bottom. Note that the
difference between schemes is becoming increasingly noticeable compared with
observations on shallow versions. Figure 13: Running Time comparison
approaching 0.8, and 0.85 train accuracy (ResNet-58): High-order schemes have
evident advantages in terms of the speed of convergence. The second-order mid-
point takes half of the time to reach 0.8, and 0.85, compared with Euler,
while the fourth and 8(9) order methods take 20–30 per cent of the time. See
details in Table.4.

##  3 Experiments

###  3.1 Datasets and implementation

In this study, we conducted most studies on the CIFAR–10 datasets [20], which
all consist of coloured natural images with 32 × 32 pixels with 50,000 images
for training and 10,000 images for testing, which has 10 target classes in
total for CIFAR-10 with 6,000 images per class.

We directly stacked blocks 64 dim after a conv(3, 64). Stochastic gradient
descent [21] under a weight decay of 1e-4 and momentum of 0.9, was used with
weight initialisation in [22] and batch norm [23]. The models were trained
with a batch size of 128, starting with an initial learning rate of 0.1, step
scheduler with [100, 150, 200, 230], factor 0.1, and 260 epochs trained.

###  3.2 Augmentations and tricks

Our focus was on the differences in behaviour for each of the stacking
strategies to verify whether the high-order stacking could provide general
guidance for neural network design, but not on pushing a certain line of
results or seeking a more considerable number. Therefore, we intentionally
compared those stacking strategies using methods that enabled fair comparisons
to be drawn.

When comparing a single block design with the corresponding deep baseline
ResNet, we used raw data and observed overfitting even under 10 layers.
Therefore, we used the auto augments’ CIFAR10 policy [24] to address this
question.

###  3.3 Robustness against learning rate

Figure.8 shows a deeper situation with different schemes. Comparing with the
case of ResNet-18 shown by Figure.4, the loss exploded due to the chain rule.
Meanwhile, higher order schemes show significant advantages, such as more
flatten loss surfaces being provided. Such a phenomenon greatly enhances the
robustness of models when different learning rates are given.

The learning rate is one of the primary hyperparameters of deep-learning
systems. Slight adjustments in the learning rate could cause vast impacts on
the whole system, large learning rates will easily lead to divergence, and
small learning rates usually fail to converge well.

Taking ResNet-86 as an example, 84 layers could consist of [42, 21, 7, 3] of
[Euler, Midpoint, RK-4, Verner’s 8(9)] blocks. Figure.9 shows how higher-order
schemes are dealing with a larger range of learning rates. ResNet-86 in the
4-th/8-th order scheme could be well trained with 2x/4x learning rate at the
point where the baseline Euler diverged

Models | Layers | Params | cifar-10
---|---|---|---
Euler X4 | 10 | 0.3M | 91.72
Midpoint X2 | 10 | 0.3M | 91.97
RK4 X1 | 10 | 0.3M | 91.51
Euler X8 | 18 | 0.59M | 92.68
Midpoint X4 | 18 | 0.59M | 92.95
RK4 X2 | 18 | 0.59M | 92.98
Euler X14 | 30 | 1.03M | 93.02
Midpoint X7 | 30 | 1.03M | 93.17
RK8 X1 Fix | 30 | 1.03M | 93.39
RK4 X3 | 26 | 0.89M | 93.41
RK8 X1 Adaptive | 30 | 1.03M | 94.17
Euler X28 | 58 | 2.07M | 92.81
Midpoint X14 | 58 | 2.07M | 94.08
RK4 X7 | 58 | 2.07M | 94.79
RK8 X2 Fix | 58 | 2.07M | 94.91
RK8 X2 Adaptive | 58 | 2.07M+2 | 95.06
Euler X42 | 86 | 3.01M | 92.73
Midpoint X21 | 86 | 3.01M | 93.87
RK4 X10 | 82 | 2.87M | 94.83
RK8 X3 Fix | 86 | 3.01M | 95.02
RK8 X3 Adaptive | 86 | 3.01M+3 | 95.15
Table 3: Classification accuracy on CIFAR-10: Various schemes with different
depth are compared. All schemes perform similar when depth smaller than 18.
After that, higher order schemes are showing more and more apparent advantages
in terms of classification accuracy.

###  3.4 Fair comparison of single block

If one compares each block design with equivalent deep baseline ResBlocks, it
is found that the improvements are stable but not apparent in the shallow
situation. However, the advantages were rapidly becoming evident in more
profound cases.

##### Euler-x2 and Midpoint-x1

Given ResNet-6, we have four layers beside the input and output layers. We
performed experiments on CIFAR-10 raw data with rare tricks; as shown in
Figure.10, loss and accuracy are stably improved.

##### Euler-x4, Midpoint-x2 and RK4-x1

If ResNet-10 was given, one could stack four baseline ResBlocks, which is the
red line in Figure.11, while Midpoint-2 and single RK-4 block benefit from
more accurate updating methods.

###  3.5 Deeper Situation with augmentation

We noticed overfitting when using ResNet shallower than 10. Therefore, we
added auto-augmentation [24] in this comparison. Without the input and output
layers, 56 layers could be divided into Baseline ResBlock X28, Midpoint x14,
RK4-Block x7, and fixed-RK8 X2.

As shown in Figure.12, baseline ResNet, represented by the red line, is
clearly improved by high-order stacking methods, using the same setting and
parameters. We noticed that the fixed-RK8 method has the potential to improve.
In the early stages it was above RK-4. In Verner’s RK-8(9) [16] design, one
should also maintain a scale factor $h$ and $error$ to adjust the step side.
Adaptive-step Verner’s RK-8(9) converges at a slightly slower rate but
outperforms all other schemes overall.

Table.3 shows more details against different depths; all schemes perform
similarly in shallow cases but show more evident gaps while getting deeper.
Moreover, there are numerous benefits to be gained from high-order, besides
the performance, such as loss landscape and robustness against
hyperparameters, Figure.13 and Table.4 show the improvement in the speed of
converging.

Schemes | Euler | Midpoint | RK-4 | Fix RK-8 | Ada RK-8
---|---|---|---|---|---
To 0.8 | 1x | 2.27x | 4.43x | 3.31x | 3.84x
To 0.85 | 1x | 2.28x | 3.92x | 3.53x | 3.03x
Accuracy | 92.71 | 94.08 | 94.79 | 94.72 | 95.06
Extra para | / | / | / | / | 2

Table 4: Ratio of time cost (ResNet-58): Note that Verner’s RK-8 costs a 1.15
increase in time compared to the others; this situation leads RK-4 to be
deemed to achieve the best performance in terms of time. Adaptive Verner’s
RK-8 appears to struggle in the mid-stage but does converge to achieve the
best overall performance in the end, with only two extra parameters needed.

###  3.6 Scaling Factors

The most effective DNNs have a relatively light stacking design between blocks
and rarely feature scaling factors on block outputs. Following the Midpoint
and RK4 scheme, the output should be compressed with a scaling factor of 0.5,
and RK4 will have 1/6 and 1/3 factors when averaging the output inside the
block.

In RK-8 [16], a constraint [0.125, 4] was given while scaling the output from
a particular layer. Although we did not further explore the adaptive stacking
design, we tried a simple test on how various scaling factors impact the same
DNNs.

Taking ResNet-58 as an example, which has two fixed Verner’s RK-8 blocks and
two extra parameters, see Table.4. Although the full Verner’s RK-8 constrains
the factor inside [0.125, 4], we did not follow the constraint when applying
the adaptive step for practical reasons. However, the results still show a
noticeable improvement, as expected.

###  3.7 Discussing

##### Stacking strategy is vital

In the past few years, some network designs have been regarded as ODE schemes
such as ResNet and PolyNet [2, 25]; some networks are designed under the
guidance of ODE schemes such as LM-ResNet and RK-Net [7, 19]. However, most of
them focus more on the subnet design, while the stacking strategy is also
critical. Block-to-block is similar to layer-to-layer, especially when some
networks will even stack 84 blocks in a single stage.

##### Three principles from High order schemes

(1). One shall apply an adaptive time step; it can easily be plugged in to any
network. In fact, several excellent studies have already proposed methods in
this manner from other perspectives. Highway networks [26], Fixup [27],
SkipInit [28], ReZero [29], and LayerScale [30] are the approaches that are
most closely related. (2). One should consider connections over long
distances, but not nested ones such as DenseNet [17] to avoid huge memory
requirements. (3). Blending from different hidden states works, but it is not
as important as the previous two points. Moreover, it requires at least 1x
extra memory.

##### Unavoidable degradation

ResNet [2] significantly alleviated the degradation problem that deeper
networks shall perform better in theoretical but worse in practice. The
shortcut does not provide extra capacity but makes models easier for training,
so does HO-ResNet. We implemented ResNet with advanced numerical methods based
on the equivalence between ResNet and Euler Forward scheme; the results in
Figure.1 show that the degradation was better handled and further alleviated
with more advanced numerical schemes; however, the degradation will still come
eventually.

##  4 Conclusion

Inspired by several existing studies [3, 7, 4, 10, 2], we assume that DNNs are
extremely complicated numerical solvers for differential equations. Most
designs are relatively low-order, while the networks are frequently deeper
than 100 or 200. We propose a general higher-order stacking strategy,
following several numerical methods, which proved that high-order numerical
methods could be adopted as general guidance for block stacking and network
design. Sufficient experiments show that high-order improvements are stable,
robust, and fully explainable in terms of math. We show that with very few
changes in stacking, ResNet has shown remarkable durability.

## References

  * [1] Q. Liao, T. A. Poggio, Bridging the gaps between residual learning, recurrent neural networks and visual cortex, CoRR abs/1604.03640 (2016).  arXiv:1604.03640.
URL http://arxiv.org/abs/1604.03640

  * [2] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, in: 2016 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2016, Las Vegas, NV, USA, June 27-30, 2016, IEEE Computer Society, 2016, pp. 770–778.  doi:10.1109/CVPR.2016.90.
URL https://doi.org/10.1109/CVPR.2016.90

  * [3] E. Weinan, A proposal on machine learning via dynamical systems, Communications in Mathematics and Statistics 5 (1) (2017) 1–11.
  * [4] R. T. Chen, Y. Rubanova, J. Bettencourt, D. K. Duvenaud, Neural ordinary differential equations, in: Advances in neural information processing systems, 2018, pp. 6571–6583.
  * [5] A. N. Gomez, M. Ren, R. Urtasun, R. B. Grosse, The reversible residual network: Backpropagation without storing activations, CoRR abs/1707.04585 (2017).  arXiv:1707.04585.
URL http://arxiv.org/abs/1707.04585

  * [6] X. Zhang, Z. Li, C. C. Loy, D. Lin, Polynet: A pursuit of structural diversity in very deep networks, in: 2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017, IEEE Computer Society, 2017, pp. 3900–3908.  doi:10.1109/CVPR.2017.415.
URL https://doi.org/10.1109/CVPR.2017.415

  * [7] Y. Lu, A. Zhong, Q. Li, B. Dong, Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations, in: J. G. Dy, A. Krause (Eds.), Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018, Vol. 80 of Proceedings of Machine Learning Research, PMLR, 2018, pp. 3282–3291.
URL http://proceedings.mlr.press/v80/lu18d.html

  * [8] L. Ruthotto, E. Haber, Deep neural networks motivated by partial differential equations, J. Math. Imaging Vis. 62 (3) (2020) 352–364.  doi:10.1007/s10851-019-00903-1.
URL https://doi.org/10.1007/s10851-019-00903-1

  * [9] M. Tan, Q. V. Le, Efficientnet: Rethinking model scaling for convolutional neural networks, in: K. Chaudhuri, R. Salakhutdinov (Eds.), Proceedings of the 36th International Conference on Machine Learning, ICML 2019, 9-15 June 2019, Long Beach, California, USA, Vol. 97 of Proceedings of Machine Learning Research, PMLR, 2019, pp. 6105–6114.
URL http://proceedings.mlr.press/v97/tan19a.html

  * [10] I. Bello, W. Fedus, X. Du, E. D. Cubuk, A. Srinivas, T.-Y. Lin, J. Shlens, B. Zoph, Revisiting resnets: Improved training and scaling strategies, arXiv preprint arXiv:2103.07579 (2021).
  * [11] X. Ding, X. Zhang, N. Ma, J. Han, G. Ding, J. Sun, Repvgg: Making vgg-style convnets great again, CoRR abs/2101.03697 (2021).  arXiv:2101.03697.
URL https://arxiv.org/abs/2101.03697

  * [12] J. Jacod, P. Protter, et al., Asymptotic error distributions for the euler method for stochastic differential equations, Annals of Probability 26 (1) (1998) 267–307.
  * [13] K. Burrage, I. Lenane, G. Lythe, Numerical methods for second-order stochastic differential equations, SIAM J. Sci. Comput. 29 (1) (2007) 245–264.  doi:10.1137/050646032.
URL https://doi.org/10.1137/050646032

  * [14] J. C. Butcher, On the implementation of implicit runge-kutta methods, BIT Numerical Mathematics 16 (3) (1976) 237–240.
  * [15] L. F. Shampine, Some practical runge-kutta formulas, Mathematics of computation 46 (173) (1986) 135–150.
  * [16] J. H. Verner, Explicit runge–kutta methods with estimates of the local truncation error, SIAM Journal on Numerical Analysis 15 (4) (1978) 772–790.
  * [17] G. Huang, Z. Liu, L. van der Maaten, K. Q. Weinberger, Densely connected convolutional networks, in: 2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017, IEEE Computer Society, 2017, pp. 2261–2269.  doi:10.1109/CVPR.2017.243.
URL https://doi.org/10.1109/CVPR.2017.243

  * [18] Y. Lee, J. Hwang, S. Lee, Y. Bae, J. Park, An energy and gpu-computation efficient backbone network for real-time object detection, in: IEEE Conference on Computer Vision and Pattern Recognition Workshops, CVPR Workshops 2019, Long Beach, CA, USA, June 16-20, 2019, Computer Vision Foundation / IEEE, 2019, pp. 752–760.  doi:10.1109/CVPRW.2019.00103.
URL
http://openaccess.thecvf.com/content_CVPRW_2019/html/CEFRL/Lee_An_Energy_and_GPU-
Computation_Efficient_Backbone_Network_for_Real-
Time_Object_CVPRW_2019_paper.html

  * [19] M. Zhu, C. Fu, Convolutional neural networks combined with runge-kutta methods, CoRR abs/1802.08831 (2018).  arXiv:1802.08831.
URL http://arxiv.org/abs/1802.08831

  * [20] A. Krizhevsky, G. Hinton, et al., Learning multiple layers of features from tiny images (2009).
  * [21] L. Bottou, Large-scale machine learning with stochastic gradient descent, in: Y. Lechevallier, G. Saporta (Eds.), 19th International Conference on Computational Statistics, COMPSTAT 2010, Paris, France, August 22-27, 2010 - Keynote, Invited and Contributed Papers, Physica-Verlag, 2010, pp. 177–186.  doi:10.1007/978-3-7908-2604-3\\_16.
URL https://doi.org/10.1007/978-3-7908-2604-3_16

  * [22] K. He, X. Zhang, S. Ren, J. Sun, Delving deep into rectifiers: Surpassing human-level performance on imagenet classification, in: 2015 IEEE International Conference on Computer Vision, ICCV 2015, Santiago, Chile, December 7-13, 2015, IEEE Computer Society, 2015, pp. 1026–1034.  doi:10.1109/ICCV.2015.123.
URL https://doi.org/10.1109/ICCV.2015.123

  * [23] S. Ioffe, C. Szegedy, Batch normalization: Accelerating deep network training by reducing internal covariate shift, in: F. R. Bach, D. M. Blei (Eds.), Proceedings of the 32nd International Conference on Machine Learning, ICML 2015, Lille, France, 6-11 July 2015, Vol. 37 of JMLR Workshop and Conference Proceedings, JMLR.org, 2015, pp. 448–456.
URL http://proceedings.mlr.press/v37/ioffe15.html

  * [24] E. D. Cubuk, B. Zoph, D. Mané, V. Vasudevan, Q. V. Le, Autoaugment: Learning augmentation policies from data, CoRR abs/1805.09501 (2018).  arXiv:1805.09501.
URL http://arxiv.org/abs/1805.09501

  * [25] W. Liu, A. Rabinovich, A. C. Berg, Parsenet: Looking wider to see better, CoRR abs/1506.04579 (2015).  arXiv:1506.04579.
URL http://arxiv.org/abs/1506.04579

  * [26] R. K. Srivastava, K. Greff, J. Schmidhuber, Highway networks, CoRR abs/1505.00387 (2015).  arXiv:1505.00387.
URL http://arxiv.org/abs/1505.00387

  * [27] H. Zhang, Y. N. Dauphin, T. Ma, Fixup initialization: Residual learning without normalization, in: 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019, OpenReview.net, 2019.
URL https://openreview.net/forum?id=H1gsz30cKX

  * [28] S. De, S. L. Smith, Batch normalization biases deep residual networks towards shallow paths, CoRR abs/2002.10444 (2020).  arXiv:2002.10444.
URL https://arxiv.org/abs/2002.10444

  * [29] T. Bachlechner, B. P. Majumder, H. H. Mao, G. W. Cottrell, J. J. McAuley, Rezero is all you need: Fast convergence at large depth, CoRR abs/2003.04887 (2020).  arXiv:2003.04887.
URL https://arxiv.org/abs/2003.04887

  * [30] H. Touvron, M. Cord, A. Sablayrolles, G. Synnaeve, H. Jégou, Going deeper with image transformers, arXiv preprint arXiv:2103.17239 (2021).


#  A novel directly energy-preserving method for charged particle dynamics

Yexin Lia,b, Ping Jiangc, Haochen Lia,b,∗ a School of Science, Beijing
University of Posts and Telecommunications, Beijing 100876, China
b Key Laboratory of Mathematics and Information Networks(Beijing University of
Posts and Telecommunications), Ministry of Education, China
c Peking University Chongqing Research Institute of Big Data, Chongqing
400031, China

###### Abstract

In this paper, we apply the coordinate increment discrete gradient (CIDG)
method to solve the Lorentz force system which can be written as a non-
canonical Hamiltonian system. Then we can obtain a new energy-preserving
CIDG-I method for the system. The CIDG-I method can combine with its adjoint
method CIDG-II which is also a energy-preserving method to form a new method,
namely CIDG-C method. The CIDG-C method is symmetrical and can conserve the
Hamiltonian energy directly and exactly. With comparison to the well-used
Boris method, numerical experiments indicate that the CIDG-C method holds
advantage over the Boris method in terms of energy-conserving.

The Fig. 5(b) in the original paper [31] contains an error. We submit the
correct Fig. 5(b) and an errata in this paper, which is described in remark
4.1.

###### keywords:

Lorentz force system; Hamiltonian system; Energy-preserving; Discrete gradient
methods.

∗Corresponding author.
E-mail address<EMAIL_ADDRESS>
##  1 Introduction

With the continuous breakthrough the theory of scientific computing and the
explosive growth of its application across diverse scientific discipline, the
numerical solution of the Hamilton system has gradually become the focus of
attention of mathematicians and physicists. For a system, a numerical method
which can conserve one or more geometric properties is called geometric
numerical integration method, or structure-preserving method [1, 2, 3, 4],
such as symplectic methods, volume-preserving methods, energy-preserving
methods, and so forth. Geometric numerical integration method has become an
essential tool to solve the Hamiltonian systems. And the method has been
effectively applied in many equations, including but not limited to Hénon-
Heiles system, nonlinear Schrödinger equation, Maxwell’s equation. [5, 6, 7,
8, 9, 10].

Because the collisions between the charged particles occur infrequently in hot
plasmas, many important phenomenons in hot plasmas can be simplified and
considered as a charged particle motion. Therefore, the motion of a charged
particle under the influence of electromagnetic fields is one of the most
fundamental process in hot plasmas [11, 12]. Long-term numerical calculation
on the trajectories of a charged particle have been widely used to study its
dynamical behaviors. However, the non-geometric method such as the standard
4th order Runge-Kutta method is not accurate in the long term simulation. For
a charged particle in the electromagnetic field, its dynamics is governed by
Lorentz force system, which can be written as a non-canonical Hamiltonian
system [13, 14]. Therefore, the solution of the Lorentz force system can be
derived by applying the geometric numerical integration method on its non-
canonical Hamiltonian system.

A considerable amount of literature on geometric numerical integration method
has been published, including the volume-preserving algorithm [15, 16], the
Boris method which can conserve phase space volume [12, 17, 18, 19], the
variational symplectic method [20, 21, 22]. The Boris method, which is the
$de$ $facto$ standard for advancing a charged particle in an electromagnetic
field and conserves phase space volume exactly. Albeit the Boris algorithm has
proven itself highly effective for charged particle dynamics in an
electromagnetic field, but it does not conserve energy and a numerical drift
may be observed inevitably [23]. It is well-known that the most essential
feature of a Hamiltonian system is the Hamiltonian energy. So We hope to
construct an energy-preserving or approximate energy-preserving algorithm to
solve the Lorentz force system numerically, and obtain long-term correct
results without numerical drift.

Based on numerical quadrature formula, the Bool discrete line integral method
[24] and the high-order energy-conserving line integral methods [25] can
conserve the energy exactly for the polynomial Hamiltonian, and conserve the
energy of the non-polynomial Hamiltonian up to the round-off error. The
invariant energy quadratization method [26] can conserve the modified energy,
which is the approximation of the original Hamiltonian energy. However, there
are few numerical methods that can directly and exactly conserve the original
Hamiltonian energy. If we want to construct an accurate energy-preserving
algorithm for the Lorentz force system, we need to focus on the algorithms
that preserve the energy of the Hamiltonian system.

Over the past few decades, a number of methods have been proposed to preserve
the energy of the Hamiltonian system. In 1996, one of the most famous energy-
preserving methods which can construct integral preserving methods for
ordinary differential equations was developed by O. Gonzalez, and the method
is known as the discrete gradient method [27]. Subsequently, the discrete
variational method for nonlinear wave equation was proposed by Matsuo [28].
E.Faou proposed that there is an energy-preserving B-series method in the
Hamiltonian system. Moreover, Quispel and McLaren derived the averaged vector
field method [29]. For the Lorentz force system in the form of the Hamiltonian
system, this study presents a novel directly energy-preserving method.

Fistly, the Lorentz force system can be write as a non-canonical Hamiltonian
form. Then we apply the coordinate increment discrete gradient method [30] to
the Hamiltonian system to get a new energy-preserving method. The new method
and its adjoint method are both energy-preserving methods. From this two
methods we can derive a new coordinate increment discrete gradient composition
(CIDG-C) method. A key advantage of the CIDG-C method is that it is
symmetrical, while simultaneously conserving the Hamiltonian energy exactly
without numerical quadrature formula.

This paper is organized as follows. In Section. 2, We introduce the dynamics
of charged particles in the electromagnetic field, which is written as a non-
canonical Hamiltonian system. By using the coordinate increment discrete
gradient method to solve the Hamiltonian system in Section. 3, we obtain a new
energy-preserving method and its adjoint method. Then the two methods are
combined into a new method, namely CIDG-C method. Numerical experiments are
presented in Section. 4 to compare the theoretical results of the new CIDG-C
method with the Boris method [18]. Section. 5 concludes this paper.

##  2 Non-canonical Hamiltonian form of the system

In this section, we firstly review the Non-canonical Hamiltonian form of the
Lorentz force system [13, 14, 15]. For a charged particle with position
$\textbf{x}(t)=(x(t),y(t),z(t))^{T}\in\mathbb{R}^{3}$ in an electromagnetic
field, its dynamics is described by the following Newton-Lorentz equation

$\frac{d^{2}\textbf{x}}{dt^{2}}=\frac{d\textbf{x}}{dt}\times\textbf{B}(\textbf{x})-\nabla
U(\textbf{x}),$ (2.1)

where $-\nabla U(\textbf{x}):=\textbf{E}(\textbf{x})$ is an electric field,
and
$\textbf{B}(\textbf{x})=(B_{1}(\textbf{x}),B_{2}(\textbf{x}),B_{3}(\textbf{x}))^{T}$
is a magnetic field.

Let $\textbf{v}=(v_{1},v_{2},v_{3})^{T}=\frac{d\textbf{x}}{dt}$ and

$\textbf{S}(\textbf{x})=\left(\begin{array}[]{ccc}0&B_{3}(\textbf{x})&-B_{2}(\textbf{x})\\\\[8.5359pt]
-B_{3}(\textbf{x})&0&B_{1}(\textbf{x})\\\\[8.5359pt]
B_{2}(\textbf{x})&-B_{1}(\textbf{x})&0\end{array}\right)=-\textbf{S}^{T}(\textbf{x}),$

in (2.1), then we can obtain the following ODEs

$\displaystyle\frac{d\textbf{x}}{dt}=\textbf{v},$ (2.2)
$\displaystyle\frac{d\textbf{v}}{dt}=\textbf{S}(\textbf{x})\textbf{v}-\nabla
U(\textbf{x}),$ (2.3)
$\displaystyle(\textbf{x}(0),\textbf{v}(0))=(\textbf{x}_{0},\textbf{v}_{0})\in\Omega\subseteq\mathbb{R}^{3}\times\mathbb{R}^{3}.$
(2.4)

Denote
$\textbf{z}=(\textbf{x}^{T},\textbf{v}^{T})^{T}=(x,y,z,v_{1},v_{2},v_{3})^{T}$,
the ODEs (2.2)-(2.4) can be written as a non-canonical Hamiltonian system

$\frac{d\textbf{z}}{dt}=K(\textbf{z})\nabla H(\textbf{z}),$ (2.5)

where

$K(\textbf{z})=\left(\begin{array}[]{cc}0&I\\\\[8.5359pt]
-I&\textbf{S}(\textbf{x})\end{array}\right)$

is a skew-symmetrical matrix, and
$H(\textbf{z})=\frac{1}{2}\textbf{v}^{T}\textbf{v}+U(\textbf{x})$ is the
Hamiltonian energy of the system, which is conserved by

$\frac{d}{dt}H(\textbf{z})=\nabla
H^{T}(\textbf{z})\frac{d\textbf{z}}{dt}=\nabla
H^{T}(\textbf{z})K(\textbf{z})\nabla H(\textbf{z})=0.$ (2.6)

##  3 Coordinate increment discrete gradient methods

Notice that the Hamiltonian energy of the system (2.5) can be exactly
preserved by the coordinate increment discrete gradient method [30]

$\frac{\textbf{z}^{n+1}-\textbf{z}^{n}}{h}=K(\frac{\textbf{z}^{n+1}+\textbf{z}^{n}}{2})\bar{\nabla}H(\textbf{z}^{n+1},\textbf{z}^{n}),$
(3.1)

i.e.,

$\frac{\textbf{z}^{n+1}-\textbf{z}^{n}}{h}=K(\frac{\textbf{z}^{n+1}+\textbf{z}^{n}}{2})\left(\begin{array}[]{c}\frac{H(z^{n+1}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{1}-z^{n}_{1}}\\\
\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})-H(z^{n+1}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{2}-z^{n}_{2}}\\\
\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})-H(z^{n+1}_{1},z^{n+1}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{3}-z^{n}_{3}}\\\
\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n}_{5},z^{n}_{6})-H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{4}-z^{n}_{4}}\\\
\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n}_{6})-H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{5}-z^{n}_{5}}\\\
\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})-H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n}_{6})}{z^{n+1}_{6}-z^{n}_{6}}\end{array}\right),$

where $h$ is the step size, and we consider the coordinate increment discrete
gradient as

$\bar{\nabla}H(\bar{\textbf{z}},\textbf{z})=\left(\begin{array}[]{c}\frac{H(\bar{z}_{1},z_{2},z_{3},z_{4},...,z_{m})-H(z_{1},z_{2},z_{3},z_{4},...,z_{m})}{\bar{z}_{1}-z_{1}}\\\
\frac{H(\bar{z}_{1},\bar{z}_{2},z_{3},z_{4},...,z_{m})-H(\bar{z}_{1},z_{2},z_{3},z_{4},...,z_{m})}{\bar{z}_{2}-z_{2}}\\\
\frac{H(\bar{z}_{1},\bar{z}_{2},\bar{z}_{3},z_{4},...,z_{m})-H(\bar{z}_{1},\bar{z}_{2},z_{3},z_{4},...,z_{m})}{\bar{z}_{3}-z_{3}}\\\
\vdots\\\
\frac{H(\bar{z}_{1},\bar{z}_{2},...,\bar{z}_{n-1},\bar{z}_{m})-H(\bar{z}_{1},\bar{z}_{2},...,\bar{z}_{m-1},z_{m})}{\bar{z}_{m}-z_{m}}\end{array}\right),$
(3.2)

with the notation $\bar{\textbf{z}}=\textbf{z}^{n+1}$,
$\textbf{z}=\textbf{z}^{n}$ and $m=6$. We denote the method (3.1) as
coordinate increment discrete gradient - I method (CIDG-I).

Then we consider the adjoint method of CIDG-I (3.1), which we denote as the
coordinate increment discrete gradient - II method (CIDG-II) of system (2.5)

$\frac{\textbf{z}^{n+1}-\textbf{z}^{n}}{h}=K(\frac{\textbf{z}^{n+1}+\textbf{z}^{n}}{2})\bar{\nabla}H(\textbf{z}^{n},\textbf{z}^{n+1}),$
(3.3)

i.e.,

$\frac{\textbf{z}^{n+1}-\textbf{z}^{n}}{h}=K(\frac{\textbf{z}^{n+1}+\textbf{z}^{n}}{2})\left(\begin{array}[]{c}\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})-H(z^{n}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})}{z^{n+1}_{1}-z^{n}_{1}}\\\
\frac{H(z^{n}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})}{z^{n+1}_{2}-z^{n}_{2}}\\\
\frac{H(z^{n}_{1},z^{n}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})}{z^{n+1}_{3}-z^{n}_{3}}\\\
\frac{H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n+1}_{5},z^{n+1}_{6})}{z^{n+1}_{4}-z^{n}_{4}}\\\
\frac{H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n+1}_{5},z^{n+1}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n+1}_{6})}{z^{n+1}_{5}-z^{n}_{5}}\\\
\frac{H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n+1}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{6}-z^{n}_{6}}\end{array}\right).$

Finally the CIDG-I method $\Phi_{h}$ (3.1) and its adjoint method CIDG-II
$\Phi^{*}_{h}$ (3.3) can be combined [6] into the following coordinate
increment discrete gradient composition (CIDG-C) method

$\Psi_{h}=\Phi_{\frac{h}{2}}\circ\Phi^{*}_{\frac{h}{2}},$ (3.4)

i.e.,

$\displaystyle\frac{\textbf{z}-\textbf{z}^{n}}{\frac{1}{2}h}=K(\frac{\textbf{z}+\textbf{z}^{n}}{2})\bar{\nabla}H(\textbf{z}^{n},\textbf{z}),$
(3.5)
$\displaystyle\frac{\textbf{z}^{n+1}-\textbf{z}}{\frac{1}{2}h}=K(\frac{\textbf{z}^{n+1}+\textbf{z}}{2})\bar{\nabla}H(\textbf{z}^{n+1},\textbf{z}),$
(3.6)

$n=0,1,2,...,N.$

######  Theorem 3.1.

CIDG-I, CIDG-II and CIDG-C methods can all conserve the Hamiltonian energy of
the system (2.5). Moreover, the CIDG-C method is symmetrical and has order 2.

###### Proof.

From the CIDG-I method (3.1), we obtain

$\displaystyle 0=$ $\displaystyle
h\bar{\nabla}H(\textbf{z}^{n+1},\textbf{z}^{n})^{T}K(\frac{\textbf{z}^{n+1}+\textbf{z}^{n}}{2})\bar{\nabla}H(\textbf{z}^{n+1},\textbf{z}^{n})=\bar{\nabla}H(\textbf{z}^{n+1},\textbf{z}^{n})^{T}(\textbf{z}^{n+1}-\textbf{z}^{n})$
$\displaystyle=$
$\displaystyle\frac{H(z^{n+1}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{1}-z^{n}_{1}}(z^{n+1}_{1}-z^{n}_{1})+...$
$\displaystyle+\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})-H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n}_{6})}{z^{n+1}_{6}-z^{n}_{6}}(z^{n+1}_{6}-z^{n}_{6})$
$\displaystyle=$
$\displaystyle[H(z^{n+1}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})]+[H(z^{n+1}_{1},z^{n+1}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})$
$\displaystyle-H(z^{n+1}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})]+...+[H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})$
$\displaystyle-H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n}_{6})]$
$\displaystyle=$ $\displaystyle
H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6}),$

i.e., the energy conservation of the method is obtained by

$H(\textbf{z}^{n+1})=H(\textbf{z}^{n}),\quad n=0,1,2,...,N.$

On the other hand, the energy conservation of the CIDG-II method (3.3) can
also be obtained by

$\displaystyle 0=$ $\displaystyle
h\bar{\nabla}H(\textbf{z}^{n},\textbf{z}^{n+1})^{T}K(\frac{\textbf{z}^{n}+\textbf{z}^{n+1}}{2})\bar{\nabla}H(\textbf{z}^{n},\textbf{z}^{n+1})=\bar{\nabla}H(\textbf{z}^{n},\textbf{z}^{n+1})^{T}(\textbf{z}^{n+1}-\textbf{z}^{n})$
$\displaystyle=$
$\displaystyle\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})-H(z^{n}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})}{z^{n+1}_{1}-z^{n}_{1}}(z^{n+1}_{1}-z^{n}_{1})+...$
$\displaystyle+\frac{H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n+1}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{6}-z^{n}_{6}}(z^{n+1}_{6}-z^{n}_{6})$
$\displaystyle=$ $\displaystyle
H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n+1}_{4},z^{n+1}_{5},z^{n+1}_{6})-H(z^{n}_{1},z^{n}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})=H(\textbf{z}^{n+1})-H(\textbf{z}^{n}).$

Because CIDG-I and CIDG-II are both energy-preserving methods, the CIDG-C
method (3.4) is also obviously energy-preserved.

Finally, take $n+1=n$, and $h=-h$ in the CIDG-C method (3.5)-(3.6), then we
have

$\displaystyle\frac{\textbf{z}-\textbf{z}^{n+1}}{-\frac{1}{2}h}=K(\frac{\textbf{z}+\textbf{z}^{n+1}}{2})\bar{\nabla}H(\textbf{z}^{n+1},\textbf{z}),$
$\displaystyle\frac{\textbf{z}^{n}-\textbf{z}}{-\frac{1}{2}h}=K(\frac{\textbf{z}^{n}+\textbf{z}}{2})\bar{\nabla}H(\textbf{z}^{n},\textbf{z}),$

which is just equal to (3.5)-(3.6). Therefore, the CIDG-C method is
symmetrical and has order 2 [6].∎

######  Remark 3.1.

In practical calculations, in order to avoid that $z^{n+1}_{i}-z^{n}_{i}$ is
used as the denominator and equals to $0$, we take (for instance)
$\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})-H(z^{n+1}_{1},z^{n+1}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{3}-z^{n}_{3}}$
as $\lim\limits_{\epsilon\to
0}\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3}+\epsilon,z^{n}_{4},z^{n}_{5},z^{n}_{6})-H(z^{n+1}_{1},z^{n+1}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{3}+\epsilon-z^{n}_{3}}$,
which is equal to
$\frac{H(z^{n+1}_{1},z^{n+1}_{2},z^{n+1}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})-H(z^{n+1}_{1},z^{n+1}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})}{z^{n+1}_{3}-z^{n}_{3}}$
when $(z^{n+1}_{3}-z^{n}_{3})\neq 0$
and equal to $\frac{\partial H}{\partial
z_{3}}(z^{n+1}_{1},z^{n+1}_{2},z^{n}_{3},z^{n}_{4},z^{n}_{5},z^{n}_{6})$ when
$(z^{n+1}_{3}-z^{n}_{3})=0$. Note that since $(z^{n+1}_{3}-z^{n}_{3})=0$, the
preservation of the Hamiltonian is guaranteed.

######  Remark 3.2.

Applying Eqs. (3.5)-(3.6) into the Eqs. (2.2)-(2.4), we can obtain the
following equivalent form of the CIDG-C method

$\displaystyle\frac{\textbf{x}-\textbf{x}^{n}}{\tau}=\frac{\textbf{v}+\textbf{v}^{n}}{2},$
$\displaystyle\frac{\textbf{v}-\textbf{v}^{n}}{\tau}=S(\frac{\textbf{x}+\textbf{x}^{n}}{2})\frac{\textbf{v}+\textbf{v}^{n}}{2}-\left(\begin{array}[]{c}\frac{U(x,y^{n},z^{n})-U(x^{n},y^{n},z^{n})}{x-x^{n}}\\\
\frac{U(x,y,z^{n})-U(x,y^{n},z^{n})}{y-y^{n}}\\\
\frac{U(x,y,z)-U(x,y,z^{n})}{z-z^{n}}\end{array}\right),$
$\displaystyle\frac{\textbf{x}^{n+1}-\textbf{x}}{\tau}=\frac{\textbf{v}^{n+1}+\textbf{v}}{2},$
$\displaystyle\frac{\textbf{v}^{n+1}-\textbf{v}}{\tau}=S(\frac{\textbf{x}^{n+1}+\textbf{x}}{2})\frac{\textbf{v}^{n+1}+\textbf{v}}{2}-\left(\begin{array}[]{c}\frac{U(x^{n+1},y^{n+1},z^{n+1})-U(x,y^{n+1},z^{n+1})}{x^{n+1}-x}\\\
\frac{U(x,y^{n+1},z^{n+1})-U(x,y,z^{n+1})}{y^{n+1}-y}\\\
\frac{U(x,y,z^{n+1})-U(x,y,z)}{z^{n+1}-z}\end{array}\right),$

where $\tau=\frac{h}{2}$, $n=0,1,2,...,N.$

##  4 Numerical experiments

In this section, we numerically test the CIDG-C method (3.5)-(3.6)and compare
the results with the effective Boris method [18], the Boole discrete line
integral (BDLI) method [24] and the line integral method with parameters (2,2)
(which we denote as LIM(2,2)) [25]. All numerical tests have been done on a
2.4 GHz Intel core i5 computer with 16 GB of memory, running Octave 8.2.0.

###  4.1 2D dynamics in a static electromagnetic field

Firstly, we consider the 2D dynamics of the charged particle in a static, non-
uniform electromagnetic field which is a significant application in the study
of the single particle motion and the guiding center dynamics. In this case
the Lorentz force system is a non-polynomial Hamiltonian system. Theoretical
analysis shows that the analytic orbit of the charged particle is a spiraling
circle which have a constant radius. The large circle corresponds to the
$\nabla\cdot\textbf{B}$ drift and the $\textbf{E}\times\textbf{B}$ drift of
the guiding center, and the small circle is a fast-scale gyromotion [12, 15].

The static, non-uniform electromagnetic field is taken as

$\textbf{B}=\nabla\times\textbf{A}=R\textbf{e}_{z},\quad\textbf{E}=-\nabla
U=\frac{10^{-2}}{R^{3}}(x\textbf{e}_{x}+y\textbf{e}_{y}),$ (4.1)

where the potentials are chosen to be
$\textbf{A}=\frac{R^{2}}{3}R\textbf{e}_{\xi}$, $U=\frac{10^{-2}}{R}$ in
cylindrical coordinates $(R,\xi,z)$ with $R=\sqrt{x^{2}+y^{2}}$. In this
example, we take the physical quantities $q=1$, $m=1$, which are normalized by
the magnetic field B and the electric field E. The Hamiltonian energy of the
system (2.5)

$H=\frac{1}{2}\textbf{v}^{T}\textbf{v}+\frac{0.01}{\sqrt{x^{2}+y^{2}}}$ (4.2)

and the angular momentum

$p=R^{2}\frac{d}{dt}\xi+\frac{1}{3}R^{3}$ (4.3)

are constant [15]. As the given electromagnetic field (4.1) changes slowly
with respect to the spatial period of the motion, the magnetic moment
$\mu=\frac{\textbf{v}_{\bot}^{2}}{2R}$ is an adiabatic invariant, where
$\textbf{v}_{\bot}$ is the component of v perpendicular to B.

We consider the initial position $\textbf{x}_{0}=(0,1,0)^{T}$ with the initial
velocity $\textbf{v}_{0}=(0.1,0.01,0)^{T}$, and the step size is $h=\pi/10$
which is one-twentieth of the characteristic gyro-period $2\pi$.

(a) (b)

Figure 1: The Boris method is applied to the simple 2D dynamics with step
$h=\pi/10$. (a) The orbit in the 100-st turn; (b) Errors of the angular
momentum $p_{\xi}$, the magnetic moment $\mu$ and the energy $H$ for
$t\in[0,5\times 10^{5}h]$.

(a) (b)

Figure 2: The BDLI method is applied to the simple 2D dynamics with step
$h=\pi/10$. (a) The orbit in the 100-st turn; (b) Errors of the angular
momentum $p_{\xi}$, the magnetic moment $\mu$ and the energy $H$ for
$t\in[0,5\times 10^{5}h]$.

We apply the Boris methods, the BDLI method, the LIM(2,2) method and the
CIDG-C method with the same step and compare the results through numerical
simulation. Fig. 2(a) shows the numerical orbits of the Boris method in the
100-st turn. Fig. 2(b) shows the error of the angular momentum $p_{\xi}$, the
magnetic moment $\mu$ and the energy $H$ for $t\in[0,5\times 10^{5}h]$. The
numerical results demonstrate that the errors of the three invariants are
bounded for a long integration time.

Fig. 2(a) and Fig. 4(a) shows the numerical orbits of the BDLI method and the
LIM(2,2) method in the 100-st turn. Fig. 2(b) and Fig. 4(b) shows the error of
the invariants of this two methods for $t\in[0,5\times 10^{5}h]$. The BDLI
method use Boole’s rule to calculate the integrand numerically and the
LIM(2,2) method use Guass-Legendre formulas with two points. They are
symmetrical and have order 2. The energy can be conserved approximatively and
the error of the energy will become smaller until round-off error if we use
higher order numerical integration formula.

(a) (b)

Figure 3: The LIM(2,2) method is applied to the simple 2D dynamics with step
$h=\pi/10$. (a) The orbit in the 100-st turn; (b) Errors of the angular
momentum $p_{\xi}$, the magnetic moment $\mu$ and the energy $H$ for
$t\in[0,5\times 10^{5}h]$.

(a) (b)

Figure 4: The CIDG-C method is applied to the simple 2D dynamics with step
$h=\pi/10$. (a) The orbit in the 100-st turn; (b) Errors of the angular
momentum $p_{\xi}$, the magnetic moment $\mu$ and the energy $H$ for
$t\in[0,5\times 10^{5}h]$.

Fig. 4(a) shows the numerical orbits of the CIDG-C method in the 100-st turn.
Fig. 4(b) shows the error of the invariants for $t\in[0,5\times 10^{5}h]$.
From the figure we can know that the error of the angular momentum $p_{\xi}$
and the magnetic moment $\mu$ are also bounded for a long integration time.
What is more, the CIDG-C method can exactly conserve the energy of the system.

And we compare the CPU times of the four methods from $t=0$ to $t=2\times
10^{4}h$ with the same step $h=\pi/10$. The Boris method takes the least
amount of time. The Boris method is essentially explicit but it does not
preserve the energy as well as the other methods. The LIM(2,2) method is 1.56
times faster than the CIDG-C method and 1.89 times faster than the BDLI
method. The BDLI method use Boole’s rule to calculate the integrand
numerically with five points and the LIM(2,2) method use Guass-Legendre
formulas with two points. The BDLI method use more points than the LIM(2,2)
method, so the BDLI method takes more time.Table 1

Table 1: CPU times of the four methods from $t=0$ to $t=2\times 10^{4}h$, with the same step $h=\pi/10$. Boris method | BDLI method | LIM(2,2) method | CIDG-C method
---|---|---|---
3.047 s | 91.656 s | 48.438 s | 75.531 s

###  4.2 Energy behavior test

Secondly, we consider the system (2.2)-(2.4) with

$\displaystyle U(\textbf{x})=x^{3}-y^{3}+\frac{1}{5}x^{4}+y^{4}+z^{4},$ (4.4)
$\displaystyle\textbf{B}(\textbf{x})=(0,0,\sqrt{x^{2}+y^{2}})^{T},$ (4.5)
$\displaystyle\textbf{x}(0)=(0,1,0.1)^{T},\textbf{v}(0)=(0.09,0.55,0.3)^{T},$
(4.6)

for which in [23, 25] it has been proven that the Boris method exhibits a
$O(th^{2})$ drift in the energy. Fig. 5(a) shows the energy drift of the Boris
method with the step size $h=10^{-2}$ over the interval $[0,3\times 10^{4}]$.
The CIDG-C method can conserve the energy with the same step size and there is
not energy drift, which is shown in Fig. 5(b). The LIM(2,2) method and the
BDLI method can also conserve the energy with the same step size, which is
shown in Fig. 5(c) and Fig. 5(d).

(a) (b)
(c) (d)

Figure 5: Hamiltonian error solving problem(4.4)-(4.6) with step $h=10^{-2}$.
(a) The Boris method; (b)The CIDG-C method ;(c)The LIM(2,2) method ;(d)The
BDLI method.

######  Remark 4.1.

Errata to A novel directly energy-preserving method for charged particle
dynamics [31]: The Fig. 5(b) in the original paper contains an error, as the
magnitude of the error on the vertical axis is incorrect. This is because we
accidentally uploaded the wrong image when submitting the final LaTeX file to
the editor. The correct image is included as Fig. 5(b) in this article.
Actually, we are more concerned that the CIDG-C method did not has an energy
drift in this example, which is why the error in the Fig. 5(b) was not
detected earlier. We apologize for our mistake.

###  4.3 2D dynamics in an axisymmetric tokamak geometry

Finally, we consider the motion of a charged particle in 2-dimensional
axisymmetric tokamak geometry without inductive electric field, which is also
an important application of the guiding center dynamics [12, 15].

The magnetic field in the toroidal coordinates $(r,\theta,\xi)$ is expressed
as

$\textbf{B}=\frac{B_{0}r}{qR}\textbf{e}_{\theta}+\frac{B_{0}R_{0}}{R}\textbf{e}_{\xi},$
(4.7)

where $B_{0}=1$, $R_{0}=1$, and $q=2$ are constant with their usual meanings.
The corresponding vector potential A is chosen to be

$\textbf{A}=\frac{z}{2R}\textbf{e}_{R}+\frac{(1-R)^{2}+z^{2}}{4R}\textbf{e}_{\xi}+\frac{\ln
R}{2}\textbf{e}_{z}.$ (4.8)

In this example, the physical quantities are normalized as usual. The energy
of the system is $H=\frac{1}{2}\textbf{v}^{T}\textbf{v}$. With the conserved
quantities, the solution orbit projected on $(R,z)$ space forms a closed
orbit. In order to use the CIDG-C method, B (4.7) is transformed to the
Cartesian coordinates $(x,y,z)$ which yields

$\textbf{B}=-\frac{2y+xz}{2R^{2}}\textbf{e}_{x}+\frac{2x-yz}{2R^{2}}\textbf{e}_{y}+\frac{R-1}{2R}\textbf{e}_{z}.$
(4.9)

(a) (b)
(c)

Figure 6: Numerical solution of the CIDG-C method with step $h=\pi/10$ for
$t\in[0,5\times 10^{5}h]$. (a) Transit orbit; (b) Banana orbit; (c) Energy
preservation.

We firstly consider the initial position $\textbf{x}_{0}=(1.05,0,0)^{T}$ and
the initial velocity $\textbf{v}_{0}=(0,2\times 4.816\times
10^{-4},2.059\times 10^{-3})^{T}$, the orbit projected on $(R,z)$ space is a
transit orbit, and it will turn to a banana orbit when the initial velocity is
changed to $\textbf{v}_{0}=(0,4.816\times 10^{-4},2.059\times 10^{-3})^{T}$.
We apply the CIDG-C method and use the step size $h=\pi/10$ which is the
$1/20$ of the gyro-period $2\pi$. The simulation over $5\times 10^{5}$ steps
is shown in Fig. 6. Fig. 6(a) shows the transit orbit solution, Fig. 6(b)
shows the banana orbit solution and Fig. 6(c) shows the energy preservation.
The CIDG-C method gives correct orbits for a very long integration time, and
also maintain the energy precisely $H$ .

##  5 Conclusions

The motion of single-particle, which satisfies the Lorentz force system, has a
major role in plasmas. The Hamiltonian energy is one of the most remarkable
features to characterize the system. In this paper, we apply the coordinate
increment discrete gradient method to the Hamiltonian system, then we derive
the CIDG-I method. The new CIDG-I method and its adjoint method are both
energy-preserving methods, and can be combined into a new CIDG-C method. The
CIDG-C method is symmetrical and can exactly conserve the Hamiltonian energy
without numerical quadrature formula. Numerical results have shown that the
CIDG-C method surpasses the Boris method in long time simulation.

Acknowledgements
This work is supported by the National Natural Science Foundation of China
(No. 12101072) and the Natural Science Foundation of Chongqing (CSTB2022NSCQ-
MSX0890).

Reference

## References

  * [1] R. Dodd, J. Eibeck, J. Gibbon, H. Morris, Solitons and Nonlinear Wave Equation, Academic Press, 1982.
  * [2] K. Feng, Difference schemes for Hamiltonian formalism and symplectic geometry, J. Comput. Math., 4 (1986) 279-289.
  * [3] K. Feng, M. Qin, Symplectic Geometric Algorithms for Hamiltonian Systems, Springer-Verlag/Zhejiang Science and Technology Publishing House, Berlin/Hangzhou, 2010.
  * [4] L. Brugnano, F. Iavernaro, Line Integral Methods for Conservative Problems, CRC Press, Boca Raton, FL, 2016.
  * [5] J. Chen, M. Qin, Multi-symplectic Fourier pseudospectral method for the nonlinear Schrödinger equation, Electron. Trans. Numer. Anal., 12 (2001) 193-204.
  * [6] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd edition, Springer-Verlag, Berlin, 2006.
  * [7] E. Hairer, S. N$\phi$rsett, G. Wanner, Solving Ordinary Differential Equations. I. Nonstiff Problems, 2nd edition, Springer, Berlin, 1993.
  * [8] L. Kong, J. Hong, J. Zang, Splitting multi-symplectic integrators for Maxwell’s equation, J. Comput. Phys., 229 (2010) 4259-4278.
  * [9] Y. Wang, Q. Li, Y. Song, Two new simple multi-symplectic schemes for the nonlinear Schrödinger equation, Chinese Phys. Lett., 25 (2008) 1538-1540.
  * [10] Q. Xu, S. Song, Y. Chen, A semi-explicit multi-symplectic splitting scheme for a 3-coupled nonlinear Schrödinger equation, Comput. Phys. Comm. 185 (2014) 1255-1264.
  * [11] P. Bellan, Fundamentals of Plasma Physics, 1st edition, Cambridge University Press, 2008.
  * [12] H. Qin, S. Zhang, J. Xiao, J. Liu, Y. Sun, W. Tang, Why is Boris algorithm so good ?, Phys. Plasmas, 20 (2013) 084503.
  * [13] R. Littlejohn, Hamiltonian formulation of guiding center motion, Phys. Fluids, 24 (1981) 1730-1749.
  * [14] P. Morrison, J. Greene, Noncanonical Hamiltonian density formulation of hydrodynamics and ideal magnetohydrodynamics, Phys. Rev. Lett., 45 (1980) 790-794.
  * [15] Y. He, Y. Sun, J. Liu and H. Qin, Volume-preserving algorithms for charged particle dynamics, J. Comput. Phys., 281 (2015), 135-147.
  * [16] R. Zhang, J. Liu, H. Qin, Y. Wang, Y. He, Y. Sun, Volume-preserving algorithm for secular relativistic dynamics of charged particles, Phys. Plasmas, 22 (2015) 044501.
  * [17] C. Birdsall, A. Langdon, Plasma Physics via Computer Simulation, Ser. Plasma Phys., Taylor $\&$ Francis, 2005.
  * [18] J. Boris, Proceedings of the Fourth Conference on numerical simulation of plasmas, Naval Research Laboratory, Washington D.C., 1970, 3-67.
  * [19] P. Stoltz, J. Cary, G. Penn, J. Wurtele, Efficiency of a Boris like integration scheme with spatial stepping, Phys. Rev. Spec. Top., Accel. Beams 5 (2002) 094001.
  * [20] R. Littlejohn, Variational principles of guiding centre motion, J. Plasma Phys., 29 (1983) 111-125.
  * [21] H. Qin, X. Guan, Variational symplectic integrator for long-time simulations of the guiding-center motion of charged particles in general magnetic fields, Phys. Rev. Lett., 100 (2008) 035006.
  * [22] H. Qin, X. Guan, W. Tang, Variational symplectic algorithm for guiding center dynamics and its application in tokamak geometry, Phys. Plasmas, 16 (2009) 042510.
  * [23] E. Hairer, C. Lubich, Energy behavior of the Boris method for charged-particle dynamics, BIT, 58 (2018) 969-979.
  * [24] H. Li, Y. Wang, A discrete line integral method of order two for the Lorentz force system, Appl. Math. Comput., 291 (2016) 207-212.
  * [25] L. Brugnano, J. I. Montijano, L. Randez, High-order energy-conserving Line Integral Methods for charged particle dynamics, J. Comput. Phys., 396 (2019) 209-227.
  * [26] H. Li, Q. Hong, An efficient energy-preserving algorithm for the Lorentz force system, Appl. Math. Comput., 358 (2019) 161-168.
  * [27] O. Gonzalez, Time integration and discrete Hamiltonian systems, J. Nonlinear Sci., 6 (1996) 449-467.
  * [28] T. Matsuo, High-order schemes for conservative or dissipative systems, J. Comput. Appl. Math., 152 (2003) 305-317.
  * [29] G. Quispel, D. McLaren, A new class of energy-preserving numerical integration methods, J. Phys. A: Math. Theor., 41 (2008) 045206.
  * [30] T. Itoh, K. Abe, Hamiltonian-conserving discrete canonical equations based on variational difference quotients, J. Comput. Phys., 77 (1988) 85-102.
  * [31] Y. Li, P. Jiang, H. Li, A novel directly energy-preserving method for charged particle dynamics, J. Comput. Appl. Math., 446 (2024) 115869.