SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

1

Enhancing approximation abilities of neural
networks by training derivatives

arXiv:1712.04473v1 [cs.NE] 12 Dec 2017

V.I. Avrutskiy

Abstract—Method for increasing precision of feedforward
networks is presented. With the aid of it they can serve as a
better tool for describing smooth functions. Namely, it is shown
that when training uses derivatives of target function up to the
fourth order, approximation can be nearly machine precise. It is
demonstrated in a number of cases: 2D function approximation,
training autoencoder to compress 3D spiral into 1D, and solving
2D boundary value problem for Poisson equation with nonlinear
source. In the first case cost function in addition to squared
difference between output and target contains squared differences
between their derivatives with respect to input variables. Training
autoencoder is similar, but differentiation is done with respect
to parameter that generates the spiral. Supplied with derivatives
up to the fourth the method is found to be 30-200 times more
accurate than regular training provided networks are of sufficient
size and depth. Solving PDE is more practical since higher
derivatives are not calculated beforehand, but information about
them is extracted from the equation itself. Classical approach is to
put perceptron in place of unknown function, choose the cost as
squared residual and to minimize it with respect to weights. This
would ensure that equation holds within some margin of error.
Additional terms used in cost function are squared derivatives of
the residual with respect to independent variables. Supplied with
terms up to the second order the method is found to be 5 times
more accurate. Efficient GPU version of algorithm is proposed.
Index Terms—Neural networks, high order derivatives, partial
differential equations, function approximation.

I. I NTRODUCTION

N

EURAL networks can be used as universal approximators [9], [11], [21], [22], [25] for both high and low
dimensional tasks. In high dimensional cases they successfully overcome the curse of dimensionality [3], [4], thus, are
excellent remedy for problems like voice recognition [12],
[17] and pattern classification [7], [23], [38]. Low dimensional
applications are not that famous since many alternatives are
available. For example, when solving differential equations
[24] (usually 1D [28], [29] or 2D [27], sometimes 3D [26])
networks inevitably have to compete with other methods like
finite differences where function is described by its values on
a set of points and in each point, in theory, those methods
can be accurate within machine precision. This quality of
approximation is not readily achievable for neural networks.
Widely developed techniques for training them [10], [16], [30],
[36], [37] are mostly focused on problems like classification,
which are not very sensitive to actual output since small
disturbances can hardly produce a shift in class. For the
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

case of direct function approximation any deviation of output
decreases accuracy.
To make neural networks more efficient tool for approximation of low dimensional smooth functions and to push their
precision to machine level this paper proposes an extended
version of backpropagation. The procedure is based on more
general cost function that includes along with output values
derivatives of those with respect to certain scalars, usually
input components. Essentially, it is a generalization of algorithms used and described in [13], [14], [20], [35], but for
any order of derivatives. Extended versions of backpropagation
described in said papers contain only low order terms and
benefits of including those are very modest.
In real world applications supplying training procedure
with high order derivatives can be a serious challenge. From
this perspective, method of solving differential equations can
benefit much more from the proposed technique, since all extra
terms of cost function can be calculated by trivial differentiation of the equation itself. Currently no papers on solving
PDE with neural networks contain a universal algorithm for
calculating gradient of a cost with respect to weights that could
be used with or without mentioned enhancements. Although,
in mathematical sense presented algorithm is an equivalent
to automatic differentiation procedure [18], it is written in
an optimal form for GPUs, thus can be very useful for
computationally demanding tasks such as solving equations
on very large grids. Even for classical training where standard
neural network libraries should be quite efficient, the proposed
code is about 3 times faster for networks used in this paper.
For the case where many derivatives are to be calculated, naive
implementations of automatic differentiation would probably
be much more slower.
II. A LGORITHM
This section is focused on implementing a gradient algorithm that can be used to minimize a more general form of
cost function for a neural network. The named function now
depends not only on network’s output, but also on derivatives
of output with respect to input. As previously stated this is
useful in at least two different cases: if some derivatives of
the target function are known, one can minimize a sum of
multiple cost terms - one for its values and one for each known
derivative. Another scenario is solving differential equations.
Here, neural network’s inputs are considered as independent
variables and output as a function, for which one writes down
a differential equation with all terms moved to either side
of it, and then square of that expression is used as a cost

© 2017 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

