SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

1

Avoiding overfitting of multilayer perceptrons by
training derivatives

arXiv:1802.10301v1 [] 28 Feb 2018

V.I. Avrutskiy

Abstract—Resistance to overfitting is observed for neural
networks trained with extended backpropagation algorithm. In
addition to target values, its cost function uses derivatives of
those up to the 4th order. For common applications of neural
networks, high order derivatives are not readily available, so
simpler cases are considered: training network to approximate
analytical function inside 2D and 5D domains and solving Poisson
equation inside a 2D circle. For function approximation, the cost
is a sum of squared differences between output and target as
well as their derivatives with respect to the input. Differential
equations are usually solved by putting a multilayer perceptron
in place of unknown function and training its weights, so that
equation holds within some margin of error. Commonly used
cost is the equation’s residual squared. Added terms are squared
derivatives of said residual with respect to the independent
variables. To investigate overfitting, the cost is minimized for
points of regular grids with various spacing, and its root mean is
compared with its value on much denser test set. Fully connected
perceptrons with six hidden layers and 2 · 104 , 1 · 106 and 5 · 106
weights in total are trained with Rprop until cost changes by
less than 10% for last 1000 epochs, or when the 10000th epoch is
reached. Training the network with 5 · 106 weights to represent
simple 2D function using 10 points with 8 extra derivatives in
each produces cost test to train ratio of 1.5, whereas for classical
backpropagation in comparable conditions this ratio is 2 · 104 .
Index Terms—Neural networks, overfitting, partial differential
equations, high order derivatives, function approximation.

I. I NTRODUCTION
VERFITTING is a common fault of neural networks
that occurs in many different applications [21], [12], [7],
[11]. Consider a case of supervised learning: network N is to
be trained to represent vector function F : X → Y known
only by a finite set of its arguments xa ⊂ X and corresponding
outputs ya ⊂ Y . It is a common practice to subdivide xa ,ya
into at least two sets: training xb ,yb and test xc ,yc . Whatever
algorithm is used for minimization, it is not run for elements of
xc ,yc . This is done to later observe network’s performance on
previously unknown data. Overfitting can be broadly described
as different behavior of cost function E = ||N (~x) − ~y || for
~x ∈ xb , ~y ∈ yb and ~x ∈ xc , ~y ∈ yc . For example, its average
values on those two sets can be orders of magnitude away.
In broad terms, the ratio between cost function averaged on
test and train set shows how much attention training procedure
pays to the actual input rather than to the function behind it.
However, this relation between averages is not important all
by itself, since the absolute values of Ec on test set xc ,yc are
the goal of training procedure.

O

V.I. Avrutskiy is with the Department of Aeromechanics and Flight Engineering of Moscow Institute of Physics and Technology, Institutsky lane 9,
Dolgoprudny, Moscow region, 141700, e-mail: avrutsky@phystech.edu

There are numerous ways to tackle with overfitting. It is
possible to directly address negative effects for Ec by introducing a special set of patterns xcv ,ycv called cross validation
set. Its elements are not fed to minimization algorithm as well,
but during the training, cost function for them is observed.
When it starts to increase, while the cost on xb ,yb continues
to lower, one can draw a conclusion that negative effect of
overfitting started to overpower positive effects of the training
itself. If one choses to stop, the decision effectively puts
a lower bond on approximation error. One might also say
that before Ec started to increase, any presence of overfitting
was not practically significant. Many ways to handle such
situations were developed [17], [13], [4]. For example, neural
network can be pruned [15], or some special statistical rule for
updating weights can be used [6], [20], [22]. Generally, the
more approximation abilities a neural network has the more it
is prone to overfitting.
If random noise is present and the network is overtrained,
i.e. backpropagation algorithm was run long enough to start
paying too much attention to the training set, it almost inevitably leads to worsening of performance on the test set,
since noise on the one is not correlated with noise on the
other. However, if data is noiseless, then even having an
extremely low error on training points can bring no negative
effects to cost on the test set. This paper is using various
derivatives of target function and, therefore, is unable to take
noise into consideration, since even the slightest variations
of neighboring points can have profound consequences to
numerically calculated derivatives even of the first, not to
mention higher orders. Vice versa, if derivatives up to say the
4th order are known then nearby points can be calculated by
Taylor approximation virtually without any noise. Therefore,
in the scope of this paper overfitting is defined not as a
worsening of performance on the test set, which was not
observed, but rather as an ability of network to pay more
attention to the training set. Thus, an absence of overfitting can
be expressed as an inability to distinguish between statistical
properties of the cost function on the test and training sets.
II. E XTENDED BACKPROPAGATION
Algorithms for supervised training of neural networks with
cost functions that include derivatives were described in a
number of publications [3], [9]. Most notable are double
backpropagation [2] and tangent propagation [19]. They both
consider image classification problem and their cost functions
in addition to squared difference between output and target
classes include the first order derivatives of that difference.

