CC BY
20
УДК 519.688
А. В. Гасников1'2'3'8, П. Е. Двуреченский1'2'4, Д. И. Камзолов1'2, Ю. Е. Нестеров3'5, В. Г. Спокойный1'2'3'4, П. И. Стецюк6, А. Л. Суворикова1'2'7,
А. В. Чернов8
Лаборатория структурных методов анализа данных в предсказательном моделировании при
МФТИ(ГУ) (ПреМоЛаб) 2Институт проблем передачи информации РАН 3Национальный исследовательский университет Высшая школа экономики 4Weierstrass Institute for Applied Analysis and Stochastics, Berlin 5Universite catholique de Louvain, CORE, Belgium 6Институт кибернетики им. В.М. Глушкова НАН Украины 7The International Research Training Group 1792, Berlin 8Московский физико-технический институт (государственный университет)
Поиск равновесий в многостадийных транспортных
моделях
В работе предлагается оргинальный способ поиска равновесий в многостадийных моделях транспортных потоков. В основе подхода лежит сочетание метода Синхорна и универсального градиентного метода.
Ключевые слова: седловая задача, энтропия, метод Синхорна, универсальный метод, неточный оракул.
A. V.Gasnikov1'2'3'8 P. E.Dvurechensky1'2'4 D. I.Kamzolov1'2 Y. E.Nesterov3'5 V. G.Spokoiny1'2'3'4 P. I.Stetsyuk6 A. L.Suvorikova1'2'7 A. V.Chernov8
1 Research laboratory in Predictive Modeling and Optimization at PhysTech (PreMoLab) 2Institute for Information transmission problems RAS 3Higher School of Economics National Research University 4Weierstrass Institute for Applied Analysis and Stochastics, Berlin 5Universite catholique de Louvain, CORE, Belgium 6V.M. Glushkov Institute of Cybernetics of NAS of Ukraine 7The International Research Training Group 1792, Berlin 8Moscow Institute of Physics and Technology (State University)
Searching for equilibriums in multistage transport models
We propose an original way to find an equilibria in multistage traffic flow models.The basis of our approach is the combination of the balancing method and the universal gradient method.
Key words: saddle problem, entropy, balancing method, Synchorn method, universal method, inexact oracle.
1. Введение
Поиск (стохастических) равновесий в многостадийных моделях транспортных потоков приводит к решению следующей седловой задачи с правильной (выпукло-вогнутой) структурой [1,2]:
{п п
У^ г^ 1п г^ + ^ с^ (х) г^ + д (х) 1,3=1 г,3=1
3 = 1 г=1
г,]=1,...,п
+ 9 (яП , (1)
п
max max < (A, L) + (ß, W) — / exp( —c^ (x) — 1 + \ + ßj )
hj=1
где д (х) и сц (х) — вогнутые гладкие функции (если ищутся не стохастические равновесия, то Су (х) могут быть негладкими), Q С Кт - выпуклое ограниченное замкнутое множество простой структуры. Легко понять, что система балансовых ограничений в (1)
п п
либо несовместна ^ Ь = ^ , либо вырождена (имеет не полный ранг). В последнем
г=1 ]=1
случае это приводит к тому, что двойственные переменные ( Л, ¡) определены с точностью до произвольной постоянной а:
( Л + ае,д — ае), е = (1,..., 1) е Мп
п п
Без ограничения общности положим ^ Ь = ^ = 1. Тогда задачу (1) можно перепи-
г=1 ]=1
сать следующим образом:
{п п
У^ 1п г^ + сгз (х) + 9 (х) } =
г, 3=1 г,з=1
3 = 1 г=1
~^0,г,]=1,...,п
= тах тах < \ Л.,Ь) + {¡1, Ш) — У^ ехщ — с^ (х) — 1 + Лг + ¡А + д(х)
{п
{Л, Ь) + {¡,Ш) + 1 — ехр(—Сц (х) + Лг + ¡^ + д(х)
, =1
= — т1п( / (х) — 9(х)), (2)
где выпуклая функция ( х) определяется следующим образом:
{п п
У] 1п + сч (х) г.
, =1 , =1
3 = 1 г=1
~^0,1,]=1,...,щ
=—ттах | {л, Ь+{¡,Ш )+1 — Е ехр(—(х)+Лг+^. (3)
Кроме того, задачу (1) можно представить так:
п
1п + с^ (х) г^ + д (х)
п п
п т1Г1 гпах < ^ 1п + Ы (х)
=ь1, Ё ^з = №з \г,]=1 г,3=1
3=1 г=1
п,п
гц>0,г^=1,...,п; £ х^ = 1
1,3 = 1
= тах тах < { Л, Ь) + {¡, Ш) — 1п I ехр( — с^ (х) + Лг + ¡Л ) + д(х)
— Й>(^(х) — д(х)), (4)
где выпуклая функция / (х) определяется как
п п
! (х) = - п шт I 1п + ^ (х)
оу ш ¿-гу ^ / у ^г]
£ г^=L^,YJzij = \1,3=1 1,3=1
j=l ¿=1
П,П
гц^0,г^=1,...,п; £ г^ = 1
1,3 = 1
шах I- 1п ( ^ ехр^-с^ (х) + А* + ^) | + (Л, Ь) + (^, Ш) \ . (5)
\1,3 = 1
Поскольку мы добавили в ограничения условие ^ = 1, являющееся следствием балан-
¿,3 = 1
совых уравнений, то это привело к тому, что двойственные переменные (Л, ^) определены с точностью до двух произвольных постоянных С\, С(Л + С\е, ц, + С^е).
В данной работе показано, как можно решать задачи (2) и (4) с помощью методов первого порядка. Мы ограничимся изучением только полноградиентных методов для этих задач, т.е. не будем рассматривать, например, рандомизацию при вычислении градиента. Планируется отдельно исследовать вопрос о возможности ускорения вычислений за счет введения рандомизации для внешней задачи. Мы предполагаем, что это может уменьшить общую трудоемкость метода только в случае, когда вспомогательная задача вычисления Ус^ (х) тоже достаточно сложная. Следует отметить, что для большинства приложений умение вычислять Усу (х) для конкретной пары (г^) позволяет одновременно рассчитать и все Усу (х), ] = 1,... ,п, без дополнительных затрат. Кроме того, планируется отдельно исследовать вопросы о том, какие подходы и в какой степени могут быть распараллелены. Оказывается, что целесообразность рандомизации тесно связанна с именно возможностью использовать параллельные вычисления.
Структура статьи следующая. В п. 2 рассматриваются алгоритмы приближенного вычисления субградиента функций / (х) и / (х). В п. 3 мы приводим алгоритм для решения задач (2) и (4). В основе оригинального подхода п. 3 лежит сочетание метода Синхорна для решения внутренней задачи оптимизации по двойственным множителям и универсального метода с неточным оракулом для внешней задачи.
2. Приближенное вычисление субградиента целевой функции
Заметим, что расчет градиента У/ (х) (в ряде транспортных приложений [1] вогнутые функции Су (х) - негладкие, тогда вместо градиентов стоит понимать суперградиенты Су (х) и субградиент f (х)) осуществляется по следующей формуле (Демьянова-Данскина-Рубинова, см., например, [1,3]):
п
У/(х) = - Е <3 (Х)УСгЗ (х) , (6)
1,3 = 1
где г*, (х) - единственное решение задачи (3). Аналогично
У/(х) = - Е 4 (х)Усгз (х), (7)
1,3 = 1
где г*^ (х) - единственное решение задачи (5).
{2 п,п 1 г € М™ : Хц = 1 > - единственное решение
г,3 = 1 )
ззадачи (3), тогда справедливо неравенство
Ц** (Х1) - г* (х2)\\1 < Ус (Х1) - с . (8)
Пусть г*(х) € Хп2 (1) = < г € М™ : = 1 ? - единственное решение задачи (5), тогда
{ г,3=1 )
справедливо неравенство (8).
п
Доказательство. Обозначим й(х) = ^ 1п х^. Из условий оптимальности в задаче (3) получаем
(V й{х(хг)) + с(хг), х(х2) — х(х\)) ^ 0, (V й(г(х2)) + с(х2), г(х\) — г(х2)) ^ 0.
Складывая эти неравенства и используя сильную выпуклость энтропии й(х) в 1-норме с константой 1 на симплексе Хп2 (1), получаем
Ус (х1) — с (х2)||те ||г(зс2) — г(х1)\\-1 ^ ^ (с(х\) — с(х2), г(х2) — г(х\)) ^ ^ (V й(г(хг)) — Vй(г(х2)), г(хг) — >
> ||г(х2) — %(х\) ||2 .
Отсюда получаем неравенство (8). Второе утверждение леммы доказывается аналогично. И
Обозначим с(х) = {сц(х))п^ € М™. Тогда V / (х) = -г* (х), где
х*(х) — единственное решение задачи (3), это матрица йс(х)/йхТ, состоящая из строк (Vс^(х))Т, г,] = 1, ...,п. Пусть в пространстве векторов х выбрана некоторая норма || ■ ||х. Через || ■ ||х,* обозначим сопряженную к ней. С учетом леммы 1, выполнено
(хг) — V/ ЫЩ* =
(йс(хг)\Т (<с(х2)\Т *. ,
х(хг) Чг(х2)
<
<
Т
х") (**(хг) —г *(х2))
+
<
( йс (хг)\
Т
|| Z*(хl) — г* (х2)Цг +
( /йс(хг)\Т («с(х2)ХГ\ «хТ ) V <хТ ) )
I / йс (хг)\Т (<Кс(х2) \
*( х2)
<
Т)
*( х2)
<
<
/ йс (хг)\
Т
Ус (хг) — с(х2)+
г^-х,*
(
( йс (хг)\Т (йс (х2)\Т\
V)
*( х2)
Здесь || ■ ||г^х,* обозначает норму линейного оператора, согласованную с 1-нормой в пространстве прообразов и нормой || ■ ||х,* в пространстве образов. Для любого х € Q выполнено включение г*(х) € Хп2(1), т.е. каждая компонента вектора г*(х) принадлежит отрезку [0,1]. Значит, если каждая компонента вектора с(х) зависит от х Липшицевым образом и ее субградиент ограничен, а также множество Q ограничено, то субградиент функции /(х) имеет ограниченную вариацию. Если субградиент функции — д(х) также имеет ограниченную вариацию субградиента, то вся целевая функция задачи (2) имеет ограниченную вариацию субградиента. Если дополнительно известно, что градиент каждой компоненты вектора с(х) и градиент функции — д(х) удовлетворяют условию Липшица, то градиент функции /(х) также будет удовлетворять условию Липшица. Однако полученная оценка константы, ограничивающей вариацию субградиента, и константы Липшица градиента получатся слишком завышенными. Тем не менее ниже предложен метод, которому не требуются оценки этих констант. Достаточно знать, что хотя бы одна из них является конечной. Аналогичные рассуждения приводят к тем же заключениям относительно субградиента функции ¡(х).
х,
х,
х,
х,
Опишем теперь алгоритм приближенного вычисления субградиента и значения функции f(x). Проводя рассуждения, аналогичные лемме 2 работы [4] с использованием теоремы Синхорна (1967), получаем, что существуют единственные с точностью до умножения на константу векторы u,v £ R+ такие, что для матрицы Z*(x) £ Rraxra, являющейся решением задачи (3), справедливо равенство
Z *(x) = diag (u)exp(-C (x))diag (v),
где diag(u), diag(v) - матрицы с диагональными элементами, равными соответствующим элементам векторов u, v, exp(-C(x)) означает поэлементную экспоненту от матрицы - C( x). При этом векторы u, могут быть вычислены с помощью итерационного процесса (метод Синхорна [4,5]):
1) U0 = exp(-C(x))e , V° = (exp(-C(x)))1e ,
2) Uk+l = exp(-C(x))vk , Vk+1 = (exp(-C(C)))1«fc+1 '
Здесь используется покомпонентное деление вектора на вектор.
Отметим, что каждая итерация этого процесса требует 0(п2) арифметических операций (а.о.). Пусть с помощью этого процесса получены приближения un, vn. Определим приближенное решение задачи (3) следующим образом:
Z(x) = diag (un _i)exp(—C (x))diag ( vN-1).
Следуя работе [6], будем использовать проективную метрику p(■, ■) для пары векторов a,b £ R++:
p(a, b) = ln (max —^ ) .
V ik akbij
Для матрицы A £ R++ra определим числа
о/лл aikaji (sJ-d(A) — Л Ц(А) = max-— и A) = V—^- .
v ' i,j,k,i a-jkau v ' ^ л/ЦА) + 1)
Тогда, согласно оценке (3.12) работы [6], если верно
1
х
1 — j(exp(—C (x)))
p^diag(un-1)exp(—C(x))vN-1,Lj + p^diag(un)exp(—C(x))vN-1,W^ ^ 5, (9)
то выполняется
е**^6) ^ Щ) ^ е**®'
ехр(-6) < < ехр(5) , з = I, ...,п. (10)
Здесь г* - элементы матрицы Z*(х), ¿^(х) - элементы матрицы ¿(х). Кроме того: р(^Иад(ь^-1)ехр(-С(х))vN+ р(^Иад(ь^)ехр(-С(х))vN^
< (-! (exp(—C (x))))N ^ p(exp(—C (xx))e,L) + p^diag (uo)exp(-C (x))e,W то есть метод имеет геометрическую скорость сходимости.
Из неравенства (10), так как г*(х) € [0,1] для всех г= 1, ...,п, следует, что |х*(х) — г^(х)1 ^ г*(х)(ехр(5) — 1) ^ 25 Уг= 1, ...,п,
при 5 ^ 1.
Пусть получено приближенное решение задачи (3):
¿(х) = йгад (им-г)ехр(—С (х))Мад (-г)
такое, что справедливо неравенство (9). Определим приближенное значение функции /(х) и приближенный градиент этой функции следующим образом:
Тогда
/ (х) = ¿з(х) 1п ¿г з(х) (х) ¿iз(х),
, =г , =г
п
^/(х) = — ^ ¿3 (х^(х).
, =г
| ¿(х) — ¡(х)1 =
п п п п
^2 ¿И(х) 1п¿И(х) — ^ (х) ¿3(х) г*з(х) 1пг*з(х) °Аз (х) г*з(х)
, =г , =г , =г
, =г
<
<
^2 (г*з(х) — ¿г.}(х)) 1пЧ¿(х)
, =г
+
^2 г*з(х)(1п г*з(х) — 1п ¿г.з(х))
, =г
+
+
^2 сгз(х) (
г* ( х) — ¿г ( х))
, =г
<
< 25 ^2 11п*гз(х)1 +п25 + 25 ^ |сц (х)1 (11)
, =г , =г
Эта оценка дает критерий остановки метода Синхорна для получения заданной точности вычисления значения функции ( х). Аналогично можно получить оценку
п п
|¿(х) — ¡(х)1 < 25 ^ 11пг*(х)1 + п25 + 25 ^ |Сг3 (х)|.
, =г , =г
Из этой оценки можно получить оценку на общую трудоемкость вычисления значения функции ( х) с заданной точностью. Также имеем
V I (х) — V } (х)
<
^ Vсг ](х) ¿г з(х) + ^ V°гз(х)
г* ( х)
, =г , =г
п
25 £ 1Сгз (х)Их',„
, =г
<
где || ■ ||х',* означает выбранную норму, которая может не совпадать с нормой || ■ ||х,*, выбранной для измерения скорости изменения субградиента функции /(х).
Теперь мы можем построить алгоритм, позволяющий приближенно вычислить значение f (х) функции ¡'(х) и ее субградиент V f (х), удовлетворяющие условиям
Ц~(х) — Дх)| < ^ и V ¡(х) — V / (х)
< ¿2.
х ,
х ,
х ,
Алгоритм 1. «Вычисление значения и производной целевой функции»
Data: Матрица C(х), векторы L,W е Хп(1), точности , 52 > 0 Result: / (х) ,Vf (х)
1) Вычислим матрицу А ^ exp(-С(х)).
2) Вычислим число j(А).
п п
3) Вычислим Ci ^52 Iсч (х)\,С2 ^52 \Wcij (х)Уж/,*.
i, 3=1 i, 3=1
4) Вычислим Uo ^ —, V0 ^ —^.
5) Вычислим a ^ р(Ае, L) + p(diag(u0)Ае, W).
6) foreach k ^ 0 do
а) Вычислим вектор а ^ Аик.
б) Вычислим ик+х ^ %, ук+1 ^ Ат^к+1.
в) Вычислим матрицу Z ^ сИад(и,к)А(1шд(Ук).
г) Вычислим число Р ^ 1-1{А) (р(^гад(ик)а, V) + р(сИад(ик+\)а, Ш)).
if ¿i ^ ¡3\ Е I %з(х)\ + п2 + 2Ci\ и 62 ^ 2f3C2 then
п п п
¡(х) <--52 zij ln zij - Е cij (х) Zij, и V f (х) <--52 zijV(х)
i, 3 = 1 i, 3 = 1 i, 3 = 1
end
else
k ^k + 1
Переходим к шагу a) end end
Algorithm 1: Вычисление значения и производной целевой функции Лемма 2. Общая трудоемкость приведенного алгоритма составляет
О (п4 + тп2 + n2N) ,
где
N = max <
1 -j(A)
ln
a E | ln4(х)| +n2 + 2C1 \i,i=1
§1(1 - j(A))
ln
2aC2
1 - 1(A) 82(1 -j(A))
Доказательство. Из геометрической скорости алгоритма следует, что после N шагов 6 справедливо неравенство
Р " 1 - 1(А) •
Это означает, что после N итераций неравенства, проверяемые на шаге 6г, выполнятся. Шаги 1, 4, 5, 6а, 6б, 6в, 6г требуют 0(п2) а.о. Шаг 2 требует 0(п4) а.о., но может быть
1
1
параллелизован. Шаг 3 требует 0(тп2) а.о., где т - размерность вектора х. Наконец, последний шаг вычисления результата работы алгоритма требует 0(тП) а.о. Используя эти оценки, получаем результат леммы.■
Аналогичные результаты можно получить для приближенного вычисления функции / (х) и ее субградиента.
К сожалению, критерий останова алгоритма 1 имеет недостатки с практической точки зрения. Это вызвано двумя причинами. Во-первых, критерий остановки требует вычисления числа 7(ехр(—С(х))), что требует порядка п4 арифметических операций и при высокой размерности задачи трудно выполнимо. Во-вторых, число 7(ехр(-С(х))) может оказаться очень близким к единице для матриц большой размерности, что приведет к необходимости делать очень много итераций для достижения требуемой точности. При этом теоретическая оценка скорости сходимости метода Синхорна является довольно сильно завышенной, как показывает практика. То есть можно ожидать, что эмпирически приближенное значение целевой функции и ее субградиента можно найти с требуемой точностью гораздо быстрее, чем это обеспечивается теоретической оценкой. Планируется посвятить отдельную публикацию аспектам практической реализации предложенного алгоритма поиска значения целевой функции и ее субградиента.
3. Универсальный метод с неточным оракулом
Как было указано в предыдущем пункте, функция р (х) = ¡(х) — д(х) принадлежит классу функций с субградиентом, удовлетворяющим условию Гельдера с параметрами
Ми > 0,и е [0,1]:
||У<р (х) — V <р (у)\\х* < МУУх — у\\1,х, у еЯ, из которого следует, что
р(у) < р(х) + (V? (х), у — х) + ||х — у\\Х+, х, у еЯ.
Кроме того, существует алгоритм, который для заданных ¿1,52 > 0 позволяет вычислить приближенно значение функции р (х) и ее субградиент Vр (х), удовлетворяющие неравенствам
1ф(х) — р(х)1 ^ ¿1 и V<р (х) — V<р(х)
< ¿2.
хх ,*
М,^
Предположим, что множество Я ограничено и имеет диаметр И в норме \\- \\х'. Отметим, что норма, выбранная для измерения константы в условии Гельдера, и норма, выбранная для измерения диаметра множества Я, могут не совпадать, что дает дополнительную гибкость. Проводя рассуждения, аналогичные п. 4.1.3.е [7] и лемме 2 [8], получим следующую лемму.
Лемма 3. Для любого 5 > 0 их е Я мы можем вычислить такие ф$ (х), ^¡р (х), что для любого
Г1 — V 6~
М ^---
[- + »
и любого у е Я справедливо неравенство
0 ^р(у) — [ф& (х) + фцр (х) ,у — х)] < М\\х — у\\Х + (12)
Доказательство. Для этого нужно вызвать алгоритм 1 с параметрами ¿1 = 74,52 = ^ЖБ и положить
Ф& (х) = / (х) — 9(х) — ¿1 = ](х) — д(х) — —,
(х) = ^¡(х) —Vg(х),
где f (х) , V f (х) - выход алгоритма 1.И
Из определения сразу следует, что для любого х е Q справедливы неравенства
о < <р(х) -ф6 (х) < 4.
В силу ограниченности множества Q, с использованием леммы 2 трудоемкость метода для получения (х), V$(p (х)может быть равномерно по х е Q оценена как
(Q
п4 + тп2 + п2 6\ -г-
где 0\,02 - некоторые константы. Следуя [8], введем сильно выпуклую с константой 1 на множестве Q прокс-функцию (1(х). Мы предполагаем, что эта функция достигает минимум на Q в точке хо и что й(хо) = 0. Определим также «расстояние» Брегмана:
С(х, у) = d(y) - d(х) - (У<1(х),у - х).
Приведем модификацию универсального быстрого градиентного метода, позволяющую работать с неточным оракулом (х), V$(p (х).
Алгоритм 2 «Универсальный быстрый градиентный метод для задачи (2)».
Data: Число М0 > 0, точность е > 0 Result: Точка у к
Определим ро(х) ^ ¡;(хо, у), Уо ^ хо, Ао ^ 0 foreach к ^ 0 do
1) Найдем Zk = argminxeQ ipk(х).
2) Найдем минимальное г к ^ 0 такое, что для числа ак+i, ik > 0, найденного из уравнения ak+lik = 2ikMk (Ак + ak+i,ik) и используемого в определениях
Ак+1,ik = Ак + ak+i,ik , Тк,ik = дк+, ^к+1,ik = £тк,ik хк+1,ik = Тк,ik Zk + (1 — Тк,ik )Ук,
хк+1,ik = arg {(Zk, у) + ak+i,ik (VSk+i,ik^ (хк+1,г k) ,У)},
Ук+i,ik = Tk,ik:tk+I,ik + (1 - Tk,ik)yk, выполнено неравенство
V&k+i,ik (Vk+i,ik ) < V&k+I,ik (хк+1,г k ) + {^&k+i,ikV (хк+1,ik ) , yk+i,ik - хк+1^ ) +
+2k-1Мк | | Ук+1,г k -хк+1,г kl IX + ^
3) Положим хк+1 = хк+1,г к, у к+1 = У к+1, i к, ak+i = ak+i,i к, тк = rk,i к. Определим Ak+i = Ак + ak+i, Mk+i = 2гк-1Мк, 5k+i = Sk+i,iк и fk+l{x) = фк (х) + (1к+1[фёк+1 (Хк+i) + (Уёк+1ф (Хк+1), х -Хк+i)}-
end
Algorithm 2: Универсальный быстрый градиентный метод для задачи (2) Теорема 1. Алгоритм 2 корректен. Кроме того, для любого k ^ 0 выполнено
Ак (ф(Ук) - J ) ^ф*к = ттфк(х),
V ¿J xEQ
где
Причем
Рк (х) = (хг) + (хг) ,х - +((хо,х)
i=0
Ак >
1-ei-vkl+3v
22+4v 6i-2M2
Значит, для всех k ' 1 верно
V(Ук) — V* <
22+4v 61-uM2
el-vkl+3v
l + v £
C(xo,x*) + 2.
Если известно число R2 > 0, такое, что £(х0,х*) ^ R2, то при Ак ' выполнено Ук) - V* < е.
Доказательство является понятной модификацией доказательства теоремы 3 работы [8]. ■
Опишем теперь, как можно восстанавливать решение двойственной задачи по информации, генерируемой алгоритмом 2. Для этого предположим, что Q = {х £ Rm : х ' х}. В этой ситуации удобно выбрать d(x) = | ||х — х||2. Тогда х0 = X и ((х, у) = | ||х — у||2. На каждой итерации алгоритма 2 необходимо решить задачу вида
mm{£(х, у) + a(g, у)} = min max { 1 ||х — х||2 + a(g, у) + (s,х — х) 1
xeQ х s'o [2 J
при заданных х, д. Эта задача легко решается в явном виде. В частности, для поиска точки к необходимо решить задачу
( 1 к Zk = arg min < —£(хо,х) + V -1 \(p&i (х*) + (V&iV (х*) ,х — х*
х'х Ак *—' Ак
{ к г=0
Эта задача эквивалентна следующей:
( 1 к аг mmmax < А^^(хо,х) + At^(хг) + (V siV (хг) ,х — х*)] + (s,х — х)
Х s' У к г=0 к
к Ак
к г=0 к
Пусть к - оптимальное значение переменной в последней задаче. Тогда в силу условий
оптимальности в этой задаче имеем
1 к а
—(хо, гк) + У^ф (хг) - вк = 0,
Ак г-^ Ак
=о
(8 к, X - гк) = 0.
Используя эти равенства и определение неточного оракула (12), получаем
Ф (х) + (вк ,Х - х) >
к
' AT [Vöi (хг) + (VsiV (хг) ,х — ^ + ^ — х*)] + (8к,х — ^ + ^ — х) =
=о к
к а 1 = У^ [Фк(Хг) + (У.ф (Хг) , гк - Хг)] + — (У(хо, гк), гк -х) >
Ак Ак
=о
к а 1
=о Ак Ак
= ~А~кФ*к - А^^^
Предположим, что £(х0,х*) ^ Я2, тогда исходная задача эквивалентна задаче шахшт |ф(х) + (в,х - х) : £(х0, х) ^ В,2} = тах'Фп^).
l
Тогда из полученной оценки на р (х) + (Sk ,х — х) и теоремы 1 получим
1 * R ^ ,
АкРк — Ä >р(^ 2 Ак
1 R £ R2
фк(Sk) = тт{р(х) + (s,х — х) : C(х0, х) ^ R2} > -ррk — R > Р(Ук) — ~ — R.
Отсюда, используя слабую двойственность, имеем
£ R2
о < Р(Ук) — фк(Sк) < 2 + ~Ак
Отметим, что предложенный подход легко распространяется на задачи композитной оптимизации, возникающие при поиске стохастических равновесий в транспортных моделях, см. [9].
4. Заключительные замечания
В приложениях часто возникают задачи, имеющие следующий вид (см., например, [10,11])
f(х) = Ф(х, у (х)) ^ min, (13)
х
при этом у (х) и Vy (х) могут быть получены из решения отдельной подзадачи лишь приближенно. Довольно типично, что существует способ, который выдает е-приближенное решение за время, зависящее от е, логарифмическим образом ~ ln (£-1Л). В данной работе намечен общий способ решения таких задач. Его наиболее важной отличительной чертой является адаптивность (самонастраиваемость), т.е. методу на вход не надо подавать никаких констант Липшица (более того, метод будет работать и в негладком случае). Метод сам настраивается локально на оптимальную гладкость функции. Это свойство метода делает его привилегированным, поскольку в реальных приложениях, чтобы что-то знать о свойствах f (х), нужно что-то знать о свойствах зависимости у (х), а это часто не доступно по постановке задачи, или при попытке оценить приводит к сильно завышенным оценкам.
В задаче (13) важно уметь эффективно пересчитывать значения у(х), а не рассчитывать их каждый раз заново (на каждой итерации внешнего цикла). Поясним сказанное. Предположим, что мы уже как-то посчитали, скажем, ( х), решив, например, с какой-то точностью соответствующую задачу оптимизации. Тогда для вычисления у (х + Ах) (на следующей итерации внешнего цикла) у нас будет хорошее начальное приближение ( х). А, как известно (см., например, [11]), расстояние от точки старта до решения (не в случае сходимости метода со скоростью геометрической прогрессии или быстрей) существенным образом определяет время работы алгоритма. Тем не менее известные нам приложения (см., [10,11]) пока как раз всецело соответствуют сходимости процедуры поиска у (х) со скоростью геометрической прогрессии. Связано это с тем, что если расчет у (х) с точностью е осуществляется за О (ln (е-^ операций, то для внешней задачи можно выбирать самый быстрый метод (а стало быть, и самый требовательный к точности). Как правило, такое сочетание оказывается недоминируемым. В частности, для рассматриваемой нами задачи планируется в отдельной статье сравнить описанный здесь подход с возможными альтернативными. Ведь, как известно [12], выпукло-вогнутые седловые задачи, к которым, безусловно, относятся задачи (2), (4), можно по-разному решать, и способ, описанный в п. 3, требует дополнительного сопоставительного анализа, которому и планируется посвятить отдельную статью. Здесь же ограничимся ссылкой на пример 4 и последующий текст из работы [11] и общим тезисом, который пока неплохо подтверждался на практике:
Если есть возможность в задаче оптимизации (в седловой задаче) явно прооптими-зировать по части переменных (или как-то эффективно это сделать с хорошей точностью), то, как правило, это и надо сделать и строить итерационный метод исходя из этого.
В реальных транспортных приложениях [1,2,9,13] достаточно сложным является расчет dj (х) и их градиентов (особенно при поиске стохастических равновесий). Тем не менее
эти задачи имеют вполне четкую привязку к решению некоторых задач на графах типа поиска кратчайших путей. Так же, как и в предыдущем абзаце для внутренней задачи, для внешней задачи можно не рассчитывать с^ (х) и их градиенты каждый раз заново, а пересчитывать; также можно допускать неточность в их вычислении (и ненулевую вероятность ошибки), надеясь на ускорение (благо метод работает в концепции неточного оракула, природа которой не принципиальна, см. [7]). Также здесь оказываются полезными идеи БАД (быстрого автоматического дифференцирования) [14,15], которые позволяют практически за то же время, что занимает вычисление самих функций, вычислять их градиенты.
Исследование А.В. Гасникова, П.Е. Двуреченского, В.Г. Спокойного и А.Л. Сувориковой в части 2 выполнено в ИППИ РАН за счет гранта Российского научного фонда (проект № 14-50-00150), исследование А.В. Гасникова и П.Е. Двуреченского в части 3 выполнено при поддержке гранта РФФИ 15-31-20571-мол_а_вед, исследование Ю.Е. Нестерова в части 3 выполнено при поддержке гранта РФФИ 13-01-12007-офи_м.
Литература
1. Гасников А.В., Дорн Ю.В., Нестеров Ю.Е, Шпирко С.В. О трехстадийной версии модели стационарной динамики транспортных потоков // Математическое моделирование. 2014. Т. 26:6. C. 34-70. arXiv:1405.7630
2. Гасников А.В. Об эффективной вычислимости конкурентных равновесий в транспортно-экономических моделях // Математическое моделирование. 2015. Т. 27, № 12. С. 121-136. arXiv:1410.3123
3. Гасников А.В., Гасникова Е.В., Нестеров Ю.Е., Чернов А.В. Об эффективных численных методах решения задач энтропийно-линейного программирования // ЖВМ и МФ. 2016. Т. 56, № 4. (принята к печати) arXiv:1410.7719
4. Cuturi M. Sinkhorn Distances: Lightspeed Computation of Optimal Transport // NIPS, 2013.
5. Benamou J.D., Carlier G., Cuturi M., Nenna L., Peyre G. Iterative Bregman Projections for Regularized Transportation Problems. e-print, 2015. arXiv:1412.5154 (to appear in SISC).
6. Franklin J., Lorenz J. On the scaling of multidimensional matrices // Linear Algebra and its applications. 1989. V. 114. P. 717-735.
7. Devolder O. Exactness, inexactness and stochasticity in first-order methods for large-scale convex optimization. CORE UCL, PhD thesis, March 2013. http://www.ecore.be/DPs/dp_1327057920.pdf
8. Nesterov Yu. Universal gradient methods for convex optimization problems // Mathematical Programming Series A. 2015. V. 152, I. 1, P. 381-404.
9. Гасников А.В., Гасникова Е.В., Двуреченский П.Е., Ершов Е.И., Лагуновская А.А. Поиск стохастических равновесий в транспортных моделях равновесного распределения потоков // Труды МФТИ, 2015. Т. 7, № 4. C. 114-128. arXiv:1505.07492
10. Bogolubsky L., Dvurechensky P., Gasnikov A., Gusev G., Nesterov Yu., Raigorodskii A, Tikhonov A., Zhukovskii M. Learning supervised PageRank with gradient-free optimization methods. e-print, 2014. arxiv:1411.4282
11. Гасников А.В., Двуреченский П.Е., Нестеров Ю.Е. Стохастические градиентные методы с неточным оракулом. e-print, 2015. arXiv:1411.4218
12. Nemirovski A. Lectures on modern convex optimization analysis, algorithms, and engineering applications.
Philadelphia: SIAM, 2013. http://www2.isye.gatech.edu/~nemirovs/Lect_ModConvOpt.pdf
13. Гасников А.В., Двуреченский П.Е., Дорн Ю.В., Максимов Ю.В. Численные методы поиска равновесного распределения потоков в моделях Бэкмана и стабильной динамики // Математическое моделирование. 2016. Т. 28. (принята к печати) arXiv:1506.00293
14. Ким К., Нестеров Ю., Скоков В., Черкасский Б. Эффективные алгоритмы для дифференцирования и задачи экстремали // Экономика и математические методы. 1984. Т. 20. С. 309-318.
15. Евтушенко Ю.Г. Оптимизация и быстрое автоматическое дифференцирование. М.: ВЦ РАН, 2013.
References
1. Gasnikov A.V., Dorn Y.V., Nesterov Y.E., Shpirko S.V. On the three-stage version of stable dynamic model. Math modeling. 2014. V. 26:6. P. 34-70. arXiv:1405.7630 (in Russian).
2. Gasnikov A.V. On the effective computability of competitive equilibrium in the transport-economic models. Math modeling. 2015. V. 27, N 12. Р. 121-136. arXiv:1410.3123 (in Russian).
3. Gasnikov A.V., Gasnikova E.V., Nesterov Y.E., Chernov A.V. On effective numerical methods for solving linear-entropy programming // JCM and MP. 2016. V. 56, N 4. (in print) arXiv:1410.7719 (in Russian).
4. Cuturi M. Sinkhorn Distances: Lightspeed Computation of Optimal Transport. NIPS, 2013.
5. Benamou J.D., Carlier G., Cuturi M., Nenna L., Peyre G. Iterative Bregman Projections for Regularized Transportation Problems. e-print, 2015. arXiv:1412.5154 (to appear in SISC).
6. Franklin J., Lorenz J. On the scaling of multidimensional matrices. Linear Algebra and its applications. 1989. V. 114. P. 717-735.
7. Devolder O. Exactness, inexactness and stochasticity in first-order methods for large-scale convex optimization. CORE UCL, PhD thesis, March 2013. http://www.ecore.be/DPs/dp_1327057920.pdf
8. Nesterov Yu. Universal gradient methods for convex optimization problems. Mathematical Programming Series A. 2015. V. 152, I. 1, P. 381-404.
9. Gasnikov A.V., Gasnikova E.V., Dvurechensky P.E., Ershov E.I., Lagunovskaya A.A. Search for the stochastic equilibria in the transport models of equilibrium flow distribution. Trudy MIPT. 2015. V. 7, N 4. P. 114-128. arXiv:1505.07492 (in Russian).
10. Bogolubsky L., Dvurechensky P., Gasnikov A., Gusev G., Nesterov Yu., Raigorodskii A, Tikhonov A., Zhukovskii M. Learning supervised PageRank with gradient-free optimization methods. e-print, 2014. arxiv:1411.4282
11. Gasnikov A.V., Dvurechensky P.E., Nesterov Y.E. Stochastic gradient method with inexact oracle. e-print, 2015. arXiv:1411.4218 (in Russian).
12. Nemirovski A. Lectures on modern convex optimization analysis, algorithms, and engineering applications. Philadelphia: SIAM, 2013. http://www2.isye.gatech.edu/~nemirovs/Lect_ModConvOpt.pdf
13. Gasnikov A.V., Dvurechensky P.E., Dorn Y.V., Maximov Y.V. Numerical methods for finding an equilibrium flow distribution in Beckmann model and stable dynamics. Math modelling. 2016. V. 28. (in print) arXiv:1506.00293 (in Russian).
14. Kim K. Nesterov Yu, Skokov V., Cherkassky B. Effective algorithms for the derivation and extremal problem. Economics and Mathematical Methods. 1984. V. 20. P. 309-318.
15. Evtushenko Y.G. Optimization and fast automatic differentiation. M.: CC RAS, 2013. (in Russian).
Поступила в редакцию 10.11.2015.