function. Neural network that satisfies the equation along with
boundary conditions with necessary precision can be viewed
as a numerical solution. Described further is an algorithm
for calculating gradient of such generalized cost function for
feedforward network with any number of hidden layers and
derivatives of any order. It is presented without derivation,
that can be found in [2].
A. Notation
Quantity that each neuron obtains from the previous layer
before applying its nonlinear sigmoid mapping will be referred
to as neuron activity. Values of neuron activities for layers are
gathered in matrices with the letter z like z κα . Index α runs
through a number of input patterns and κ through neurons
of that layer. The absolute values of indices like |α| and |κ|
denote the total number of enumerated objects. Here they are
referring to the number of input patterns and the number of
neurons in layer respectively. Different layers are denoted by
different Greek letters which also appear in weights matrices
connecting those layers: W θκ would be a matrix that is used
to pass from the layer with neurons denoted by κ to the layer
with those denoted by θ. In addition to neurons activities z κα ,
each layer has derivatives of those activities with respect to
certain variables denoted by Latin letters a, b, c, d . . . and those
derivatives are matrices of the same size |κ| × |α|. They are
denoted as follows:
∂2
∂ κα
κα
z = zaκα , 2 z κα = zaa
∂a
∂a
and so on. Initializing them at the input layer and propagating
forward is required to calculate cost function. The last layer
neurons are considered to be linear and designated by ω. The
cost function includes derivatives of z ωα and can be written
as
X
ωα
E=
e(z ωα , zaωα , zaa
, . . .)
α

where e is a scalar function defined on space of output
vector and all its derivatives. The letter for the input layer
neurons is β. In formulas where Latin variables like a, b
are used explicitly their number is limited to four however
generalization is straightforward.
Before writing forward and backward propagation formulas
some more notations are to be made. Let’s fix the order of
variables, derivatives with respect to which are propagated,
for example: {a, b, c, d}. Any derivative now can be written
as a vector s of integers, where each component is the
order of derivative with respect to the corresponding variable:
s = {sa , sb , sc , sd }. Difference two such vectors is defined
element-wisely: s − r = {si − ri }.PThe absolute value is
defined as sum of components |s| = si
i

B. Forward pass
βα
The first step is to initialize z βα , zaβα , zaa
, zbβα , . . . at the
βα
input layer of the neural network. z is a conventional input
matrix, of which one have to take derivatives with respect
to a, b and so on. Let’s say a is the first component of the
input vector, then zaβα is a matrix with the first row elements

2

equal to 1 and all others equal to 0: zaβα = δβ−1 (Kronecker
βα
βα
delta). All higher derivatives that include a like zab
, zaa
are zero. Therefore, initializing derivatives with respect to the
input of neural network is trivial. If b is some variable that
all components of the input vector depend on, for example,
input matrix is z βα = f (α, β) sin b, then zbβα = f (α, β) cos b,
and other higher derivatives will be non zero and should be
calculated accordingly.
After input layer matrix is generated and its derivatives are
calculated, they are to be propagated forward. Consider two
generic successive layers with indices κ and θ. The expression
to propagate from κ to θ is as follows:
zsθα = tθ + W θκ ×

∂ |s|
σ(z κα )
∂s

(1)

Here tθ is a threshold vector added to each column of matrix
product to the right of it. Scalar sigmoid function σ is applied
independently to each element of z κα . Values themselves
(s = ~0, ∂ |s| /∂s ≡ 1) are of course propagated by simple
matrix multiplication and nonlinear sigmoid:
z θα = tθ + W θκ × σ(z κα )
The first order derivative propagation formula looks the same
for all variables. Let’s write it for a:
∂
∂z θα
≡ zaθα = W θκ ×
σ(z κα ) = W θκ × (σ 0 (z κα ) · zaκα )
∂a
∂a
Expression in parentheses is an element-wise product. In
similar formulas of forward pass sigmoid’s argument as well
as “·” sign will be omitted, so it is written as σ 0 zaκα . All
∂
derivatives like ∂a
σ(z κα ) are calculated using the chain rule:
derivative of σ must be multiplied by derivative of the inner
∂ κα
function ∂a
z = zaκα . This leads to derivative of next layer
θα
za being dependent on derivative of previous zaκα .
The second order derivative propagation:
∂ 2 z θα
∂2
θα
= W θκ ×
≡ zab
σ(z κα ) =
∂a∂b
∂a∂b
κα
= W θκ × (σ 00 zaκα zbκα + σ 0 zab
)
For non mixed second derivatives like ∂ 2 /∂a2 one can
2
κα
simplify the last part to (σ 00 (zaκα ) + σ 0 zaa
). The same
holds for further formulas - all Latin variables are considered
distinct, and in case they are not, it is useful to simplify expressions first to avoid unnecessary arithmetic and/or memory
operations. Any mixed second derivative depends on three
terms of previous layer: said second derivative plus both first
order ones. Non mixed second derivatives depend only on two
terms: first and second derivatives with respect to that variable.
The third order:
h
θα
= W θκ × σ 000 zaκα zbκα zcκα +
zabc
i
κα κα
κα κα
κα κα
κα
+ σ 00 · (zab
zc + zac
zb + zbc
za ) + σ 0 zabc

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

The forth order:
h
θα
zabcd
=W θκ × σ IV zaκα zbκα zcκα zdκα +

D. Computational efficiency

