CC BY
26
Реализация 3-х проходного метода быстрого автоматического дифференцирования
А.Ю. Горчаков
Аннотация— В статье рассматривается реализация 3-х проходного метода быстрого автоматического дифференцирования. Добавление к стандартному 2-х проходному методу еще одного прохода приводит к существенному уменьшению потребления оперативной памяти, при незначительном (25%) увеличении количества вычислений. Приведены оценки производительности и требований к объему оперативной памяти 2-х и 3-х проходного методов. Проведен численный эксперимент на примере обратной задачи восстановления начальных данных для уравнения переноса.
Ключевые слова— оптимальные управление, быстрое автоматическое дифференцирование, уравнение переноса, обратная задача восстановления начальных данных.
управления на сетках большой размерности.
В данной статье рассматривается методика уменьшения используемого объема оперативной памяти, за счет незначительного увеличения объема вычислений. Методика реализована автором в модифицированной версии пакета быстрого автоматического дифференцирования Adept [8]. Рассмотрение проводится на примере обратной задачи восстановления начальных данных для уравнения переноса. Необходимо отметить, что идея уменьшения требований к объему оперативной памяти не является абсолютно новой, в частности похожий метод (называемый revolve) можно найти в пакете ADOL-C [9].
I. Введение
При численном решении задач оптимизации часто используются градиентные методы минимизации функций, которые требуют умения вычислять не только значения оптимизируемых функций, но и их градиентов. При этом крайне важно вычислять точные значения градиентов за минимально возможное время.
Вычисление градиентов методом конечных разностей иногда связано с большими трудностями, и занимает достаточно много времени в случае задач большой размерности.
Более перспективным представляется использование методологии быстрого автоматического
дифференцирования (обобщенная БАД-методология, см. [1]-[4]). Но, при очевидных преимуществах БАД-методологии, таких как, точное вычисление градиента, независимость времени расчета от размерности вектора управлений, имеется и определенный недостаток, а именно - высокие требования к объему оперативной памяти. Указанный недостаток начинает широко проявляться при решении задач оптимального
Работа выполнена при выполнена при финансовой поддержке РФФИ, проект 16-07-00458 A, программы Президиума РАН I.5 "Интеллектуальные информационные технологии и системы". А.Ю. Горчаков - ведущий математик Вычислительного центра им. А.А. Дородницина Федерального исследовательского центра «Информатика и управление» Российской академии наук. andrgor 12@gmail. com
II. МЕТОДОЛОГИЯ БЫСТРОГО АВТОМАТИЧЕСКОГО ДИФФЕРЕНЦИРОВАНИЯ
Пусть г Е Яп и и Е Яг векторы. Дифференцируемые функции и) и Ф(2, и) определяют отображения
Яп х Яг ^ Я1, Ф: Яп х Яг ^ Яг. Вектора г и и удовлетворяют следующей системе из п нелинейных скалярных уравнений Ф(г, и) = 0. Если матрица Ф^ (г, и) невырожденна, то сложная функция ^(и) = и)
дифференцируема и ее градиент относительно переменных и вычисляется по формуле: йП
— = Ши^(и),и) + ФТиШ), и) • р, (1)
где вектор р Е Яп находится из решения линейной алгебраической системы:
Щ!(г, и) + Ф1 (г, и) •р = 0, (2)
Здесь и далее индекс Т обозначает транспонирование, нижние индексы 7, и обозначают частные производные
д№ д№ т функций по векторам г и и: Ши = — , И^ = —, Ф'и =
дФТ дФт _ .
——, Ф1 = ——, так же будем обозначать 1-е, 1-е
ди дг
компоненты векторов г и и как г^, , и¿, щ.
Предположим, что каждый компонент вектора последовательно выражается только через компоненты вектора и и предыдущие г^, т.е. 1</</., тогда формулы (1)-(2) можно записать в следующем виде:
п
VI = и)Р] + И^(г,и) (3)
]=1+1
(4)
dD. sr „
du = Wu. (z, u) + Ф1.(z, u) • P]
1 ]=i
Системы (3)-(4) называют явными [1], и возникают они, например, в задачах оптимального управления, если интегрирование проводится по явным схемам и по некоторым неявным схемам.
Для примера рассмотрим подробно дифференцирование скалярной функции двух переменных:
f(ui, U2) = щ sin(u1) + Ui*U2, (5)
в точке и1 = 0.3, и2 = 0.5.
Вычисление значения этой функции представляем, как последовательность формул, содержащих
арифметические действия и основные элементарные функции:
i = 1 Zl = sin(u1) sin(0.3) = 0.296
i = 2 z2 = u2*z1 0.5 * 0.296 = 0.148
i = 3 z3 = u1*u2 0.3 * 0.5 = 0.15
i = 4 Z4 = z2 + z3 = f(u1, u2) 0.148 + 0.15 = 0.298
Переходя от одной формулы к другой в сторону увеличения индекса i, вычисляем значение функции. Далее в обратном порядке определяем компоненты вектора p и все частные производные функции по независимым переменным u.
df df
Р4 =L°, Рз = Р2 = Р1 = ^ = ^ =0
i = 4 Рз = Рз + Р4 = 1.0
Р2 = Р2 + Р4 = 1.0 df df
i = 3 —- = —— + и2*р3 =0.5 * 1.0 = 0.5
ди1 ди1
df df
= -L- + и1*р3 = 0.3 * 1.0 = 0.3
ди2 ди2 1 нз
df df i = 2 — = — + Z1*P2 = 0.596 ди2 ди2
p1 = p1 + u2*p2 = 0.5 * 1.0 = 0.5 df df
i = 1 —- = —- + cos(u1) *p1 =0.978
ди1 du1
Можно представить процесс вычислений в виде информационного графа (см. рис.1), в таком случае говорят, что сначала расчеты ведутся слева-направо, а затем справа-налево.
Рис.1 графическое представление процесса вычислений В работе [10] дана оценка временной сложности вычисления градиента Тд/Т0 < 3, где Т0 полное время, требуемое для вычисления значения функции, Тд дополнительное время, требуемое для вычисления всех частных производных функции; оценка дана без учета времени, необходимого для работы с памятью компьютера. Напомним, что в случае численного дифференцирования для расчета градиента потребуется
времени не менее (1 + п)Т0, где n число необходимых нам частных производных.
Необходимо отметить, то что при вычислениях требуется сохранять в памяти компьютера значения компонентов векторов z и р. Для многих задач оптимального управления, решаемых конечно-разностными методами, это означает, что необходимый объем памяти прямо пропорционален количеству узлов расчетной сетки. Например, для расчетной сетки с 4.8 * 1011 узлами и вычислениями с двойной точностью (8 байт на число) потребовалось бы порядка 7 терабайт памяти.
III. Постановка задачи восстановления начальных
ДАННЫХ
Рассмотрим задачу восстановления начальных данных для одномерного уравнения переноса пассивной примеси в жидкости движущейся с постоянным значением скорости: дф дц>
— + с-^ =0, с = 1, (х, О EQ, dt ах
ц>(х, 0) = w1(x), р(0,0 = W2 (О,
х £ [0,1], t £ (0,1].
(6)
(7)
(8)
Здесь х - пространственная координата; t - время; с -скорость движения жидкости; ф(х^) - концентрация примеси в точке x в момент времени ^ прямоугольник Q = (0,1)х(0,1); w2(x) - заданная функция, определяющая концентрацию примеси в точке х=0; ш1(х) - искомая функция, определяющая концентрацию примеси в начальный момент времени. Пусть на прямой х+t=1 задана функция ср (х, £), где (х, Ь 6 Q). Ее можно рассматривать как концентрацию примеси, полученную в результате некоторых экспериментальных исследований. Поставим следующую задачу: требуется определить оптимальное управление и,(х) = м^Дх) и соответствующее решение (х, £) системы (6)-(8), при котором интегральный функционал W(u) = 1/01[^(х, 1 -х) - ф(х, 1 - х)]2 ах (9)
достигал бы минимально возможного значения. Задачи, подобные этой, обычно решаются численно с помощью некоторого метода спуска, который требует знания градиента функционала (9).
Перейдем к дискретному варианту системы (6)-(8). Для этого покроем прямоугольник Q регулярной расчетной сеткой с шагами к=13 и т=1/Ы, где N -количество интервалов по времени, а 3 - количество интервалов по оси х (для удобства расчетов возьмем N кратным 3, т.е. N=^3). Используем хорошо известную схему «уголок», аппроксимирующую уравнение (6) на трехточечном шаблоне с первым порядком точности:
<1-П + c<Pj-(Pj-i = 0 ,
= Wlj = Uj,
Ф0 = W2n,
h
(10) (11) (12)
W(u)=2 ^(VÏ^-VÏ-^ )2, (13)
i=i
где tf = (p(xj, tn), ф7- = <p(xj, tn),Xj = jh, tn = пт, = Wij = Uj = u(xj),
= = Ш2 (¡п), 0 К] <],0 <п < N.
Требуется найти вектор и, при котором сложная дискретизированная функция Ш(и) принимает минимальное значение.
1У.ОПИСАНИЕ МЕТОДИКИ ВЫЧИСЛЕНИЯ ГРАДИЕНТА 3-Х ПРОХОДНЫМ МЕТОДОМ
С целью снижения требований к объему памяти, предлагается использовать 3-х проходный метод автоматического дифференцирования. Идея метода основана на том, что для «явных» по времени разностных схем вычисления можно остановить, сохранив один временной слой Фп - значения ^",1 < ]' < } (или несколько временных слоев, в зависимости от используемой разностной схемы), а затем продолжить их в любой момент используя сохраненные значения. Иначе говоря, для таких разностных схем можно делать «контрольные точки». Используя эту особенность, разобьем расчетную сетку по оси времени на М-групп, с К=Ы/М временных слоев в каждой группе и будем хранить в памяти только последний временной слой в каждой группе и временные слои, относящиеся к текущей группе. Приведем описание метода.
МЕТОД 3Х
1. Разбить расчетную сетку по оси времени на М групп, с К=Ы/М временных слоев в группе.
2. В цикле по п от 1 до (М-1)*К выполнить
a. вычислить Фп,
b. для п =К,2К,...,(М-2)*Ксохранить Фп в памяти
3. Для группы М вычислить Фп и частные
производные
dW
>дф
(М-1)К i
-1 <]<]
4.
В цикле по m от M-1 до 1 выполнить a. вычислить Фп,
b. вычислить
dW,
'дер)
(m-i)K
-1 <]<]
5. как градиент функционала W(u) вернуть
dW,
о ,1 <]<]
/dtf
Графически этот метод можно представить следующим образом:
ООО*
ООО*
I
• ••О
(
Здесь черные стрелки представляют собой классические проходы слева-направо и справа-налево, а белые стрелки это дополнительный проход, при котором вычисляются и запоминаются промежуточные временные слои Фп. Белые шарики - временные слои рассчитанные, но не хранящиеся в памяти; черные -рассчитанные и сохраненные; серые - рассчитанные ранее.
Легко показать, что временная сложность вычисления градиента для 3-х проходного метода Тд/Т0 < 4. При этом в оперативной памяти необходимо хранить значения <р™ в J*(M-2+K) и значения компонентов вектора p в J*K узлах расчетной сетки. При оптимальном выборе M требования по памяти уменьшаются с J*N до ] * 2 * V2 * N. Для расчетной сетки 16000x48000 вместо 11.4 Гб потребуется 75.6 Мб оперативной памяти.
Метод 3Х реализован как дополнительный программный модуль в пакете быстрого автоматического дифференцирования Adept функцией класса Stack:
compute_adjoint3X(double *x, double *g, int N, int M, adouble (*F_k)(const adouble *x, const int n1, const int n2, adouble *PHI) ), где x - вектор независимых переменных, g - вектор частных производных функции, N - количество интервалов по времени, M - количество групп, F_k - ссылка на пользовательскую функцию, которая принимая как параметры вектор x, n1 и n2 индексы начала и конца группы, возвращает значение функционала и PHI - вектор Фп.
V. Численный эксперимент
Задача (6)-(9) решалась при начальном условии:
fx-0,2\2 fx-0.4 \
w-tix) = е f 0.07 ) + е f о .о7 )
и w2 (t) = 0
(14)
Разностная схема и вычисление функционала были записаны на языке C и продифференцированы с помощью программного пакета Adept [8]. Вычисления производились на персональном компьютере с процессором Intel Core i7-4770, 3.4 GHz, 8Gb оперативной памяти. Значение функционала W(u) и его частных производных вычислялось в точке
u(t) = ср^), (15)
на трех расчетных сетках различной размерности, стандартным методом и методом 3X. Графики wi(х), и(х), приведены на рис.3
Рис.2 графическое представление метода 3X.
12
А А
А А
/ W \
0 2 0 / v \
J V
0 2 0.6 1 0 0.4 0.8
Рис. 3 Графики функций (14) - синяя линия и (15) -красная линия
Значение M=N/600 (т.е. 600 временных слоев в группе). Значения частных производных функционала вычисленных обоими методами совпали друг с другом с машинной точностью. Для расчетной сетки 1 оба метода используют только оперативную память, при этом метод 3X работает немного медленнее чем стандартный метод. Для расчетной сетки 2 стандартный метод БАД начинает использовать более медленную виртуальную память компьютера, соответственно наблюдается уменьшение скорости его работы. При переходе к расчетной сетке 3 скорость работы стандартного метода БАД уже значительно меньше скорости работы метода 3X.
Необходимо отметить, что программный пакет Adept сохраняет в оперативной памяти не только значения в узлах расчетной сетки, но и дополнительные внутренние переменные, поэтому объем используемой оперативной памяти выше расчетной.
VI. ЗАКЛЮЧЕНИЕ
В работе рассмотрена реализация 3-х проходного метода быстрого автоматического дифференцирования для задач оптимального управления. Метод реализован
как дополнительный программный модуль в пакете быстрого автоматического дифференцирования Adept. По сравнению со стандартным методом БАД существенно снижены требования к оперативной памяти компьютера (вместо линейной зависимости объема оперативной от количества интервалов по времени N, получена зависимость 2V2W), при незначительном (на 25%) увеличении объема вычислений. Особенно актуально использование метода при расчетах на GPU (графические ускорители), так как в отличие от CPU, которому можно предоставить 1Тб и более оперативной памяти, для видеокарт и профессиональных графических ускорителей существует ограничение 16/32Гб.
Метод применим для задач оптимального управления любой размерности по пространственным координатам. Единственное требование, предъявляемое к разностным схемам - возможность делать «контрольные точки». Метод совместим с большинством современных пакетов автоматического дифференцирования. В случае необходимости дополнительного снижения требований к оперативной памяти можно, добавляя вспомогательные проходы, получать методы 4X, 5X и т.д.
VII. БИБЛИОГРАФИЯ
[1] Айда-Заде К.Р., Евтушенко Ю.Г. Быстрое автоматическое дифференцирование // Математическое моделирование, 1989. V. 1, pp. 121-139.
[2] Евтушенко Ю.Г., Мазурик В.П. Математическое обеспечение оптимизации. М: Знание. 1989.
[3] Evtushenko Y.G. Automatic differentiation viewed from optimal control theory. Proceedings of the Workshop on Automatic Differentiation of Algorithms: Theory, Implementation and Application. Philadelphia: SIAM, 1991.
[4] Griewank A, Corliss G.F., eds. Automatic Differentiation of Algoritms. Theory, Implementation and Application. Philadelphia: SIAM, 1991.
[5] Горчаков А.Ю., Посыпкин М.А. Эффективность методов локального поиска в задаче минимизации энергии плоского кристалла. Международный научный журнал «современные информационные технологии и ИТ-образование». - 2017. - V. 13. - N. 2. - pp. 97-102.
[6] Евтушенко, Ю. Г., Лурье, С. А., Посыпкин, М. А., Соляев, Ю. О. (2016). Применение методов оптимизации для поиска равновесных состояний двумерных кристаллов. Журнал вычислительной математики и математической физики, 56(12), 2032-2041.
[7] Turkin A., Thu A. Benchmarking Python Tools for Automatic Differentiation. International Journal of Open Information Technologies. -2016. - V. 4. - N. 9.
[8] Hogan R. J. Fast reverse-mode automatic differentiation using expression templates in C++. ACM Transactions on Mathematical Software (TOMS). - 2014. - V. 40. - N. 4. - pp. 26.
[9] Griewank A., Walther A. Algorithm 799: revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation. ACM Transactions on Mathematical Software (TOMS). 2000. V. 26. N. 1. pp. 19-45.
[10] Евтушенко Ю.Г. Оптимизация и быстрое автоматическое дифференцирование. Научное издание ВЦ РАН, 2013.
Таблица 1. Сравнение стандартного метода БАД и 3X
№№ Расчетная Использованная Время
сетка (JxN) память расчетов
стандарт/3Х стандарт/ЗХ
(Мб) (сек)
1 4000x12000 2197.8/110.8 4.1/5.8
2 8000x24000 8790.5/222.8 53.3/29.6
3 16000x48000 35158.4/450.6 1429.2/127.0
The implementation of a 3-pass method for fast automatic differentiation
Andrei Y. Gorchakov
Abstract— in this paper, we consider the implementation of the three-pass method of fast automatic differentiation. Adding one more pass to the standard 2-pass method leads to a significant decrease in the consumption of RAM, with an insignificant (25%) increase in the number of calculations. Estimates of performance and requirements for the amount of RAM of 2 and 3-pass methods are given. A numerical experiment is performed on the example of the inverse problem of reconstructing the initial data for the transport equation.
Keywords— optimal control, fast automatic differentiation, transport equation, inverse problem of recovering initial data.