© 2018 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or
reuse of any copyrighted component of this work in other works.

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

Double backpropagation calculates them with respect to values
of individual input pixels, so derivatives of output can be
viewed as a slope of network’s reaction to subtle distortions.
Derivatives of target with respect to distortions are zero since
they do not change class, thus, training tries to minimize
derivatives of output with respect to input for each meaningful
pattern. Tangent propagation is more elaborate and calculates
derivatives along special directions of input space. For each
image I, a direction is another pseudo-image J that can be
multiplied by small factor α and added to I, so that I+αJ simulates an infinitesimal geometrical distortion of I, for example,
a rotation through small angle. Such distortion do not change
class either, so derivative of output with respect to α is trained
to be zero. Described methods are rarely used in practice,
which might be due to their limited benefits, an increase in
training time and a tricky process of calculating derivatives.
This study does not consider complicated mappings like image
classifiers, and tasks are limited to somewhat artificial, where
functions are from space C∞ , and their derivatives are easily
computable. This allows to closely investigate effects of using
high order terms on perceptron training.
Another case when derivatives appear in cost function, is
solving differential equations [8], [10], [14], [18]. Classical
approach is to construct a neural network that maps each vector of independent variables to solution. Consider a boundary
value problem inside a region Γ for function u(x, y) written
as:
U (x, y, u, ux , uy , uxx , ...) = 0
u|∂Γ = f
The first step is to make a substitution using relation:
u(x, y) = v(x, y) · φ(x, y) + f

(1)

where φ(x, y) is a known function that vanishes on the
boundary ∂Γ and is non zero inside Γ, and v(x, y) is to be
constructed by a neural network [10]. The term f is now a
smooth continuation of the boundary condition into Γ that is
supposed to exist. The equation is written as:
V (x, y, v, vx , vy , vxx , ...) = 0
and the boundary condition is simply:
v|∂Γ < ∞
A network with two inputs, one output and a suitable number
of hidden layers and neurons is then created. Its weights are
initialized and trained to minimize cost function E - a measure
of a local cost function e = V 2 on Γ calculated using a
finite set of grid points. If it reaches small enough values,
v(x, y) · φ(x, y) + f can be considered as an approximate
numerical solution, and as long as the procedure converges
to a finite v, the boundary condition is satisfied. An extension
of classical backpropagation is required due to the presence in
e terms like vxx , which are derivatives of the neural network’s
output with respect to input. Various types of differential equations can be solved using this approach [9], [23]. The smaller
effect of overfitting is, the fewer points are required to capture
the behavior of function e on Γ for proper discretization of E
and, therefore, to solve the equation in the entire region.

2

Study [1] found significant improvements from using high
order derivatives in learning process. For classical supervised
training, the cost, which is equal to squared difference between
output and target, was added with squared differences of
their derivatives with respect to input up to the 4th order.
In the study target function was chosen as a piece of 2D
Fourier series and all derivatives were readily available. Such
modifications were found to enhance precision up to ∼ 100
times and allowed to use three times less points for training,
but required to increase the number of neurons and layers
as well. Another obvious cost of such enhancements was that
derivatives of target must be known beforehand. For the case of
solving partial differential equations the situation is different.
The expression V behind the cost function e is simply all
terms of the equation put together, therefore, it is possible to
introduce extra derivatives simply by applying operators like
∂/∂x, ∂/∂y to V and adding squares of the results to the local
cost function. For partial differential equations inclusion of
extra derivatives up to the second order allowed to use 8 times
less grid points and obtain 5 times smaller error. Paper [1],
however, was mostly focused on precision and implemented a
cascade procedure when training is started with all available
derivatives, and then higher order terms are turned off as
the process continues. A contradiction between the necessity
to increase network size and decrease in minimum number
of points that prevent overfitting leads to this study. Since
absolute precision is not the goal, the process of altering cost
functions is omitted and instead the one with maximum number of available derivatives is used. Overfitting properties of
extended training are investigated for both cases - supervised
learning and solving PDE.
III. M ETHOD
A. Direct function approximation
Neural network is trained to represent a scalar function
f : Rn → R1 defined inside a region Γ ∈ Rn . Components
of the input vector ~x are denoted by xi . Network’s output is
denoted by N = N (~x) = N (x1 , ..., xn ). To gather a cost
function the following terms are used:
2
X  ∂k
∂k
N − kf
(2)
ek =
∂xki
∂xi
i
Here i can run through all or some of the input variables. If
k = 0, all terms are the same and the sum can be omitted
turning the expression into classical backpropagation cost:
(N − f )2 . The next order term k = 1 is used in double
backpropagation. The cost of order s for one pattern (i.e. local)
is defined as follows:
s
X
elocal
=
ek
(3)
s
k=0

