АВТОМАТИЧЕСКОЕ ДИФФЕРЕНЦИРОВАНИЕ ФУНКЦИЙ, ВЫРАЖЕННЫХ ПРОГРАММНЫМ КОДОМ
ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА
УДК 519.688
К. К. СЕМЕНОВ
АВТОМАТИЧЕСКОЕ ДИФФЕРЕНЦИРОВАНИЕ ФУНКЦИЙ, ВЫРАЖЕННЫХ ПРОГРАММНЫМ КОДОМ
Алгебра дуальных чисел рассматривается как инструмент автоматического дифференцирования функций, выраженных программным кодом. Проанализированы специфические свойства дуальных чисел, необходимые для данного применения.
Ключевые слова: автоматическое дифференцирование, дуальные числа.
Введение. В практике программирования часто возникают задачи, когда необходимо
помимо вычисления значений некоторой функции f ( x) вычислять значения ее производной
df ( x)
dx
.
В
частности,
в
современной
метрологии
существует
задача
автоматической
оценки
интервала значений погрешности результатов вычислений y = f ( x) с помощью программ
обработки, исходные данные x для которых являются результатами прямых измерений. Классическим способом учета трансформированной погрешности является оценка на основе
производной f ′( x) . Специфика задачи состоит в том, что функция f ( x) реализуется в виде
программы вычислений и соответственно задана своим исходным кодом.
Традиционные способы оценки производной f ′( x) на практике сопряжены с трудно-
стями и имеют существенные недостатки. Реализация вычислений f ′( x) в виде отдельной
программной процедуры является избыточным решением и требует предварительного анали-
за функции f ( x) . Использование метода конечных разностей требует обоснованного выбора
значения приращения аргумента х, поскольку задача численного дифференцирования является некорректной и ее необходимо решать методом регуляризации [1].
Известен способ автоматического вычисления вместе с функцией f ( x0 ) при некотором
значении ее аргумента х = x0 значения производной этой функции f ′( x0 ) в той же точке [2].
Данный метод основан на алгебре дуальных чисел [3] и носит название „автоматическое
дифференцирование“, подчеркивающее одновременность вычислений значений f ( x0 ) и
f ′( x0 ) на основе только лишь исходного кода функции f ( x) . К сожалению, ощущается не-
хватка отечественной литературы, посвященной данному методу. Настоящая статья призвана частично восполнить данный недостаток.
Дуальные числа. Рассмотрим множество чисел z, таких, что z = x + ε∆ , где x ∈ R,
∆ ∈ R , а ε ≠ 0 — такая инфинитезимальная единица, что для нее выполняется точное равен-
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
Автоматическое дифференцирование функций, выраженных программным кодом
35
ство ε2 = 0 . Возможны математические операции сложения и умножения ε с вещественными числами. Числа z называют дуальными. Их составляющую x принято называть действительной частью z и обозначать Re z , а ∆ — инфинитезимальной частью z и обозначать Inf z .
Рассмотрим разложение в ряд Тейлора гладкой функции f ( x) одного действительного
аргумента:
f
(x)
=
f
(
x0
)
+
∞
∑
n=1
1 n!
d
n f ( x0
dxn
)
(
x
−
x0 )n .
Выполним его аналитическое продолжение на множество дуальных чисел: заменим в
выражении разложения f ( x) в ряд действительный аргумент х на дуальный z = x + ε∆ . То-
гда значение функции
f
( x + ε∆) =
f
∑(
x0
)
+
∞ n=1
1 n!
d
n
f ( x0
dxn
)
(
x
−
x0
+
ε∆
)n
,
согласно биному Ньютона, равно
∑ ∑f
( x0 ) +
∞ 1 d n f ( x0
n=1 n! dxn
)
n
Cnm
m=0
(x
−
x0 )n−m
εm∆m ,
что с учетом свойства ε2 = 0 есть
( )∑f
( x0
)
+
∞ n=1
1 n!
d n f ( x0
dxn
)
( x − x0 )n + nε∆ ( x − x0 )n−1
=
∑ ∑=
⎛ ⎜⎝⎜
f
(
x0
)
+
∞ n=1
1 n!
d
n f ( x0
dxn
)
(
x
−
x0
)n
⎞ ⎟⎠⎟
+
ε∆
∞ m=1
(m
1 −
1)!
d
m f ( x0
dxm
)
(
x
−
x0 )m−1
=
f
(x) +
ε∆
df (x) dx
.
Получаем, что результат вычислений является дуальным числом, одна составляющая
которого есть точное значение функции f ( x) , а другая — точное значение производной этой
же функции при том же значении аргумента х. Таким образом, возможно использовать дуальные числа для автоматического вычисле-
ния значения производной f ′( x) точно вплоть до ошибки округления. Для того чтобы вос-
пользоваться этим приемом, на практике необходимо создать новый тип данных. Поскольку
дуальные
числа
z=
x + ε∆
суть пары
вещественных
чисел
⎛ ⎜ ⎝
x ∆
⎞ ⎟ ⎠
,
то
можно
на
языке
С++
вве-
сти их следующим образом.
тип данных, описывающий дуальное число.
struct dual {
double real;
// действительная часть.
double inf;
// инфинитезимальная часть.
}
Для того чтобы реализовать автоматическое дифференцирование функции f ( x) , пред-
ставленной своим исходным кодом, необходимо перейти в программе от вычислений с действительными данными (double) к вычислениям с дуальными числами (dual).
Вычисления, производимые в любой программе математической обработки, являют собой вычисления значений сложных функций, например, для вычисления значения
( ( ) )y = f g1 g2 (...gn ( x)...), х , х необходимо начать с вычисления значения yn = gn ( x) ,
продолжить вычислениями yn−1 = gn−1 ( yn , x) , yn−2 = gn−2 ( yn−1, x) и т.д. Как известно,
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
36 К. К. Семенов
( ( ) )производная сложной функции y = f g1 g2 (...gn ( x)...), х , х по аргументу х вычисляется
по цепному правилу, т.е.
dgi
( yi+1,
dx
x)
=
∂gi
( yi+1,
∂yi+1
x)
⋅
dgi+1
( yi+2
dx
,
x)
+
∂gi
(
yi +1 , ∂x
x)
для i = 1, 2, …, n. Следовательно, можно вычислять значение производной
df ( x)
dx
в последо-
( ( ) )вательности вычисления значения функции y = f g1 g2 (...gn ( x)...), х , х , т.е. начиная с
dgn (
dx
x
)
,
продолжая
вычислением
и т.д. вплоть до
dgn−1 ( yn
dx
,
x)
=
∂gn−1 ( yn ,
∂yn
x)
⋅
dgn ( x)
dx
+
∂gn−1 ( yn ,
∂x
x)
df ( x1 )
dx1
=
∂g1 ( y2
∂y2
,
x)
⋅
dg2 (
dx
x)
+
∂g1
( y2 ,
∂x
x)
.
Таким образом, можно выполнить автоматическое дифференцирование одновременно с вычислениями основной функции.
На практике математические преобразования, реализуемые компьютерными программами, сводятся к тому, что все функции gi , вообще говоря, являют собой последовательность
элементарных арифметических действий из набора {+, −, ×, /} и элементарных функций —
математических примитивов из стандартных библиотек. Можно констатировать, что приведенное рассуждение о вычислении производной сложной функции справедливо и для функции, реализованной программой. Таким образом, для работы с дуальными числами необходимо задать основные арифметические действия с ними и указать, как производить вычисления элементарных математических функций при дуальном аргументе.
С математической точки зрения эта задача равносильна определению алгебры дуальных
чисел Λ = R2 , +, × как соответствующей структуры над множеством пар действительных
чисел R2 с операциями сложения и умножения.
Программно-ориентированное задание алгебры дуальных чисел. Рассмотрим мно-
жество пар действительных чисел R2 и зададим для двух произвольных его элементов
z1
=
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
и
z2
=
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
арифметические
операции
следующим образом.
Сложение будем выполнять по правилу
z1
+
z2
=
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
+
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
x1 y1
+ +
x2 y2
⎞ ⎟
.
⎠
(1)
Операция обладает свойствами коммутативности и ассоциативности. Действительно,
выполнены соответствующие соотношения z1 + z2 = z2 + z1 и (z1 + z2 ) + z3 = z1 + (z2 + z3 ) .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
Автоматическое дифференцирование функций, выраженных программным кодом
37
Существует
нейтральный
элемент
операции
сложения
∅
=
⎛ ⎜ ⎝
0⎞ 0 ⎟⎠
,
справедливо
равенство
z
+
∅
=
⎛ ⎜ ⎝
x y
⎞ ⎟ ⎠
+
⎛ ⎜ ⎝
0 0
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
x y
⎞ ⎟ ⎠
=
z
.
Вычитание — обратная сложению операция. При выводе правил для компонентов ее
результата
требуется
решить
уравнение
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
+
⎛ ⎜ ⎝
x3 y3
⎞ ⎟ ⎠
относительно
⎛ ⎜ ⎝
x3 y3
⎞ ⎟ ⎠
.
Получаем
систему
уравнений
⎧ ⎨ ⎩
x1 y1
= =
x2 y2
+ +
x3 y3
,
ее
решением
является
⎧ ⎨ ⎩
x3 y3
= =
x1 y1
− −
x2 y2
.
Следовательно,
z1
−
z2
=
⎛ ⎜ ⎝
x1 ⎞
y1
⎟ ⎠
−
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
x1 y1
− −
x2 y2
⎞ ⎟ ⎠
.
(2)
Умножение: введем операцию умножения как
z1
× z2
=
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
y1
⋅
x1 ⋅ x2 x2 + y2
⋅
x1
⎞ ⎟ ⎠
.
(3)
Операция умножения также коммутативна и ассоциативна, поскольку z1 × z2 = z2 × z1 и
( z1
× z2
)× z3
=
z1
× (z2
×
z3 ) .
Для
нее
существует
нейтральный
элемент
1=
⎛ ⎜ ⎝
1⎞ 0 ⎠⎟
.
Действитель-
но,
z
×1
=
⎛ ⎜ ⎝
x y
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
1 0
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
y
x ⋅1 ⋅1+ 0⋅
⎞
х
⎟ ⎠
=
⎛ ⎜ ⎝
x⎞
y
⎟ ⎠
=
z.
Деление — обратная умножению операция. При выводе правил для компонентов ее ре-
зультата
требуется
решить
уравнение
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
x3 y3
⎞ ⎟ ⎠
относительно
⎛ ⎜ ⎝
x3 y3
⎞ ⎟ ⎠
.
Получаем
сис-
тему уравнений
⎧
⎨ ⎩
y1
x1 = x2 x3 = y2 x3 + y3 x2
⇒
⎧
⎪⎪ ⎨
⎪ ⎩⎪
y3
x2
x3 =
= y1
x1 x2
− y2
x1 x2
⇒
⎧
⎪⎪
⎨
⎪ ⎪⎩
y3
x3
=
x1 x2
=
y1x2 − y2 x1 x22
.
В качестве вывода получаем, что при x2 ≠ 0
z1 z2
=
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
⎛ x1 ⎞
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
=
⎜ ⎜ ⎜ ⎜ ⎝
y1 x2
x2 − x22
y2
x1
⎟
⎟ ⎟
.
⎟
⎠
(4)
Имеет место также дистрибутивный закон: z1 × (z2 + z3 ) = z1 × z2 + z1 × z3 .
Для таким образом введенной алгебры Λ = R2 , +, × дуальных чисел справедливо сле-
дующее свойство:
⎛ ⎜ ⎝
0 y1
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
0 y2
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
0 0
⎞ ⎟ ⎠
,
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
38 К. К. Семенов
т.е. для квадрата чисто инфинитезимального числа z всегда имеет место свойство z2 = 0. Таким образом, из свойств арифметических операций (1)—(4) для инфинитезимальной едини-
цы ε вытекает ее свойство ε2 = 0 . В силу того что справедливы соотношения
Re (z1 ± z2 ) = Re (z1) ± Re (z2 ),
Re (z1 × z2 ) = Re(z1) Re(z2 ),
Re (z1 / z2 ) = Re (z1 ) / Re (z2 ),
алгебра Λ описывает такое расширение множества действительных чисел на двумерную
плоскость, которое включает в себя поле действительных чисел как составную часть. Дейст-
вия с инфинитезимальными составляющими дуальных чисел не влияют на значения их дей-
ствительных составляющих.
Результат
вычисления
элементарной
функции
g (z)
от
дуального
аргумента
z
=
⎛ ⎜ ⎝
x⎞
y
⎟ ⎠
равен
g
(
z
)
=
⎛ ⎜ ⎝
g ( x) yg′( x)
⎞ ⎟ ⎠
,
что следует из рассмотренного выше аналитического продолжение ряда Тейлора для функции
g ( x) . Указанное свойство задает правило вычисления элементарных функций от дуального
аргумента, которое необходимо использовать для перегрузки математических примитивов.
процедура для вычислений значения
тригонометрического синуса
dual sin (dual x)
{
dual z;
z.real = sin(x.real);
z.inf = x.inf * cos(x.real);
return (z); }
процедура для вычислений значения гиперболического синуса
dual sinh (dual x) { dual z; z.real = sinh(x.real); z.inf = x.inf * cosh(x.real); return (z); }
Рассмотрим пример. Пусть требуется вычислить значение функции
( )f ( x) = sin ( x) exp −x2 + x в точке x0 . В соответствии с представленными выше правилами
произведем последовательно вычисления, заменив действительный аргумент x = x0 дуаль-
ным
z0
=
⎛ x0
⎜ ⎝
1
⎞ ⎟
.
⎠
f
(z0
)
=
sin
⎡⎛ ⎢⎜ ⎣⎝
x0 1
⎞⎤ ⎟⎥ ⎠⎦
×
exp
⎡ ⎢−
⎛ ⎜
⎣⎝
x0 1
⎞ ⎟
×
⎛ ⎜
⎠⎝
x0 1
⎞⎤ ⎟⎥ ⎠⎦
+
⎛ ⎜ ⎝
x0 1
⎞ ⎟ ⎠
=
⎛ ⎜⎝1
sin (
⋅ cos
x0 ) ( x0
⎞
)
⎟ ⎠
×
exp
⎡⎛ ⎢⎢⎣⎜⎝⎜
− x02 −2x0
⎞⎤ ⎟⎟⎠⎥⎥⎦
+
⎛ ⎜ ⎝
x0 1
⎞ ⎟ ⎠
=
( ( ) ) ( ) ( ) ( )=
⎛ ⎜ ⎝
sin ( cos (
x0 x0
) )
⎞ ⎟ ⎠
×
⎛ ⎜ ⎜⎜⎝
exp −x02 −2x0 exp −
x02
⎞ ⎟ ⎟⎟⎠
+
⎛ ⎜ ⎝
x0 1
⎞ ⎟ ⎠
=
⎛ ⎜ ⎜⎜⎝
−2
x0
sin
(
x0
sin ( ) exp
x0 −
) exp
x02 +
− x02
cos (
x0
)
exp
− x02
⎞ ⎟ ⎟⎟⎠
+
⎛ ⎜ ⎝
x0 1
⎞ ⎟ ⎠
=
( ) ( ) ( )=
⎛ ⎜ ⎜⎜⎝
−2
x0
sin
(
x0
sin ( x0 ) exp −x02 ) exp −x02 + cos
(
+ x0
x0
) exp
− x02
⎞
+
⎟ 1⎟⎟⎠
.
( ) ( )Действительно, f ′( x) = −2x exp −x2 sin ( x) + cos ( x) exp −x2 + 1 , что в точности сов-
падает с результатом вычислений при помощи дуальных чисел.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
Автоматическое дифференцирование функций, выраженных программным кодом
39
Если,
например,
x0
=
0,
то получим
z0
=
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
и
sin
z0
=
sin
⎡⎛ 0 ⎞⎤
⎢⎜ ⎣⎝
1
⎟⎥ ⎠⎦
=
⎛ sin 0 ⎞
⎜⎝1⋅
cos
0
⎟ ⎠
=
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
,
( )exp −z02
=
exp
⎡ ⎢− ⎣
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
0 1
⎞⎤ ⎟⎥ ⎠⎦
=
exp
⎡⎛ ⎢⎜ ⎣⎝
0 0
⎞⎤ ⎟⎥ ⎠⎦
=
⎛ ⎜ ⎝
exp 0 0 ⋅ exp 0
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
1⎞ 0 ⎠⎟
.
Следовательно,
( )f (z0 ) = sin z0 × exp
− z02
+
z0
=
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
1 0
⎞ ⎟ ⎠
+
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
=
⎛0⎞
⎜ ⎝
1
⎟ ⎠
+
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
0⎞ 2 ⎟⎠
,
откуда получаем, что f (0) = 0 и f ′(0) = 2 .
Рассмотренный метод давно и с успехом применяется в самых различных приложениях. Представленное выше описание соответствует так называемому прямому ходу (forward mode) автоматического дифференцирования, когда значение производной вычисляется параллельно вычислению значения функции при движении от начала алгоритма к его концу. Метод легко обобщается на случай одновременного вычисления всех частных производных одновременно с вычислением значения функции.
Из всего сказанного следует, что метод автоматического дифференцирования характеризуется следующими свойствами:
— автоматическое дифференцирование функций, выраженных программным кодом, носит аналитический характер и позволяет получить точное (вплоть до ошибок округления) значение производной;
— использование подхода не влияет на сходимость задействованных в выполняемой программе итерационных методов и не искажает конечного результата вычислений;
— автоматическое дифференцирование позволяет выполнить программу вычислений лишь единожды в отличие от приближенных методик численного дифференцирования на основе конечных разностей;
— использование автоматического дифференцирования позволяет не изменять порядок вычислений в программах.
Заключение. Метод автоматического дифференцирования функций, выраженных программным кодом, обладает рядом полезных для метрологической практики свойств, что обусловливает возможность его применения в метрологических целях [4, 5]. В частности, с помощью метода можно автоматически получать оценки характеристик трансформированной погрешности результатов вычислений, выполняемых в измерительных системах с цифровой обработкой сигналов [6].
СПИСОК ЛИТЕРАТУРЫ
1. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. М.: Наука, Главная редакция физикоматематической литературы, 1979.
2. Corliss G., Faure С., Griewank A., Hascoлt L., Naumann U. Automatic Differentiation Bibliography // Automatic Differentiation of Algorithms: From Simulation to Optimization. Springer, 2002. P. 383—425.
3. Яглом И. М. Комплексные числа и их применение в геометрии. М.: Государственное изд-во физикоматематической литературы, 1963.
4. Hall B. D. Calculating measurements uncertainty using automatic differentiation // Measurement Science and Technology. 2002. Vol. 13, N 4. P. 421—427.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
40 К. Г. Коротков, Д. В. Орлов, Е. Н. Величко
5. Luca Mari. A computational system for uncertainty propagation of measurement results // Measurement. 2009. Vol. 42, is. 6. P. 844—855.
6. Семенов К. К., Солопченко Г. Н. Теоретические предпосылки реализации метрологического автосопровождения программ обработки результатов измерений // Измерительная техника. 2010. № 6. С. 9—14.
Константин Константинович Семенов
Сведения об авторе — аспирант; Санкт-Петербургский государственный политехни-
ческий университет, кафедра измерительных информационных технологий; E-mail: semenov.k.k@gmail.com
Рекомендована кафедрой измерительных информационных технологий
Поступила в редакцию 23.06.11 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
УДК 519.688
К. К. СЕМЕНОВ
АВТОМАТИЧЕСКОЕ ДИФФЕРЕНЦИРОВАНИЕ ФУНКЦИЙ, ВЫРАЖЕННЫХ ПРОГРАММНЫМ КОДОМ
Алгебра дуальных чисел рассматривается как инструмент автоматического дифференцирования функций, выраженных программным кодом. Проанализированы специфические свойства дуальных чисел, необходимые для данного применения.
Ключевые слова: автоматическое дифференцирование, дуальные числа.
Введение. В практике программирования часто возникают задачи, когда необходимо
помимо вычисления значений некоторой функции f ( x) вычислять значения ее производной
df ( x)
dx
.
В
частности,
в
современной
метрологии
существует
задача
автоматической
оценки
интервала значений погрешности результатов вычислений y = f ( x) с помощью программ
обработки, исходные данные x для которых являются результатами прямых измерений. Классическим способом учета трансформированной погрешности является оценка на основе
производной f ′( x) . Специфика задачи состоит в том, что функция f ( x) реализуется в виде
программы вычислений и соответственно задана своим исходным кодом.
Традиционные способы оценки производной f ′( x) на практике сопряжены с трудно-
стями и имеют существенные недостатки. Реализация вычислений f ′( x) в виде отдельной
программной процедуры является избыточным решением и требует предварительного анали-
за функции f ( x) . Использование метода конечных разностей требует обоснованного выбора
значения приращения аргумента х, поскольку задача численного дифференцирования является некорректной и ее необходимо решать методом регуляризации [1].
Известен способ автоматического вычисления вместе с функцией f ( x0 ) при некотором
значении ее аргумента х = x0 значения производной этой функции f ′( x0 ) в той же точке [2].
Данный метод основан на алгебре дуальных чисел [3] и носит название „автоматическое
дифференцирование“, подчеркивающее одновременность вычислений значений f ( x0 ) и
f ′( x0 ) на основе только лишь исходного кода функции f ( x) . К сожалению, ощущается не-
хватка отечественной литературы, посвященной данному методу. Настоящая статья призвана частично восполнить данный недостаток.
Дуальные числа. Рассмотрим множество чисел z, таких, что z = x + ε∆ , где x ∈ R,
∆ ∈ R , а ε ≠ 0 — такая инфинитезимальная единица, что для нее выполняется точное равен-
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
Автоматическое дифференцирование функций, выраженных программным кодом
35
ство ε2 = 0 . Возможны математические операции сложения и умножения ε с вещественными числами. Числа z называют дуальными. Их составляющую x принято называть действительной частью z и обозначать Re z , а ∆ — инфинитезимальной частью z и обозначать Inf z .
Рассмотрим разложение в ряд Тейлора гладкой функции f ( x) одного действительного
аргумента:
f
(x)
=
f
(
x0
)
+
∞
∑
n=1
1 n!
d
n f ( x0
dxn
)
(
x
−
x0 )n .
Выполним его аналитическое продолжение на множество дуальных чисел: заменим в
выражении разложения f ( x) в ряд действительный аргумент х на дуальный z = x + ε∆ . То-
гда значение функции
f
( x + ε∆) =
f
∑(
x0
)
+
∞ n=1
1 n!
d
n
f ( x0
dxn
)
(
x
−
x0
+
ε∆
)n
,
согласно биному Ньютона, равно
∑ ∑f
( x0 ) +
∞ 1 d n f ( x0
n=1 n! dxn
)
n
Cnm
m=0
(x
−
x0 )n−m
εm∆m ,
что с учетом свойства ε2 = 0 есть
( )∑f
( x0
)
+
∞ n=1
1 n!
d n f ( x0
dxn
)
( x − x0 )n + nε∆ ( x − x0 )n−1
=
∑ ∑=
⎛ ⎜⎝⎜
f
(
x0
)
+
∞ n=1
1 n!
d
n f ( x0
dxn
)
(
x
−
x0
)n
⎞ ⎟⎠⎟
+
ε∆
∞ m=1
(m
1 −
1)!
d
m f ( x0
dxm
)
(
x
−
x0 )m−1
=
f
(x) +
ε∆
df (x) dx
.
Получаем, что результат вычислений является дуальным числом, одна составляющая
которого есть точное значение функции f ( x) , а другая — точное значение производной этой
же функции при том же значении аргумента х. Таким образом, возможно использовать дуальные числа для автоматического вычисле-
ния значения производной f ′( x) точно вплоть до ошибки округления. Для того чтобы вос-
пользоваться этим приемом, на практике необходимо создать новый тип данных. Поскольку
дуальные
числа
z=
x + ε∆
суть пары
вещественных
чисел
⎛ ⎜ ⎝
x ∆
⎞ ⎟ ⎠
,
то
можно
на
языке
С++
вве-
сти их следующим образом.
тип данных, описывающий дуальное число.
struct dual {
double real;
// действительная часть.
double inf;
// инфинитезимальная часть.
}
Для того чтобы реализовать автоматическое дифференцирование функции f ( x) , пред-
ставленной своим исходным кодом, необходимо перейти в программе от вычислений с действительными данными (double) к вычислениям с дуальными числами (dual).
Вычисления, производимые в любой программе математической обработки, являют собой вычисления значений сложных функций, например, для вычисления значения
( ( ) )y = f g1 g2 (...gn ( x)...), х , х необходимо начать с вычисления значения yn = gn ( x) ,
продолжить вычислениями yn−1 = gn−1 ( yn , x) , yn−2 = gn−2 ( yn−1, x) и т.д. Как известно,
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
36 К. К. Семенов
( ( ) )производная сложной функции y = f g1 g2 (...gn ( x)...), х , х по аргументу х вычисляется
по цепному правилу, т.е.
dgi
( yi+1,
dx
x)
=
∂gi
( yi+1,
∂yi+1
x)
⋅
dgi+1
( yi+2
dx
,
x)
+
∂gi
(
yi +1 , ∂x
x)
для i = 1, 2, …, n. Следовательно, можно вычислять значение производной
df ( x)
dx
в последо-
( ( ) )вательности вычисления значения функции y = f g1 g2 (...gn ( x)...), х , х , т.е. начиная с
dgn (
dx
x
)
,
продолжая
вычислением
и т.д. вплоть до
dgn−1 ( yn
dx
,
x)
=
∂gn−1 ( yn ,
∂yn
x)
⋅
dgn ( x)
dx
+
∂gn−1 ( yn ,
∂x
x)
df ( x1 )
dx1
=
∂g1 ( y2
∂y2
,
x)
⋅
dg2 (
dx
x)
+
∂g1
( y2 ,
∂x
x)
.
Таким образом, можно выполнить автоматическое дифференцирование одновременно с вычислениями основной функции.
На практике математические преобразования, реализуемые компьютерными программами, сводятся к тому, что все функции gi , вообще говоря, являют собой последовательность
элементарных арифметических действий из набора {+, −, ×, /} и элементарных функций —
математических примитивов из стандартных библиотек. Можно констатировать, что приведенное рассуждение о вычислении производной сложной функции справедливо и для функции, реализованной программой. Таким образом, для работы с дуальными числами необходимо задать основные арифметические действия с ними и указать, как производить вычисления элементарных математических функций при дуальном аргументе.
С математической точки зрения эта задача равносильна определению алгебры дуальных
чисел Λ = R2 , +, × как соответствующей структуры над множеством пар действительных
чисел R2 с операциями сложения и умножения.
Программно-ориентированное задание алгебры дуальных чисел. Рассмотрим мно-
жество пар действительных чисел R2 и зададим для двух произвольных его элементов
z1
=
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
и
z2
=
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
арифметические
операции
следующим образом.
Сложение будем выполнять по правилу
z1
+
z2
=
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
+
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
x1 y1
+ +
x2 y2
⎞ ⎟
.
⎠
(1)
Операция обладает свойствами коммутативности и ассоциативности. Действительно,
выполнены соответствующие соотношения z1 + z2 = z2 + z1 и (z1 + z2 ) + z3 = z1 + (z2 + z3 ) .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
Автоматическое дифференцирование функций, выраженных программным кодом
37
Существует
нейтральный
элемент
операции
сложения
∅
=
⎛ ⎜ ⎝
0⎞ 0 ⎟⎠
,
справедливо
равенство
z
+
∅
=
⎛ ⎜ ⎝
x y
⎞ ⎟ ⎠
+
⎛ ⎜ ⎝
0 0
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
x y
⎞ ⎟ ⎠
=
z
.
Вычитание — обратная сложению операция. При выводе правил для компонентов ее
результата
требуется
решить
уравнение
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
+
⎛ ⎜ ⎝
x3 y3
⎞ ⎟ ⎠
относительно
⎛ ⎜ ⎝
x3 y3
⎞ ⎟ ⎠
.
Получаем
систему
уравнений
⎧ ⎨ ⎩
x1 y1
= =
x2 y2
+ +
x3 y3
,
ее
решением
является
⎧ ⎨ ⎩
x3 y3
= =
x1 y1
− −
x2 y2
.
Следовательно,
z1
−
z2
=
⎛ ⎜ ⎝
x1 ⎞
y1
⎟ ⎠
−
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
x1 y1
− −
x2 y2
⎞ ⎟ ⎠
.
(2)
Умножение: введем операцию умножения как
z1
× z2
=
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
y1
⋅
x1 ⋅ x2 x2 + y2
⋅
x1
⎞ ⎟ ⎠
.
(3)
Операция умножения также коммутативна и ассоциативна, поскольку z1 × z2 = z2 × z1 и
( z1
× z2
)× z3
=
z1
× (z2
×
z3 ) .
Для
нее
существует
нейтральный
элемент
1=
⎛ ⎜ ⎝
1⎞ 0 ⎠⎟
.
Действитель-
но,
z
×1
=
⎛ ⎜ ⎝
x y
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
1 0
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
y
x ⋅1 ⋅1+ 0⋅
⎞
х
⎟ ⎠
=
⎛ ⎜ ⎝
x⎞
y
⎟ ⎠
=
z.
Деление — обратная умножению операция. При выводе правил для компонентов ее ре-
зультата
требуется
решить
уравнение
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
x3 y3
⎞ ⎟ ⎠
относительно
⎛ ⎜ ⎝
x3 y3
⎞ ⎟ ⎠
.
Получаем
сис-
тему уравнений
⎧
⎨ ⎩
y1
x1 = x2 x3 = y2 x3 + y3 x2
⇒
⎧
⎪⎪ ⎨
⎪ ⎩⎪
y3
x2
x3 =
= y1
x1 x2
− y2
x1 x2
⇒
⎧
⎪⎪
⎨
⎪ ⎪⎩
y3
x3
=
x1 x2
=
y1x2 − y2 x1 x22
.
В качестве вывода получаем, что при x2 ≠ 0
z1 z2
=
⎛ ⎜ ⎝
x1 y1
⎞ ⎟ ⎠
⎛ x1 ⎞
⎛ ⎜ ⎝
x2 y2
⎞ ⎟ ⎠
=
⎜ ⎜ ⎜ ⎜ ⎝
y1 x2
x2 − x22
y2
x1
⎟
⎟ ⎟
.
⎟
⎠
(4)
Имеет место также дистрибутивный закон: z1 × (z2 + z3 ) = z1 × z2 + z1 × z3 .
Для таким образом введенной алгебры Λ = R2 , +, × дуальных чисел справедливо сле-
дующее свойство:
⎛ ⎜ ⎝
0 y1
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
0 y2
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
0 0
⎞ ⎟ ⎠
,
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
38 К. К. Семенов
т.е. для квадрата чисто инфинитезимального числа z всегда имеет место свойство z2 = 0. Таким образом, из свойств арифметических операций (1)—(4) для инфинитезимальной едини-
цы ε вытекает ее свойство ε2 = 0 . В силу того что справедливы соотношения
Re (z1 ± z2 ) = Re (z1) ± Re (z2 ),
Re (z1 × z2 ) = Re(z1) Re(z2 ),
Re (z1 / z2 ) = Re (z1 ) / Re (z2 ),
алгебра Λ описывает такое расширение множества действительных чисел на двумерную
плоскость, которое включает в себя поле действительных чисел как составную часть. Дейст-
вия с инфинитезимальными составляющими дуальных чисел не влияют на значения их дей-
ствительных составляющих.
Результат
вычисления
элементарной
функции
g (z)
от
дуального
аргумента
z
=
⎛ ⎜ ⎝
x⎞
y
⎟ ⎠
равен
g
(
z
)
=
⎛ ⎜ ⎝
g ( x) yg′( x)
⎞ ⎟ ⎠
,
что следует из рассмотренного выше аналитического продолжение ряда Тейлора для функции
g ( x) . Указанное свойство задает правило вычисления элементарных функций от дуального
аргумента, которое необходимо использовать для перегрузки математических примитивов.
процедура для вычислений значения
тригонометрического синуса
dual sin (dual x)
{
dual z;
z.real = sin(x.real);
z.inf = x.inf * cos(x.real);
return (z); }
процедура для вычислений значения гиперболического синуса
dual sinh (dual x) { dual z; z.real = sinh(x.real); z.inf = x.inf * cosh(x.real); return (z); }
Рассмотрим пример. Пусть требуется вычислить значение функции
( )f ( x) = sin ( x) exp −x2 + x в точке x0 . В соответствии с представленными выше правилами
произведем последовательно вычисления, заменив действительный аргумент x = x0 дуаль-
ным
z0
=
⎛ x0
⎜ ⎝
1
⎞ ⎟
.
⎠
f
(z0
)
=
sin
⎡⎛ ⎢⎜ ⎣⎝
x0 1
⎞⎤ ⎟⎥ ⎠⎦
×
exp
⎡ ⎢−
⎛ ⎜
⎣⎝
x0 1
⎞ ⎟
×
⎛ ⎜
⎠⎝
x0 1
⎞⎤ ⎟⎥ ⎠⎦
+
⎛ ⎜ ⎝
x0 1
⎞ ⎟ ⎠
=
⎛ ⎜⎝1
sin (
⋅ cos
x0 ) ( x0
⎞
)
⎟ ⎠
×
exp
⎡⎛ ⎢⎢⎣⎜⎝⎜
− x02 −2x0
⎞⎤ ⎟⎟⎠⎥⎥⎦
+
⎛ ⎜ ⎝
x0 1
⎞ ⎟ ⎠
=
( ( ) ) ( ) ( ) ( )=
⎛ ⎜ ⎝
sin ( cos (
x0 x0
) )
⎞ ⎟ ⎠
×
⎛ ⎜ ⎜⎜⎝
exp −x02 −2x0 exp −
x02
⎞ ⎟ ⎟⎟⎠
+
⎛ ⎜ ⎝
x0 1
⎞ ⎟ ⎠
=
⎛ ⎜ ⎜⎜⎝
−2
x0
sin
(
x0
sin ( ) exp
x0 −
) exp
x02 +
− x02
cos (
x0
)
exp
− x02
⎞ ⎟ ⎟⎟⎠
+
⎛ ⎜ ⎝
x0 1
⎞ ⎟ ⎠
=
( ) ( ) ( )=
⎛ ⎜ ⎜⎜⎝
−2
x0
sin
(
x0
sin ( x0 ) exp −x02 ) exp −x02 + cos
(
+ x0
x0
) exp
− x02
⎞
+
⎟ 1⎟⎟⎠
.
( ) ( )Действительно, f ′( x) = −2x exp −x2 sin ( x) + cos ( x) exp −x2 + 1 , что в точности сов-
падает с результатом вычислений при помощи дуальных чисел.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
Автоматическое дифференцирование функций, выраженных программным кодом
39
Если,
например,
x0
=
0,
то получим
z0
=
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
и
sin
z0
=
sin
⎡⎛ 0 ⎞⎤
⎢⎜ ⎣⎝
1
⎟⎥ ⎠⎦
=
⎛ sin 0 ⎞
⎜⎝1⋅
cos
0
⎟ ⎠
=
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
,
( )exp −z02
=
exp
⎡ ⎢− ⎣
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
0 1
⎞⎤ ⎟⎥ ⎠⎦
=
exp
⎡⎛ ⎢⎜ ⎣⎝
0 0
⎞⎤ ⎟⎥ ⎠⎦
=
⎛ ⎜ ⎝
exp 0 0 ⋅ exp 0
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
1⎞ 0 ⎠⎟
.
Следовательно,
( )f (z0 ) = sin z0 × exp
− z02
+
z0
=
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
×
⎛ ⎜ ⎝
1 0
⎞ ⎟ ⎠
+
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
=
⎛0⎞
⎜ ⎝
1
⎟ ⎠
+
⎛ ⎜ ⎝
0 1
⎞ ⎟ ⎠
=
⎛ ⎜ ⎝
0⎞ 2 ⎟⎠
,
откуда получаем, что f (0) = 0 и f ′(0) = 2 .
Рассмотренный метод давно и с успехом применяется в самых различных приложениях. Представленное выше описание соответствует так называемому прямому ходу (forward mode) автоматического дифференцирования, когда значение производной вычисляется параллельно вычислению значения функции при движении от начала алгоритма к его концу. Метод легко обобщается на случай одновременного вычисления всех частных производных одновременно с вычислением значения функции.
Из всего сказанного следует, что метод автоматического дифференцирования характеризуется следующими свойствами:
— автоматическое дифференцирование функций, выраженных программным кодом, носит аналитический характер и позволяет получить точное (вплоть до ошибок округления) значение производной;
— использование подхода не влияет на сходимость задействованных в выполняемой программе итерационных методов и не искажает конечного результата вычислений;
— автоматическое дифференцирование позволяет выполнить программу вычислений лишь единожды в отличие от приближенных методик численного дифференцирования на основе конечных разностей;
— использование автоматического дифференцирования позволяет не изменять порядок вычислений в программах.
Заключение. Метод автоматического дифференцирования функций, выраженных программным кодом, обладает рядом полезных для метрологической практики свойств, что обусловливает возможность его применения в метрологических целях [4, 5]. В частности, с помощью метода можно автоматически получать оценки характеристик трансформированной погрешности результатов вычислений, выполняемых в измерительных системах с цифровой обработкой сигналов [6].
СПИСОК ЛИТЕРАТУРЫ
1. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. М.: Наука, Главная редакция физикоматематической литературы, 1979.
2. Corliss G., Faure С., Griewank A., Hascoлt L., Naumann U. Automatic Differentiation Bibliography // Automatic Differentiation of Algorithms: From Simulation to Optimization. Springer, 2002. P. 383—425.
3. Яглом И. М. Комплексные числа и их применение в геометрии. М.: Государственное изд-во физикоматематической литературы, 1963.
4. Hall B. D. Calculating measurements uncertainty using automatic differentiation // Measurement Science and Technology. 2002. Vol. 13, N 4. P. 421—427.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12
40 К. Г. Коротков, Д. В. Орлов, Е. Н. Величко
5. Luca Mari. A computational system for uncertainty propagation of measurement results // Measurement. 2009. Vol. 42, is. 6. P. 844—855.
6. Семенов К. К., Солопченко Г. Н. Теоретические предпосылки реализации метрологического автосопровождения программ обработки результатов измерений // Измерительная техника. 2010. № 6. С. 9—14.
Константин Константинович Семенов
Сведения об авторе — аспирант; Санкт-Петербургский государственный политехни-
ческий университет, кафедра измерительных информационных технологий; E-mail: semenov.k.k@gmail.com
Рекомендована кафедрой измерительных информационных технологий
Поступила в редакцию 23.06.11 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 12