κα κα κα
κα κα κα
κα κα κα
+ σ 000 (zab
zc zd + zac
zb zd + zad
zb zc +
κα κα κα
κα κα κα
κα κα κα
+ zbc
za zd + zbd
za zc + zcd
za zb )+
κα
κα
κα
κα
+ σ 00 (zaκα zbcd
+ zbκα zacd
+ zcκα zabd
+ zdκα zabc
+
i
κα κα
κα κα
κα κα
κα
+ zab
zcd + zac
zbd + zad
zbc ) + σ 0 zabcd

In general, the derivative with vector s uses those and only
those r derivatives from the previous layer, for which s − r
has no negative components. For example, if cost function
would use only on ∂ 4 /∂a4 , one still needs to calculate ∂ 3 /∂a3 ,
∂ 2 /∂a2 , ∂/∂a and values themselves for all layers.
C. Backward pass
After obtaining the final layer matrices z ωα , zaωα , zbωα . . .
one can calculate the cost function:
X
E=
e(z ωα , zaωα , zbωα . . .)
α

and more importantly, its derivatives with respect to elements
of all used matrices:
∂E
∂E
∂E
, ωα , ωα . . .
ωα
∂z
∂za ∂zb
which are in turn matrices of the same size |ω| × |α|. They
are propagated backward using the following relation (now
one moves from θ to κ):


X
∂E
∂ |s−r| 0 κα
∂E
r
θκ T
=
C
·
σ
(z
)·
[W
]
×
(2)
s
∂zrκα
∂(s − r)
∂zsθα
s
where sum is taken through all s that were used for
forward propagation and for which s − r has no negative
components. Differential operator is applied only to σ 0 . Csr
is a combinatorial coefficient:
Csr =

Y
i

si !
ri !(si − ri )!

