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

min

E 7 =i ^ i₃ = L, E ? =i ^ i₃ = W j ^x ij > 0, i, j = 1,..., n

П max< y^Q I i,j = 1

X ij ln X ij +

П f Cij (VWij + g(y)^ = i,j=1

= max max yeQ х.ргу—

{

( A, L) + (ц, W}- ^f exp( -_{С <7}- (y) - 1 + A i + i,j=1

ц ) ⁺ g⁽y)},

где g(y) и C ij (y) >  0 — вогнутые гладкие функции (если ищутся не стохастические равновесия, то C ij (y) могут быть негладкими), Q - множество простой структуры, например,

Q = {y Е R ⁿ : y > у}.

Легко понять, что система балансовых ограничений в (1) либо несовместна п        п

X ^ і = 52 W j , либо вырождена (имеет не полный ранг). В последнем случае это приво- і=1        з=1

дит к тому, что двойственные переменные (А, м) определены с точностью до произвольной постоянной С :

(А + Се,ц-Се), е = (1,..., 1).

^ '- 1 у п

Задачу (1) также можно переписать следующим образом (не ограничивая общности, п        п считаем 52 ^ =52 Wj = 1):

і=1      з=1

min п             п

£ xгj = L^, £ a^j = Wj j = 1             г = 1

п xtj > 0, £ x^j = 1,i,j = 1, ..., n

п max< yGQ I .^ і,3 = 1

Хц ln Х із +

п

^ С із ⁽У) ^х із ⁺ 9^ = і,3=1

п

= ^maQ A ^maK -{^(А ,^ ^{+ (} ^ ^,W ^ - ^lnl ^ ^ехР (^—С із Ы ^{+ А} і ⁺ M )] ⁺ 9⁽у^ =

У Q ,ц                                  і,з = 1

= - mmf (у), v^Q где выпуклая функция /(у) определяется как п               п

^{/ (}У) =

max п            п

£ x^j = L£, £ x^j = Wj j = 1            г = 1

п,п xгj > 0, Й xгj = 1,i,j = 1,..., n

{ ^— ^ Х із ^ln Х із ^— ^ ^С із ⁽у)Х із ^— 9⁽у)} =

• і,з = 1                і,з = 1

п

_A ^miⁿ_n{^ln [^ exp^{( —С} із ⁽у) ^{+ А} і ⁺^ j )]^-(^A, L ) ^- ( м, ^W) ^- 9⁽у)}.            ⁽³⁾

,^           і,з=1

п

Поскольку мы добавили в ограничения условие X Х із = 1, являющееся следствием і,з=1

балансовых уравнений, то это привело к тому, что двойственные переменные (А, м) определены с точностью до двух произвольных постоянных С \ , С ^ : (А + С \ е, м + С ^ е).

Заметим, что расчет градиента V /(у) (в ряде транспортных приложений вогнутые функции С із (у) — негладкие, тогда вместо градиентов стоит понимать суперградиенты С із (у) и субградиент / (у)) осуществляется по следующей формуле (Демьянова-Данскина-Рубинова, см., например, [2, 6]):

І2 exp( — С із ( у ) + А * + M * ) V c із (у)

V /(у) = — ^і-,^з=-п --V9(у) =

X ^exP ( ^{— С} із ⁽у^{) + А} * ⁺м*)

і,з=1

п

= ^— ^ ^Х із ^(А*, M*^{) Vc} із ⁽у) ^{— V} 9⁽у^{) ,                                (4)}

где (А*,м*) — решение задачи (3), не важно какое именно, градиент V/(у) от выбора С\, Сц (см. выше) не зависит. В данной статье мы (так же, как ив [1]) ограничимся изучением только полноградиентных методов для задачи (2), т.е. не будем рассматривать, например, рандомизацию при вычислении градиента по формуле (4). Планируется отдельно исследовать вопрос о возможности ускорения вычислений за счет введения рандомизации для внешней задачи. На данный момент нам представляется (см. формулу (13) в п. 3), что это может принести дивиденды только в случае, когда вспомогательная задача расчета V Cij (у) достаточно сложная. Тут требуется много оговорок, в частности, в большинстве приложений умение рассчитывать Vcij (у) для конкретной пары (i,j) без дополнительных затрат позволяет заодно рассчитать и все Vc^j(у), j = 1,..., п. Также отдельно планируется исследовать вопрос о том, какие подходы и насколько хорошо допускают распараллеливание. Вопрос о целесообразности рандомизации оказывается завязанным и на вопрос о возможности распараллеливания.

Структура статьи следующая. В п. 2 мы рассматриваем популярную в последнее время (в связи с большим числом приложений) задачу вычисления барицентра Вассерштейна различных вероятностных мер [7]. Эта задача оказывается тесно связанной с задачей (1). Мы разбираем в статье этот пример, потому что он хорошо проясняет возможные альтернативы предлагаемому нами основному подходу решения задач (1), (2), изложенному в п. 3. В основе подхода п. 3 (см. также [1, 8]) лежит сочетание метода балансировки для решения внутренней задачи оптимизации по двойственным множителям и универсального градиентного метода с неточным оракулом для внешней задачи (2). В работе [1] мы сконцентрировались на обобщении универсального метода на случай неточного оракула. В данной работе акцентируем внимание на получении условий на шум оракула для внешней задачи, исходя из оценок точности решения внутренней задачи. Таким образом, данная статья дополняет работу [1] в теоретическом плане. Практическим экспериментам планируется посвятить отдельную работу.