A network is trained on a set of M patterns, so the total cost
is:
M
1 X local
e
(~xα )
Es =
M α=1 s
Results for s = 0 will be compared with the ones for s = 4.

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

B. Solving differential equation
Boundary value problem inside 2-dimensional region Γ is
considered:
U (x1 , x2 , u, ux1 , ux1 x1 , ...) = 0
u|∂Γ = f
According to (1) the function is substituted, and the equation
is written for v:
V (x1 , x2 , v, vx1 , vx1 x1 ...) = 0
v|∂Γ < ∞
which is to be represented by a neural network v = N .
Boundary condition is satisfied automatically. Similarly to (2)
terms of different order are written as:
2
2  k
X
∂
V
ek =
∂xki
i=1
And one can omit the sum for k = 0. Classical method of
solving PDE involves only e0 = V 2 . The local cost function
of the order s is similar:
s
X
ek
(4)
elocal
=
s
k=0

Used grids are very close to regular, therefore, the total cost
on Γ can be discretized
function which is a measure of elocal
s
on M points as:
Es =

µ (elocal
)
s
~
x∈Γ

M
1 X local
'
e
(~xα )
M α=1 s

Results for s = 0 will be compared with the ones for s = 4.
IV. R ESULTS
Analytical functions are approximated inside 2D box and
5D unit sphere. Boundary value problem is solved inside 2D
unit circle. Regions are denoted as Γ and their boundaries
as ∂Γ. For all cases the cost is minimized on grids that are
generated in similar manner and comprised of two parts internal and surface. The internal part is a regular grid, which
is generated in three steps. At first, a Cartesian grid with
spacing λ inside sufficient volume is created. Two random
vectors in input space are then chosen, and the grid is rotated
from one to another. Finally, it is shifted along each direction
by a random value from interval [− λ4 , λ4 ], and after that
all points outside of Γ are thrown away. The surface part
contains equidistant points from ∂Γ with distance τ = λ unless
otherwise stated.
Training is based on RProp [16] procedure with parameters
η+ = 1.2, η− = 0.5. Weights are forced to stay in [−20, 20]
interval, no min/max bonds for steps are imposed. Initial step
∆0 is set to 2·10−4 . Due to the absence of a minimal step a lot
of them are reduced to zero after certain number of epochs. To
tackle that, after each 1000 epochs all zero steps are restored
to 10−6 . The precision is single. Weights
are initialized [5]
√
with random values from range ±2/ s, where s is a number
of senders. Thresholds are initialized in range ±0.1. All of

3

the layers are fully connected and have non linear sigmoid
activation function:
1
σ(x) =
1 + exp(−x)
except for the input and output layers which are linear. To
make sure that backpropagation is run long enough to produce
overfitting, weak rules for stopping are implemented. During
both classical and extended
p training modes, the root mean
he0 itrain is tracked and after its
of e0 for grid points
best value has changed by less than 10% for the last 1000
epochs, or when the 10000th epoch is reached, the training is
stopped. The
p weights from the best iteration are saved, and
for them he0 itest is calculated on a much finer grid. All
results are linked not to the number of grid points, but to
the number of parameters N that were used for training. The
value N coincides with the number of grid points M , when
only values of function are used in the cost, and is equal to
(p + 1)M , if p extra derivatives in each point are trained.
Plotting results against N also equalizes the total number of
arithmetic operations per one epoch provided networks are
of the same architecture. For each training grid, method and
network the process is repeated 5 times, and the average
values are presented. Space above and below each plot is filled
according to the maximum and minimum values obtained in
those 5 runs.
A. Direct function approximation
1) 2D function: networks are trained to represent the following expression inside [−1, 1]2 box:


