CC BY
10
УДК 519.688
А. С. Аникин1, А. В. Гасникое2'3, А.Ю. Горнов1
Институт динамики систем и теории управления им. В.М. Матросова СО РАН 2Институт проблем передачи информации им. А.А. Харкевича РАН 3Московский физико-технический институт (государственный университет)
О неускоренных эффективных методах решения разреженных задач квадратичной оптимизации
В работе исследуется класс разреженных задач квадратичной оптимизации и его обобщения, на которых будут эффективно работать специальные модификации обычного (неускоренного) прямого градиентного метода с 11-нормой и метода условного градиента (Франк-Вульфа).
Ключевые слова: huge-scale оптимизация, разреженность, 11-норма, метод Франк-Вульфа.
A.S. Anikin1, A. V. Gasnikov2'3, A. Yu. Gornov1
1Matrosov Institute for System Dynamics and Control Theory of SB RAS 2Kharkevich Institute for Information Transmission Problems of RAS 3Moscow Institute of Physics and Technology (State University)
Nonaccelerated efficient numerical methods for sparse quadratic optimization problems and their generalizations
We investigate the primal gradient method with 11-norm and the conditional gradient method (both are nonaccelerated methods).We show that these methods can outperform well known accelerated approaches for some classes of sparse quadratic problems. Moreover, we discuss some generalizations.
Key words: huge-scale optimization, sparsity, l1-norm, Frank-Wolfe method.
1. Введение
В работе [1] (см. также [2]) были предложены два различных способа решения задачи минимизации (на единичном симплексе) квадратичного функционала специального вида, связанного с задачей PageRank [3]. Первый способ базировался на обычном (неускоренном) прямом градиентном методе в 11-норме [4]. Второй способ базировался на методе условного градиента [5] (Франк-Вульф). Общей особенностью обоих способов является возможность учитывать разреженность задачи в большей степени, чем известные альтернативные способы решения отмеченной задачи. За счёт такого свойства итоговые оценки времени работы упомянутых методов в ряде случаев получаются лучше, чем оценки для ускоренных методов и их различных вариантов [2, 6]. Естественно было задаться вопросами: на какие классы задач и на какие простые множества (входящие в ограничения) можно перенести разработанные в [1] алгоритмы? Настоящая статья имеет целью частично ответить на поставленные вопросы.
В п. 2 мы показываем (во многом следуя изначальным идеям Ю.Е. Нестерова [7,8]), что прямой градиентный метод в 11-норме из [1] можно распространить на общие задачи квадратичной минимизации. При этом можно перейти от задач на симплексе к задачам на всем пространстве (это даже заметно упростит сам метод и доказательство оценок). Но наиболее интересно то, что удаётся распространить этот метод на общие аффинно-сепарабельные задачи с сепарабельными композитами [6]. Такие задачи часто возникают при моделировании компьютерных и транспортных сетей [9-11], а также в анализе данных [12]. Основными
конкурирующими методами здесь являются различные варианты (прямые, двойственные) ускоренных покомпонентных методов [6,7,9,13,14], которые чаще интерпретируются как различные варианты методов рандомизации суммы [6, 12, 15-22]. Таким образом, в п. 2 предлагается новый подход к этим задачам.
В п. 3 мы распространяем метод Франк-Вульфа из [1] на более общий класс задач квадратичной минимизации. К сожалению, выйти за пределы описанного в п. 3 класса задач нам не удалось. Тем не менее стоит отметить интересную (саму по себе) конструкцию, позволившую перенести метод с единичного симплекса на неотрицательный ортант. Стоит также отметить, что эта конструкция существенным образом опиралась на ряд новых идей, успешно использовавшихся и в других контекстах.
В п. 4 мы предлагаем новый рандомизированный метод (дополнительно предполагая сильную выпуклость функционала задачи), который работает быстрее неускоренного покомпонентного метода, и позволяет для ряда задач в большей степени воспользоваться их разреженностью, чем позволяют это сделать ускоренные покомпонентные методы.
«За кадром» описываемых далее сюжетов стоит задача разработки эффективных численных методов решения системы линейных уравнений Ах = b в пространстве огромной размерности (от десятков миллионов переменных и выше). При этом качество решения оценивается согласно \\Ах — Ь\\2. Полученные в работе результаты позволяют надеяться, что предложенные в пп. 2-4 методы будут доминировать остальные в случае, когда (приводим наиболее важные условия):
1) матрица А такова, что в каждом столбце и каждой строке не более s « л/п отличных от нуля элементов;
2) решение х* системы Ах = b таково, что « \\ж*\\2, или точнее: всегда имеет место неравенство
Г\\2 <Mi \\®*\\2 ,
которое (в нашем случае) можно уточнить следующим образом:
\\®*\\2 <Ml «^\И\2 .
Сравнительному анализу предложенных в данной работе методов с альтернативными методами поиска решения системы Ах = b (в частности, с популярными сейчас ускоренными покомпонентными методами [6]) посвящено приложение.
2. Прямой неускоренный градиентный метод в И-норме
Рассматривается следующая задача:
f (х) = 1 (Ах,х) — (b,x) ^ min, (1)
2 хей"
где квадратная матрица А предполагается симметричной неотрицательно определённой. Также предполагается, что в каждом столбце и каждой строке матрицы А не более s « п элементов отлично от нуля. Мы хотим решить задачу (1) с точностью е:
f (XN) — f (Х*) <
Заметим, что тогда
\\AxN — Ъ\\2 < л/2Amax(А)£.
Эту задачу предлагается решать обычным градиентным методом, но не в евклидовой норме, а в 11-норме [4] или (что то же самое для данной задачи) неускоренным покомпонентным
методом с выбором максимальной компоненты [7] (если гк определяется не единственным образом, то можно выбрать любого представителя):
^ = xfc + argmin |/ (xk) + lv/(xk),h) + L \\h\\j) = ( У i df(^) ] = V (2)
heR l \ / z xi - x dx.k , 1 = 1 ,
■h
где i* = argmax
df(xk)
L = max lAij |. Точку старта x0 итерационного процесса вы-
i'j=1'...,n
3=1,...,п
берем таким образом, чтобы лишь одна из компонент была отлична нуля. Любая итерация такого метода (за исключением самой первой - первая требует О(п) арифметических операций) может быть осуществлена за О^^2п). Действительно, хк+1 и хк отличаются в одной компоненте. Следовательно,
V /(хк+1) = Ахк+1 - Ь = Ахк - Ь + 1 Ае= V/(х*) + 1 Ц!^ Аег.,
где
ei = (0,..., 0,1, 0,..., 0)
Таким образом, по условию задачи (матрица A разрежена) V/(xfc+1) отличается от V/(хк) не более чем в s компонентах. Следовательно, ik+1 можно пересчитать за 0(slog2n), храня массив Цд/(хк)/dxj|}™=1 в виде специальной структуры данных (кучи) для поддержания максимального элемента [1,8].
Согласно [4] необходимое число итераций
2 max I Aij R
i'j=1'...,n
N = ———-,
где
Rj = sup{\\x -x*\\j : f(x) < f(x0)} .
Таким образом, общее число арифметических операций (время работы метода) можно оценить следующим образом:
( max I Aij I Rj n + s log2 n ——-
Рассмотрим теперь задачу вида
m n
/(x) = Y1/k (Abx) + Y19k (xk) ^ ¿ir^ k=1 k=1 X
где все функции (скалярного аргумента) Д, д^ имеют равномерно ограниченные числом L первые две производные,1 а матрица A = [Aj,..., Am]T такова, что в каждом ее столбце не больше sm ^ m ненулевых элементов, а в каждой строке - не больше sn ^ n. Описанный выше подход (формула (2)) позволяет решить задачу за время
( L max llA^I^ R?\
„ i=1,...,n 112
О n + snsm log2 n-
I e )
где A^ - i-й столбец матрицы A.
1g'k можно равномерно ограничивать числом L max llA^II .
г=1,...,п 11 11 2
Собственно, задача из [1] является частным случаем приведённой выше общей постановки (в смысле выбора функционала):
1 п 1нл„ 1|2 , 7 „Л2 ^
/(х) = -||Ах||2 + ^ ,тш,
2 2 —* (х,е)=1
где
, , (у, если у > 0,
(у)+ = { 0, У> ,
0, если у < 0.
В принципе, можно перенести результаты, полученные в этом пункте на случай, когда оптимизация происходит на симплексе. Это частично (но не полностью [1]) решает основную проблему данного метода: большое значение R2, даже когда ||ж*||2 небольшое. Однако мы не будем здесь этого делать.
3. Метод условного градиента (Франк—Вульфа)
Рассмотрим теперь следующую задачу2:
f{x) = 1 (Ax, х) - (^ х) ^ min, (3)
где квадратная матрица A предполагается симметричной неотрицательно определённой. Также предполагаем, что в каждом столбце и каждой строке матрицы A не более s ^ п элементов отлично от нуля, и в векторе b не более s элементов отлично от нуля.
Предположим сначала, что мы знаем такой R, что решение задачи (3) удовлетворяет условию
х* е£га^) = |х е R+ : ¿х* < R j .
Выберем одну из вершин Sn(R) и возьмём точку старта х1 в этой вершине. Далее будем действовать по индукции, шаг которой имеет следующий вид. Решаем задачу
(v f (хк ),у min . (4)
\ / yeSn(R)
Введём (если iк определяется не единственным образом, то можно выбрать любого представителя)
df (хк)
г к = argmm—-.
г=1,...,га дх1
Обозначим решение задачи (4) через
= ( R • ак, если df (хк)/дхгк < 0, Ук = \ 0, если df (хк)/дх%к > 0.
Положим
2
хк+1 = (1 -7к )хк + 1кук, 1к = —, к = 1,2,... Имеет место следующая оценка [5,23,24]:
2 L R2
f (хм) - }(х*) < ¡(хм) - k=max.N{f (хк) + (V/(хк),ук - х*)} < , (5)
2Отличие от задачи (1) в том, что рассматривается неотрицательный ортант вместо всего пространства -раздутием исходного пространства в два раза к такому ограничению можно прийти из задачи оптимизации на всем пространстве.
где
R2 = max ||w — ж||2 , Lv = max (h,Ah), 1 <р < ж. р x,yesn(R) р р 1Ш1Р<1
С учётом того, что оптимизация происходит на Sn(R), мы выбираем р = 1. Несложно показать, что этот выбор оптимален. В результате получим, что
Таким образом, чтобы
достаточно сделать
R? = 4R2, Li = max |. (6)
i,j=1,...,n
f(xN) — f(x*) < г,
N = 8maxi,i=i)...,ra A |R2
итераций.
Сделав один раз дополнительные вычисления стоимостью О(п), можно так организовать процедуру пересчёта V / (хк) и вычисление гк, что каждая итерация будет иметь сложность О (в 1о§2 п). Для этого вводим
к-1
Рк = П(1 - 7г), гк = хк/¡Зк, % = 7к/рк+1.
г=1
Тогда итерационный процесс можно переписать следующим образом:
гк+1 = гк + ~кук.
Пересчитывать Агк+1 — Ь при известном значении Агк можно за О(в). Далее, задачу
ъ к+1 = argmin д / (хк+1)/дХг
1=1,...,П
можно переписать как
гk+i = argmin ( ¡Агk+1] — -р— ) . i=1,...,n \L Рк+1/
Поиск ik+1 можно осуществить за 0(s log2n) (см. п. 2). Определив в конечном итоге zN, мы можем определить xN, затратив дополнительно не более
0(n + ln N) = О(п)
арифметических операций. Таким образом, итоговая оценка сложности описанного метода будет
( max I Aij |RfN n + s log2 n ——-
где R2 = ||x*||2 (следует сравнить с аналогичной формулой из п. 2).
Отметим, что при этом мы можем пересчитывать (это не надо делать, если известно значение f(x*), например, для задачи (3) естественно рассматривать постановки с f(x*) = 0)
f(xk+1) + (V f (xk+1), yk+1 — xk+^
также за О(,в log2n). Следуя [23], можно немного упростить приведённые рассуждения за счёт небольшого увеличения числа итераций. А именно, можно использовать оценки (при этом следует полагать 7к = 2(к + 2)-1)
¡(Хк) — Дх*) < (V /(хк),хк — ук) ,
_min (V /(xk),xk - yk^ <
7 LiRj
к=1,...,м \ " К 7' а/ - N + 2 Вернёмся теперь к предположению, что изначально известно К. На практике, как правило, даже если мы можем как-то оценить К, то оценка получается слишком завышенной. Используя формулы (5), (6). мы можем воспользоваться следующей процедурой рестартов по параметру К.
Выберем сначала К = К0 = 1. Делаем предписанное этому К число итераций и проверяем (без дополнительных затрат) критерий останова (все необходимые вычисления уже были сделаны по ходу итерационного процесса)
2
8 max I Aij |R
/(xN) - /(x") + (V /(xk), Ук -xk)} < ^N'^ 1-< e.
Если он выполняется, то мы угадали и получили решение. Если не выполняется, то полагаем R := \R (х > 1), и повторяем рассуждения. Остановимся поподробнее на вопросе оптимального выбора параметра х [25]. Обозначим R* = \\x* \\1. Пусть
R0Xr-1 < R* < R0Xr.
Общее число итераций будет пропорционально
х2г+2 1 х4 / п* \ 2
2 , „ 4 , „ 2r X - 1 ^ X I R \
1 + х2 + х4 + х2Г = 2 , <
X2 -1 - X2 - 1\К0; •
Выберем х = \/2, исходя из минимизации правой части этого неравенства. При этом общее число итераций возрастёт не более чем в четыре раза по сравнению со случаем, когда значение К известно заранее.
Описанный выше подход распространяется и на задачи
1 (х) = 1 ||Ах -ь II2 ^ .
Матрица А такова, что в каждом ее столбце не больше 8т ^ т ненулевых элементов, а в каждой строке - не больше 8п ^ п. В векторе Ь не более 8т элементов отлично от нуля. Описанный выше подход (на базе метода Франк-Вульфа) позволяет решить задачу за время
О
( max 1^112 R2\
i=1,...,n 112 1
n + SnSm log2 n-
V /
где К2 = Ух* || 1 (следует сравнить с аналогичной формулой из п. 2, полученной для более общего класса задач).3
В связи с написанным выше в этом пункте, заметим, что задача может быть не разрежена, но свойство разреженности появляется в решении при использовании метода Франк-Вульфа, что также может заметно сокращать объем вычислений. Интересные примеры имеются в работах [26] (п. 3.3) и [27].
В последнее время методы условного градиента переживают бурное развитие в связи с многочисленными приложениями к задачам машинного обучения. В связи с этим появились интересные обобщения классического метода Франк-Вульфа. Отметим, например, работы [28, 29]. Интересно было бы понять: возможно ли перенести (а если возможно, то как именно и в какой степени) результаты этого пункта на более общий класс задач (чем класс задач с квадратичным функционалом) и на более общий класс методов?
3Здесь так же, как и в задаче (3), требуются рестарты по неизвестному параметру Нх, который явно
используется в методе в качестве размера симплекса, фигурирующего при решении вспомогательной задачи ЛП на каждой итерации. Однако, как было продемонстрировано выше, все это приводит к увеличению числа итераций не более чем в четыре раза.
4. Неускоренный рандомизированный градиентный спуск в сильно выпуклом случае
Рассмотрим для большей наглядности снова задачу (1). Будем считать, что f(x) -ß-сильно выпуклый функционал в 2-норме, т.е. Amin(A) > ß > 0. Сочетая методы [30] и [31], введём следующий рандомизированный метод:
_x^_
2
ß ■ (k + 1)
IV(xfc )||isign
fdf(xk) \
\dxi(xk) J
^i(xk),
x1 _ 0,
где
ei _ (0,...,0,1,0,...,0),
i(xk) _ i с вероятностью ц^.^ Тогда4 [30]
df (xk) dxi
i _ 1,
,n.
N
E
f(yk) - f(x*) <
fc=i
£_£e ||v/(xfc)||!
ß ■ (N + 1)
<
2 max E
k=1,...,N
|V f (xk )|
ß ■ (N + 1)
E
V/(yk)
<
4L max E
k=1,...,N
|V f (xk )|
ß ■ (N + 1)
где L _ Amax(^4), а
N
yk _
N ■ (N + 1)
k
k=1
На базе описанного метода построим новый метод, который «следит» за последовательностью к ■ хк), пересчитывая (с некоторой частотой)
V
Метод «дожидается» момента5 N _ 0(nL/ß), когда
||Vf (yN)\1 < 2 ||Vf (x1)!! ,
и перезапускается с x1 :_ yN. Можно показать, что необходимое число таких перезапусков для достижения точности е (по функции) будет ~ log2(e-1). Более того, подобно [6], можно так организовать описанную выше процедуру, чтобы попутно получить и оценки вероятностей больших отклонений (детали мы вынуждены здесь опустить).
Предположим теперь, что V f (x + hei) отличается от V f (x) (при произвольных x, h и ei) не более чем в s компонентах ( s ^ n).6 Для задачи (1) это имеет место (можно
4Второе неравенство может быть достаточно грубым [6].
5Это довольно грубая оценка на N, поскольку использующееся при ее получении правое неравенство
1 < ||V/(у
Ч/IV f (yN"12
< n,
2<
может быть сильно завышенным.
6Это условие можно обобщить с сохранением оценок, например, на случай, когда f (х) = д(х) + h(x), и существует такая (скалярная) функция a(x,he-i) > 0, что Vf (х+hei) = Vg(x + hei) + Vh(x + hei) отличается от a(x,hei) • Vg(x) + Vh(x) не более чем в s компонентах. Для этого приходится перезаписать исходный метод: отличие в том, что теперь вводится двойная рандомизация (рандомизация согласно вектору Vg(x) и независимая рандомизация согласно вектору Vh(x)).
1
2
1
2
1
2
2
2
к
x
2
1
также говорить и о задаче (3), с очевидной модификацией описанного метода - все оценки сохраняются). Можно так организовать процедуру выбора момента N (с сохранением свойства N = О(пЬ//л)), что амортизационная (средняя) сложность итерация метода (с учетом затрат на проверку условий остановки на каждом перезапуске) будет О(slog2n) [1,3]. Казалось бы, что мы ничего не выиграли по сравнению с обычными (неускоренным) покомпонентными методами [6]. Количество итераций и стоимость одной итерации одинаковы для обоих методов (в категориях О() с точностью до логарифмических множителей). В действительности, ожидается, что предложенный нами алгоритм будет работать заметно быстрее (неускоренного) покомпонентного метода, поскольку, как уже отмечалось, при получении этих оценок мы пару раз использовали потенциально довольно грубые неравенства. С другой стороны, пока мы говорили только о задаче (1) (в разреженном случае), для которой эффективно работают ускоренные покомпонентные методы [6] с такой же оценкой стоимости одной итерации, но заметно лучшей оценкой для числа итераций N = О(пл/Щ) (с точностью до логарифмического множителя). В связи с этим может показаться, что предложенный в этом пункте метод теоретически полностью доминируем ускоренными покомпонентными методами. На самом деле, это не совсем так. Хорошо известно, см., например, [6], что существующие сейчас всевозможные модификации ускоренных покомпонентных методов, которые могут сполна учитывать разреженность задачи (что проявляется в оценке стоимости итерации О (в)), применимы лишь к специальному классу задач [6]. Для общих задач стоимость одной итерации ускоренного покомпонентного метода будет О(п) независимо от разреженности задачи. Например, на данный момент неизвестны такие модификации ускоренных покомпонентных методов, которые бы позволяли сполна учитывать разреженность в задаче [6, 25, 31]
(т \
^ ехр (А1х)\ — Ъ) + | ||х||2 ^ тт
В то время, как описанный выше метод, распространим и на эту задачу с оценкой амортизационной стоимости одной итерации О(втвп ^2п).
5. Приложение
В этом приложении мы напомним некоторые полезные факты о поиске решения (псевдорешения) системы линейных уравнений Ах = Ь [32,33]. К близким задачам и методам приводит изучение проекции точки на аффинное множество Ах = Ь [9,31].
Хорошо известно, что задача поиска такого вектора х, что А х = , может быть эффективно полиномиально решена даже в концепции битовой сложности (например, простейшим методом Гаусса или более современными методами внутренней точки). Однако требуется, учитывая специфику задачи, так подобрать метод, чтобы можно было искать решения систем огромных размеров. Несмотря на более чем двухвековую историю, эта область математики (эффективные численные методы решения систем линейных уравнений) до сих пор бурно развивается. О чем, например, говорит недавний доклад Д. Спильмана на международном математическом конгрессе [35].
Пусть наблюдается вектор
Ь = Ах + е,
где матрица А - известна, ек € N(0, а2) - независимые одинаково распределенные случайные величины, к = 1,...,т (ненаблюдаемые). Оптимальная оценка7 неизвестного вектора параметров х определяется решением задачи
Пх) = 2IIА х — 6||2 . (7)
7То есть несмещенная и с равномерно (по х) минимальной дисперсией. При такой (статистической) интерпретации также стоит предполагать, что п < т.
Введём псевдообратную матрицу (Мура-Пенроуза) А* = (АТА)-1АТ, если столбцы матрицы А линейно независимы (первый случай), и А^ = АТ(ААТ)-1, если строки матрицы А линейно независимы (второй случай). В первом случае х* = А^Ь - единственное решение задачи (7), во втором случае х* = А^Ь - решение задачи (7) с наименьшим значением 2-нормы.8 Вектор х* называют псевдорешением задачи Ах = Ь.
Поиск решения системы Ах = Ь был сведён выше к решению задачи (7). Для задачи выпуклой оптимизации (7) существуют различные эффективные численные методы. Например, метод сопряжённых градиентов, сходящийся со скоростью (стоимость одной итерации этого метода равна О (ей) - числу ненулевых элементов в матрице А)
1 \\Ахи - Ъ\\1 = /(хм) - ¡(х*) < г,
2
N = О
|уАтах(АТА)||х*1^ ,
где N < п. На практике при не очень больших размерах матрицы А (тах{т, п} < 106) эффективно работают квазиньютоновские методы типа Ь-БРОБ [34]. А при совсем небольших размерах (тах{т,п} < 103-104) - и методы внутренней точки [33]. Однако нам интересны ситуации, когда тах{т,п} ^ 106. В таких случаях эффективно использовать прямые и двойственные ускоренные покомпонентные методы для задачи (7). Эти методы дают следующие оценки [6] (одна итерация у обоих методов в среднем требует О( ) арифметических
операций)9:
2ЦАх* -ь 112
< (прямой метод) ,
, 1 */2 = п £ ЦА^Ь, ^ = ||х*||2, Е [||Ахм - Ь||2] < е (двойственный метод),
о^п^Щ , Ц/2 = т р |А||2, Ку = ||л,
где у* - решение (если решение не единственно, то можно считать, что у* - решение с наименьшей 2-нормой) «двойственной» задачи
- 2 \\аТУ\\2 + <ь, у)
Напомним, что (Ак)Т - к-я строка матрицы А, а - к-й столбец.
Из написанного выше не очевидно, что для получения решения х* необходимо исходить именно из решения задачи оптимизации (7) и, как следствие, исходить из определяемого этой задачей критерия точности решения (2-норма невязки). Например, в настоящей работе мы исходили в основном из задачи (1), которая хотя и похожа по свойствам на задачу (7), но все же приводит к отличному от (7) критерию точности решения. Вместо задачи (7)
8В приложениях во втором случае чаще ищут не псевдорешение, а вводят в функционал регуляриза-тор [9], например, исходя из байесовского подхода (дополнительного априорного вероятностного предположения относительно вектора неизвестных параметров) или исходя из желания получить наиболее разреженное решение. В таком случае многое (но не все) из того, что написано в данной статье удаётся сохранить (должным образом модифицировав). При этом появляются и дополнительные новые возможности за счёт появления новых степеней свободы в выборе регуляризатора [27].
9Далее для наглядности мы дополнительно предполагаем, что f (х*) = 0. Если f (х*) = f * > 0, то из того, что /(хм) — f * < е будет следовать — 6||2 — /* < е//* (это простое, но полезное наблюдение
было сделано Дмитрием Островским (ЩША, Гренобль)).
в определённых ситуациях вполне осмысленно бывает рассматривать даже негладкую (но выпуклую) задачу
/(х) = ||Аж ^ т|п .
Для такой задачи можно предложить специальным образом рандомизированный вариант метода зеркального спуска [36], который гарантирует выполнение (считаем /(х*) = 0)
Е [||Ахм -Ь||те]
за время
/ тах ЦАкЩд2'
„ I . к=\,...,т
О 8тп + 8т 1о§2 т-^-
если в каждом столбце матрицы А не больше 8т < т ненулевых элементов. Отметим [36], что обычный (не рандомизированный) метод зеркального спуска гарантирует выполнение
||Ажм -Ь||те
за время
/ тах ЦАк ||2Д2' О I т + п + 8п8т log2 т-^2-
если в каждом столбце матрицы А не больше 8т < т ненулевых элементов, а в каждой строке - не больше < п.
Выше мы обсуждали так называемые вариационные численные методы, которые сводили поиск решения системы Ах = Ь к решению задачи выпуклой оптимизации. Однако имеется большой класс итерационных численных методов (типа метода простой итерации), которые исходят из перезаписи системы Ах = Ь в эквивалентном виде. Например, в виде ж = Ах + Ь, где А = I — А, и организации вычислений по формуле10 хк+1 = Ахк + Ь [39]. Для ряда таких методов также удаётся получить оценки скорости сходимости, причём не только для нормы невязки в системе, но и непосредственно для оценки невязки в самом решении ||хм — х*||. К сожалению, такие оценки получаются весьма пессимистичными, причём эти оценки типично реализуются на практике (закон Мерфи). Поясним сказанное примером, восходящим к работе Красносельского-Крейна [40]. Пусть
хк+1 = Ахк + ь,
где А = I — А - положительно определённая матрица п х п с собственными числами 0 < Х1 < ... < Хп < 1, а х0 - выбирается равновероятно из В^(0) - 2-шара радиуса Д с центром в 0. Обозначим через
5к = хк+1 -хк = Ахк + Ь-хк = Ь- Ахк
контролируемую невязку. В качестве критерия останова итерационного процесса выберем момент N первого попадания 5м € В£(0). Заметим, что
||х* -х*||2 =
(' -А)
1 5к
А-15к
10Интересно заметить, что критерием сходимости такого итерационного процесса является: Атах(А) < 1 [32]. Однако, если при этом ||А||2,2 = а"тах(^4) = Атах(АТА) > 1 (к счастью, такие ситуации не типичны [37]), то существуют такие х0 (причём на них довольно легко попасть), стартуя с которых, итерационная процедура на первых итерациях приводит к резкому росту нормы невязки — х*Ц2, т.е. к наличию «горба» у такой последовательности. Высота этого горба может быть (в зависимости от конкретной задачи) сколь угодно большой [32,38]. Поскольку все вычисления проходят с конечной длиной мантиссы, то такой резкий рост может приводить к вычислительной неустойчивости процедуры, что неоднократно наблюдалось в численных экспериментах [32]. Если сттах(А) < 1, то таких проблем нет, и последовательность — х*Ц2 мажорируется геометрической прогрессией с основанием <гтах(А).
Имеет место следующий результат. При Д ^ те вероятность выполнения неравенств
0.999< Цхм -х*||2 < ^
1 - Хп 1 — Хп
стремится к единице. Поскольку типично, что при больших п число Хп близко к 1, то вряд ли можно рассчитывать на то, что ||хм - х*||2 мало, если установлена только малость ЦАхм -Ь||2.
Отметим, что оценки скорости сходимости многих итерационных методов (типа простой итерации) исходят из первого метода Ляпунова, т.е. из спектральных свойств матрицы А. При плохих спектральных свойствах вводится регуляризация. По-сути, большинство таких методов можно проинтерпретировать как численный метод решения соответствующей вариационной переформулировки задачи. Регуляризация при этом имеет чёткий смысл (восходящий к пионерским работам А.Н. Тихонова) - сделать функционал сильно выпуклым и использовать его сильную выпуклость при выборе метода и получении оценки скорости сходимости [41]. Однако вариационный подход позволяет сполна использовать и второй метод Ляпунова, что даёт возможность распространить всю мощь современных численных методов выпуклой оптимизации на решение задачи А х = . В частности, использовать рандомизацию [36] или/и свойства неравнозначности компонент в решении х* и разреженность А, которые использовались при разработке методов из настоящей работы.
Ю.Е. Нестеровым недавно было отмечено, что возможные прорывы в разработке новых эффективных численных методов (в том числе применительно к решению системы Ах = Ь) можно ожидать от рандомизированных квазиньютоновских методов. К сожалению, здесь имеются большие проблемы, унаследованные от обычных (не рандомизированных) квазиньютоновских методов, с получением оценок скоростей сходимости, адекватных реальной (наблюдаемой на практике) скорости сходимости (все известные сейчас оценки никак не объясняют быстроту сходимости этих методов на практике). Тем не менее недавно появилось несколько интересных работ Гаверса-Ричтарика в этом направлении [42-44].
6. Заключение
Авторы выражают благодарность Ю.Е. Нестерову за ряд ценных замечаний и внимательное отношение к работе.
Исследование в части 2 выполнено при поддержке гранта РФФИ 14-01-00722-а, исследование в части 3 - при поддержке гранта РФФИ 15-31-20571-мол_а_вед, исследование в части 4 выполнено в ИППИ РАН за счёт гранта Российского научного фонда (проект № 14-50-00150).
Литература
1. Аникин А.С., Гасников А.В., Горнов А.Ю., Камзолов Д.И., Максимов Ю.В., Нестеров Ю. Е. Эффективные численные методы решения задачи PageRank для дважды разреженных матриц // Труды МФТИ. 2015. Т. 7, № 4. С. 74-94. arXiv:1508.07607
2. Гасников А.В., Двуреченский П.Е., Нестеров Ю.Е. Стохастические градиентные методы с неточным оракулом // Труды МФТИ. 2016. — Т. 8, № 1. С. 41-91. arxiv:1411.4218
3. Гасников А.В., Дмитриев Д.Ю. Об эффективных рандомизированных алгоритмах поиска вектора PageRank // ЖВМ и МФ. 2015. Т. 54, № 3. C. 355-371. arXiv:1410.3120
4. Allen-Zhu Z, Orecchia L. Linear coupling: An ultimate unification of gradient and mirror descent // e-print. 2014. arXiv:1407.1537
5. Jaggi M. Revisiting Frank-Wolfe: Projection-free sparse convex optimization // Proceedings of the 30th International Conference on Machine Learning, Atlanta, Georgia, USA, 2013. https://sites.google.com/site/frankwolfegreedytutorial
6. Гасников А.В., Двуреченский П.Е., Усманова И.Н. О нетривиальности быстрых (ускоренных) рандомизированных методов // Труды МФТИ. 2016. Т. 8, № 2 (в печати) arXiv:1508.02182
7. Nesterov Y.E. Efficiency of coordinate descent methods on large scale optimization problem // SIAM Journal on Optimization. 2012. V. 22, N 2. P. 341-362.
http://www1.se.cuhk.edu.hk/~sqma/SEEM5121_Spring2015/Nesterov-CD-2012.pdf
8. Nesterov Y. E. Subgradient methods for huge-scale optimization problems // CORE Discussion Paper 2012/2. 2012.
http://www.optimization-online.org/DB_FILE/2012/02/3339.pdf
9. 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 // e-print. 2015. arXiv:1508.00858
10. Гасников А.В., Гасникова Е.В., Двуреченский П.Е., Ершов Е.И., Лагуновская А.А. Поиск стохастических равновесий в транспортных моделях равновесного распределения потоков // Труды МФТИ. 2015. Т. 7, № 4. С. 114-128. arXiv:1505.07492
11. Гасников А.В., Двуреченский П.Е., Дорн Ю.В., Максимов Ю.В. Численные методы поиска равновесного распределения потоков в модели Бэкмана и модели стабильной динамики // Математическое моделирование. 2016. Т. 28. (в печати) arXiv:1506.00293
12. Johnson B., Zhang T. Accelerating stochastic gradient descent using predictive variance reduction // In Advances in Neural Information Processing Systems (NIPS). 2013. http://www.stat.rutgers.edu/home/tzhang/pubs.html
13. Fercoq O., Richtarik P. Accelerated, parallel and proximal coordinate descent // e-print. 2013. arXiv:1312.5799
14. Qu Z, Richtarik P. Coordinate Descent with Arbitrary Sampling I: Algorithms and Complexity // e-print. 2014. arXiv:1412.8060
15. Le Roux N., Schmidt M., Bach F. A stochastic gradient method with an exponential convergence rate for strongly-convex optimization with finite training sets //In Advances in Neural Information Processing Systems (NIPS). 2012. arXiv:1202.6258
16. Konecny J., Richtarik P. Semi-Stochastic gradient descent methods // e-print. 2013. arXiv:1312.1666
17. Konecny J., Liu J., Richtarik P., Takac M. Mini-batch semi-stochastic gradient descent in the proximal setting // e-print. 2015. arXiv:1504.04407
18. Agarwal A., Bottou L. A lower bound for the optimization of finite sums // e-print. 2014. arXiv:1410.0723
19. Shalev-Shwartz S., Zhang T. Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization // In Proceedings of the 31th International Conference on Machine Learning, ICML 2014, Beijing, China, 21-26 June 2014. P. 64-72. arXiv:1309.2375
20. Zheng Q., Richtarik P., Zhang T. Randomized dual coordinate ascent with arbitrary sampling // e-print. 2014. arXiv:1411.5873
21. Lan G., Zhou Y. An optimal randomized incremental gradient methods // Technical Report, Department of Industrial and Systems Engineering, University of Florida, July 7, 2015,
updated in October, 2015.
http://www.ise.ufl.edu/glan/files/2015/10/OptRandom10-18.pdf
22. Lin Q., Lu Z, Xiao L. An Accelerated Randomized Proximal Coordinate Gradient Method and its Application to Regularized Empirical Risk Minimization // SIAM J. Optim. 2015. V. 25, N 4. P. 2244-2273.
http://research.microsoft.com/pubs/228730/spdc_paper.pdf http://research.microsoft.com/pubs/258430/APCG_ERM_2015.pdf
23. Harchaoui Z., Juditsky A., Nemirovski A. Conditional gradient algorithms for norm-regularized smooth convex optimization // Math. Program. Ser. B. 2015. V. 152. P. 75-112. http://www2.isye.gatech.edu/~nemirovs/ccg_revised_apr02.pdf
24. Nesterov Yu. Complexity bounds for primal-dual methods minimizing the model of objective function // CORE Discussion Papers. 2015/03.
https://www.uclouvain.be/cps/ucl/doc/core/documents/coredp2015_3web.pdf
25. Гасников А.В., Гасникова Е.В., Нестеров Ю.Е., Чернов А.В. Об эффективных численных методах решения задач энтропийно-линейного программирования // ЖВМ и МФ. 2016. Т. 56, № 4. С. 523-534. arXiv:1410.7719
26. Bubeck S. Convex optimization: algorithms and complexity //In Foundations and Trends in Machine Learning. 2015. V. 8, N 3-4. P. 231-357. arXiv:1405.4980
27. Nesterov Yu., Nemirovski A. On first order algorithms for l\ / nuclear norm minimization // Acta Numerica. 2013. V. 22. P. 509-575.
http://www2.isye.gatech.edu/~nemirovs/ActaFinal_2013.pdf
28. Lan G., Zhou Y. Conditional gradient sliding for convex optimization // SIAM J. Optim. 2015. http://www.ise.ufl.edu/glan/files/2015/09/CGS08-31.pdf
29. Lacost-Julien S., Jaggi M. On the Global Linear Convergence of Frank-Wolfe Optimization Variants // e-print. 2015. arXiv:1511.05932
30. Lacost-Julien S., Schmidt M., Bach F. A simpler approach to obtaining 0(1/1) convergence rate for the projected stochastic subgradient method // e-print. 2012. arXiv:1212.2002
31. Аникин А.С., Двуреченский П.Е., Гасников А.В., Тюрин А.И., Чернов А.В. Двойственные подходы к задачам минимизации простых сильно выпуклых функционалов при аффинных ограничениях // ЖВМ и МФ. 2016. Т. 56. (подана) arXiv:1602.01686
32. Рябенький В.С. Введение в вычислительную математику. М.: Физматлит, 2008.
33. Nemirovski A. Introduction to linear optimization // ISYE 661. Lecture Notes in Georgian Institute of Technology, 2012.
http://www2.isye.gatech.edu/~nemirovs/OPTI_LectureNotes2015.pdf
34. Nocedal J., Wright S.J. Numerical Optimization. Springer, 2006.
35. Spielman D. Algorithms, graph theory, and linear equations in Laplacian matrices // Proc. of the International Congress of Mathematicians. Hyderabad, India, 2010. P. 1-23. http://www.cs.yale.edu/homes/spielman/PAPERS/icm10post.pdf
36. Аникин А.С., Гасников А.В., Горнов А.Ю. Рандомизация и разреженность в задачах huge-scale оптимизации на примере работы метода зеркального спуска // Труды МФТИ. 2016. Т. 8, № 1. С. 11-24. arXiv:1602.00594
37. Опойцев В.И. Устойчивые системы большой размерности // Автомат. и телемех. 1986. № 6. С. 43-49.
38. Поляк Б.Т., Тремба А.А., Хлебников М.В., Щербаков П.С., Смирнов Г.В. Большие отклонения в линейных системах при ненулевых начальных условиях // Автомат. и телемех. 2015. № 6. C. 18-41.
39. Saad Y. Iterative methods for sparse linear systems. SIAM, 2003. http://www-users.cs.umn.edu/~saad/books.html
40. Красносельский М.А., Крейн С.Г. Замечание о распределении ошибок при решении системы линейных уравнений при помощи итерационного процесса // УМН. 1952. Т. 7, № 4(50). С. 157-161.
41. Бакушинский А.Б., Кокурин М.Ю. Итерационные методы решения некорректных операторных уравнений с гладкими операторами. М.: УРСС, 2002.
42. Gower R., Richtarik P. Randomized Iterative Methods for Linear Systems // SIAM. J. Matrix Anal. & Appl. 2015. V. 36(4). P. 1660-1690. arXiv:1506.03296
43. Gower R., Richtarik P. Stochastic Dual Ascent for Solving Linear Systems // e-print. 2015. arXiv:1512.06890
44. Gower R., Richtarik P. Randomized Quasi-Newton Updates are Linearly Convergent Matrix Inversion Algorithms // e-print. 2015. arXiv:1602.01768
References
1. Anikin A.S., Gasnikov A.V., Gornov A.Yu., Kamzolov D.I., Maksimov Yu.V., Nesterov Yu.E. Effective numerical methods for huge-scale linear systems with double-sparsity and applications to PageRank. Proceedings of MIPT. 2015. V. 7, N 4. P. 74-94. arXiv:1508.07607
2. Gasnikov A.V., Dvurechensky P.E., Nesterov Yu.E. Stochastic gradient methods with inexact oracle. Proceedings of MIPT. 2016. V. 8, N 1. P. 41-91. arxiv:1411.4218
3. Gasnikov A.V., Dmitriev D.Yu. On efficient randomized algorithms for finding the PageRank vector. Zh. Vychisl. Mat. Mat. Fiz. 2015. V. 54, N 3. P. 355-371. arXiv:1410.3120
4. Allen-Zhu Z., Orecchia L. Linear coupling: An ultimate unification of gradient and mirror descent. e-print. 2014. arXiv:1407.1537
5. Jaggi M. Revisiting Frank-Wolfe: Projection-free sparse convex optimization. Proceedings of the 30th International Conference on Machine Learning, Atlanta, Georgia, USA, 2013. https://sites.google.com/site/frankwolfegreedytutorial
6. Gasnikov A.V., Dvurechensky P.E., Usmanova I.N. About accelerated randomized methods. Proceedings of MIPT. 2016. V. 8, N 2. (in print) arXiv:1508.02182
7. Nesterov Y.E. Efficiency of coordinate descent methods on large scale optimization problem. SIAM Journal on Optimization. 2012. V. 22, N 2. P. 341-362.
http://www1.se.cuhk.edu.hk/~sqma/SEEM5121_Spring2015/Nesterov-CD-2012.pdf
8. Nesterov Y. E. Subgradient methods for huge-scale optimization problems. CORE Discussion Paper 2012/2. 2012.
http://www.optimization-online.org/DB_FILE/2012/02/3339.pdf
9. 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. e-print. 2015. arXiv:1508.00858
10. Gasnikov A.V., Gasnikova E.V., Dvurechensky P.E., Ершов E.I., Lagunovskaya A.A. Search for the stochastic equilibria in the transport models of equilibrium flow distribution. Proceedings of MIPT. 2015. V. 7, N 4. P. 114-128. arXiv:1505.07492
11. Gasnikov A.V., Dvurechensky P.E., Dorn Yu.V., Maximov Yu.V. Searching equillibriums in Beckmann's and Nesterov-de Palma's models. Mathematical Models and Computer
Simulations. 2016. V. 28. (in print) arXiv:1506.00293
12. Johnson B., Zhang T. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in Neural Information Processing Systems (NIPS). 2013. http://www.stat.rutgers.edu/home/tzhang/pubs.html
13. Fercoq O, Richtarik P. Accelerated, parallel and proximal coordinate descent. e-print. 2013. arXiv:1312.5799
14. Qu Z., Richtarik P. Coordinate Descent with Arbitrary Sampling I: Algorithms and Complexity. e-print. 2014. arXiv:1412.8060
15. Le Roux N., Schmidt M., Bach F. A stochastic gradient method with an exponential convergence rate for strongly-convex optimization with finite training sets. In Advances in Neural Information Processing Systems (NIPS). 2012. arXiv:1202.6258
16. Konecny J., Richtarik P. Semi-Stochastic gradient descent methods. e-print. 2013. arXiv:1312.1666
17. Konecny J., Liu J., Richtarik P., Takac M. Mini-batch semi-stochastic gradient descent in the proximal setting. e-print. — 2015. arXiv:1504.04407
18. Agarwal A, Bottou L. A lower bound for the optimization of finite sums. e-print. 2014. arXiv:1410.0723
19. Shalev-Shwartz S., Zhang T. Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization. In Proceedings of the 31th International Conference on Machine Learning, ICML 2014, Beijing, China, 21-26 June 2014. P. 64-72. arXiv:1309.2375
20. Zheng Q., Richtarik P., Zhang T. Randomized dual coordinate ascent with arbitrary sampling. e-print. 2014. arXiv:1411.5873
21. Lan G., Zhou Y. An optimal randomized incremental gradient methods. Technical Report, Department of Industrial and Systems Engineering, University of Florida, July 7, 2015, updated in October, 2015.
http://www.ise.ufl.edu/glan/files/2015/10/DptRandom10-18.pdf
22. Lin Q., Lu Z, Xiao L. An Accelerated Randomized Proximal Coordinate Gradient Method and its Application to Regularized Empirical Risk Minimization . SIAM J. Optim. 2015. V. 25, N 4. P. 2244-2273.
http://research.microsoft.com/pubs/228730/spdc_paper.pdf http://research.microsoft.com/pubs/258430/APCG_ERM_2015.pdf
23. Harchaoui Z., Juditsky A., Nemirovski A. Conditional gradient algorithms for norm-regularized smooth convex optimization. Math. Program. Ser. B. 2015. V. 152. P. 75-112. http://www2.isye.gatech.edu/~nemirovs/ccg_revised_apr02.pdf
24. Nesterov Yu. Complexity bounds for primal-dual methods minimizing the model of objective function. CORE Discussion Papers. 2015/03.
https://www.uclouvain.be/cps/ucl/doc/core/documents/coredp2015_3web.pdf
25. Gasnikov A.V., Gasnikova E.V., Nesterov Yu.E., Chernov A.V. Entropy linear programming. Zh. Vychisl. Mat. Mat. Fiz. 2016. V. 56, N 4. P. 523-534.
arXiv:1410.7719
26. Bubeck S. Convex optimization: algorithms and complexity. In Foundations and Trends in Machine Learning. 2015. V. 8, N 3-4. P. 231-357. arXiv:1405.4980
27. Nesterov Yu., Nemirovski A. On first order algorithms for l\ / nuclear norm minimization. Acta Numerica. 2013. V. 22. P. 509-575.
http://www2.isye.gatech.edu/~nemirovs/ActaFinal_2013.pdf
28. Lan G., Zhou Y. Conditional gradient sliding for convex optimization. SIAM J. Optim. 2015. http://www.ise.ufl.edu/glan/files/2015/09/CGS08-31.pdf
29. Lacost-Julien S., Jaggi M. On the Global Linear Convergence of Frank-Wolfe Optimization Variants. e-print. 2015. arXiv:1511.05932
30. Lacost-Julien S., Schmidt M., Bach F. A simpler approach to obtaining 0(1/1) convergence rate for the projected stochastic subgradient method. e-print. 2012. arXiv:1212.2002
31. Anikin A.S., Gasnikov A.V., Turin A.I., Chernov A.V. Dual approaches to the strongly convex simple function minimization problem under affine restrictions. Zh. Vychisl. Mat. Mat. Fiz. 2016. V. 56. (submitted) arXiv:1602.01686
32. Ryabenkij V.S. Introduction to computational mathematics. Moscow: Fizmatlit, 2008.
33. Nemirovski A. Introduction to linear optimization. ISYE 661. Lecture Notes in Georgian Institute of Technology, 2012.
http://www2.isye.gatech.edu/~nemirovs/0PTI_LectureNotes2015.pdf
34. Nocedal J., Wright S.J. Numerical Optimization. Springer, 2006.
35. Spielman D. Algorithms, graph theory, and linear equations in Laplacian matrices. Proc. of the International Congress of Mathematicians. Hyderabad, India, 2010. P. 1-23. http://www.cs.yale.edu/homes/spielman/PAPERS/icm10post.pdf
36. Anikin A.S., Gasnikov A.V., Gornov A.Yu. Randomization and sparsity in huge-scale optimization on the Mirror Descent example. Proceedings of MIPT. 2016. V. 8, N 1. P. 1124. arXiv:1602.00594
37. Opoitsev, V.I. High-dimensional stable systems . Autom. Remote Control. 1986. V. 47, N 6. P. 768-774.
38. Polyak B.T., Tremba A.A., Khlebnikov M.V., Shcherbakov P.S., Smirnov G.V. Large deviations in linear control systems with nonzero initial conditions. Autom. Remote Control. 2015. V. 76, N 6. P. 957-976.
39. Saad Y. Iterative methods for sparse linear systems. SIAM, 2003. http://www-users.cs.umn.edu/~saad/books.html
40. Krasnoselskii M.A., Krein, S.G. Remark on the distribution of errors in the solution of a system of linear equations by means of an iterative process. Uspehi Matem. Nauk. 1952. V. 7, N 4(50). P. 157-161. (in Russian)
41. Bakushinsky A.B., Kokurin M.Yu. Iterative methods for ill-posed operator equations with smooth operators. Moscow: URSS, 2002.
42. Gower R., Richtarik P. Randomized Iterative Methods for Linear Systems. SIAM. J. Matrix Anal. & Appl. 2015. V. 36(4). P. 1660-1690. arXiv:1506.03296
43. Gower R., Richtarik P. Stochastic Dual Ascent for Solving Linear Systems. e-print. 2015. arXiv:1512.06890
44. Gower R., Richtarik P. Randomized Quasi-Newton Updates are Linearly Convergent Matrix Inversion Algorithms . e-print. 2015. arXiv:1602.01768
Поступила в редакцию 23.03.2016
