CC BY
45
УДК 519.688
А. В. Гасников1'2, П. Е. Двуреченский2'4, В. Г. Спокойный1'2'3'4, П. И. Стецюк5,
А. Л. Суворикова1'2'6
1 Московский физико-технический институт (государственный университет)
2Институт проблем передачи информации РАН 3Национальный исследовательский университет Высшая школа экономики 4Weierstrass Institute for Applied Analysis and Stochastics, Berlin 5Институт кибернетики им. В.М. Глушкова НАН Украины 6The International Research Training Group 1792, Berlin
Суперпозиция метода балансировки и универсального градиентного метода для поиска энтропийно-регуляризованного барицентра
Т| и г* г*
Вассерштейна и равновесий в многостадийных моделях
транспортных потоков
Представлен обзор современных численных методов поиска барицентра Вассерштейна конечного семейства вероятностных мер с одинаковым конечным носителем. Такие задачи в последнее время стали очень популярны в связи с всевозможными приложениями к сравнительному анализу изображений, в частности, к обнаружению разладок в ряде изображений. Например, подобные задачи возникают при изучении деятельности головного мозга.
В основном мы исходили из цикла работ М. Си1ип с соавторами. Общая идея этих работ - найти барицентр вероятностных мер согласно энтропийно-регуляризованному расстоянию Вассерштейна. Такое (регуляризованное) расстояние можно заметно эффективнее посчитать, чем исходное расстояние Вассерштейна. В одной из работ отмеченного цикла [8] содержалась идея сочетания метода Синхорна (балансировки) для решения внутренней задачи (расчета соответствующих регуляризованных расстояний и их субградиентов) и быстрого градиентного метода для решения внешней задачи (поиск барицентра). К сожалению, в описанном авторами виде метод оказался не пригодным для использования на практике (не было также никаких теоретических гарантий его сходимости). В [1] показано, как можно доработать данный метод (в частности, доказана сходимость предложенной модификации). Однако мы были сконцентрированы на другом приложении (к поиску равновесий в многостадийных транспортных моделях).
В данной работе мы рассматриваем оба отмеченных приложения. Главным результатом работы является разработка в общем случае (не только для этих двух приложений) концепции суперпозиции методов, когда мы можем выделить в исходной задаче часть переменных, по которым задача эффективно решается внутреннем методом (но допускается, что лишь приближенно) при замороженных остальных переменных. А по оставшейся группе переменных запускается внешний метод, на каждой итерации которого требуется запускать внутренний метод. В статье получен частичный ответ на довольно общий вопрос: как оптимально сочетать эти методы (внутренней и внешний), т.е. насколько точно надо решать на каждой итерации внешнего метода внутреннюю задачу, чтобы минимизировать общее время работы метода при заданной точности решения, которую хотим получить?
Ключевые слова: седловая задача, энтропия, метод балансировки, метод Синхорна, универсальный метод, неточный оракул, задача Монжа-Канторовича, расстояние Вассерштейна, барицентр Вассерштейна.
A. V. Gasnikov1'2, P.E. Dvurechensky2'4, V. G. Spokoiny1'2'3'4, P.I. Stetsyuk5,
A. L. Suvorikova1'2'6
1 Moscow Institute of Physics and Technology (State University) 2 Institute for Information Transmission Problems RAS 3Higher School of Economics National Research University 4Weierstrass Institute for Applied Analysis and Stochastics, Berlin 5V. M. Glushkov Institute of Cybernetics of NAS of Ukraine 6The International Research Training Group 1792, Berlin
Superposition of the balancing algorithm and the universal gradient method for search of the regularized Wasserstein brycenter and equilibria in multistage transport models
The paper presents a survey of modern computational algorithms the mean (in the sense of 2-Wasserstein distance) of a finite set of empirical probability measures. Each measure is assumed to have the same support. The problem is quite popular due to numerous applications to comparative image analysis, for example the detection of structural changes in a set of fMRI images.
The work is inspired by a number of papers written by M. Cuturi and his colleagues. They propose to use regularized 2-Wasserstein distance instead of the original one. The new approach significantly improves computational efficiency. The paper [8] proposes to combine Sinkhorn algorithm with the fast gradient method. The first algorithm computes the regularized distance and its subgradient, while the second one is used for barycenter search. However, this approach uses some predefined parameters that cannot be selected in practice. Moreover, there is no theoretical guarantee that it converges. The paper [1] proposes a modification of the algorithm, which solves the above problems.
The work [1] is mainly focused on searching equilibria in the multistage traffic model, whereas this paper investigates its application to the search for Wasserstein barycenters. The key result is a general concept that allows constructing the superposition of algorithms.
It is used if the assumption holds. The set of variables can be divided into two groups (A) and (B), s.t. by fixing (A) one can obtain (an inexact) solution by (B). Thus, the initial problem can be reduced to successive optimization by (A) and (B). The present work aims to answer the following questions: how to combine the steps and how accurate one should be when optimizing by (B) so that to achieve the predefined accuracy?
Key words: saddle point problem, entropy, balancing algorithm, Sinkhorn's algorithm, universal gradient method, inexact oracle, Monge-Kantorovich problem, Wasserstein distance, Wasserstein barycenter.
1. Введение
Данная статья представляет собой продолжение недавно опубликованной работы [1]. Для удобства мы повторим во введении постановку задачи и основные обозначения.
Поиск (стохастических) равновесий в многостадийных моделях транспортных потоков приводит к решению следующей седловой задачи с правильной (выпукло-вогнутой) структурой [2] - [5]:
п п
min maxi Y^ хц ln x^ + V^ Cij(у)хц +
£?=i = Li, E?=1 = yeQl J J
Xij > 0,i,j = l,...,n i,3 = 1 i,3 =1
п
= Удл^Д L} + w) - S exp(-cü(y) - 1 + Ai + »3) + 9(y)}, (1)
i, 3 = 1
где g(y) и Cij(y) ^ 0 - вогнутые гладкие функции (если ищутся не стохастические равновесия, то с^(у) могут быть негладкими), Q - множество простой структуры, например,
Q = {у е Rn : у > у}.
Легко понять, что система балансовых ограничений в (1) либо несовместна
п п
^г = Е , либо вырождена (имеет не полный ранг). В последнем случае это приво-
г=1 з=1
дит к тому, что двойственные переменные (Л, определены с точностью до произвольной постоянной С:
(Л + Се,^ - Се), е = (1,..., 1).
п
Задачу (1) также можно переписать следующим образом (не ограничивая общности,
п п
считаем ^ Ьг = ^ = 1):
г=1 з=1
шт
Ё = Lí, Ё х^ =
> 0, £ =1,1,3 = 1,...,п
шах уея I
п п
{ ^ ХИ 1п ХИ + ^ СИ (У)хИ +
1,3 = 1
1,3 = 1
шах шах уея А/еж™
V) + {ц., W) - ^ ехр(-Сгз (у) + Хг + ^) г,з=1
= - ш^п / (У),
+
где выпуклая функция /(у) определяется как
(2)
f (у) = ^ шах
Ё = Lí, Ё = 3=1 ¿=1
п п
{- ^ Хг 3 1п хг з - ^ Сг з(у)3
, =1 , =1
11п^У2 еХР(-Сгз(У) + ^г + »3) -{Х,Ь)-{^,№)
А,/л€Кп I. I *—' J
, =1
(3)
Поскольку мы добавили в ограничения условие Е Хгз = 1, являющееся следствием
, =1
балансовых уравнений, то это привело к тому, что двойственные переменные (Л, у) определены с точностью до двух произвольных постоянных С а, С/: (Л + С а е, Ц- + С/е).
Заметим, что расчет градиента V/(у) (в ряде транспортных приложений вогнутые функции Сгз(у) - негладкие, тогда вместо градиентов стоит понимать суперградиенты сг з(у) и субградиент f (у)) осуществляется по следующей формуле (Демьянова-Данскина-Рубинова, см., например, [2,6]):
п
Е ехр(-Сгз(у) + X* + Vciз(y)
V/(у) = - ^^--Vg(y) =
Е ех^-Сгз(у) + X* +
, =1
^ хгз(\*, сгз(у) -Vg(y),
, =1
(4)
где (Х*,^*) - решение задачи (3), не важно какое именно, градиент V/(у) от выбора Са, С/ (см. выше) не зависит. В данной статье мы (так же, как и в [1]) ограничимся изучением только полноградиентных методов для задачи (2), т.е. не будем рассматривать, например, рандомизацию при вычислении градиента по формуле (4). Планируется отдельно исследовать вопрос о возможности ускорения вычислений за счет введения рандомизации для
'Х4
внешней задачи. На данный момент нам представляется (см. формулу (13) в п. 3), что это может принести дивиденды только в случае, когда вспомогательная задача расчета достаточно сложная. Тут требуется много оговорок, в частности, в большинстве приложений умение рассчитывать ^с^(у) для конкретной пары (г,]) без дополнительных затрат позволяет заодно рассчитать и все Ус^(у), ] = 1, ...,п. Также отдельно планируется исследовать вопрос о том, какие подходы и насколько хорошо допускают распараллеливание. Вопрос о целесообразности рандомизации оказывается завязанным и на вопрос о возможности распараллеливания.
Структура статьи следующая. В п. 2 мы рассматриваем популярную в последнее время (в связи с большим числом приложений) задачу вычисления барицентра Вассерштейна различных вероятностных мер [7]. Эта задача оказывается тесно связанной с задачей (1). Мы разбираем в статье этот пример, потому что он хорошо проясняет возможные альтернативы предлагаемому нами основному подходу решения задач (1), (2), изложенному в п. 3. В основе подхода п. 3 (см. также [1,8]) лежит сочетание метода балансировки для решения внутренней задачи оптимизации по двойственным множителям и универсального градиентного метода с неточным оракулом для внешней задачи (2). В работе [1] мы сконцентрировались на обобщении универсального метода на случай неточного оракула. В данной работе акцентируем внимание на получении условий на шум оракула для внешней задачи, исходя из оценок точности решения внутренней задачи. Таким образом, данная статья дополняет работу [1] в теоретическом плане. Практическим экспериментам планируется посвятить отдельную работу.
2. Поиск барицентра Вассерштейна
К похожей на (1) задаче приводит поиск барицентра Монжа-Канторовича (в западной литературе чаще говорят барицентра Вассерштейна [7, 8])1. Изложим вкратце постановку задачи [8] - [10]. Вводится энтропийно регуляризованное транспортное расстояние (см. рис. 1 в [11]) с матрицей НсуН™^, сформированной из квадратов попарных расстояний °гз = от носителя г меры Ь до носителя ] меры Ш (Ь.Ш € 5п(1)):
п п А(Ь.Ш) = шт {7 ^ Xу 1пХгз + ^ с^х^
XI3 = Ь1, хг3 = ^3 i. 7=1 i. 7=1
3 = 1 г=1 ™ ™
XIз ^ 0,г, ] = 1, ..., п
«Л (А. Ь) + (,, IV) - 7 ± е*р( — Л =
к ^=1 )
шах< (X. V) — 7 У^ Шп 1п
Ае к« 1 ^
3 = 1
£ ехр( —)]}
(5)
Н*ш (А)
Определим при Ь € 5п(1) функцию Нщ(Ь) = А(Ь.Ш). Эта гладкая на Ь € 5п(1) функция с градиентом (см. утверждение 3 [10]):
1
(Ь) = Л*.
хСтрого говоря, мы будем искать барицентр вероятностных мер не согласно настоящему (негладкому) расстоянию Вассерштейна (на наш взгляд, исторически более правильно это расстояние называть расстоянием Монжа-Канторовича-Добрушина), как это можно было подумать из названия, а согласно энтропийно-регуляризованному расстоянию Вассерштейна [8].
где Л* - единственное решение (5), удовлетворяющее условию2 {X*, е) = 0. Отсюда следует, что
HwW = тах л (К L) — Hw
Lesn (1) I
п 1 п ,
= 7 £ ln[i £exp()
j=1 j i=1 '
Теперь можно перейти к изложению основной конструкции. Задача поиска барицентра Вассерштейна3 записывается следующим образом:
т
THwk (L) ^ min . (6)
Les"(i)
К сожалению, в такой формулировке мы не можем оценить константу Липшица градиента функционала (6), явно входящую в большинство современных быстрых (ускоренных) численных методов. Однако оказывается (см. п. 3), что существуют быстрые методы, которым для работы не требуется такая информация (константа Липшица градиента). Перепишем задачу (6), следуя п. 3 работы [10], следующим образом:
У] Hwk (Lk)
k=1
^ max
Li = Lm\ A\L! e S„(1) -^m — 1 = Am , Lm e s„(i)
т—1 т—1
^ ьЩ^ ^ ,Ьк) - (Ьк >} + {- ^ ^ ,Ьт) - М ^ , к=1 к=1
т—1 т—1
£ Нк (Хк) + (- £ Хк) ^ шш ^ (7)
к=1 к=1 '"''
Ь* = (\к) для любого к = 1, ....,т — 1, где Ь* - единственное решение задачи (6),
(Лк}т=—1 - единственное решение задачи (7). Важное свойство функционала задачи (7) -равномерная ограниченность константы Липшица градиента (следует из [15]). Задача безусловной минимизации (7) может быть эффективно решена различными способами (в зависимости от того, насколько велики п и т). В частности, для больших п и т, неплохо с задачей справляются различные модификации метода сопряженных градиентов и быстрых градиентных методов [10]. Структура задачи (7) позволяет эффективно использовать покомпонентные методы (см., например, [16,17]), которые к тому же хорошо параллелятся для данной задачи. Задача (7) хорошо также решается с помощью распределенных вычислений [18].
В приложениях к поиску разладки требуется много раз перерешивать задачу (7), которую для симметричности перепишем следующим образом:
£ и
* (Xk) ^ min
Lwk (
k=1 £ хк=о
fc=i
2Решая задачу (6) каким-нибудь прокс-методом с КЬ прокс-структурой [12], легко понять, что от того, как именно выбирать Л*, задаваемое с точностью до сдвига всех компонент на одно и то же произвольное число, метод зависеть не будет. Единственное, для чего имеет смысл стремиться к выполнению этого нормирующего условия, так это для лучшей практической обработки (меньшее накопление ошибок округления из-за конечной длины мантиссы) экспоненциального взвешивания компонент градиента, возникающего на каждом шаге итерационного процесса при выборе КЬ прокс-структуры.
3К сожалению, пока не так много известно о статистической обоснованности использования расстояния Вассерштейна. Другими словами, хотелось бы иметь связь барицентра Вассерштейна с оценками максимального правдоподобия, ну или хотя бы с состоятельными оценками для соответствующих схем экспериментов. Пока установлена только связь с состоятельными оценками [13,14].
немного смещая окно, т.е. заменяя каждый раз несколько первых слагаемых в сумме
HW..,Hwr)
на столько же новых (которые, как ожидается, близки к HW (\т)). В таком случае предлагается в итерационном процессе стартовать при сдвиге окошка с того, на чем остановились на прошлом положении окошка, экстраполируя \т на вновь пришедшие слагаемые. Ясно, что для новой задачи набор
( \r+l \т— 1 \т \т \т \
И* , Л* , Л* , Л* , Л*
с которого стартуем, уже не будет оптимальным, однако мы вправе надеяться на его близость к оптимальному набору, что существенно сокращает число последующих итераций. Интересной, особенно в данном контексте, представляется возможность использования (и интерпретации) распределенных вычислений [18].
Может показаться, что подход, сводящий поиск решения задачи (6) к задаче (7), не доминируем, поскольку, в отличие от задачи (6), в задаче (7) мы можем явно выписать функционал и по простым формулам, рассчитать градиент, который к тому же имеет равномерно ограниченную константу Липшица. С одной стороны, это, действительно, преимущество, но получено оно дорогой ценой - ценой раздутия пространства, в котором происходит оптимизация почти в т, раз. И это раздутие скажется не только на сложности одной итерации. Для задачи (6) осуществление одной итерации будет еще более дорогим ввиду необходимости на каждой итерации решать т отдельных подзадач расчета VHwk (L). Скажется это, прежде всего, на числе необходимых итераций. В следующем разделе будет отмечено, что расчет VHwk (L) с помощью метода балансировки не намного сложнее расчета (\k). При этом задача (6) решается в пространстве намного меньшей размерности, и мы вправе ожидать, что необходимое число итераций может быть намного меньше, чем для задачи (7). Кроме того, задача (6) решается на компакте, т.е. размер решения (если быть точным, то расстояние от точки старта до решения), входящий в оценку необходимого числа итераций, заведомо ограничен размером симплекса. Задача (7) - задача безусловной оптимизации, причем без свойства сильной выпуклости функционала. Размер ее решения может быть большим, и входит он в оценки необходимого числа итераций так же, как и для задачи (6), к сожалению, степенным образом (для быстрых (ускоренных) методов можно ожидать линейной зависимости необходимого числа итераций от этого размера). Наконец, для постановок задач об обнаружении разладки (см. выше) также ожидается, что использовать близость решений прямых задач (6) при смещении окошка удастся намного лучше, чем близость в решении двойственных задач (7). В итоге выгода от подхода, связанного с переходом к задаче (7), уже не столь очевидна и требует отдельного и более аккуратного исследования с решающей ролью численных экспериментов.
В ряде задач требуется искать параметрический барицентр Вассерштейна. В таком случае в одном из вариантов постановки предполагают наличие параметрической зависимости L(0) £ Sn(1), в £ в, где размерность вектора параметров dim в ^ п. К сожалению, в этом случае нельзя гарантировать с помощью стандартных приемов [19, с. 86] выпуклости задачи
т
£ HWk (ВД) ^ min, (8)
k=1
за исключением случая, когда L(9) = АО + Ь, а в - выпуклое множество. В этом случае конструкция (7) видоизменяется следующим образом:
т
- £ HWk iLk) ^ k=1
i = Lm 1 Am , Lm—i £ Sn
Lm = Ae + bi X, Lm e sn(i),e e в
max
Li = Lm i \\ Li e Sn(1)
т—1 т—1
£ Н^ (Ак) + H*Wm{- £ \к — л)+ ^Ав + 6> ^ min х , (9)
к=1 к=1 '"¿'е е '
L* = (Х'к) для любого к = 1,....,т — 1. Как следствие, нет никаких гарантий, что
изложенная выше конструкция, связанная с переходом к двойственной задаче (7) и восстановлением решения прямой задачи (6) по явным формулам через двойственные множители, в общем случае будет работать хотя бы для поиска локальных решений.4 Другими словами, необходимо искать глобальный оптимум задачи (8) исходя из работы с прямой задачей (8). Один из вариантов того, как это можно делать, будет описан в следующем пункте.5
Однако при другом варианте постановки (более предпочтительном) можно задавать зависимость L(9) с помощью аффинных равенств и выпуклых неравенств:
min .
Ае + BL = с g(e, L) < 0
l е s„(i); е е е
Многие параметрические зависимости представимы в таком виде [21]. В частности, отметим полиэдральные представления Фурье-Моцкина [21], возникающие, например, в робастной оптимизации:
т
У^ HWk (L) ^ min .
к=1 Le{beS„(l): Зв:Ав+ВЬ^с]
Можно переписать задачи (7), (9) и на эти случаи, причем сделать это корректно в том смысле, что правомочность подхода полностью сохранится. При этом принципиально ничего из сказанного выше не поменяется. Подробнее об этом планируется написать в отдельной работе.
В действительности, в приложениях наиболее интересен случай, когда ищется барицентр именно расстояний Вассерштейна,6 а не энтропийно-регуляризованных расстояний
[8] - [11]. Другими словами, интересно изучать предельное поведение j ^ 0+ (см. п. 3.1
[9], утверждение 1 [10], п. 3 и конец п. 4 [11]). К сожалению, методы из пп. 2, 3 оказываются весьма чувствительными к этому предельному переходу. Для метода из этого раздела константа Липшица градиента в задаче (7) будет расти как 7-1, соответственно, число итераций будет увеличиваться (при использовании быстрых (ускоренных) методов) как j-1/2. Еще более плохое поведение (см. [6]) можно ожидать от метода балансировки, использующегося в подходе из п. 3. Планируется в отдельной публикации исследовать вопрос о том, как следует действовать при малых j > 0. В частности, в вырожденном случае j = 0. По-видимому, в этом случае поможет философия искусственного сглаживания 7 [15], в которой искусственно введенная энтропийная регуляризация уже задается с четко заданным
4Впрочем, есть результаты (см. формулу (8) п. 3 § 2 главы 8 [20]) о локальной сходимости обычного градиентного спуска для задачи (9) при некоторых дополнительных предположениях.
5При этом правомочность подхода п. 3 для постановки задачи (8) имеет место при дополнительном предположении, что метод стартует из выпуклой окрестности точки минимума с небольшим запасом, допускающим возможность по ходу итерационного процесса оказаться дальше от решения, чем в начальный момент.
6Численные методы поиска «честного» барицентра Вассерштейна вероятностных мер в основном строятся на том, что когда меры заданы на прямой, существуют эффективные способы решения задачи поиска барицентра [7]. Далее проектируют (считают преобразования Радона) меры на случайные прямые и решают одномерные задачи. По их решениям восстанавливают решение исходной задачи [22,23]. В отличие от других подходов, здесь существенно используется структура матрицы затрат сц = l2j (в дискретном случае). Интересно было бы исследовать вопрос о применимости этого подхода к постановкам задач о разладках, в которых требуется много раз пересчитывать барицентр (см. выше). Также интересно было бы сравнить описанные подходы с остальными. Этому планируется посвятить отдельную публикацию.
7Это сглаживание правильно называть двойственным сглаживанием, поскольку для того чтобы добиться гладкости в прямой негладкой задаче, которая имеет лежандровское (седловое) представление [15,21], в это представление, которое также можно понимать как двойственное, вводят аддитивным образом с небольшим коэффициентом сильно выпуклый (вогнутый) функционал. Этот функционал и обеспечивает
Y^Hwk (L) ^ к=1
коэффициентом регуляризации ^ > 0, зависящим от итоговой точности по функции, с которой требуется решить задачу. Другой способ - использовать менее чувствительные (чем метод балансировки) способы решения двойственной задачи, например, при небольших значениях п ожидается, что лучше сработает г-алгоритм Н. З. Шора и некоторые его обобщения [24,25]. В данной работе мы фиксируем ^ > 0 и далее уже не будем возвращаться к подобного рода вопросам.
В заключение отметим, что поиск барицентра Вассерштейна в случае т = 1 может быть осуществлен явно: Ь = Ш. Обоснование этого частного результата представляется довольно полезным для понимания основной конструкции этого раздела.
3. Универсальный метод с неточным оракулом
Из п. 2 следует, что внутренняя задача максимизации по (Л, р) может быть явно решена по р при фиксированном Л, и наоборот (это верно для задач (1) и (2) и приводит к одним и тем же формулам). Собственно, таким образом, получается метод балансировки расчета матрицы корреспонденций по энтропийной модели, см., например, [6] (тесно связанный с методом Синхорна [9,11]), как метод простой итерации для явно выписываемых условий экстремума (принципа Ферма): Л = Л(^), р = М(А). Метод балансировки имеет вид ([А]о = [р]о = 0):
В этих формулах «—1» в экспоненте для метода (2) (в отличие от метода (1)) можно не писать, поскольку двойственные множители определяются неоднозначным образом с большим произволом для задачи (2) (см. выше), достаточным для справедливости этого замечания.
Оператор (X, р) ^ (Л(р), М(А)) является сжимающим в метрике Биркгофа-Гильберта р [26]. Это означает, что после N ~ 1п(ст-1) итераций метода балансировки можно получить такие (\м), что { {Х*(у), р*(у))} - двумерное аффинное множество решений (см. п. 1):
Причем на практике наблюдается очень быстрая сходимость, т.е. коэффициент пропорциональности не большой [6]. Таким образом, мы можем приближенно решить внутреннюю задачу.8
Далее предлагается воспользоваться прямодвойственным (эта важно, поскольку нужно восстанавливать двойственные переменные) универсальным методом [27] для решения
гладкость (а еще точнее липшицевость градиента) в прямой задаче. В нашем случае мы исходим из задачи о перемещении масс (Монж-Канторович), являющейся задачей ЛП. Однако для удобства вычисления расстояний Вассерштейна мы перешли к двойственной задаче. Мы хотим сделать гладкой двойственную задачу, потому что именно с ней в дальнейшем и идет работа. Для этого двойственное сглаживание (в нашем случае энтропийное) применяется к двойственной задаче для двойственной задачи к транспортной задаче, т.е. применяется просто к транспортной задаче.
8Заметим, что в пп. 3.1, 3.2 работы [9] предлагается за счет раздутия прямого пространства с помощью обобщения описанного метода балансировки Брэгмана (метода проекций Брэгмана) решать задачу поиска барицентра напрямую, т.е. отпадает необходимость в решении внешней задачи. Плата за это достаточно большая - увеличение размера прямого пространства в т раз, но метод при этом хорошо параллелится.
или
(10)
внешней задачи оптимизации по у (имеется видеопрезентация с описанием этого метода [28]). К сожалению, в формулировке (1) (в отличие от формулировки (2)) кроме того, что внешняя задача гладкая (при условии гладкости с^(у) [29,30]), больше ничего о ней сказать нельзя (константа Липшица градиента не ограничена). Также не понятна гладкость задачи (6). Поэтому и по ряду других причин, о которых будет сказано далее, было отдано предпочтение универсальному методу, оптимально адаптивно настраивающемуся на гладкость функционала f (у) на текущем участке пребывания итерационного процесса.9 Однако нам потребуется использовать этот метод в варианте с неточным оракулом, выдающим градиент [31]. Напомним (см. п. 1), что мы решаем задачу 2, представимую в виде (здесь в тах представлении х = х,
Q = {Xij ^ 0,i,j = 1,...,п : ^ Xij = Li, ^ Xij = Wj^,
3=1 i=l
2n\.
а в min представлении x = (A, ß), Q = R
f (y) = max Ф(ж, у) = min Ф(х, у) ^ min у £ Q.
Далее везде будем считать, что у £ Q.
Определение 1 (см. главу 4 [32]). (5, Ь)-оракул выдает (на запрос, в котором указывается только одна точка у) такие [F(у), G(y)), что и для любых у, у' £ Q
0 < f (у') - F(y)-(G(y),y' - у) < |У - у\\2 + S. Из определения 1 сразу следует, что для любого х £ Q
F(х) < f (х) < F (х) + S
и для любых х,у £ Q
f (У) > f (У) -(G(x),y - х)- 5.
Из последнего свойства получаем, что определение (ö, Ь)-оракула можно понимать как обобщение на гладкие задачи классического понятия негладкой выпуклой оптимизации: ¿-субградиента (см. п. 5 § 1 главы 5 [20]). В приводимом далее утверждении в первой его части следует сохранить обозначения для задачи (2), (3) и следует обозначить х = А, у = L для задачи (5), (6); а во второй части утверждения следует обозначить х = (X,ß), у = у для задачи (2), (3). Таким образом, на задачу (2), (3) можно посмотреть с двух разных ракурсов, однако второй ракурс менее привлекателен ввиду необходимости рассмотрения ограниченных множеств Q, что в интересующих нас приложениях место не имеет.
9Бытует мнение, с которым столкнулись и авторы данной статьи, что любой универсальный метод должен чем-то платить за свою универсальность, и в этой связи возникает много вопросов, в частности: насколько дорога эта плата? В принципе, в статье [27] довольно подробно проясняется этот момент. Тем не менее мы повторим здесь соображения из [27]. Действительно, плата за универсальность есть. Универсальный метод из работы [27] может сделать где-то в 4 раза больше обращений к оракулу для задачи с более менее одинаковой константой Липшица градиента во всей области (где довелось пройти итерационному процессу), по сравнению с обычным быстрым градиентным методом [15]. Тем не менее замечательная особенность универсального метода не только в том, что он настраивается на гладкость задачи и применим к любым задачам, но и в том, что (в отличие от подавляющего большинства методов) этот метод локально настраивается на гладкость функционала. И для сильно неоднородных функционалов типично, что универсальный метод делает заметно меньше итераций, чем, скажем, быстрый градиентный метод (плата за универсальность уже учтена в отмеченном выше потенциально возможном увеличении числа итераций в 4 раза в худшем случае). Примеры, поясняющие сказанное, имеются в работе [27].
Утверждение 1. Если ф(у) = тах^(х,у), где ^(х,у) - выпуклая по у и вогнутая по х
хея
функция, и найден такой х € , что
ф(у) — Щх,у) < 5,
то субградиент ду Ф(х,у) есть ¿-субградиент функции ф(у) в точке у. Если
<^(у) = тш Ф(х, у), где Ф(х,у) - выпуклая по совокупности переменных функция, и хе<С
найден такой х € , что
тах(Фх(ж,у),х — г) ^ 5, С)
то
Ф(х,у) — ¡р(у) < 6
и субградиент Фу(ос, у) = дуФ(х,у) есть ¿-субградиент функции (р(у) в точке у.
Доказательство. Ограничимся доказательством только второй части этого утверждения. Доказательство первой части см. на с. 124 (лемма 13) книги [20]. Из выпуклости Ф(х,у) по совокупности переменных имеем
Ф(х', у') ^ Ф(х, у) + (Фх(х, у),х' — х) + {Фу(х, у), у' — у). (11)
Определим зависимость х(у) из соотношения
1р(у) = тщФ(х,у) = Ф(х(у),у). хес)
Заметим, что Ф(х,у) ^ <р(у). Положим в (11) х' = х(у'), х = х. Тогда
р(у') = Ф(х', у') ^ Ф(х, у) + (Фх(х, у),х' — х) + {Фу(х, у), у' — у) ^ ^ ^(у) + {Фх(х, у), х(у') — х) + (Фу(х, у), у' — у) ^
> ^(У) + Ф(х,у),у' — у) — 5.
В последней формуле мы использовали, что
(Фх(х,у),х — х(у') < 6.
В свою очередь из выпуклости Ф(х,у) по х (для всех допустимых у) имеем
Ф(х,у) — Ф(х(у'),у) < (Фх(х,у),х — х(у')).
Беря в этой формуле у' = у и воспользовавшись определением х(у), получаем, что
Ф(х, у) — <р(у) < (Фх(х, у),х — х(у)).
Однако не хочется довольствоваться возможностью находить только ¿-субградиент (из утверждения 1 эта возможность очевидна), поскольку в определенных ситуациях явно можно рассчитывать на некоторую гладкость итоговой (внешней) задачи (2). Понятие (5,Ь)-оракула в некотором смысле налагает наиболее слабые условия на возможные неточности в вычислении функции и градиента, при которых можно рассчитывать, что скорость сходимости метода, учитывающего гладкость (липшицевость градиента функционала) задачи, сильно не пострадает (см. теорему 1 ниже).
На первый взгляд может показаться, что применимость описанной концепции (5, V)-оракула к задаче (1) следует из следующего результата (см. п. 4.2.2 [32]).
Утверждение 2. Пусть подзадача энтропийно-линейного программирования (ЭЛП) в (2) решена (по функции) с точностью 5, т.е. найден такой х(с), удовлетворяющей балансовым ограничениям, что
п п
Хц(с)Ы Х^(с) + сцхц(с)
1,3 = 1 1,3 = 1
п п
mm
3 = 1
^ ^ ^ Xi j In Xij + ^ ^ C-i jXij J* ^
= i,3 = l i,3 = l
Xij ^ 0, i, j = 1, ..., n
Тогда для функции
f (с) = - ^ mill j ^ Xij ln Xij + ^ CijXijJ
/ - -^zj — ^г > А-/ ^zj
з=1 г=1
хгj ^ 0,i,j = 1,..., п
Wj
i, 3 = 1 i ,3 = 1
набор
(п п пп \
^ Хц(с) lnХц(с) + ^ CijXij(C)^Xij(C)^ п = I i, 3=1 i, 3=1 iJ /
является (5, 2 max Ci j)-оракулом.
i, ]=1,...,п
К сожалению, большинство методов (в том числе метод балансировки) не удовлетворяют одному пункту утверждения 2, а именно, они выдают вектор х, который лишь приближенно удовлетворяет балансовым ограничениям (в утверждении требование точного удовлетворения балансовых ограничений является существенным и не может быть как-то равнозначно релаксировано). Связанно это с тем, что для задачи ЭЛП, когда ограничений намного меньше числа прямых переменных, обычно решается двойственная задача, по которой восстанавливается решение прямой задачи [6,33]. Как следствие, приобретается невязка и в ограничениях. Собственно, в представлении градиента функционала по формуле (4) имеются два способа. Первый через двойственные множители (Х,р), второй через решение прямой задачи х. Функционал прямой задачи сильно выпуклый по х, поскольку энтропия 1-сильно выпуклая функция в 1-норме [15]. Поэтому сходимость в решении прямой задачи по функции обеспечивает сходимость и по аргументу, что и означает возможность определения c хорошей точностью градиента по формуле (4) через х. Другая ситуация возникает, если смотреть на двойственную задачу к задаче ЭЛП (в приводимом далее утверждении следует обозначить х = (Х,р), у = у для задачи (2), (3)).
Утверждение 3. Пусть <р(у) = тщФ(ж,у), где Ф(х,у) - такая достаточно гладкая, вы-
xeQ
пуклая по совокупности переменных функция, что10
\\ЧФ(х',у') - УФ(ж,у)||2 < L||(х',у') - (х,у)11 2.
Пусть для произвольного у £ Q (считаем, что множество Q содержит внутри себя шар радиуса более л/25/L) можно найти такой х = х(у) £ Q, что
тах(Фх(ж,у),х — z) ^ 5.
zeQ
Тогда
Ф(х,у) - у(у) < S,
10 Это утверждение имеет достаточно простую геометрическую интерпретацию. Проекция надграфика выпуклой функции будет выпуклым множеством, то есть, в свою очередь, надграфиком некоторой выпуклой функции. Кривизна границы у полученного при проектировании множества будет не больше, чем была у исходного множества. Это следует из того, что проектирование — сжимающий оператор.
\\Vfiy') — У^(у)\\2 < Ц\у' — у\\2,
и ^Ф(х,у) — 25, Фу(х,у)^ будет (65, 2Ь)-оракулом для <^(у) на выпуклом множестве, полученном из множества Q отступанием от границы во внутрь Q на расстояние л/25/Ь (по условию это множество не пусто).
Доказательство. По условию задачи имеем при всех допустимых значениях аргументов Ф:
Ап
( Ф XX Ф ху \ = 8ир А Фхх Ф ху А
у Фух Фуу у |М2<1\ ' Фух Фуу /
М < ь.
(12)
Заметим, что также по условию при всех допустимых значениях аргументов Ф: У 0, Фхх У 0, Фуу У 0
Фхх Фху Фух Фуу
Ф ух = Су, Ф хх = ^
Фуу = ФТуу.
Для упрощения последующих рассуждений (в частности, чтобы не работать с псевдообратными матрицами) будем, считать, что матрица Фхх У 0 положительно определена (исходя из условий, гарантировать можно лишь неотрицательную определенность). Также будем считать (в интересующих нас приложениях к задачам (2), (6) это имеет место), что зависимость х(у), определяемая из соотношения
1р(у) = шщФ(х,у) = Ф( х(у),у)
однозначным образом, и удовлетворяет соотношению
Фх{х(у),у) = 0,
у
из которого имеем
* ^л „Л дхШ , ^
-ху
Ф хх( х(у),у)
ду
+ Фху (х(у),у) = 0,
т.е.
\дх/ду \ \ = \ \ дх^дуу \ \ = — Ф-ХФху.
Поскольку !^(у) = Ф(х(у),у), то
ууу = \ \ дх/ду \ \ Т Фхх \ \ дх/ду \ \ + \ \ дх/ду \ \ Т Фху + Фух \ \ дх/ду \ \ + Фуу
= Ф уу — ФухФ-хФху.
С учетом этой формулы и из формулы дополнения по Шуру [34], получаем
Фхх Фху
Фух Ф у у
Ех 0 ФухФ-х Еу
Фх х
0Ф
0
у у
Ех Ф-хФху 0 Еу
где Ех, Еу - единичные матрицы соответствующих размеров. Поскольку
Ех 0 ФухФ-х Еу
Ех Ф-хФху
Еу
Т
и эти матрицы полного ранга, то из (12) имеем, что
8Ир (Н, 1£уу к) = \тах(^уу) ^
2<1
^ Шах{\тах(Фхх)Атах(<Руу)} =
( Фхх Фху \
^ Фух Фуу )
^ ь.
0
Таким образом, установлено, что
г 2
Р(У) < Р(х) + {Чф),у - х) + 2\\У - х||2. Согласно утверждению 1
Р(У) > Р(х) + {Фу (х, У-х) - &
Далее проведем рассуждения аналогично рассуждениям на с. 107 (и немного отлично от с. 115) диссертации [32]. Вычитая из первого неравенства второе, получим
г.. и2
(Фу(Х,у) -Ур(х),у-х) < 2\\У-XI2 +
Положим £ > 0:
= Фу(х,у) -Ур(х)
ЦФу(Х,У) -Чф)\\2 ,
получим
и т / \ „ , , м Ы 6 ||фу(х, У) ФЦ.2 < у + -.
Минимизируя правую часть неравенства по £ > 0, при получаем
II Фу {х, у) -V р(х) ||2 < ^25Ь. Отсюда и из утверждения 1 находим
г.. и2
р(у) < р(х) + (Ур(х),у -х) + 2\У -х||2 <
< р(х) + (Фу(х, у), у -х) + л/25Ъ\\у-х\\2 + г\.У -х||2 ^
< Ф{х,у) - 25 + (Фу(х, у), у -х) + л/25Ъ\\у - х\\2 + Ы ..У -х||2 + 26 <
2
< Ф{х, у) - 26 + (Фу(х, у), у -х) + Ы\\у - х^2 + 65. С учетом того, что (см. утверждение 1)
Р(У) > Р(х) + (Фу(х, У),У -х)
> Фу(х, У) - 25 + (Фу(х, у), у -х),
из определения 1 получаем доказываемое утверждение.
Это утверждение позволяет установить гладкость задачи (2), (3) (но не (5), (6)). Таким образом, для (5), (6) необходимость использования универсального метода для внешней задачи является отражением надежды сходиться быстрее, чем в негладком случае, в то время как для (2), (3) использование универсального метода для внешней задачи является скорее отражением желания настраиваться на правильную константу Липшица градиента. Можно, конечно, пытаться использовать приведенные выше формулы, однако из способа рассуждений (см., например, доказательство утверждения 3) видно, что полученная таким образом константа Липшица может оказаться завышенной.
К сожалению, практическое применение утверждения 3 натыкается на следующие сложности:
1) необходимости отступать от границы множества во внутрь на л/25/Ь,
2) необходимости рассмотрения ситуации (см. доказательство утверждения 3):
\9хг/д Уз \ = -Ф- Ф
ху,
3) необходимости предположения о компактности множества Q, иначе невозможно будет добиться выполнения условия
шах(Фх(ж,у),х — г) ^ 5. гея
Первая сложность на практике разрешима за счет возможности доопределения функционала задачи с сохранением всех свойств на л/2£/Ь-окрестность множества (заметим, что доопределение часто не требуется, поскольку функционал и так задан «с запасом»). Например, для рассматриваемых нами транспортных приложений с Q = {у : у ^ у} это возможно [2] - [5]. Сложность 2 часто вообще не возникает (разве что оговорка о существовании Ф-1 впрочем, приведенные выше рассуждения можно провести, сохранив все результаты в идентичном виде, так, что эта оговорка будет не нужна), поскольку Q совпадает со всем (двойственным) пространством. А вот сложность 3, действительно, портит дело. К сожалению, простых теоретически обоснованных способов борьбы с этой сложностью мы пока не знаем. Тем не менее полезно заметить, что в действительности нужно гарантировать выполнение (см. доказательство утверждения 1)
{Ф х (х(у) ,у),х(у) — х(у')) < 5,
где точки у и у' близки, поскольку возникают на соседних итерациях внешнего метода. С учетом ожидаемой «близости» х = х(у) и х(у), мы можем заменить в этом критерии настоящее множество Q, которое, как правило, совпадает со всем пространством, на шар конечного радиуса. Более детальные исследования (для задачи (2), (3)) и практические эксперименты показывают, что для выполнения приведенного выше условия достаточно обеспечить для внутреннего итерационного процесса {хк} ^ х(у) условия
I | Фх(Хк ,у) |1 2 | | Хк |1 2 < 5/2, | | Фх(Хк ,у) |1 2 < 5.
Соответствующее Хк = (Хк, Ц-к) порождает нужное х(у) = Хк. С учетом специфики рассматриваемой нами задачи (2), (3), имеем следующий критерий (возвращаемся к обозначениям
(2), (3)):
\\Ах(Хк,ц.к) — Ь\\2 ||(Ак,ц-к)\\2 ^ 5/2, \\Ах(Хк,ц.к) — Ь\\2 ^ 5,
где х(Хк, ц-к) определяется в формуле (4), а введённая линейная система балансовых уравнений Ах = Ь есть общая запись аффинных (транспортных) ограничений:
п п
^Жу = = Шу, г,] = 1, ...,п.
3=1 г=1
В связи со сказанным выше заметим, что (это следует из оценки (10)) метод балансировки обеспечивает сходимость и по аргументу, что для других методов (без введения регуляризации) решения двойственной задачи, вообще говоря, нельзя гарантировать. Это свойство наряду с линейной скоростью сходимости метода (со скоростью геометрической прогрессии) позволяет надеяться, что выбранный критерий является достаточно точным (точнее, не слишком грубым).
Принципиально важно для гладкого случая (сг¿(у) - функции с липшицевым градиентом), как это будет следовать из дальнейших оценок (см. теорему 1), не просто уметь решать двойственную задачу, т.е. находить (А, так, чтобы была сходимость по аргументу, а делать это так, чтобы сложность решения задачи зависела от точности ее решения логарифмическим образом. Выше мы отмечали, что это имеет место для метода балансировки. Также это имеет место и для быстрых методов, примененных к регуляризованной двойственной задачи. При фиксации параметра регуляризации, исходя из итоговой желаемой точности, быстрые градиентные методы (для сильно выпуклых функций) решают
регуляризованную двойственную задачу так, что зависимость сложности от точности ее решения - логарифмическая.
Хочется, чтобы при решении внешней задачи в (2), т.е. задачи
mmf(у), у^Я
можно было не задумываться ни о какой гладкости. Если она есть, то метод бы это хорошо учитывал, не требуя знания констант Липшица градиента (это намного более существенно для возможности применять описанный подход к поиску барицентра Вассерштейна вероятностных мер, см. п. 2), если ее нет, то метод также работал бы оптимальным (для негладкого случая) образом. Именно таким свойством и обладает универсальный метод [27], работающий и в концепции неточного оракула [31] (см. определение 1).
Заметим [27], что можно погрузить задачу с гёльдеровым градиентом (г £ [0,1])
/(у') -V/(у)\\* \
(в том числе и негладкую задачу с ограниченной нормой разности субградиентов при V = 0) в класс гладких задач с оракулом, характеризующимся точностью и
г = и
Гу (1 -V)
25(1+ г)
1-у
1+у
Это позволяет даже в случае, когда можно рассчитывать только на ¿-субградиент11 (с ограниченной нормой субградиента (разности субградиентов), причем, какой именно константой ограниченной, методу знать не обязательно), все равно работать в концепции (5,Ь)-оракула.
Итак, у нас есть внешняя задача (2)
ттf (у), у£ Я
для которой обращение к (5, Г)-оракулу за значением функции и градиента стоит ~ 1п($-1). Насколько быстро мы можем решить такую задачу, т.е. при каком N(е) можно гарантировать, что
ДУм(е)) - т^п1(У) < е? Ответ можно получить из следующего результата.
Теорема 1 (см. [1,21,27,31]). Существует однопараметрическое семейство универсальных градиентных методов (параметр р £ [0,1]), не получающих на вход, кроме р, больше никаких параметров (в частности, не использующих значения и К - «расстояние» от точки старта до решения, - априорно не известное), которое приводит к следующей оценке на требуемое число итераций:
N„(6) = О М
2
1+2ру+у
если 5 < 0(еN„(8)~р)
11На ¿-субградиент всегда можно рассчитывать согласно утверждению 1. Причем, как уже отмечалось раньше, для получения ¿-субградиента не нужна сходимость по аргументу для вспомогательной задачи.
Из теоремы 1 можно заключить, что если мы рассчитываем на некоторую гладкость /(у), то стоит выбирать значение параметра р = 1, при этом общие трудозатраты машинного времени будут
0(Ч(е)(Т 1п(е"1)+ Г)), (13)
где Т - время вычисления (суб-)градиента функционала (в основном это вычисления {^Су(у)}'ПуП=1 [29,30]), Т - время решения вспомогательной задачи методом балансировки с относительной точностью 1%. Численные эксперименты показывают, что на одном современном ноутбуке при п ~ 102 время Т ~ 1 с. [31], что сопоставимо с временем Т для таких п [29].
Выгода от описанной выше конструкции по сравнению с обычным способом решения исходной задачи минимизации (2), (3) сразу по совокупности всех переменных (см., например, [5]) заключается в гарантированном не увеличении константы Липшица градиента в оценке необходимого числа итераций (см. утверждение 3) и ожидаемое уменьшение в этой же оценке расстояния от точки старта до (неизвестного априори) решения. Выгода здесь вполне может достигать одного порядка и более. При этом можно ожидать лишь незначительного увеличения стоимости одной итерации. Причем стоит иметь в виду, что при оптимизации сразу по всем переменным требуется рассчитывать градиент функционала по большему числу переменных, чем в описанном выше подходе, что также играет нам на пользу. В конечном итоге, сокращение числа итераций заметно превалирует над небольшим увеличением стоимости одной итерации.
Что касается задач (5), (6), то описанный выше подход представляется естественным и не имеющим альтернатив в рассматриваемом классе. Альтернативные методы, с которыми можно сравнивать, мы упоминали по ходу статьи, но все они были предложены из принципиально других подходов. Сравнению (практическому) всех этих методов планируется посвятить отдельную публикацию.
Резюмируем ключевой результат этого раздела (и всей статьи) следующим образом.
Для решения задач типа (2), (6) или (8) предлагается использовать универсальный метод из работы [27] (а точнее его модификацию из [31]). Если рассчитываем на гладкость12 f (у), то полагаем в методе р = 1. Если на гладкость рассчитывать не приходится13, то полагаем р = 0. В обоих случаях, кроме априорной подсказки относительно параметра р, методу больше ничего от нас знать не надо!
Авторы выражают глубокую признательность Ю. Е. Нестерову за серию продуктивных обсуждений.
Исследование А. В. Гасникова, П.Е. Двуреченского, В. Г. Спокойного и А. Л. Сувориковой в части 2 выполнено в ИППИ РАН за счет гранта Российского научного фонда (проект № 14-50-00150); исследование П.Е. Двуреченского в части 3 выполнено при поддержке гранта РФФИ 15-31-20571-мол_а_вед; исследование А. В. Гасникова в части 3 выполнено при поддержке гранта РФФИ 15-31-70001 мол_а_мос.
Литература
1. Гасников А.В., Двуреченский П.Е., Камзолов Д.И., Нестеров Ю.Е., Спокойный В.Г., Стецюк П.И., Суворикова А.Л., Чернов А.В. Поиск равновесий в многостадийных транспортных моделях // Труды МФТИ. 2015. Т. 7, № 4. С. 143-155.
12В этом случае как раз существенна логарифмическая сложность приближенного вычисления градиента и значения функции /(у) от точности, обеспеченная методом балансировки.
13В этом случае точность решения вспомогательной задачи расчета ¿-субградиента можно завязать на желаемую точность решения задачи (2) е (или (6) или (8)) по формуле 6 = О(е) (см. теорему 1 при и = 0), с константой порядка 1.
2. Гасников А.В., Дорн Ю.В., Нестеров Ю.Е, Шпирко С.В. О трехстадийной версии модели стационарной динамики транспортных потоков // Математическое моделирование. 2014. Т. 26:6. C. 34-70. arXiv:1405.7630
3. Гасников А.В. Об эффективной вычислимости конкурентных равновесий в транспортно-экономических моделях // Математическое моделирование. 2015. Т. 27, № 12. С.121-136. arXiv:1410.3123
4. Бабичева Т.С., Гасников А.В., Лагуновская А.А., Мендель М.А. Двухстадийная модель равновесного распределения транспортных потоков // Труды МФТИ 2015. Т. 7, № 3. С. 31-41. https://mipt.ru/upload/medialibrary/971/31-41.pdf
5. Гасников А.В., Гасникова Е.В., Мациевский С.В., Усик И.В. О связи моделей дискретного выбора с разномасштабными по времени популяционными играми загрузок // Труды МФТИ. 2015. Т. 7, № 4. С. 129-142. arXiv:1511.02390
6. Гасников А.В., Гасникова Е.В., Нестеров Ю.Е., Чернов А.В. Об эффективных численных методах решения задач энтропийно-линейного программирования // ЖВМ и МФ. 2016. Т. 56, № 4. С. 523-534. arXiv:1410.7719
7. Agueh M., Carlier G. Barycenters in the Wasserstein space // SIAM J. Math. Anal. 2011. V. 43, N 2. P. 904-924.
8. Cuturi M., Doucet A. Fast Computation of Wasserstein Barycenters // ICML. 2014. http://www.iip.ist.i.kyoto-u.ac.jp/member/cuturi/
9. Benamou J.D., Carlier G., Cuturi M., Nenna L., Peyre G. Iterative Bregman Projections for Regularized Transportation Problems // e-print, 2015. (to appear in SISC) arXiv:1412.5154
10. Cuturi M., Peyre G., Rolet A. A Smoothed Dual Formulation for Variational Wasserstein Problems // e-print, 2015. arXiv:1503.02533
11. Cuturi M. Sinkhorn Distances: Lightspeed Computation of Optimal Transport // NIPS. 2013.
12. Немировский А.С., Юдин Д.Б. Сложность задач и эффективность методов оптимизации. М.: Наука, 1979. http://www2.isye.gatech.edu/simnemirovs/Lect_EMCO.pdf
13. Boissard E., Le Gouic T., Loubes J.-M. Distribution's Template Estimate with Wasserstein Metrics // e-print, 2013. (to be published in Bernoulli) arXiv:1111.5927
14. Bigot J., Klein T. Consistent estimation of a population barycenter in the Wasserstein space // e-print, 2015. arXiv:1212.2562
15. Nesterov Y. Smooth minimization of nonsmooth function // Math. Program. Ser. A. 2005. V. 103, N 1. P. 127-152.
16. Fercoq O., Richtarik P. Accelerated, Parallel and Proximal Coordinate Descent // e-print, 2013. arXiv:1312.5799
17. Qu Z., Richtarik P. Coordinate Descent with Arbitrary Sampling I: Algorithms and Complexity // e-print, 2014. arXiv:1412.8060
18. Boyd S., Parikh N., Chu E., Peleato B., Eckstein J. Distributed optimization and statistical learning via the alternating direction method of multipliers // Foundations and Trends in Machine Learning. 2011. V. 3(1). P. 1-122. http://stanford.edu/boyd/papers.html
19. Boyd S., Vandenberghe L. Convex optimization. Cambridge University Press, 2004. http://stanford.edu/boyd/cvxbook/
20. Поляк Б.Т. Введение в оптимизацию. М.: Наука, 1983.
21. Nemirovski A. Lectures on modern convex optimization analysis, algorithms, and engineering applications. Philadelphia: SIAM, 2013.
http://www2.isye.gatech.edu/simnemirovs/Lect_ModConvOpt.pdf
22. Rabin J., Peyer G., Delon J. Bernot M. Wasserstein barycenter and its applications to texture mixing // LNCS. 2011. Proc. SSVM'11. Springer. V. 6667. P. 435-446. https://hal.archives-ouvertes.fr/hal-00476064/document
23. Bonnel N., Pfister H. Sliced Wasserstein barycenter of multiple densities // Harvard Technical Report. 2013. TR-02-13.
ftp://ftp.deas.harvard.edu/techreports/tr-02-13.pdf
24. Стецюк П.И. Методы эллипсоидов и r-алгоритмы. Кишинев: Эврика, 2014.
25. Стецюк П.И., Гасников А.В. NLP-программы и r-алгоритм в задаче энтропийно-линейного программирования // Теория оптимальных решений. Киев: Институт кибернетики им. В.М.Глушкова НАН Украины, 2015. С. 73-78.
26. Franklin J., Lorenz J. On the scaling of multidimensional matrices // Linear Algebra and its applications. 1989. V. 114. P. 717-735.
27. Nesterov Yu. Universal gradient methods for convex optimization problems // CORE Discussion Paper 2013/63. 2013; Math. Program. Ser. A. 2015. V. 152. P. 381-404. https://www.uclouvain.be/cps/ucl/doc/core/documents/coredp2013_26web.pdf
28. Nesterov Yu. http://www.youtube.com/watch?v=Fm9h92pcbvg
29. Гасников А.В., Двуреченский П.Е., Дорн Ю.В., Максимов Ю.В. Численные методы поиска равновесного распределения потоков в моделях Бэкмана и стабильной динамики // Математическое моделирование. 2016. Т. 28, № 10. С. 40-64. arXiv:1506.00293
30. Гасников А.В., Гасникова Е.В., Ершов Е.И., Двуреченский П.Е., Лагуновская А.А. Поиск стохастических равновесий в моделях равновесного распределения потоков // Труды МФТИ. 2015. Т. 7, № 4. С. 114-128. arXiv:1505.07492
31. Гасников А.В, Двуреченский П.Е., Камзолов Д.И. Градиентные и прямые методы с неточным оракулом для задач стохастической оптимизации // Динамика систем и процессы управления. Труды Международной конференции, посвященной 90-летию со дня рождения академика Н.Н. Красовского. Екатеринбург, 15-20 сентября 2014. Издательство: Институт математики и механики УрО РАН им. Н.Н. Красовского (Екатеринбург). 2014. С. 111-117. arXiv:1502.06259
32. Devolder O. Exactness, inexactness and stochasticity in first-order methods for large-scale convex optimization. CORE UCL, PhD thesis, March 2013.
33. Anikin A., Dvurechensky P., Gasnikov A, Golov A., Gornov A., Maximov Yu., Mendel M., Spokoiny V. Modern efficient numerical approaches to regularized regression problems in application to traffic demands matrix calculation from link loads // Proceedings of International inference ITAS-2015. Russia, Sochi, September, 2015. arXiv:1508.00858
34. Zhang F. The Schur Complement and Its Applications. Springer, 2005.
References
1. Gasnikov A., Dvurechensky P., Kamzolov D., Nesterov Y., Spokoiny V., Stetsyuk P., Suvorikova A, Chernov A. Search for the stochastic equilibria in the transport models of equilibrium flow distribution. Proceedings of MIPT. 2015. V. 7, N 4. P. 143-155. (in Russian)
2. Gasnikov A., Dorn Yu., Nesterov Yu., Shpirko S. On the three-stage version of stable dynamic model. Matem. Mod. 2014. V. 26:6. P. 34-70. arXiv:1405.7630 (in Russian)
3. Gasnikov A. About reduction of searching competetive equillibrium to the minimax problem in application to different network problems. Mathematical modelling. 2015. V. 27, N 12. P. 121-136. (in Russian) arXiv:1410.3123
4. Babicheva T., Gasnikov A., Lagunovskaya A, Mendel M. Two-stage model of equilibrium distributions of traffic flows. 2015. V. 7, N 3. P. 31-41. (in Russian) https://mipt.ru/upload/medialibrary/971/31-41.pdf
5. Gasnikov A., Gasnikova E., Matsievsky S., Usik I. Searching of equilibriums in hierarchical congestion population games. Proceedings of MIPT. 2015. V. 7, N 4. P. 129-142. (in Russian) arXiv:1511.02390
6. Gasnikov A., Gasnikova E., Nesterov Y., Chernov A. Entropy linear programming. Computational Mathematics and Mathematical Physics. 2016. V. 56, N 4. P. 523-534. (In Russian). arXiv:1410.7719
7. Agueh M., Carlier G. Barycenters in the Wasserstein space. SIAM J. Math. Anal. 2011. V. 43. N 2. P. 904-924.
8. Cuturi M., Doucet A. Fast Computation of Wasserstein Barycenters. ICML, 2014. http://www.iip.ist.i.kyoto-u.ac.jp/member/cuturi/
9. Benamou J.D., Carlier G., Cuturi M., Nenna L., Peyre G. Iterative Bregman Projections for Regularized Transportation Problems. e-print, 2015. (to appear in SISC). arXiv:1412.5154
10. Cuturi M., Peyre G., Rolet A. A Smoothed Dual Formulation for Variational Wasserstein Problems. e-print, 2015. arXiv:1503.02533
11. Cuturi M. Sinkhorn Distances: Lightspeed Computation of Optimal Transport. NIPS. 2013.
12. Nemirovsky A., Yudin D. Problem Complexity and Method Efficiency in Optimization. SIAM Review 27.2 (1985): 264.
13. Boissard E., Le Gouic T., Loubes J.-M. Distribution's template estimate with Wasserstein metrics. e-print, 2013. (to be published in Bernoulli). arXiv:1111.5927
14. Bigot J., Klein T. Consistent estimation of a population barycenter in the Wasserstein space. e-print, 2015. arXiv:1212.2562
15. Nesterov Y. Smooth minimization of nonsmooth function. Math. Program. Ser. A. 2005. V. 103. N 1. P. 127-152.
16. Fercoq O., Richtarik P. Accelerated, Parallel and Proximal Coordinate Descent. e-print, 2013. arXiv:1312.5799
17. Qu Z., Richtarik P. Coordinate Descent with Arbitrary Sampling I: Algorithms and Complexity. e-print, 2014. arXiv:1412.8060
18. Boyd S., Parikh N., Chu E., Peleato B., Eckstein J. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning. 2011. V. 3(1). P. 1-122. http://stanford.edu/boyd/papers.html
19. Boyd S., Vandenberghe L. Convex optimization. Cambridge University Press, 2004. http://stanford.edu/boyd/cvxbook/
20. Polyak B. Introduction to optimization. New York: Optimization Software, 1987.
21. Nemirovski A. Lectures on modern convex optimization analysis, algorithms, and engineering applications. Philadelphia: SIAM, 2013.
http://www2.isye.gatech.edu/simnemirovs/Lect_ModConvOpt.pdf
22. Rabin J., Peyer G., Delon J. Bernot M. Wasserstein barycenter and its applications to texture mixing. LNCS. 2011. Proc. SSVM'11. Springer. V. 6667. P. 435-446. https://hal.archives-ouvertes.fr/hal-00476064/document
23. Bonnel N., Pfister H. Sliced Wasserstein barycenter of multiple densities. Harvard Technical Report. 2013. TR-02-13.
ftp://ftp.deas.harvard.edu/techreports/tr-02-13.pdf
24. Stetsyuk P. Ellipsoid methods and r-algorithms. Kishinev: Evrika. 2014. (in Russian)
25. Stetsyuk P., Gasnikov A. NLP-programms and r-algorithms in linear-entropic programming. Optimal Decision Theory. Kiev: Glushkov Institute of Cybernetics NASU, 2015. P. 73-78.
26. Franklin J., Lorenz J. On the scaling of multidimensional matrices. Linear Algebra and its applications. 1989. V. 114. P. 717-735.
27. Nesterov Yu. Universal gradient methods for convex optimization problems. CORE Discussion Paper 2013/63. 2013; Math. Program. Ser. A. 2015. V. 152 P. 381-404. https://www.uclouvain.be/cps/ucl/doc/core/documents/coredp2013_26web.pdf
28. Nesterov Yu. http://www.youtube.com/watch?v=Fm9h92pcbvg
29. Gasnikov A., Dvurechensky P., Dorn Y., Maksimov Y. Searching equillibriums in Beckmann's and Nesterov-de Palma's models. Mathematical modelling. 2016. V. 28, N 10. P. 40-60. (in Russian). arXiv:1506.00293
30. Gasnikov A., Gasnikova E., Ershov E., Dvurechensky P., Lagunovskaya A. Search for the stochastic equilibria in the transport models of equilibrium flow distribution. Proceedings of MIPT. 2015. V. 7, N 4. C. 114-128. arXiv:1505.07492
31. Gasnikov A., Dvurechensky P., Kamzolov D. Gradient and direct methods with inexact oracle in stochastic optimization. Proceedings of the International Conference «Systems Dynamics and Control Processes» (SDCP'2014) dedicated to the 90th Anniversary of Academician Nikolay Nikolayevich Krasovskii, P. 111-117. Ekaterinburg, 2014. (in Russian) arXiv:1502.06259
32. Devolder O. Exactness, inexactness and stochasticity in first-order methods for large-scale convex optimization. CORE UCL, PhD thesis, March 2013.
33. Anikin A., Dvurechensky P., Gasnikov A., Golov A., Gornov A., Maximov Yu., Mendel M., Spokoiny V. Modern efficient numerical approaches to regularized regression problems in application to traffic demands matrix calculation from link loads. Proceedings of International Conference ITAS-2015. Russia, Sochi, September, 2015. arXiv:1508.00858
34. Zhang F. The Schur Complement and Its Applications. Springer, 2005.
Поступила в редакцию 09.03.2016