Propagating those matrices backward is required to calculate
gradient of the cost function with respect to weights and
thresholds:
X ∂E
∂E
∂ |s|
=
σ(z κα )]T
(3)
×
[
θα
∂Wθκ
∂z
∂s
s
s
X ∂E
∂E
=
k
∂t
∂z κα
α

3

(4)

In (3) sum is taken through all used derivatives. Expression for
gradient of E with respect to the weights Wθκ between layers
κ and θ includes matrices like like zaκα calculated during the
forward pass on layer κ - they emerge from the term
∂ |s|
σ(z κα )
∂s
as well as matrices ∂E/∂zsθα obtained on layer θ during the
backward pass.

ωα
Consider an example when cost function depends on zabcd
.
That requires one to initialize and propagate 16 following
matrices:

z κα
κα
zbc

zaκα
κα
zbd

zbκα
κα
zcd

zcκα
κα
zabc

zdκα
κα
zabd

κα
zab

κα
zac

κα
zad

κα
zacd

κα
zbcd

κα
zabcd

This can be readily done using (1) and is approximately
16 times longer than classical forward pass as one have to
compute 16 different matrix products. Applying nonlinear
sigmoids and performing element wise operations on matrices
require a small fraction of the total time.
Backpropagation is a bit more nuanced. Let’s consider layer
κα
κ and the derivative of E with respect to matrix zab
. To
κα
obtain ∂E/∂zab one should use the matrices of layer θ. As
previously described, s − r for r = {1, 1, 0, 0} can be non
negative for:
s = {1, 1, 0, 0}, {1, 1, 1, 0}, {1, 1, 0, 1}, {1, 1, 1, 1}
thus, according to (2), four of those matrices will be used
(sigmoid’s argument is omitted again):




∂ 0
∂E
∂E
∂E
0
θκ T
θκ T
] × θα + σ · [W ] × θα
κα = σ · [W
∂zab
∂c
∂zab
∂zabc




2
∂ 0
∂E
∂
∂E
θκ T
0
θκ T
+ σ · [W ] × θα +
σ · [W ] × θα
∂d
∂c∂d
∂zabd
∂zabcd
2

∂
σ 0 are calculated using the chain rule just the
Parts like ∂c∂d
same way as in forward propagation:

∂ 2 0 κα
κα
σ (z ) = σ 000 zcκα zdκα + σ 00 zcd
∂c∂d
Note that now instead of σ, formulas contain its derivative
κα
spawned here by the chain rule
σ 0 . All terms zcκα , zdκα , zcd
were already obtained during the forward pass. Expression
κα
∂E/∂zab
that one has to calculate due to the presence of
derivative ∂/∂a∂b contains four matrix products, unlike the
case of forward pass when any extra derivative was propagated
by only one product, therefore, higher computational demand
might be expected, especially for low order derivatives like
∂E/∂z κα that have 16 matrix products. This is not quite true
because whenever the matrix from the next layer ∂E/∂zsθα
is encountered, it always multiplied from the left by the
same term [W θκ ]T , thus, matrix product is always the same,
and only subsequent element-wise operations are different.
Therefore, computationally expensive matrix multiplications
can be done once for every backpropagated derivative.
A trick to reduce memory complexity is worth
mentioning. Instead of calculating pure matrix products
[W θκ ]T ×∂E/∂zsθα and storing them in some cache for
further use one can put them in memory dedicated for
∂E/∂zrκα with r = s and then calculate ∂E/∂zrκα in a
proper order. Namely, as soon as evaluation of ∂E/∂zrκα is
started for any r, one have to multiply pure matrix residing
in that memory by σ 0 (z κα ), thus, it cannot be used anywhere
else. To tackle this situation one should start with r = ~0
and, thus, spoil s = ~0, which is not used anywhere else
since s − r have to be positive. Then one can pick any

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

first order derivative, for example, r = {1, 0, 0, 0} and spoil
s = {1, 0, 0, 0}, which is again not a problem, since the only
other case when {1, 0, 0, 0} − r has no negative components
is when r = ~0, but that term was already calculated. The
same way all first order derivatives can be calculated, then
all second order ones and so on. When following this pattern
no conflicts are encountered - this can be easily seen when
one writes down all formulas for backward pass.
E. Hardware efficiency.
Some technical details are required to get maximum hardware efficiency. Note, that this study uses CUDA with cuBLAS
[31], but other platforms probably behave the same way. For
all matrix multiplications cublasSgemm routine is used, and
the rest of the code is written in such a way that delay between
calls of cublasSgemm is minimal, and matrices that it operates
on are as big as possible. Since heights of matrices are fixed by
network’s configuration, they can only be as long as possible,
i.e. it is preferably to process all input patterns in one batch.
As for the delay between cublasSgemm calls, two options are
available: either memory for matrices of all derivatives in all
layers is permanently allocated in the GPU, or intermediate
data saves/loads occur in background while information is
propagated from one layer to another. Permanent allocation
is, of course, suitable for systems with enough GPU memory.
Intermediate loading and saving can be used in any system
with fast enough communications between device and the host.
More precisely, the time for matrix |α| × |κ| to load must be
smaller than the time required to multiply it from the right by
matrix |κ| × |θ|. If that’s the case, one can hide all memory
operations behind computations in the way described below.
Consider a situation when K derivatives are to be calculated
and the number of input patterns is |α|. Two CUDA streams
are created: A and B. For each stream one allocates GPU
memory bank that fits 2K matrices of the size nmax × |α|,
where nmax is number of neurons in the thickest layer. Those
banks will be referred to as MA and MB . The first half of
each bank is used for storing matrices like zs and denoted
with the superscript +. The second half is for ∂E/∂zs , and
its superscript is −. Streams are able to cast forward and
backward propagation formulas with minimum delay when
source and destination are located in banks. Weights of neural
network are always stored on device. Host memory has to fit
everything that is produced by forward propagation.
During the forward pass all functions and data transfers
operate only on the first halves of memory banks. The procedure is as follows: input matrices are initialized and put into
MA+ . Stream B uses forward propagation formulas (1) on MA+
and puts results to MB+ . As soon as it finishes calculations it
starts saving MB+ to the host. At the same time A starts to
evaluate forward pass formulas using the data from MB+ . The
result overwrites bank MA+ . This produces the third layer data.
When A completes this operation and starts to save, B is done
saving and waiting for signal from A to start calculating data
for layer the 4 and so on.
Let’s say stream A reaches the last layer n. Derivatives
of cost are calculated and put into MA− . At this point MB+

4

TABLE I
RUNNING TIMES FOR CLASSICAL BACKPROPAGATION
n

Theano

TensorFlow

CUDA C

128

17.4

10

3.1

256

19.9

11.8

5.7

512

27.6

22.2

14.9

contains the data from layer n−1st obtained during the forward
pass, thus, all matrices used in formulas (3), (4) of ∂E/∂W
for the last matrix of weights and thresholds are already stored
in banks. Said formulas are evaluated by stream B. After
finishing, B uses backward pass formulas (2) on MB+ and
MA− to propagate ∂E/∂zs from nth to n − 1st layer. Results
are stored in MB− . Meanwhile stream is A loading n − 2nd
layer data from the host into MA+ and finishes before B ends
its calculations. As soon as B is done, stream A calculates
∂E/∂W for another weights matrix using MB− and MA+ ,
propagates ∂E/∂zs backward using the same data and stores
results in MA− . Meanwhile stream B is loading n − 3rd layer
data from the host into MB+ . This process is repeated until
gradient of E with respect to all weights is calculated. No
memory conflicts are encountered, since different halves of
banks are used in turns.
CUDA allows [32] simultaneous copying and execution as
long as it is done by different streams. That means saving and
loading operations are completely hidden by calculations and
do not affect total time and efficiency. Another good practice
is to use the latest available CUDA and cuBLAS libraries.
F. Performance test
CUDA C code was written using presented suggestions.
It is tested on fully connected perceptrons with 7 layers of
the same width n. Classical backpropagation is run for 1000
epochs with 2048 input patterns processed in one batch. It is
compared against Keras 2.0.8 with backends theano [6] 0.9.0
and TensorFlow [1] 1.3.0. The system is Amazon Linux 3.3
run on EC2 p2.xlarge instance with one GK210 core of Tesla
K80 available. Driver version is 375.66, CUDA version is
8.0. Results are gathered in Table I. This paper uses neural
networks with n equal to 64 and 128. Provided TensorFlow
can scale its performance for cases when many derivatives are
being propagated, the gain is around 300%.
III. R ESULTS
For all applications RProp [33] procedure is used with
parameters: η+ = 1.2, η− = 0.5. Weights are forced to stay in
[−20, 20] interval, no min/max bonds for steps are imposed.
Initial steps ∆0 are set to 2·10−4 , unless otherwise stated. The
absence of minimal step causes many of them being reduced
to zero during the training. To tackle that ∆ are restored
back to certain values where stated. Weights matrices W θκ
p are
initialized [15], [19] with random values from range ±2/ |κ|.
Thresholds are random from ±0.1 range. All layers, but the
input and output, are nonlinear with sigmoid function:
1
σ(x) =
1 + exp(−x)

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

E4 =

5

1X
α
α
α
α
c10 (zxxxx
− fxxxx
)2 + c11 (zxxxy
− fxxxy
)2 +
N α
α
α
α
α
+ c12 (zxxyy
− fxxyy
)2 + c13 (zyyyx
− fyyyx
)2 +
α
α
+ c14 (zyyyy
− fyyyy
)2

Here and further subscripts of z and f denote derivatives.
Constants ci are used to equalize the magnitudes of derivatives
of different orders, since for terms of Fourier series like
sin 10x they grow exponentially. For example:
c1 = c2 = 0.2, c3 = 0.025
c9 = 0.0025, c13 = 0.001, c14 = 0.0005
The first training phase uses:
Etot = E4 + E3 + E2 + E1 + E0

Fig. 1. Function generated by piece of random 2D Fourier series

unless otherwise stated. All input patterns are processed in one
batch.
A. 2D function approximation
Target function is generated by 2D Fourier series:
f (x, y) =

10
10 X
X
rnk
sin nx cos ky + . . .
n
·k
n=1

and lasts for 3000 epochs. Steps ∆ are reset back to 2 · 10−4
after every 1000 epochs. The second phase discards the fourth
derivatives so that:
Etot = E3 + E2 + E1 + E0
and goes on for 4000 epochs with initial step ∆0 = 10−5 .
Without setting ∆0 low the cost function was disturbed too
much on the first step to get back to previous values. The last
phase uses only the first and second derivatives:
Etot = E2 + E1 + E0

k=1

where rnk is a set of random coefficients in range [−1, 1].
Three similar terms with other combinations of sin and cos
have different constants and are omitted for brevity. The total
number of parameters is 400. The particular realization is
shown on Fig. 1. Cost function for classical training is:
N

E0 =

1X α
(z − f α )2
N α=1

where α denotes the number of input pattern. Network is a
fully connected perceptron with the following configuration
of layers:
2, 128, 128, 128, 128, 1
Extended training is divided into phases with different cost
functions. Each of them is gathered using the following terms:
1X
E1 =
c1 (zxα − fxα )2 + c2 (zyα − fyα )2
N α
1X
α
α 2
α
α 2
E2 =
c3 (zxx
− fxx
) + c4 (zxy
− fxy
) +
N α
α
α 2
+ c5 (zyy
− fyy
)
1X
α
α
α
α
E3 =
c6 (zxxx
− fxxx
)2 + c7 (zxxy
− fxxy
)2 +
N α
α
α
α
α
+ c8 (zyyx
− fyyx
)2 + c9 (zyyy
− fyyy
)2

and 3000 iterations of weights updates are made. Initial step
∆0 = 10−5 . Each iteration from the first, second and third
phases takes approximately 15, 10 and 6 more time than one
iteration of regular backpropagation, so in order to compensate
for total running time conventional training without derivatives
is run for 100,000 epochs with resetting ∆ back to 2 · 10−4
after each 10,000.
Input for the extended algorithm is generated as a Cartesian
grid in [−1, 1]2 (this range does not coincide with the period
of f ) with spacing 1/13, a total of 729 points. In each point
along with values 14 different derivatives are calculated. To
provide classical training with the same amount of information
about the target function one might want to increase number
of points 15-fold. However, setting the grid spacing lower than
1/26 (2809 points) produced no changes in training results,
so in order to save time it is not increased any further.
Net initialization and training is repeated 100 times for
both methods. Trained networks are verified against analytical
expression on grid with spacing 1/41. Plots on Fig. 2 depict
the probability F (E) of obtaining a network with root mean
square error less than E. An average error for the extended
training is 200 times smaller than for classical one. Each
trained network has maximum error approximately 6 times
higher than the root mean square one regardless of the training
method. Apart from that their F (E) plots look very similar
and are not included.

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

0.8

0.8

0.6

0.6
F(E)

1.0

F(E)

1.0

6

0.4

0.4

0.2

E4

0.2

E4

E0
0.0

E0
0.0

-5.5

-5.0

-4.5

-4.0

-3.5

-3.0

-2.5

log E
Fig. 2. Cumulative distribution function of RMS of error E for approximation
of 2D expression. E0 is the classical training and E4 is the extended variant.

B. Autoencoder for 3D curve
Target vectors are generated by a parametric curve in 3D
space:
{x(t), y(t), z(t)} = {cos t, sin t, t/π} = f ω (t)
ω = 1, 2, 3
They are to be compressed into one dimensional representation
and restored back by fully connected perceptron [5], [8] with
the following configuration of layers:
3∗ , 64, 64, 1∗ , 64, 64, 3∗
Asterisks denote linear layers. Cost function for classical
training is:
3
N
1 X X ωα
E0 =
(z − f ωα )2
N α=1 ω=1

Extended training is divided into similar phases. To construct
cost functions the following terms are used:
1 XX
E1 =
2(ztωα − ftωα )2
N α ω
E2 =

1 XX
ωα
3(ztt
− fttωα )2
N α ω

E3 =

1 XX
ωα
ωα 2
4(zttt
− fttt
)
N α ω

E4 =

1 XX
ωα
ωα 2
5(ztttt
− ftttt
)
N α ω

Total cost functions are put together it the same way as for 2D
approximation, but their training durations are different. The
first phase lasts for 5000 epochs and after each 1000 iterations
steps ∆ are reset back to 2 · 10−4 . The second phase lasts for
4000 epochs with ∆0 = 10−5 and the last one for 2000 epochs

-6

-5

-4

-3

log E
Fig. 3. Cumulative distribution function of RMS of error E for autoencoder.
E0 is the classical training and E4 is the extended variant. Three curves are
for different components of the output vector.

with the same ∆0 = 10−5 . Even though now each iteration
takes approximately 5, 4 and 3 more time respectively, the
regular training is run for 126000 steps, as if they were as
slow as in the previous example.
Input vectors for the extended algorithm are generated as
points of curve f ω (t) with values of t taken from interval
[−2π, 2π] with spacing 1/5, 63 points in total. In each
point values and 4 derivatives of them with respect to t are
calculated. Classical training produces slight overfitting with
this spacing, so it is lowered to 1/10 (126 points). Decreasing
it further does not change approximation quality.
Training is repeated 200 times. Networks are verified on
vectors generated by t ∈ [−2π, 2π] with step 1/100 (few of
them near 2π are excluded since they lie outside of the training
set interval). Plots on Fig. 3 depict the probability F (E) of
obtaining a network with root mean squared error less than E
for each component separately. Average error for the extended
algorithm is 30 times smaller than for classical training. Each
trained network has maximum error approximately 5 times
higher than the root mean square one regardless of the training
method, plots for them are not included.
C. Boundary value problem for Poisson equation
Previous use cases of the extended backpropagation require
prior calculation of derivatives of target function in order to increase precision. For real world applications those derivatives,
especially high order ones, can be very difficult to obtain.
However, in case of solving partial differential equations,
necessary relations between higher derivatives can be acquired
during the training by using only the equation itself.
Consider a boundary value problem inside 2D region
Γ : x2 + y 2 < 1 for Poisson equation with nonlinear source:
3
uxx + uyy = u2 + u3
2

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

u|x2 +y2 =1 = −2
One can move all terms to the left side and make the following
notation:
3
ε = uxx + uyy − u2 − u3 = 0
2
To construct extra relations for cost function, derivatives of ε
will be used. For example, one can write a trivial consequence
of the previous equation:
9
εx = uxxx + uyyx − 2uux − u2 ux = 0
2
In this example, five ”new” equations are formed:
εx = 0, εy = 0, εxx = 0, εxy = 0, εyy = 0

(5)

As for boundary condition, direct inclusion of it into cost
function does not usually produce very accurate results [34].
A better choice is to use function substitution [27], after which
the boundary condition is satisfied automatically:
u = (1 − x2 − y 2 ) · n(x, y) − 2
Now n(x, y) can be found using a neural network with two
inputs and one output. The initial and additional equations (5)
are now written in terms of unknown function n. Boundary
condition for n(x, y) is:
n|x2 +y2 =1 < ∞
It holds for any MLP unless training produces increasingly
high weights (author never observed that in practice). Classical
approach for solving PDE is to minimize the following cost
function with respect to weights:
ε0 = ε2
Its performance is compared against the extended training that
uses the following extra terms:
ε1 = ε2x + ε2y
ε2 = ε2xx + ε2xy + ε2yy

7

for 1500 epochs with ∆0 = 2 · 10−5 . The last phase has 1000
epochs, ∆0 = 2 · 10−5 and the cost function is:
εtot = ε0
For this case, instead of giving classical training extra advantages, the total running times for both algorithms is to be
equalized. Grid for the extended backpropagation is generated
with spacing λ = 0.15 and contains 167 points. Classical
algorithm demonstrates heavy overfitting for this λ, so, before equalizing, the number of grid points is to be chosen.
Overfitting stops around λ = 0.05 which produces about
1300 points. The first, second and third stages of the extended
training use 15, 10 and 6 derivatives in each point respectively.
The only stage of classical training operates with 6. Using
mentioned phases lengths and ratio between numbers of grid
points, one can estimate a total number of 1700 epochs for the
regular training. Using smaller grids that produce overfitting
to increase the number of epochs does not lead to significant
changes. During the classical training after the first 1000
epochs ∆ are reset back to 2 · 10−4 .
Training for both approaches is repeated 100 times. Networks are verified against the analytical solution:
4
ua = 2
x + y2 − 3
on Cartesian grid with λ = 0.02 and about 8000 points.
Results are presented on Fig. 4. Curves are now much closer
than in previous two examples, partially because only two
extra orders of derivatives are used. Another reason of smaller
difference is network’s errors being very close to the machine
noise. If it was the only source of the error, RMS would be
around 0.36 · 10−6 = 10−6.44 .
One can notice that cost function for the last phase of
the extended training is the same as for regular training. As
unusual as it might be, this produces the ratio of cost on train
and test set εtrain /εtest of the order of 1, usually close to 1.2.
Classical training with the same cost used from the beginning
on a grid with 167 points results in said ratio being of the
order of 10, usually as high as 20.

For both cases discretization of cost’s measure µ is the same:

and lasts for 3500 epochs. After each 1000 epochs ∆ are reset
back to 2 · 10−4 . The second phase uses:

IV. C ONCLUSION
A training process that enhances approximation abilities of
fully connected feedforward neural networks was presented.
It is based on including extra derivatives in the cost function,
and was demonstrated to work well for low dimensional
cases. Using derivatives up to the fourth order, precision of
approximation for 2D analytical function was increased 200
times. With both - amount of information about the function
and running time equalized, classical training was far behind.
Number of epochs was so big that it is safe to say that
the regular training without derivatives can never achieve the
same level of accuracy and is always at least one order of
magnitude behind. The method was found to be not that useful
for neural networks with few hidden layers or with small
number of neurons in each layer. For example, in the case of
training autoencoder, network with the following configuration
of layers (asterisks denote linear activation function):

εtot = ε1 + ε0

3∗ , 64, 1∗ , 64, 3∗

N
1 X
ε(xα , yα , n(xα , yα ), ...)
E = µ (ε) '
N α=1
x,y∈Γ

where (xα , yα ) are points of regular grids. They are generated
in the following way: at first a Cartesian grid inside [−1, 1]2 is
created with spacing λ. All points that are not from Γ are then
excluded and additional points from the boundary with spacing
λ are added. Neural network has the following configuration
of layers:
2, 64, 64, 64, 64, 64, 64, 1
For the extended algorithm, training is divided into three
phases. The first one uses:
εtot = ε2 + ε1 + ε0

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

1.0

0.8

F(E)

0.6

0.4

0.2

E2
E0

0.0
-6.5

-6.0

-5.5

-5.0

-4.5

log E
Fig. 4. Cumulative distribution function of RMS of error E for the case of
solving partial differential equation. E0 is the classical procedure and E2
includes derivatives of residual up to the 2nd order.

has failed to perform at least one order of magnitude better
than the classical training. For deeper networks the gap between algorithms becomes quite substantial: after including
two hidden layers of width 64, it performed about 40 times
better than regular training, which showed almost no change
in precision. However, for real neural network applications
like classification one can imagine hard times calculating high
order derivatives. From this point of view, solving partial
differential equations can benefit much more from the proposed enhancements, as all information about extra derivatives
can be obtained via simple differentiation of the equation
itself. In presented example of solving PDE, the precision of
classical method was not too bad, however it required a grid
with spacing three times smaller. Even though the extended
training achieved 5 times smaller error using the same running
time, more important result was lowering the number of grid
points. If the same three-fold increase in grid spacing can
be demonstrated for higher dimensional PDEs, the extended
training would be even more advantageous.
ACKNOWLEDGMENT
Author owes a great debt of gratitude to his scientific
advisor to E.A. Dorotheyev. Special thanks for invaluable
support are due to Y.N. Sviridenko, A.M. Gaifullin, I.A. Avrutskaya and I.V. Avrutskiy without whom this work would be
impossible.
R EFERENCES
[1] Martı́n Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng
Chen, Craig Citro, Greg S Corrado, Andy Davis, Jeffrey Dean, Matthieu
Devin, et al. Tensorflow: Large-scale machine learning on heterogeneous
distributed systems. arXiv preprint arXiv:1603.04467, 2016.
[2] V.I. Avrutskiy. Backpropagation generalized for output derivatives. arXiv
preprint, 2017.
[3] Andrew R Barron. Universal approximation bounds for superpositions
of a sigmoidal function. IEEE Transactions on Information theory,
39(3):930–945, 1993.

8

[4] Andrew R Barron. Approximation and estimation bounds for artificial
neural networks. Machine Learning, 14(1):115–133, 1994.
[5] Yoshua Bengio et al. Learning deep architectures for ai. Foundations
and trends® in Machine Learning, 2(1):1–127, 2009.
[6] James Bergstra, Olivier Breuleux, Pascal Lamblin, Razvan Pascanu,
Olivier Delalleau, Guillaume Desjardins, Ian Goodfellow, Arnaud Bergeron, Yoshua Bengio, and Pack Kaelbling. Theano: Deep learning on
gpus with python. 2011.
[7] Christopher M Bishop. Neural networks for pattern recognition. Oxford
university press, 1995.
[8] Hervé Bourlard and Yves Kamp. Auto-association by multilayer
perceptrons and singular value decomposition. Biological cybernetics,
59(4):291–294, 1988.
[9] Pierre Cardaliaguet and Guillaume Euvrard. Approximation of a
function and its derivative with a neural network. Neural Networks,
5(2):207–220, 1992.
[10] Jia-Ren Chang and Yong-Sheng Chen. Batch-normalized maxout network in network. arXiv preprint arXiv:1511.02583, 2015.
[11] George Cybenko. Approximation by superpositions of a sigmoidal
function. Mathematics of Control, Signals, and Systems (MCSS),
2(4):303–314, 1989.
[12] Li Deng and Xiao Li. Machine learning paradigms for speech recognition: An overview. IEEE Transactions on Audio, Speech, and Language
Processing, 21(5):1060–1089, 2013.
[13] Harris Drucker and Yann Le Cun. Improving generalization performance
using double backpropagation. IEEE Transactions on Neural Networks,
3(6):991–997, 1992.
[14] 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.
[15] Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the Thirteenth
International Conference on Artificial Intelligence and Statistics, pages
249–256, 2010.
[16] Ian J Goodfellow, David Warde-Farley, Mehdi Mirza, Aaron Courville,
and Yoshua Bengio. Maxout networks. arXiv preprint arXiv:1302.4389,
2013.
[17] Alex Graves, Santiago Fernández, Faustino Gomez, and Jürgen Schmidhuber. Connectionist temporal classification: labelling unsegmented
sequence data with recurrent neural networks. In Proceedings of the 23rd
international conference on Machine learning, pages 369–376. ACM,
2006.
[18] Andreas Griewank and Andrea Walther. Evaluating derivatives: principles and techniques of algorithmic differentiation. SIAM, 2008.
[19] 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.
[20] Shouling He, Konrad Reif, and Rolf Unbehauen. Multilayer neural
networks for solving a class of partial differential equations. Neural
networks, 13(3):385–396, 2000.
[21] Kurt Hornik. Approximation capabilities of multilayer feedforward
networks. Neural networks, 4(2):251–257, 1991.
[22] Kurt Hornik, Maxwell Stinchcombe, and Halbert White. Multilayer
feedforward networks are universal approximators. Neural networks,
2(5):359–366, 1989.
[23] 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.
[24] 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.
[25] Vera Kurkova. Kolmogorov’s theorem and multilayer neural networks.
Neural networks, 5(3):501–506, 1992.
[26] IE Lagaris, A Likas, and DI Fotiadis. Artificial neural network methods
in quantum mechanics. Computer Physics Communications, 104(1-3):1–
14, 1997.
[27] 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.
[28] 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.
[29] AJ Meade and Alvaro A Fernandez. The numerical solution of
linear ordinary differential equations by feedforward neural networks.
Mathematical and Computer Modelling, 19(12):1–25, 1994.

SUBMITTED TO IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

[30] Dmytro Mishkin and Jiri Matas. All you need is a good init. arXiv
preprint arXiv:1511.06422, 2015.
[31] CUDA Nvidia. Cublas library. NVIDIA Corporation, Santa Clara,
California, 15(27):31, 2008.
[32] CUDA Nvidia. Programming guide, 2010.
[33] 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.
[34] 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.
[35] 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.
[36] Nitish Srivastava, Geoffrey E Hinton, Alex Krizhevsky, Ilya Sutskever,
and Ruslan Salakhutdinov. Dropout: a simple way to prevent neural
networks from overfitting. Journal of machine learning research,
15(1):1929–1958, 2014.
[37] Li Wan, Matthew Zeiler, Sixin Zhang, Yann Le Cun, and Rob Fergus.
Regularization of neural networks using dropconnect. In International
Conference on Machine Learning, pages 1058–1066, 2013.
[38] Guoqiang Peter Zhang. Neural networks for classification: a survey.
IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews), 30(4):451–462, 2000.

9

