Автоматическое дифференцирование функций, выраженных программным кодом Текст научной статьи по специальности «Математика»

ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА

УДК 519.688

К. К. Семенов

АВТОМАТИЧЕСКОЕ ДИФФЕРЕНЦИРОВАНИЕ ФУНКЦИЙ, ВЫРАЖЕННЫХ ПРОГРАММНЫМ КОДОМ

Алгебра дуальных чисел рассматривается как инструмент автоматического дифференцирования функций, выраженных программным кодом. Проанализированы специфические свойства дуальных чисел, необходимые для данного применения.

Ключевые слова: автоматическое дифференцирование, дуальные числа.

Введение. В практике программирования часто возникают задачи, когда необходимо помимо вычисления значений некоторой функции / (х) вычислять значения ее производной

#(х) В й й

-. В частности, в современной метрологии существует задача автоматической оценки

дх

интервала значений погрешности результатов вычислений у = / (х) с помощью программ

обработки, исходные данные х для которых являются результатами прямых измерений. Классическим способом учета трансформированной погрешности является оценка на основе производной /'(х) . Специфика задачи состоит в том, что функция / (х) реализуется в виде программы вычислений и соответственно задана своим исходным кодом.

Традиционные способы оценки производной /'(х) на практике сопряжены с трудностями и имеют существенные недостатки. Реализация вычислений /'(х) в виде отдельной программной процедуры является избыточным решением и требует предварительного анализа функции / (х) . Использование метода конечных разностей требует обоснованного выбора

значения приращения аргумента х, поскольку задача численного дифференцирования является некорректной и ее необходимо решать методом регуляризации [1].

Известен способ автоматического вычисления вместе с функцией / (х0 ) при некотором

значении ее аргумента х = хо значения производной этой функции /'(хо ) в той же точке [2]. Данный метод основан на алгебре дуальных чисел [3] и носит название „автоматическое дифференцирование", подчеркивающее одновременность вычислений значений / (хо) и

/'(хо ) на основе только лишь исходного кода функции /(х) . К сожалению, ощущается нехватка отечественной литературы, посвященной данному методу. Настоящая статья призвана частично восполнить данный недостаток.

Дуальные числа. Рассмотрим множество чисел г, таких, что г = х + вА, где х е Я, Ае Я, а в ^ 0 — такая инфинитезимальная единица, что для нее выполняется точное равен-

2

ство s = 0 . Возможны математические операции сложения и умножения s с вещественными числами. Числа z называют дуальными. Их составляющую х принято называть действительной частью z и обозначать Re z , а А — инфинитезимальной частью z и обозначать Inf z .

Рассмотрим разложение в ряд Тейлора гладкой функции f (х) одного действительного аргумента:

f (х ) = f (Х0 ) + t 1 ^^ (х - х. у.

n=1n! dK

Выполним его аналитическое продолжение на множество дуальных чисел: заменим в выражении разложения f (х) в ряд действительный аргумент х на дуальный z = х + sA . Тогда значение функции

f (х + sA) = f (х. ) + £ 1 ^^ (х - х. +sA)n,

n=1 П! dx

согласно биному Ньютона, равно

f (х. ) + t 1dfl ZC (х - х. )n-m sm Am,

n=1П! dc m=0

что с учетом свойства s2 = 0 есть

f (х0) + t 1 ((х - х0 )n + nsA(х - х0 )n-1) =

f ( ) » 1 dnf (*0)( )n J + » 1 dmf (х0)( )m-1 ,( ) A df (х) f (х0) + >--(х-х0) + sAt--(х-х0) = f (х) + sA 7.

n=1 n! dхn v J m=1(m-1)! Ckm v V dr

Получаем, что результат вычислений является дуальным числом, одна составляющая которого есть точное значение функции f (х), а другая — точное значение производной этой

же функции при том же значении аргумента х.

Таким образом, возможно использовать дуальные числа для автоматического вычисления значения производной f '(х) точно вплоть до ошибки округления. Для того чтобы воспользоваться этим приемом, на практике необходимо создать новый тип данных. Поскольку

дуальные числа г = х + вА суть пары вещественных чисел сти их следующим образом.

( x >

vAy

то можно на языке С++ вве-

тип данных, описывающий дуальное число.

struct dual {

double real; // действительная часть.

double inf; } // инфинитезимальная часть.

Для того чтобы реализовать автоматическое дифференцирование функции f (х), представленной своим исходным кодом, необходимо перейти в программе от вычислений с действительными данными (double) к вычислениям с дуальными числами (dual).

Вычисления, производимые в любой программе математической обработки, являют собой вычисления значений сложных функций, например, для вычисления значения

y = f ( (g2 (...gn (х)...), x), x) необходимо начать с вычисления значения yn = gn (х), продолжить вычислениями yn-1 = gn-1 (yn, х) , yn-2 = gn-2 (yn-1, х) и т.д. Как известно,

производная сложной функции у = /( (^ (..gn (х) . ), х), х) по аргументу х вычисляется по цепному правилу, т.е.

^ (+1, х) = ^ (у,.+1, х) ^+1 (у.+2, х) + ^ (у{+1, х)

дх

ду.

+1

дх

дх

для . = 1, 2, ..., п. Следовательно, можно вычислять значение производной

д/ (х )

дх

в последо-

вательности вычисления значения

^п (х)

функции у = у ( ^ ( g2 (... gn (х)...), х), х) ,

т. е. начиная с

дх

, продолжая вычислением

Фп-1 (уп , х) = ^п-1 (уп , х) ^п ( х) + д^п-1 (уп , х)

дх

дуп

дх

дх

и т.д. вплоть до

(х1) = ^1 (у2, х) ^2 (х) + (у2, х)

дх1

ду2

дх

дх

Таким образом, можно выполнить автоматическое дифференцирование одновременно с вычислениями основной функции.

На практике математические преобразования, реализуемые компьютерными программами, сводятся к тому, что все функции gi, вообще говоря, являют собой последовательность

элементарных арифметических действий из набора {+, -, х, /} и элементарных функций —

математических примитивов из стандартных библиотек. Можно констатировать, что приведенное рассуждение о вычислении производной сложной функции справедливо и для функции, реализованной программой. Таким образом, для работы с дуальными числами необходимо задать основные арифметические действия с ними и указать, как производить вычисления элементарных математических функций при дуальном аргументе.

С математической точки зрения эта задача равносильна определению алгебры дуальных

чисел Л = ^Я2, +, х^ как соответствующей структуры над множеством пар действительных

чисел Я2 с операциями сложения и умножения.

Программно-ориентированное задание алгебры дуальных чисел. Рассмотрим множество пар действительных чисел Я 2 и зададим для двух произвольных его элементов

арифметические операции следующим образом.

г =

Г х ^

и г^ =

Г х2 ^

V у у

V у2 у

Сложение будем выполнять по правилу

хх

г1 + г 2 =

V л.

V у2.

х + х^

V у + у2 У

(1)

Операция обладает свойствами коммутативности и ассоциативности. Действительно, выполнены соответствующие соотношения г1 + г2 = г2 + г1 и (г1 + г2 ) + г3 = г1 +(г2 + гз).

Существует нейтральный элемент операции сложения 0 =

( 0 ^ V 0 у

справедливо равенство

г + 0 =

Г х Л Г0Л Г х Л

V у у

V 0 У

= г.

V у у

Вычитание — обратная сложению операция. При выводе правил для компонентов ее

результата требуется решить уравнение

х

V у у

Гх2 ]-

V у2 у

х

V уз

х

относительно

у

V уз у

. Получаем

I х — х2 + хз I хз — х х2

систему уравнений \ , ее решением является ■< . Следовательно,

[у = у2 + уЗ IуЗ = у - у2

г — г^ —

Г X Л

V у1 у

V у2 у

Умножение: введем операцию умножения как

г1 х г 2 =

V у1.

г

V у2 у

х1 — х^

V у1 — У2 у

х1 ' х2

V у ' х2 + у2 ' х1 у

(2)

(з)

Операция умножения также коммутативна и ассоциативна, поскольку г1 х г2 — г2 х г1 и

Г1 Л

(г1 х г2 ) х гз — г1 х (г2 х гз ) . Для нее существует нейтральный элемент 1 — . Действитель-

V 0 у

х

но, г х 1 —

V у у

Г1Л Г х '1 Л

V 0 у

у'1 + 0' х

х

V у у

— г.

Деление — обратная умножению операция. При выводе правил для компонентов ее ре-

зультата требуется решить уравнение тему уравнений

Г х1 Г х2 Л Г хз Л

V у у I у2 у I уз,

относительно

Г хз Л V уз у

. Получаем сис-

х1 — х2 хз [у — у2 хз + уз х2

х

х—

1

уз х2 — у1 — у2-

х2

х — ■

12

уз —

у1 х2 — у2 х1

В качестве вывода получаем, что при х2 ^ 0

(

(хл /гул

V у1 л

V у2 у

х1 х2

у1 х2 — у2 х1

v

х22

(4)

Имеет место также дистрибутивный закон: г1 х (г 2 + гз ) — г1 х г 2 + г1 х гз. Для таким образом введенной алгебры Л — ^Я2, +, х^ дуальных чисел справедливо следующее свойство:

Г 0 Л г 0 Л Г 0 Л

х

V у2 У

V у у

V 0 У

х

2

х

г

1

г

2

т.е. для квадрата чисто инфинитезимального числа г всегда имеет место свойство г2 = 0. Таким образом, из свойств арифметических операций (1)—(4) для инфинитезимальной единицы в вытекает ее свойство в2 = 0 .

В силу того что справедливы соотношения

Яе (г1 ± г2 ) = Яе ( г1 )± Яе ( г 2 ),

Яе (г1 х г 2 ) = Яе (г1 )Яе (г 2 ),

Яе (г1 / г 2 ) = Яе (г1 )/Яе (г 2 ),

алгебра Л описывает такое расширение множества действительных чисел на двумерную плоскость, которое включает в себя поле действительных чисел как составную часть. Действия с инфинитезимальными составляющими дуальных чисел не влияют на значения их действительных составляющих.

Результат вычисления элементарной функции g (z) от дуального аргумента z =

f g (x) ^

f x >

V У J

равен

g (z ) =

yg '(x )

что следует из рассмотренного выше аналитического продолжение ряда Тейлора для функции g (х). Указанное свойство задает правило вычисления элементарных функций от дуального аргумента, которое необходимо использовать для перегрузки математических примитивов.

процедура для вычислений значения тригонометрического синуса процедура для вычислений значения гиперболического синуса

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 (х) = sin(х)exp(-х2) + x в точке x0 . В соответствии с представленными выше правилами

произведем последовательно вычисления, заменив действительный аргумент х = Хо дуаль-

ным z0 =

f хо ^

V 1 J

f (z0 ) = sin

f

f хо ^

V 1 J

x exp

ГхЛ íx~\

V 1 J

V 1 J

íx,a f

V 1 J

sin

(хо)

x exp

f X2 ^

- x0 V-2 xo J

f xo ^

V 1 J

sin

cos

(x0) (x0)

Л

exP(-x0) -2x0 exp (-xO2)

f xo ^

V 1 J

f x0 ^

V 1 J

1•cos(x0 ) sin (x0 )exp (-xO) -2x0 sin(x0)exp(-x0) + cos(x0)exp(-x0) sin (x0 ) exp (-xÓj ) + x0 -2x0 sin(x0 )exp(-x2 ) + cos(x0 )exp(-x2) +1

Действительно, f'(x) = -2xexp(-x2)sin(x) + cos(x)exp(-x2) +1, что в точности совпадает с результатом вычислений при помощи дуальных чисел.

Если, например, х0 = 0, то получим z0 =

( 0 1 V 1,

и

sin z0 = sin

( 0 1 v 1,

( sin0 1 (01

V1 • cos0 у

exp

R)

= exp

(01 (01

— V 1У x V 1У = exp

( 01 V 0 у

V1У

( exp0 1 0 • exp0

(11 V 0 у

Следовательно,

f (z0 ) =

sin z0 x exp

H)

( 01 (11 ( 01

+ z0 =

V1У

V 0 У

V1У

(01 (01 +

v 1 у

V1У

( 01 V 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.

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 г.

технологий

УДК 004.045

К. Г. Коротков, Д. В. Орлов, Е. Н. Величко

ПРИМЕНЕНИЕ МЕТОДА ГАЗОРАЗРЯДНОЙ ВИЗУАЛИЗАЦИИ ДЛЯ АНАЛИЗА РАЗЛИЧНЫХ ЖИДКОСТЕЙ

Приведены результаты применения компьютерной программно-аппаратной системы для оценки свойств жидкостей, основанной на измерении и компьютерной обработке стимулированной электромагнитным полем электрофотонной эмиссии с поверхности жидкости, на базе хорошо известного метода газоразрядной визуализации.

Ключевые слова: структура воды, газоразрядная визуализация, анализ крови.

Введение. В последние годы большое внимание уделяется изучению структурных свойств воды и возможности наделения воды информацией. Имеется достаточно много противоречивой информации о „памяти воды" [1]. Согласно одной из квантово-механических моделей, наблюдаемые экспериментально явления обусловлены процессами формирования кластеров и клатратов преимущественно на атомах примесей [2]. Для введения этих понятий в контекст современного научного мышления прежде всего необходим набор доказательных и воспроизводимых экспериментальных фактов. Вода — сложный объект исследований, ее свойства зависят от большого числа факторов, поэтому необходимо параллельно использовать несколько независимых методик, а также требуется разрабатывать новые информативные методы исследования свойств воды.

Экспериментальные данные, полученные группой ученых под руководством А. Нильсо-на, ведущего специалиста Стенфордского центра синхротронного излучения (Stanford Synchrotron Radiation Lightsource), позволили сделать важный шаг в направлении разгадки 66 уникальных свойств воды благодаря новейшим методам изучения строения жидкостей с использованием мощного рентгеновского излучения, получаемого с помощью больших ускорителей элементарных частиц — синхротронов.

Ученые выяснили, что молекулы воды одновременно формируют не один, а два типа структур, сосуществующих вне зависимости от температуры. Один тип представляет собой формирования примерно по 100 молекул и напоминает структуру льда. Структура второго типа жидкости гораздо менее упорядочена. Увеличение температуры вплоть до точки кипения воды приводит к некоторому искажению структуры молекулярных формирований, уменьшению их количества и доминированию разупорядоченной структуры [3].