x2 
1
−x21 − x22 + 1 +2 tanh x1 · sin
−sin x1 ·cos x2
f=
2
2
A series of grids with various densities is generated. It starts
with λ = 0.073 and a grid containing 804 points. For each
next grid λ is increased, so that the total number of points
would be approximately 10% less. The last grid has λ = 1.45
and 5 points, 4 on the boundary and 1 near the middle. In each
point derivatives from the 1st to the 4th order are calculated
with respect to x1 and x2 . Extended training is using i = 1, 2
in expression (2) and s = 4 in (3), therefore, the total number
of parameters for a grid of M points is N = 9M . For classical
training s = 0 and N = M . To plot curves on the same
interval of N extended training is run only on grids from 88
to 5 points (λ from 0.24 to 1.45), and classical one is run on
grids from 804 to 43 points (λ from 0.073 to 0.37). It is worth
mentioning that calculating the target function 9M times is
exactly enough to evaluate derivatives up to the 4th order with
respect to x1 and x2 in M points using the first order finite
difference stencil. To investigate how network’s architecture
affects overfitting three different layer configurations are used:
2, 64, 64, 64, 64, 64, 64, 1
2, 512, 512, 512, 512, 512, 512, 1
2, 1024, 1024, 1024, 1024, 1024, 1024, 1
They comprise 2 · 104 , 1 · 106 and 5 · 106 weights respectively
and are referred to on figures by those numbers. After training

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

IV: 5⨯106

2⨯104

1⨯106

1.2

5⨯106

4

3

IV: 2⨯104

0.6
0.4
0.2

log (

2

0.8

〈e0〉test /

〈e0〉train )

IV: 1⨯106

log (

〈e0〉train )

IV: 5⨯106

1.0

5

〈e0〉test /

4

0.0
1
-0.2

0

200

400

600

800

N

0
400

600

800

N
Fig. 1. Approximating 2D function. Base 10 logarithm of the ratio between
root mean square errors on test and train sets is plotted against the number
of parameters used for training. Solid line is for the network with 5 · 106
weights trained by extended algorithm, and the rest curves are for networks
trained by classical algorithm.

is finished, the performance of each network is tested on a grid
with 3400 points
p
p and λ = 0.035. Base 10 logarithm of the ratio
he0 itest / he0 itrain is plotted on fig. 1. The network with
the most number of weights trained by extended algorithm
is compared with three networks trained with classical cost
function. The ratio itself and its variance, which is related
to amount of filling above and below each curve, are much
smaller for extended method. Fig. 2 shows overfitting ratio for
networks trained by extended algorithm. Plots for different
architectures are very close to each other, unlike when the
same networks are trained with classical cost function. From
about 200 parameters and more, which corresponds to 23 grid
points, all curves lie mostly below the zero, which means that
the root mean square error calculated on the dense test grid
is a bit smaller than the one calculated on the training grid.
Fig. 3 depicts the final root mean square error on the test
set. It represents an inevitably more important aspect in a
sense that the previous plots alone do not necessary mean that
training was successful. Only one curve is shown for extended
algorithm, since others lie very close to it.
2) 5D function: networks are trained to represent the following expression inside 5D unit sphere Γ : r ≤ 1

1
x4 
− sin x1 cos x2
f = (−x22 − x25 + 1) + 2 tanh x3 · sin
2
2
A series of grids is generated. Unlike 2D case, τ = 1.6λ.
The first grid has λ = 0.336 and 1579 points. Each next grid
has about 10% less points. The last grid has λ = 1.1 and
11 points, 10 of them are on the boundary and 1 near the
middle. In each point derivatives from the 1st to the 4th order
are calculated with respect to two randomly chosen directions

Fig. 2. Approximating 2D function. Base 10 logarithm of the ratio between
root mean square errors on test and train sets is plotted against the number
of grid parameters. Comparison between networks with different capacities
trained by extended algorithm, solid line is the same as on fig. 1.

0

-1

〈e0〉test

200

log

0

IV: 5⨯106

2⨯104

1⨯106

5⨯106

-2

-3

-4

-5
0

200

400

600

800

N
Fig. 3. Approximating 2D function. Base 10 logarithm of root mean square
error on the test set. The solid line is for extended algorithm, which utilizes
derivatives up to the IV order.

xp and xq out of five possible. Extended training is using
i = p, q in expression (2) and s = 4 in (3), therefore the total
number of parameters for a grid of M points is N = 9M .
Extended training is run on grids from 166 to 11 points (λ
from 0.55 to 1.1) and classical training on grids from 1579 to
98 points (λ from 0.336 to 0.62). Calculating target function
9M times is again exactly enough to evaluate derivatives with
respect to two mentioned variables. Test grid has 11000 points
and λ = 0.15. Network configurations are the same as for 2D
case, except now input layers have 5 neurons instead of 2.

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