• 2.    Поиск барицентра Вассерштейна

К похожей на (1) задаче приводит поиск барицентра Монжа–Канторовича (в западной литературе чаще говорят барицентра Вассерштейна [7, 8]) ¹ . Изложим вкратце постановку задачи [8] – [10]. Вводится энтропийно регуляризованное транспортное расстояние (см. рис. 1 в [11]) с матрицей H c ^j H „jⁿ , сформированной из квадратов попарных расстояний C ij = I 2j от носителя i меры L до носителя j меры W (L, W Е £ „ (1)):

A(L,W ) =

min п            п

£ x ^j = L ^ , £ x ^j = W j

7 = 1             г = 1

X ij >  0, г, j = 1, ..., n

П

{7 £ $ ij lnx ij i,j=1

n

⁺ £ C ij X ij^ = i,j=1

{ n

(A,L) + (M,W)- 7 £ exp(- i,j=1

^C ij ^{+ A} i ⁺ ^ j

- 1

{n„ max

Ae R ⁿ

< A ^, L )— 7 £W , ln[ wW £ exp( -2^ )] .

j=1            j i=1             ')

^..

h w (A)

Определим при L Е £ _n (1) функцию Н\ ү (L) = A(L, W). Эта гладкая на L Е £ _n (1) функция с градиентом (см. утверждение 3 [10]):

VHW (L) = A*, где А* — единственное решение (5), удовлетворяющее условию2 (А*, е) = 0. Отсюда следует, что

H W (А) = max |(A,L ) - H _w (L)} =

L e s n (1) n             ⁿ

= 7£ »', ln[± £_exp( - )].

J=1           J 1=1            '

Теперь можно перейти к изложению основной конструкции. Задача поиска барицентра Вассерштейна ³ записывается следующим образом:

m

£ ^H w _k ^{( L )}

k=i

^ min . LeS n (i)

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

Перепишем задачу (6), следуя п. 3 работы [10], следующим образом:

^- £^H W k ^{( L} k ) ^

m max

L i = L _m | A 1 ,L 1 e S _n (1)

k=1 ...

L _m - 1 = L_m| A ^m-1 ,L _m e S „ (1)

m—1

£ max    ( А ^к ,L _k ) - H w_k (L k ) ■ + max

L _k GS n (1)l                          )    L _m e 8 n (1)

k=1

m—1

I ^{( —} £ ^А ^ , ^L m ) ^{- H} W m ^{( L} m ) I

k=1

^ min , A 1 ,...,A ^m-1

m—1m

£ HWk

k=1                              k=1

L* = ^HW_к(А^к) для любого к = 1,....,m — 1, где L* - единственное решение задачи (6), {А^к}k=1¹— единственное решение задачи (7). Важное свойство функционала задачи (7) -равномерная ограниченность константы Липшица градиента (следует из [15]). Задача безусловной минимизации (7) может быть эффективно решена различными способами (в зависимости от того, насколько велики п и m). В частности, для больших п и m неплохо с задачей справляются различные модификации метода сопряженных градиентов и быстрых градиентных методов [10]. Структура задачи (7) позволяет эффективно использовать покомпонентные методы (см., например, [16, 17]), которые к тому же хорошо параллелятся для данной задачи. Задача (7) хорошо также решается с помощью распределенных вычислений [18].

В приложениях к поиску разладки требуется много раз перерешивать задачу (7), которую для симметричности перепишем следующим образом:

m

£HWk(Ак) ^ mmin к=1             Е Ak=0

k=1

немного смещая окно, т.е. заменяя каждый раз несколько первых слагаемых в сумме

HW1 (A^,..,HW„ (A^r)

на столько же новых (которые, как ожидается, близки к HW (Am)). В таком случае предлагается в итерационном процессе стартовать при сдвиге окошка с того, на чем остановились на прошлом положении окошка, экстраполируя Am на вновь пришедшие слагаемые. Ясно, что для новой задачи набор fx г+1    \ т— 1 \т \т \т \

^(A* , ..., ^A*     , ^A* , ^A* , ..., ^A* ^.),

⏟⏞

г с которого стартуем, уже не будет оптимальным, однако мы вправе надеяться на его близость к оптимальному набору, что существенно сокращает число последующих итераций. Интересной, особенно в данном контексте, представляется возможность использования (и интерпретации) распределенных вычислений [18].

Может показаться, что подход, сводящий поиск решения задачи (6) к задаче (7), не доминируем, поскольку, в отличие от задачи (6), в задаче (7) мы можем явно выписать функционал и по простым формулам, рассчитать градиент, который к тому же имеет равномерно ограниченную константу Липшица. С одной стороны, это, действительно, преимущество, но получено оно дорогой ценой – ценой раздутия пространства, в котором происходит оптимизация почти в т раз. И это раздутие скажется не только на сложности одной итерации. Для задачи (6) осуществление одной итерации будет еще более дорогим ввиду необходимости на каждой итерации решать т отдельных подзадач расчета VHw_k (+ Скажется это, прежде всего, на числе необходимых итераций. В следующем разделе будет отмечено, что расчет VHw_k +) с помощью метода балансировки не намного сложнее расчета ^HW_k(A^k). При этом задача (6) решается в пространстве намного меньшей размерности, и мы вправе ожидать, что необходимое число итераций может быть намного меньше, чем для задачи (7). Кроме того, задача (6) решается на компакте, т.е. размер решения (если быть точным, то расстояние от точки старта до решения), входящий в оценку необходимого числа итераций, заведомо ограничен размером симплекса. Задача (7) – задача безусловной оптимизации, причем без свойства сильной выпуклости функционала. Размер ее решения может быть большим, и входит он в оценки необходимого числа итераций так же, как и для задачи (6), к сожалению, степенным образом (для быстрых (ускоренных) методов можно ожидать линейной зависимости необходимого числа итераций от этого размера). Наконец, для постановок задач об обнаружении разладки (см. выше) также ожидается, что использовать близость решений прямых задач (6) при смещении окошка удастся намного лучше, чем близость в решении двойственных задач (7). В итоге выгода от подхода, связанного с переходом к задаче (7), уже не столь очевидна и требует отдельного и более аккуратного исследования с решающей ролью численных экспериментов.

В ряде задач требуется искать параметрический барицентр Вассерштейна. В таком случае в одном из вариантов постановки предполагают наличие параметрической зависимости +#) G + (1), Ө Е Ө, где размерность вектора параметров dimӨ ^ п. К сожалению, в этом случае нельзя гарантировать с помощью стандартных приемов [19, c. 86] выпуклости за- дачи

т

^^HW_k +^(Ө)) ^ mn, к=1

за исключением случая, когда ^(Ө) = АӨ + b, а Ө — выпуклое множество. В этом случае конструкция (7) видоизменяется следующим образом:

т

^-^ ^HWk (L_k) ^ к=1

max

Li = L’ A¹, Li е S„(1) • ••

L-i = Lm| A-’ , Lm—i е S„(1) Lm = A9 + ь| A, L- е S„(1), 0 е ө m—1                   m—1

^ ^HW_k (^^fc) ^{+ H}W_m -- ^ ^X - ^A) + (A,^^ + I) ^    шт- _ ,        (9)