IV: 5⨯106

2⨯104

1⨯106

5

3.0

5⨯106

IV: 5⨯106
2.5

〈e0〉test /

1.5
1.0

〈e0〉test /

4

2.0

3

log (

〈e0〉train )

5

IV: 1⨯106
IV: 2⨯104

log (

〈e0〉train )

6

2

0.5

1

0.0
0

500

1000

1500

N

0
1000

1500

N
Fig. 4. Approximating 5D function. Base 10 logarithm of the ratio between
root mean square errors on test and train sets is plotted against the total number
of parameters used for training. Solid line is for extended cost function, and
the rest are for classical.

Similarly to previous case, fig. 4 shows overfitting ratio
for the biggest network trained with extended mode, and all
networks trained with classical mode. Unlike 2D case, even
on the most dense grids overfitting is very strong for larger
networks trained without derivatives. Fig. 5 plots the same
ratio for three architectures trained with extended cost. It is
again small, barely different and has low variance between
training attempts. Finally, fig. 6 compares performance of
networks on the test set. It shows that extended mode not
only produces very close test and train errors, but also a
better approximation. Different curves for extended training
are not shown since they lie very close to each other. One
can notice that among results of classical training, the 1 · 104
network (the line marked with triangles) shows both: lower
overfitting on fig. 4 and better test cost on fig. 6. This could
make the network a preferable choice for the task. However,
when derivatives are used, any architecture demonstrates much
smaller overfitting ratio and nearly one extra order in final
precision. The only advantage of smaller network left is
training time. Despite this being a 5D task, for each input
vector it is enough to train derivatives with respect to two
directions, provided they are randomly chosen from point to
point.

Poisson equation inside 2D unit circle Γ : x21 + x22 ≤ 1 with
vanishing boundary condition is considered:

u|∂Γ = 0

0.0
-0.5

IV: 5⨯106

2⨯104

1⨯106

5⨯106

-1.0
-1.5
-2.0
-2.5
-3.0
0

500

1000

1500

N
Fig. 6. Approximating 5D function. Base 10 logarithm of root mean square
error on the test set. Solid line is for extended algorithm, which utilizes
derivatives up to the IV order.

According to (1) the following substitution can be made:
u = v(x1 , x2 ) · (1 − x21 − x22 )
where v is to be found using a neural network. To make results
comparable to previous cases, analytical solution is chosen as:

B. Solving differential equation

4u = g

Fig. 5. Approximating 5D function. Base 10 logarithm of the ratio between
root mean square errors on test and train sets is plotted against the total number
of grid parameters. Comparison between networks with different capacities
trained with extended algorithm, solid line is the same as on fig. 4.

〈e0〉test

500

log

0

(5)

ua = f (x1 , x2 ) · (1 − x21 − x22 )
where f (x1 , x2 ) is the expression used for 2D approximation.
The function ua is substituted into (5) to calculate the source
g that would produce it. This way a neural network in the

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

where i = 1, 2. Due to overlapping this produces only 24
different derivatives of v. Therefore, an equivalent classical
grid for M points in this mode is N = 5M . Extended training
is run on grids from 139 to 3 points (λ from 0.16 to 1.62), and
classical one is run on grids from 704 to 16 points (λ from
0.07 to 0.6). Test grid has 3000 points and λ = 0.033. Neural
networks have the same configurations as for 2D function
approximation. Fig. 7 compares overfitting ratio between one
network trained by extended method and all three processed
by classical one. Fig. 8 compares performance of the same
networks on the test set. Results are similar to those on fig. 1
and 3. Overfitting ratios for different architectures trained with
extended cost are nearly indistinguishable, so instead, fig. 9
depicts the maximum absolute difference between analytical
and numerical solution averaged across 5 solving attempts.
Extended mode is about one order ahead of classical method
on grids that are not too dense.
Particular details for boundary value problem are worth
mentioning. As grid spacing λ is increased, the last ”accurate”
solution (max |u − ua | < 3 · 10−5 ) is obtained by classical
method on a grid with λ = 0.09 (441 points) and by extended
one with λ = 0.29 (52 points). The last ”meaningful” solution
(max |u − ua | < 1 · 10−3 ) is obtained by classical training on
a grid with λ = 0.17 (139 points) and by extended one with
λ = 0.68 on 12 points. Those results are very close to what
one can expect from the 4th and 6th order finite difference
stencils. Solution u is based on the equation for v, in which
terms 4xi ∂v/∂xi , i = 1, 2 would be the biggest source of
stencil error ω, and for 2D Poisson equation, the value of ω is
very close to the maximum deviation from the exact solution.
The leading term of error of the first order derivative calculated
with the 4th order central stencil is λ4 /30 · ∂ 5 v/∂x5i and with
the 6th order stencil is λ6 /140 · ∂ 7 v/∂x7i . After substituting
those expressions into the error source and noting that all
derivatives of v are of the order of 1, the 4th order stencil
has ω = 1 · 10−5 for λ = 0.09 and ω = 1 · 10−4 for λ = 0.17.
The 6th order stencil has ω = 1.6 · 10−5 for λ = 0.29 and
ω = 3 · 10−3 for λ = 0.68. However, the 6th order stencil
requires at least 7 points in each direction, therefore, it can not
be used on a grid with 12 points in total. On the 52-point grid
its implementation is possible, but unlikely that precise, since
in a lot of points derivatives will have to be calculated using