k=1                          k=1                                    ”'e е Ө ’

L* = VHW_k(X*) для любого к = 1,....,m — 1. Как следствие, нет никаких гарантий, что изложенная выше конструкция, связанная с переходом к двойственной задаче (7) и восстановлением решения прямой задачи (6) по явным формулам через двойственные множители, в общем случае будет работать хотя бы для поиска локальных решений.⁴Другими словами, необходимо искать глобальный оптимум задачи (8) исходя из работы с прямой задачей (8). Один из вариантов того, как это можно делать, будет описан в следующем пункте.⁵

Однако при другом варианте постановки (более предпочтительном) можно задавать зависимость L(6) с помощью аффинных равенств и выпуклых неравенств:

m

^^HW_k^(L) ^ к=1

min .

А6 + BL = с 9(6, L) 6 0 L е S„(1); в е ө

Многие параметрические зависимости представимы в таком виде [21]. В частности, отметим полиэдральные представления Фурье–Моцкина [21], возникающие, например, в робастной оптимизации:

m

52 ^Hw_k^(L) ^   , ^min ,.

_k=1              Le{LgS„(1): 3 9: A0+BL6c^

Можно переписать задачи (7), (9) и на эти случаи, причем сделать это корректно в том смысле, что правомочность подхода полностью сохранится. При этом принципиально ничего из сказанного выше не поменяется. Подробнее об этом планируется написать в отдельной работе.

В действительности, в приложениях наиболее интересен случай, когда ищется барицентр именно расстояний Вассерштейна,6 а не энтропийно-регуляризованных расстояний [8] - [11]. Другими словами, интересно изучать предельное поведение 7 ^ 0+ (см. п. 3.1 [9], утверждение 1 [10], п. 3 и конец п. 4 [11]). К сожалению, методы из пп. 2, 3 оказываются весьма чувствительными к этому предельному переходу. Для метода из этого раздела константа Липшица градиента в задаче (7) будет расти как 7-1, соответственно, число итераций будет увеличиваться (при использовании быстрых (ускоренных) методов) как 7-1/2. Еще более плохое поведение (см. [6]) можно ожидать от метода балансировки, использующегося в подходе из п. 3. Планируется в отдельной публикации исследовать вопрос о том, как следует действовать при малых 7 > 0. В частности, в вырожденном случае 7 = 0. По-видимому, в этом случае поможет философия искусственного сглаживания 7 [15], в которой искусственно введенная энтропийная регуляризация уже задается с четко заданным коэффициентом регуляризации 7 > 0, зависящим от итоговой точности по функции, с которой требуется решить задачу. Другой способ – использовать менее чувствительные (чем метод балансировки) способы решения двойственной задачи, например, при небольших значениях п ожидается, что лучше сработает г-алгоритм Н. З. Шора и некоторые его обобщения [24,25]. В данной работе мы фиксируем 7 > 0 и далее уже не будем возвращаться к подобного рода вопросам.

В заключение отметим, что поиск барицентра Вассерштейна в случае т = 1 может быть осуществлен явно: L = W. Обоснование этого частного результата представляется довольно полезным для понимания основной конструкции этого раздела.

• 3.    Универсальный метод с неточным оракулом

Из п. 2 следует, что внутренняя задача максимизации по (А, р) может быть явно решена по р при фиксированном А, и наоборот (это верно для задач (1) и (2) и приводит к одним и тем же формулам). Собственно, таким образом, получается метод балансировки расчета матрицы корреспонденций по энтропийной модели, см., например, [6] (тесно связанный с методом Синхорна [9, 11]), как метод простой итерации для явно выписываемых условий экстремума (принципа Ферма): А = Л(р), р = М(А). Метод балансировки имеет вид ([А]о = [р]о = 0):

1 ^п

^[Аг^]&+1

^[pj h+1

^-ln[l E^exp(^-cЫ ^{- 1 + [}Pjh)],

^{— 1} + ^[Ai^]fc)

1 ⁿ

^{- ln} [^ E ^exP (^-c(у)

или

^[pj^]fe+1

1     ¹¹

^{- ln} [^j E ^exp(^-co

^{— 1} + ^[Аг^]fe+1) .

В этих формулах « —1» в экспоненте для метода (2) (в отличие от метода (1)) можно не писать, поскольку двойственные множители определяются неоднозначным образом с бoльшим произволом для задачи (2) (см. выше), достаточным для справедливости этого замечания.

Оператор (А, р) ^ (Л(р), М(А)) является сжимающим в метрике Биркгофа-Гильберта р [26]. Это означает, что после N ~ ln(^ ¹) итераций метода балансировки можно получить такие (An, pn), что { (А*(у), р*(у) } - двумерное аффинное множество решений (см. п. 1):

р((^an, Pn); {(^А*(у),Д*(у))}) 6^. ⁽¹⁰⁾

Причем на практике наблюдается очень быстрая сходимость, т.е. коэффициент пропорциональности не большой [6]. Таким образом, мы можем приближенно решить внутреннюю задачу.⁸

Далее предлагается воспользоваться прямодвойственным (эта важно, поскольку нужно восстанавливать двойственные переменные) универсальным методом [27] для решения внешней задачи оптимизации по у (имеется видеопрезентация с описанием этого метода [28]). К сожалению, в формулировке (1) (в отличие от формулировки (2)) кроме того, что внешняя задача гладкая (при условии гладкости Oj (у) [29,30]), больше ничего о ней сказать нельзя (константа Липшица градиента не ограничена). Также не понятна гладкость задачи (6). Поэтому и по ряду других причин, о которых будет сказано далее, было отдано предпочтение универсальному методу, оптимально адаптивно настраивающемуся на гладкость функционала f (у) на текущем участке пребывания итерационного процесса.9 Однако нам потребуется использовать этот метод в варианте с неточным оракулом, выдающим градиент [31]. Напомним (см. п. 1), что мы решаем задачу 2, представимую в виде (здесь в max представлении ж = ж,

Q = ^ж_гj > 0, г,у = 1,..., п

п            п

52 ^ = ^^^j = ^Wj j=1           г=1

,

а в min представлении ж = (А, ц), Q = R²ⁿ):

f (у) = max Ф(ж, у) = min Ф(ж, у) ^ min у Е Q. xEQ          xEQ

Далее везде будем считать, что у Е Q.

Определение 1 (см. главу 4 [32]). (5, Ғ)-оракул выдает (на запрос, в котором указывается только одна точка у) такие (F(у),С(у)), что и для любых у, у’ Е Q

0 6f (у') -F (у)-(О(у),у' -у) 6 ^ ^у' - у^²+ 5.

Из определения 1 cразу следует, что для любого ж Е Q

F(ж) 6 f (ж) 6 F(ж) + 5

и для любых ж, у Е Q f (у) > f (у) - Мж), у -х} - 5

Из последнего свойства получаем, что определение (5, Ғ)-оракула можно понимать как обобщение на гладкие задачи классического понятия негладкой выпуклой оптимизации: 5-субградиента (см. п. 5 § 1 главы 5 [20]). В приводимом далее утверждении в первой его части следует сохранить обозначения для задачи (2), (3) и следует обозначить ж = А, у = L для задачи (5), (6); а во второй части утверждения следует обозначить ж = (А, ц), у = у для задачи (2), (3). Таким образом, на задачу (2), (3) можно посмотреть с двух разных ракурсов, однако второй ракурс менее привлекателен ввиду необходимости рассмотрения ограниченных множеств Q, что в интересующих нас приложениях место не имеет.

Утверждение 1. Если ^(у) = тахФ(ж,у), где Ф(ж,у) - выпуклая по у и вогнутая по ж xEQ функция, и найден такой ж G Q, что

^(у) - Ф(ж, у) 6 5, то субградиент ду Ф(ж,у) есть 5-субградиент функции ф(у) в точке у. Если <р(у) = тіпФ(ж,у), где Ф(ж,у) - выпуклая по совокупности переменных функция, и xEQ найден такой ж G Q, что тах(Фх(ж,у),ж — г) 6 5,

ZEQ то

Ф(ж,у) — у(у) 6 5

и субградиент Ф_у(ж, у) = Э_уФ(ж, у) есть 5-субградиент функции у(у) в точке у.

Доказательство. Ограничимся доказательством только второй части этого утверждения. Доказательство первой части см. на с. 124 (лемма 13) книги [20]. Из выпуклости Ф(ж,у) по совокупности переменных имеем

^Ф(ж‘,у’) > ^Ф(ж,у) + (^Фх^(ж,у^),ж‘ ^{— ж}) + ^Ф,^(ж,у),у’^—у).                ⁽¹¹⁾

Определим зависимость ж(у) из соотношения

У (у) = тіпФ(ж,у) = Ф (ж(у),у).

xEQ

Заметим, что Ф(ж,у) > у(у). Положим в (11) ж’ = ж(у’), ж = ж. Тогда

У(у’) = Ф(ж‘, у’) > Ф(ж, у) + (Фх(ж, у),ж‘ — ж) + (Ф_у(ж, у), у’ — у) > > ^Ы ⁺ (^Фх^(ж, у^),ж(у'^{) — ж}) ^{+ Ф},^(ж, у), у’ ^—у) > > ^Ы ⁺ (^Фу^(ж, у),у‘ ^—у) ^{— 5}.

В последней формуле мы использовали, что

(Фх(ж,у),ж — ж(у’)) 6 5.

В свою очередь из выпуклости Ф(ж,у) по ж (для всех допустимых у) имеем

Ф(ж,у) — Ф(ж(у‘),у) 6 (Фх(ж,у),ж — ж(у‘)).

Беря в этой формуле у’ = у и воспользовавшись определением ж(у), получаем, что

Ф(ж, у) — у(у) 6 (Фх(ж, у),ж — ж(у)У

Однако не хочется довольствоваться возможностью находить только 5-субградиент (из утверждения 1 эта возможность очевидна), поскольку в определенных ситуациях явно можно рассчитывать на некоторую гладкость итоговой (внешней) задачи (2). Понятие (5, L)-оракула в некотором смысле налагает наиболее слабые условия на возможные неточности в вычислении функции и градиента, при которых можно рассчитывать, что скорость сходимости метода, учитывающего гладкость (липшицевость градиента функционала) задачи, сильно не пострадает (см. теорему 1 ниже).

На первый взгляд может показаться, что применимость описанной концепции (5, L)-оракула к задаче (1) следует из следующего результата (см. п. 4.2.2 [32]).

Утверждение 2. Пусть подзадача энтропийно-линейного программирования (ЭЛП) в (2) решена (по функции) с точностью 6, т.е. найден такой X(с), удовлетворяющей балансовым ограничениям, что

п                    п

Е Xij (c)ln X^ (с) + Е Сгз X^ (с) - г,3 = 1                        г,з = 1

min п            п

£ Xij = L^, £ Xij = Wj i=1            i=1

Xij > 0, i, j = 1, ..., n

п

{ Е ^X-з

In Xij

п

+ Е г,3=1

^Сгз *^гз

}

6 6.

Тогда для функции

7/ \ f (с) = — пп

£ Xij = Li, £ Xij = Wj j=1

Xij > 0,i,j = 1,..., n

{ Е XijlnХгз + Е счX'з^ г,з=1

п

п

набор

(п                  п                 пп \

Е Xij(с) lnXij(с) + Е сгзXij(с), ^Xij(с)}     j г,3=1                      г,з=1                        г,з= / является (6, 2 max сгз)-оракулом.

г,з=1,...,п J

К сожалению, большинство методов (в том числе метод балансировки) не удовлетворяют одному пункту утверждения 2, а именно, они выдают вектор X, который лишь приближенно удовлетворяет балансовым ограничениям (в утверждении требование точного удовлетворения балансовых ограничений является существенным и не может быть как-то равнозначно релаксировано). Связанно это с тем, что для задачи ЭЛП, когда ограничений намного меньше числа прямых переменных, обычно решается двойственная задача, по которой восстанавливается решение прямой задачи [6, 33]. Как следствие, приобретается невязка и в ограничениях. Собственно, в представлении градиента функционала по формуле (4) имеются два способа. Первый через двойственные множители (А,^), второй через решение прямой задачи X. Функционал прямой задачи сильно выпуклый по X, поскольку энтропия 1-сильно выпуклая функция в 1-норме [15]. Поэтому сходимость в решении прямой задачи по функции обеспечивает сходимость и по аргументу, что и означает возможность определения c хорошей точностью градиента по формуле (4) через X. Другая ситуация возникает, если смотреть на двойственную задачу к задаче ЭЛП (в приводимом далее утверждении следует обозначить X = (А, ^), у = у для задачи (2), (3)).

Утверждение 3. Пусть у (у) = minФ(x,у), где Ф(x, у) - такая достаточно гладкая, вы-xEQ пуклая по совокупности переменных функция, что10

||VФ^(x‘^,у‘^{) —} VФ^(x,у⁾^₂6 L^^(X‘,y‘) ^{— (X,}y)^2.

Пусть для произвольного у Е Q (считаем, что множество Q содержит внутри себя шар радиуса более 26/L) можно найти такой X = X(у) Е Q, что max(Фж(X,у),X — г) 6 6.

ZEQ

Тогда

ФО^у) ^—у⁽у) 6 ⁶,

^уСу') - Vy(y)^₂6 іД\у' - у^, и ^Ф(ж,у) — 25, Фу(ж, у)) будет (65, 2L)-оракулом для у(у) на выпуклом множестве, полученном из множества Q отступанием от границы dQ во внутрь Q на расстояние ^25/L (по условию это множество не пусто).

Доказательство. По условию задачи имеем при всех допустимых значениях аргумен- тов Ф:

^max

^Фжж ^Фжу

^Фуж ^Фуу

) = sup ( Һ, І№61\

^Фжж ^Фжу

^Фуж ^Фуу

һ^ 6 L.

Заметим, что также по условию при всех допустимых значениях аргументов Ф:

^Фжж ^Фжу

^Фуж ^Фуу

^ 0,

Фхх У 0,  Фуу ^ 0, Фуж = Ф^Ту,  Фжж = Ф^Тж,  Фуу = ФТу.

Для упрощения последующих рассуждений (в частности, чтобы не работать с псевдооб-ратными матрицами) будем, считать, что матрица Ф_ЖЖ^ 0 положительно определена (исходя из условий, гарантировать можно лишь неотрицательную определенность). Также будем считать (в интересующих нас приложениях к задачам (2), (6) это имеет место), что зависимость ж (у), определяемая из соотношения

У (у) = тіпФ(ж,у) = Ф(ж(у),у) жед однозначным образом, и удовлетворяет соотношению

Фж(ж(у),у) = 0, из которого имеем

Фжж(ж(у),у) ||^ || + Фжу (ж(У),у) = 0, т.е.

||дж(ду¹¹= ¹¹dXi/dyj¹¹= -^Ф-ж^Фжу.

Поскольку у(у) = Ф(ж(у),у), то

Ууу = дЭх(Уу\\^т Фжж||дж/ду|| + |^|дж/ду|\^ТФжу + Фуж^дж/дуЦ + Фуу =

= ^Фуу ^{— Ф}уж^Фжж ^Фжу.

С учетом этой формулы и из формулы дополнения по Шуру [34], получаем

^Фжж ^Фжу

^Фуж ^Фуу

Е_ж   0

^Фуж^Фжж Еу

Фжж 0

^{0 Ф}уу

^Еж^Фжж ^Фжу

0   Еу

где Е_ж, Еу - единичные матрицы соответствующих размеров. Поскольку

Е_ж   0

^Фуж^Фжж ^Еу

^Еж ^Фжж ^Фжу

0   ^Еу

Т

и эти матрицы полного ранга, то из (12) имеем, что sup (һ, Ууу һ = ^max (Ууу) 6 b^b261

^тах{^тах^(Фжж), ^тах^(ууу)} — ^max

(

^Фжж ^Фжу

^Фуж ^Фуу

)6L

Таким образом, установлено, что у(у) 6 у(ж) + ^у(ж), у - ж) + 11\y - ж||2.

Согласно утверждению 1

У ⁽у) > у^{(ж) +} (Ф_у(ж, у), у - ж) - 5.

Далее проведем рассуждения аналогично рассуждениям на с. 107 (и немного отлично от с. 115) диссертации [32]. Вычитая из первого неравенства второе, получим

(^Фу (^ж,у) - ^Vу⁽ж^),у - ^ж) 6 ¹1|у - ж||2 + 5.

Положим t > 0:	-   =   ф у ( ж, у) - V У( ж ) t У Х   11фу( ж, у) -V У( ж )h 2 ,
получим	11 ф у ( ж, у) - V У( ж ) 11 2 6 у + 5 .

Минимизируя правую часть неравенства по t >  0, при t = ^25/1 получаем

|^|ф у ( ^ж, у) - ^V у^{( ж )}|^| ₂ 6 2^1.

Отсюда и из утверждения 1 находим

У⁽У) 6 У^{( ж ) +} ( ^V У^{( ж ) ,} У ^-ж ) ^{+ 1} 1¹ У ^{- ж} |¹ 2 6

6 у^{( ж ) +} (^ф у ^{( ж,} у^{) ,} у ^-ж ) ⁺ V^^ ^-ж || ₂ ⁺ |||у — ^{ж |} 2 6

6 Ф(ж,у) - 25 + (Ф у (ж, у), у - ж) + V 25L||y - ж^| ^ + ^^ - ж ^{|| 2} + 25 6

6 Ф(ж,у) - 25 + (Ф у (ж, у), у - ж) + Ь^у - ж ^ 2 + 65.

С учетом того, что (см. утверждение 1)

у⁽у) >  У^{( ж ) +} (^Фу ^{( ж} , у), у ^{- ж}( ^{- 5}

> фу (ж, у) - 25 + (Фу (ж, у), у - ж), из определения 1 получаем доказываемое утверждение.

Это утверждение позволяет установить гладкость задачи (2), (3) (но не (5), (6)). Таким образом, для (5), (6) необходимость использования универсального метода для внешней задачи является отражением надежды сходиться быстрее, чем в негладком случае, в то время как для (2), (3) использование универсального метода для внешней задачи является скорее отражением желания настраиваться на правильную константу Липшица градиента. Можно, конечно, пытаться использовать приведенные выше формулы, однако из способа рассуждений (см., например, доказательство утверждения 3) видно, что полученная таким образом константа Липшица может оказаться завышенной.

К сожалению, практическое применение утверждения 3 натыкается на следующие сложности:

• 1)    необходимости отступать от границы множества Q во внутрь на ^25/1,

• 2)    необходимости рассмотрения ситуации (см. доказательство утверждения 3):

II д ^ж г /ду ^ ¹¹ = ^{- Ф -1 ф} _Ж у ,

• 3)    необходимости предположения о компактности множества Q, иначе невозможно будет добиться выполнения условия

тах{Ф _ж (х, у), х — г) 6 5. z e Q

Первая сложность на практике разрешима за сче т возможности доопределения функционала задачи с сохранением всех свойств на ^25/L-окрестность множества Q (заметим, что доопределение часто не требуется, поскольку функционал и так задан «с запасом»). Например, для рассматриваемых нами транспортных приложений с Q = { у : у >  у } это возможно [2] – [5]. Сложность 2 часто вообще не возникает (разве что оговорка о существовании Ф -1 , впрочем, приведенные выше рассуждения можно провести, сохранив все результаты в идентичном виде, так, что эта оговорка будет не нужна), поскольку Q совпадает со всем (двойственным) пространством. А вот сложность 3, действительно, портит дело. К сожалению, простых теоретически обоснованных способов борьбы с этой сложностью мы пока не знаем. Тем не менее полезно заметить, что в действительности нужно гарантировать выполнение (см. доказательство утверждения 1)

(Фж(х(у),у),х(у) — х(у'Х) 6 5, где точки у и у' близки, поскольку возникают на соседних итерациях внешнего метода. С учетом ожидаемой ««близости» X = Х(у) и х(у), мы можем заменить в этом критерии настоящее множество Q, которое, как правило, совпадает со всем пространством, на шар конечного радиуса. Более детальные исследования (для задачи (2), (3)) и практические эксперименты показывают, что для выполнения приведенного выше условия достаточно обеспечить для внутреннего итерационного процесса {хк} ^ х(у) условия

||Ф х (Я к ,у) І 2 ||Xfc H 2 6 5/2, ||Ф _ж (х к ,у) І 2 6 5.

Соответствующее х к = (Х к , ^ к ) порождает нужное х(у) = х к . С учетом специфики рассматриваемой нами задачи (2), (3), имеем следующий критерий (возвращаемся к обозначениям (2), (3)):

^Ах(Ак,Цк) — b^ ||(Afc,цк)||2 6 5/2,   ||Ax(Xfc,цк) — b^ 6 5, где х(\к, ^к) определяется в формуле (4), а введённая линейная система балансовых уравнений Ах = b есть общая запись аффинных (транспортных) ограничений:

п            п

У^ х ^з =Ь г/^ хц = W j ,   г,] = 1, ...,П.

j=1           г=1

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

Принципиально важно для гладкого случая (aj (у) - функции с липшицевым градиентом), как это будет следовать из дальнейших оценок (см. теорему 1), не просто уметь решать двойственную задачу, т.е. находить (X, ^) так, чтобы была сходимость по аргументу, а делать это так, чтобы сложность решения задачи зависела от точности ее решения логарифмическим образом. Выше мы отмечали, что это имеет место для метода балансировки. Также это имеет место и для быстрых методов, примененных к регуляризованной двойственной задачи. При фиксации параметра регуляризации, исходя из итоговой желаемой точности, быстрые градиентные методы (для сильно выпуклых функций) решают регуляризованную двойственную задачу так, что зависимость сложности от точности ее решения – логарифмическая.

Хочется, чтобы при решении внешней задачи в (2), т.е. задачи min / (у), yeQ можно было не задумываться ни о какой гладкости. Если она есть, то метод бы это хорошо учитывал, не требуя знания констант Липшица градиента (это намного более существенно для возможности применять описанный подход к поиску барицентра Вассерштейна вероятностных мер, см. п. 2), если ее нет, то метод также работал бы оптимальным (для негладкого случая) образом. Именно таким свойством и обладает универсальный метод [27], работающий и в концепции неточного оракула [31] (см. определение 1).

Заметим [27], что можно погрузить задачу с гёльдеровым градиентом (v Е [0,1])

IV(у ‘ ) — V /(у) һ * 6 L _v || у ' - у һ "

(в том числе и негладкую задачу с ограниченной нормой разности субградиентов при v = 0) в класс гладких задач с оракулом, характеризующимся точностью 6 и

L = L „

L _v ⁽¹ — v ) 26(1 + v)

1-^ і+^

Это позволяет даже в случае, когда можно рассчитывать только на 6-субградиент ¹¹ (с ограниченной нормой субградиента (разности субградиентов), причем, какой именно константой ограниченной, методу знать не обязательно), все равно работать в концепции (6, L)-оракула.

Итак, у нас есть внешняя задача (2)

min / (у), yeQ для которой обращение к (6, L)-оракулу за значением функции и градиента стоит ~ ln(6 1). Насколько быстро мы можем решить такую задачу, т.е. при каком N(е) можно гарантиро- вать, что

^{/ (}У N (е) ) ^{— тіп /} (у) ^{6 е} ? y e Q

Ответ можно получить из следующего результата.

Теорема 1 (см. [1, 21, 27, 31]) . Cуществует однопараметрическое семейство универсаль ных градиентных методов (параметр р Е [0,1]) , не получающих на вход, кроме р, больше никаких параметров (в частности, не использующих значения L _p и R - «расстояние» от точки старта до решения, – априорно не известное ) , которое приводит к следующей оценке на требуемое число итераций:

N _p ( e )=O| inf ^G[0,1]

LX+

е

1+2pv+v

,

если 6 <  О (eN _p (е) ^p ).

11 На 5-субградиент всегда можно рассчитывать согласно утверждению 1. Причем, как уже отмечалось раньше, для получения 5-субградиента не нужна сходимость по аргументу для вспомогательной задачи.

Из теоремы 1 можно заключить, что если мы рассчитываем на некоторую гладкость /(у), то стоит выбирать значение параметра р = 1, при этом общие трудозатраты машинного времени будут

о(^( _£ )(Т ln( _£ ^-1 ) + Т)), (13)

где Т - время вычисления (суб-)градиента функционала (в основном это вычисления {Vc -j (у)}^ =1 [29,30]), Т - время решения вспомогательной задачи методом балансировки с относительной точностью 1%. Численные эксперименты показывают, что на одном современном ноутбуке при п ~ 10 2 время Т ~ 1 с. [31], что сопоставимо с временем Т для таких п [29].

Выгода от описанной выше конструкции по сравнению с обычным способом решения исходной задачи минимизации (2), (3) сразу по совокупности всех переменных (см., например, [5]) заключается в гарантированном не увеличении константы Липшица градиента в оценке необходимого числа итераций (см. утверждение 3) и ожидаемое уменьшение в этой же оценке расстояния от точки старта до (неизвестного априори) решения. Выгода здесь вполне может достигать одного порядка и более. При этом можно ожидать лишь незначительного увеличения стоимости одной итерации. Причем стоит иметь в виду, что при оптимизации сразу по всем переменным требуется рассчитывать градиент функционала по большему числу переменных, чем в описанном выше подходе, что также играет нам на пользу. В конечном итоге, сокращение числа итераций заметно превалирует над небольшим увеличением стоимости одной итерации.

Что касается задач (5), (6), то описанный выше подход представляется естественным и не имеющим альтернатив в рассматриваемом классе. Альтернативные методы, с которыми можно сравнивать, мы упоминали по ходу статьи, но все они были предложены из принципиально других подходов. Сравнению (практическому) всех этих методов планируется посвятить отдельную публикацию.

Резюмируем ключевой результат этого раздела (и всей статьи) следующим образом.

Для решения задач типа (2), (6) или (8) предлагается использовать универсальный метод из работы [27] (а точнее его модификацию из [31]) . Если рассчитываем на гладкость ¹² /(у) , то полагаем в методе р = 1 . Если на гладкость рассчитывать не приходится ¹³ , то полагаем р = 0 . В обоих случаях, кроме априорной подсказки относительно параметра р, методу больше ничего от нас знать не надо!

Авторы выражают глубокую признательность Ю. Е. Нестерову за серию продуктивных обсуждений.

Исследование А. В. Гасникова, П. Е. Двуреченского, В. Г. Спокойного и А. Л. Сувориковой в части 2 выполнено в ИППИ РАН за счет гранта Российского научного фонда (проект № 14-50-00150); исследование П. Е. Двуреченского в части 3 выполнено при поддержке гранта РФФИ 15-31-20571-мол_а_вед; исследование А. В. Гасникова в части 3 выполнено при поддержке гранта РФФИ 15-31-70001 мол_а_мос.