1⨯106

5⨯106

〈e0〉test /

〈e0〉train )

6

log (

∂ ∂2 ∂3 ∂4
,
,
,
∂xi ∂x2i ∂x3i ∂x4i

2⨯104

7

5
4
3
2
1
0
0

100

200

300

400

500

600

700

N
Fig. 7. 2D Boundary value problem. Base 10 logarithm of the ratio between
root mean square of residuals on test and train grids is plotted against the
equivalent number of grid points. Solid line is for extended cost function and
the rest are for classical.

1

IV: 5⨯106

2⨯104

1⨯106

5⨯106

0

〈e0〉test

Extended training casts 8 extra operators on V :

IV: 5⨯106

-1
-2

log

process of minimizing the equation’s residual would have to
fit exactly the same function as in 2D approximation case. The
distinction from direct approximation is that instead of being
a squared difference between output and target, e0 is a square
of the equation’s residual. A series of grids is generated with
the first one having λ = 0.07 and 704 points and the last one
having λ = 1.62 and 3 points. Classical mode is using s = 0
in 4, therefore, in addition to values themselves, 4 derivatives
of v encountered in the residual V must be calculated in each
point:
∂
∂2 ∂2
∂
,
, 2, 2
∂x1 ∂x2 ∂x1 ∂x2

6

-3
-4
-5
0

100

200

300

400

500

600

700

N
Fig. 8. 2D Boundary value problem. Base 10 logarithm of root mean square
residual on the test set plotted against the equivalent number of grid points.
Solid line is for extended algorithm, which utilizes derivatives up to the IV
order.

not central but forward or backward stencil, approximation
error ω for which is 20 times higher.
C. Eliminating few trivial explanations
The simplest explanation of lower overfitting for extended
training could be an early stopping. To summarize data the
lower 3-quantiles for the number of epochs were calculated

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

0
IV,

log (max |u-ua|)

to the δ being 20% less. Usage of additional derivatives does
not restrict training to a smaller subset of weights.
Another simple explanation of overfitting ratio can be
obtained as follows: consider 1-dimensional function f approximated by a network N trained on n  1 points uniformly
distributed on [0, 1]. Error e on the test set can be approximated
as a difference between network’s value N and target f in the
middle of an interval between two neighboring training points
a and b. Its length is 1/n and from Taylor expansion in a one
can write:


0 1
1 ∼ [N − f ]
N − f + (N − f )
ea+ 2n
1 '
a+ 2n
2n a

1⨯106

1⨯106

-1

7

-2

-3

-4

-5

After dividing this expression by N − f one obtains the
overfitting ratio between e on train and test sets:

-6

(N − f )0 1
etest
∼1+
etrain
N − f 2n

0

100

200

300

400

500

600

700

N
Fig. 9. 2D Boundary value problem. Base 10 logarithm of the maximum
deviation from analytical solution calculated on the test set, N is an equivalent
number of grid points. Solid line is for extended algorithm, which utilizes
derivatives of the residual up to the IV order.

for each problem, training mode and network size. For 2D
classical approximation they are 7000, 6000 and 4000 and for
extended one 5000, 4000 and 4000 for 2 · 104 , 1 · 106 and
5 · 106 networks respectively. For 5D case they are 6000, 4000
and 3000 for classical and 6000 for all networks on extended
mode. Finally the boundary value problem has got 5000, 3000
and 2000 for classical and 4000 for all networks on extended
algorithm. The difference is not quite significant. Even if
extended training was run as long as classical whenever
possible, it would not bring any noticeable changes to the
plots.
Another concern arises from a specific property of RProp.
If a landscape of a cost function has a lot of oscillations,
then steps for many weights will be reduced to zero quite
quickly, and the procedure will not go along those directions.
If inclusion of high order derivatives creates such ”ripples”,
training could effectively be constrained to a network of much
smaller capacity, therefore, it can exhibit lower overfitting.
To investigate this, the following measurement was made: the
amount of non zero steps δ was observed for 2 · 104 network
during classical and the 4th order training for 5D function
approximation with N = 449. In both
was
p cases the procedure
−4
run
for
6000
epochs.
It
resulted
in
he
i
=
1
·
10
and
0
p train
p
−2
−3
for
classical
and
,
he
i
=
3
·
10
he
i
=
2
·
10
0
0
test
train
p
−3
he0 itest = 4 · 10
for extended training. In both cases
δ ∼ 300, except for small intervals of epochs close to the
beginning or resetting of steps. The ratio δextended /δclassical
was calculated with the following results: minimum 0.45, the
lower 3-quantile 0.8, and median 0.82. Value δ is relatively
small, so another measurement was made: an average probability of a weight to be changed within 10 epochs is estimated
as 11% for classical and 14% for extended training. It does
not seem possible that overfitting ratios of 300 and 2 are due

(6)

If one chooses to split n points into two parts, distribute half of
them uniformly and then use the rest to calculate derivatives
of target in each point using the first order finite difference
stencil, the error in between points and the ratio will be written
as:


1
ea+ n1 ' N − f + (N − f )0
n a
etest
(N − f )0 1
∼1+
etrain
N −f n

(7)

The term (N − f ) is being minimized by classical algorithm
and its derivative (N − f )0 is what is included in the first
order extended training. From numerical experiments it was
found that, if a neural network trained by classical algorithm
shows the train set precision ε it has the first derivative on
this set off by about 10ε (this factor also depends on the
network size). For extended mode, the precisions for different
trained derivatives are approximately the same , and the first
non trained derivative is similarly 10 off. Therefore, one can
estimate test to train ratio 6:
etest
10ε 1
5
∼1+
=1+
etrain
ε 2n
n
And the ratio 7:
etest
1
1
∼1+
=1+
etrain
n
n
For rather small n used in this paper, the difference between
logarithms of those expressions is about 0.2. The same analysis
expanded by higher order terms and applied to 2D and 5D
cases can more or less explain differences between curves for
classical and extended training shown on figure 10, provided
ratios (N − f ){k} /(N − f ) are obtained from experiment.
Thus, overfitting ratio can be seen as a consequence of
two reasons: extended training being simply a higher order
approximation and existence of synergy between derivatives,
which happens when inclusion of higher order terms enhances
precision of the low order ones. However, this description is
by no means complete. Note the solid line on fig. 7: from
the initial number of 139 grid points down to 11 logarithm of
the overfitting ratio for the network with 5 · 106 weights does
not exceed 0.03. Meanwhile, logarithm of the test precision

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

IV

III

II

I



6

〈e0〉test /

4
3

log (

〈e0〉train )

5

2
1
0
0

500

1000

1500

N
Fig. 10. Approximation of 5D function by the network with 106 weights.
Overfitting ratio lowers as the order of included derivatives increases. Horizontal axis is the number of parameters used for training. Curves are marked
according to the maximum order of trained derivatives.

p
log he0 itest represented by the solid line on fig. 8 increases
from −4.5 to −2.0. It means that a network with few millions
of connections have not established them in a way that would
satisfy 9 conditions in 11 training points more precisely than
anywhere else in Γ. Similar situation can be spotted in 2D and
5D cases of function approximation.
V. C ONCLUSION
Unexpectedly low overfitting was observed when derivatives
of target up to the 4th order were included in backpropagation.
Cost function values on train and test sets were very similar
whether neural network approximated a smooth analytical
expression or solved partial differential equation. Increasing
the capacity of a network from thousands up to millions of
weights barely affected those results. Test and train precisions
were less than 8% different for a network with 5 · 106 weights
up until 23 points left for 2D function approximation and
11 points left for solving 2D Poisson equation. Whereas
approximating smooth analytical functions can be regarded as
highly artificial task, solving partial differential equations on
a very sparse grids can have more practical applications.
ACKNOWLEDGMENT
Author expresses a profound gratitude to his scientific advisor E.A. Dorotheyev and gratefully indebeted to Y.N. Sviridenko, A.M. Gaifullin, I.A. Avrutskaya and I.V. Avrutskiy
without whom this work would not be possible.
R EFERENCES
[1] VI Avrutskiy. Enhancing approximation abilities of neural networks by
training derivatives. arXiv preprint arXiv:1712.04473, 2017.

8

[2] Harris Drucker and Yann Le Cun. Improving generalization performance
using double backpropagation. IEEE Transactions on Neural Networks,
3(6):991–997, 1992.
[3] Gary William Flake and Barak A Pearlmutter. Differentiating functions
of the jacobian with respect to the weights. In Advances in Neural
Information Processing Systems, pages 435–441, 2000.
[4] Babak Hassibi, David G Stork, and Gregory J Wolff. Optimal brain
surgeon and general network pruning. In Neural Networks, 1993., IEEE
International Conference on, pages 293–299. IEEE, 1993.
[5] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Delving
deep into rectifiers: Surpassing human-level performance on imagenet
classification. In Proceedings of the IEEE international conference on
computer vision, pages 1026–1034, 2015.
[6] Geoffrey E Hinton, Nitish Srivastava, Alex Krizhevsky, Ilya Sutskever,
and Ruslan R Salakhutdinov. Improving neural networks by preventing
co-adaptation of feature detectors. arXiv preprint arXiv:1207.0580,
2012.
[7] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. Imagenet
classification with deep convolutional neural networks. In Advances
in neural information processing systems, pages 1097–1105, 2012.
[8] Manoj Kumar and Neha Yadav. Multilayer perceptrons and radial
basis function neural network methods for the solution of differential
equations: a survey. Computers & Mathematics with Applications,
62(10):3796–3811, 2011.
[9] IE Lagaris, A Likas, and DI Fotiadis. Artificial neural network methods
in quantum mechanics. Computer Physics Communications, 104(1-3):1–
14, 1997.
[10] Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis. Artificial
neural networks for solving ordinary and partial differential equations.
IEEE Transactions on Neural Networks, 9(5):987–1000, 1998.
[11] Steve Lawrence and C Lee Giles. Overfitting and neural networks:
conjugate gradient and backpropagation. In Neural Networks, 2000.
IJCNN 2000, Proceedings of the IEEE-INNS-ENNS International Joint
Conference on, volume 1, pages 114–119. IEEE, 2000.
[12] Steve Lawrence, C Lee Giles, and Ah Chung Tsoi. Lessons in
neural network training: Overfitting may be harder than expected. In
AAAI/IAAI, pages 540–545. Citeseer, 1997.
[13] Yann LeCun, John S Denker, and Sara A Solla. Optimal brain damage.
In Advances in neural information processing systems, pages 598–605,
1990.
[14] A Malek and R Shekari Beidokhti. Numerical solution for high
order differential equations using a hybrid neural network optimization
method. Applied Mathematics and Computation, 183(1):260–271, 2006.
[15] Russell Reed. Pruning algorithms-a survey. IEEE transactions on Neural
Networks, 4(5):740–747, 1993.
[16] Martin Riedmiller and Heinrich Braun. A direct adaptive method
for faster backpropagation learning: The rprop algorithm. In Neural
Networks, 1993., IEEE International Conference on, pages 586–591.
IEEE, 1993.
[17] Warren S Sarle. Stopped training and other remedies for overfitting.
Computing science and statistics, pages 352–360, 1996.
[18] Yazdan Shirvany, Mohsen Hayati, and Rostam Moradian. Multilayer
perceptron neural networks with novel unsupervised training method
for numerical solution of the partial differential equations. Applied Soft
Computing, 9(1):20–29, 2009.
[19] Patrice Simard, Yann LeCun, John Denker, and Bernard Victorri. Transformation invariance in pattern recognitiontangent distance and tangent
propagation. Neural networks: tricks of the trade, pages 549–550, 1998.
[20] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever,
and Ruslan Salakhutdinov. Dropout: A simple way to prevent neural
networks from overfitting. The Journal of Machine Learning Research,
15(1):1929–1958, 2014.
[21] Igor V Tetko, David J Livingstone, and Alexander I Luik. Neural
network studies. 1. comparison of overfitting and overtraining. Journal
of chemical information and computer sciences, 35(5):826–833, 1995.
[22] Wojciech Zaremba, Ilya Sutskever, and Oriol Vinyals. Recurrent neural
network regularization. arXiv preprint arXiv:1409.2329, 2014.
[23] CJ Zúñiga-Aguilar, HM Romero-Ugalde, JF Gómez-Aguilar,
RF Escobar-Jiménez, and M Valtierra-Rodrı́guez. Solving fractional
differential equations of variable-order involving operators with mittagleffler kernel using artificial neural networks. Chaos, Solitons &
Fractals, 103:382–403, 2017.

