CC BY
20
УДК 519.85 ББК 22.18
ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ ДЛЯ ЗАДАЧ УСЛОВНОЙ ОПТИМИЗАЦИИ БОЛЬШОЙ РАЗМЕРНОСТИ С ДЕКОМПОЗИЦИЕЙ ОГРАНИЧЕНИЙ
Величко А. С.1
(Институт автоматики и процессов управления ДВО РАН, Владивосток)
Для класса задач условной оптимизации рассматриваются параллельные алгоритмы на основе метода проекции градиента и декомпозиции ограничений задачи условной оптимизации. Использование предлагаемого подхода продемонстрировано для класса задач линейного программирования большой размерности. Эффективность алгоритма показана на специальном вычислительно трудоемком наборе тестовых задач.
Ключевые слова: параллельный алгоритм, метод проекции градиента, декомпозиция, большая размерность.
Введение
В работе [5] для решения задач условной оптимизации рассматривался алгоритм двойственных отсечений, на основе которого строились его различные модификации и параллельные реализации. Разработанные алгоритмы использовались для эффективного решения задач линейного программирования большой размерности путем их декомпозиции по ограничениям, в частности при решении задач двухэтапного стохастического программирования.
Несмотря на вычислительную эффективность, данный подход ограничен в своем использовании возможностью решения
1 Андрей Сергеевич Величко, кандидат физико-математических наук, научный сотрудник (vandre@dvo.ru). 60
только структурированных двублочных задач, параллельный алгоритм предполагает только двухпроцессорный вариант выполнения. В работе [4] предложен альтернативный подход к распараллеливанию решения квадратичной оптимизационной задачи с ограничениями, основанный на двойственной постановке оптимизационной задачи и нелинейном методе Якоби.
Совершенствование численных методов является актуальной проблемой для решения задач с большим числом ограничений и/или переменных. Для задач линейного программирования классические алгоритмы типа симплекс-метода имеют неполиномиальную оценку сложности, а известные алгоритмы с полиномиальной вычислительной сложностью и линейной оценкой скорости сходимости при практическом использовании для задач большой размерности могут демонстрировать очень медленную сходимость по числу итераций и времени работы алгоритма.
В данной статье рассматриваются алгоритмы решения задач условной оптимизации на основе метода проекций градиента и декомпозиции ограничений. Этот подход является более универсальным по сравнению с построением двойственных отсечений [5], поскольку, во-первых, он применяется при решении задач условной оптимизации с ограничениями общего вида и произвольным числом подмножеств множества ограничений исходной задачи, во-вторых, - может использоваться для решения нелинейных задач и позволяет реализовать различные модификации последовательных и параллельных алгоритмов.
1. Метод проекций градиента для задач условной оптимизации
Рассмотрим экстремальную задачу
(1) min h(x),
xeV
где h( ) - выпуклая дифференцируемая функция на допустимом
множестве V = П Ci,Im = 1 , m, которое представляет собой
ie/m
пересечение конечного числа выпуклых замкнутых множеств Ci
в Rn. Будем в дальнейшем считать, что V = 0 и рассматриваемая задача (1) имеет решение.
Класс методов проекции градиента, применяемый для решения задачи (1), заключается в построении последовательности точек xs+1 = Pv(xs — Asgs), s = 0,1,..., где Pv - оператор проекции точки на множество V, вектор —gs - антиградиент функции h в точке xs в предположении конечности и дифференцируемо-сти функции h на множестве V. Вектор xs+1 является решением задачи поиска проекции точки на выпуклое замкнутое множество
и xs+1 = arg min ||x — xs + Asgs||2.
xev
Начальное приближение x0 и шаговые множители As лучше выбирать с учетом специфики задачи.
Пристальный интерес к задаче нахождения проекции точки на множество, поиска точки, принадлежащей пересечению выпуклых множеств (convex feasibility problem) и проекции на пересечение множеств (best approximation problem) появился еще в 60-х годах XX века. К наиболее известным работам, повлиявшим на развитие этого класса численных методов оптимизации, можно отнести работы Л.М. Брэгмана [2], И.И. Еремина [3], Б.Т. Поляка [9], Е.Г. Гольштейна [8], Baushcke H. [18], Censor Y. [20], Bertsekas D. [19].
Выбор определенной последовательности шаговых множителей As в конкретном классе задач может повлиять на скорость и характер сходимости метода проекций градиента. Рассматриваемый метод был сначала обоснован для шаговых множителей, зависящих от константы Липшица для градиента функции h [14]. В работе Б.Т. Поляка [15] рассматривался шаг в виде As = y(h(xs) — h(x*))/||gs||2, 0 < 7 < 2, где x* - оптимальное решение. Зачастую в литературе используется последовательность шагов вида As = a/(b + s) [12, 14], где a, b - выбираемые положительные парамеры. В работе [13] было показано, что метод проекций градиента можно использовать для более общего случая оператора Pv и последовательности шаговых множителей As.
2. Декомпозиция ограничений и блочные задачи линейного программирования
Задачи линейного программирования могут рассматриваться как частный случай задачи (1), когда множество V = {x : aix ^ qi, i = 1, M} - пересечение линейных полупространств, h(x) = cx и градиент gs функции h на множестве V равен вектору с, т.е. задача (1) принимает вид
(2) min{cx : aix ^ qi, i = 1, M}.
x
Класс задач линейного программирования и задач с ограничениями типа линейных неравенств является важным объектом для разработки эффективных численных методов оптимизации, поскольку зачастую на промежуточных шагах алгоритмов для нелинейных или недифференцируемых задач требуется решение такого класса вспомогательных задач. Рассматриваемый класс линейных задач взят лишь в качестве примера, наглядно демонстрирующего применяемую методику.
Рассмотрим два способа декомпозиции ограничений задачи (1).
Первый подход состоит в представлении множества V допустимых решений задачи (2) в виде пересечения множеств Ci = {x £ Rn : Dix ^ bi},i £ Im = {1,...,m}, содержащих различные ограничения исходной задачи. Здесь Di - матрицы коэффициентов левой части ограничений, bi - векторы.
При таком определении множества V алгоритм проекций градиента для решения задачи (2) представляет собой итерационный процесс xs+1 = Ps(xs — Лsgs), s = 0,1,..., где Ps - оператор проекции вектора xs — Л^ на множество Cis, индексы {is} последовательно принимают значения из множества Im, т.е. io = 1, ii = 2,..., im-1 = m, im = 1,... вектор градиента gs не зависит от xs и равен вектору c.
Вектор xs+1 = Ps(xs — Лsc) определяется в результате решения задачи
(3) min llx — xs + ЛscH2.
xeCi,
Второй подход к декомпозиции ограничений задачи (2) заключается в другом представлении допустимого множества V. Определим сначала векторы гг € Ега и д = г1 х ■ ■ ■ х гт € Етга. Пусть V = П Сг, где множества Сг = {д : ^ Ьг}, г € 1т
г€/т+1
и множество Ст+1 = {д : гг = г^- для всех г = j}.
В этом случае алгоритм проекций для решения задачи (4) можно представить в виде следующего итерационного процесса. На каждом шаге в = 0,1,... на первом этапе решаются задачи гг = Рг(д5 — Л^д5), где Рг - оператор проекции вектора (д5 — Л^д5) на множество Сг, г € 1т, вектор д5 = д5 х ■ ■ ■ х дт, где д| - градиент функции Л, в точке г|, он не зависит от д5 и равен вектору с. Вектор г представим в виде декартова произведения Р1 х ■ ■ ■ х Рт.
В результате решения задач первого этапа оптимальный вектор г* для конкретного г будет определяться только оптимальным значением компоненты р* и не будет зависеть от р для всех j = г, поскольку множество Сг определяется только г-й компонентой декартова произведения, составляющего вектор д5. Отметим, что векторы Рг определяются в результате решения задач (3) для = г и х = рг. В силу специфики определения множества Сг возникающие задачи (3) на первом этапе каждого шага рассматриваемого алгоритма могут решаться параллельно для разных г. На втором этапе шага в алгоритма полагаем д5+1 =
Рт+1(Р1 х ••• х рт).
Рассмотрим задачу блочного линейного программирования с непустым допустимым множеством и конечной целевой функцией на множестве допустимых решений:
т
(4) /(у*) = Ш1п£ /¿(у), /¿(у) = тт{Сг^ : Л*** ^ ^ — Вгу},
г=1
где г € 1т,у - вектор связывающих переменных; ¿г - векторы блочных переменных в подзадачах определения функций /г(у); Аг, Вг - матрицы; ^г - векторы соответствующих размерностей. Обозначим допустимое множество задачи (4) за
Рассмотрим более подробно первый способ представления
допустимого множества задачи (4) в виде пересечения множеств Сг.
Пусть Сг = {(гг х у) : ^ — Вгу} - замкнутые выпуклые множества, г € /т. По построению множеств Сг очевидно, что = V = П Сг. Определим вектор х = (21, г2,..., 2т, у)',
г€/т
где ' - операция транспонирования, тогда в обозначениях задачи
т
(1) имеем Л,(х) = ^ /¿(у) = сх, с = (с1, С2,..., Ст, 0)'.
г=1
Рассмотрим задачу (3) для некоторого г = г8. Поскольку множество Сг содержит ограничения только на векторы тогда при решении задачи (3) будут определены оптимальные значения векторов ~г, Уг и = У+1 = уг, г € /т, а = г? — Л5с для ^ € /т\{г}.
Рассмотрим второй подход к описанию допустимого множества Рт задачи (4).
Определим множества Сг = {({гг}, {иг}) : ^ —
Вгиг},г € /т и множество Ст+1 = {({гг}, {иг}) : и1 = и2 = ■ ■ ■ = ит}. Очевидно, что = V = П Сг.
г€/т+1
Обозначим X = (г1, г2,..., 2т, и1, и2,..., ит)', тогда в обо-
т
значениях задачи (1) имеем Л,(Х) = ^ /г(иг) = сХ, здесь с =
г=1
(с1, с2, . . . , ст, 0, . . . , 0)'.
Рассмотрим задачу (3) для некоторого г = г8. Поскольку множество Сг содержит ограничения только на векторы иг, тогда получим = — Л5с и = и? для ] € /т\{г}. Решение задачи (3) для множества Ст+1 приводит к определению опти-
т
т
г=1
мальных иг = тт X] и| и = 2?. Значит полагаем = и и
= 2г.
3. Разработка параллельных алгоритмов и вычислительные эксперименты
В работах [1,3, 16, 19] для задач решения систем неравенств, в том числе соответствующих необходимым и достаточным усло-
65
виям экстремума для задач линейного программирования, предлагается осуществлять проекции на отдельные ограничения типа равенств-неравенств или на систему ограничений типа равенств. Такие задачи имеет аналитическое решение, в результате получающиеся алгоритмы являются привлекательными с точки зрения простоты реализации на ЭВМ. Это направление использования проекционных итерационных методов берет свое начало с работы С. Качмажа 1937 г. (перевод на английский язык см. в [22]).
В работе Д.В. Долгого и Е.А. Нурминского [10] алгоритмы проекций рассматривались для специального представления допустимого множества задачи (1) в виде внешне заданных полиэдров, а ускоренный вариант метода предполагает представление допустимого множества в виде пересечения двух множеств, которые динамически меняются.
В данной статье, в отличие от приведенных выше работ, во-первых, предлагается рассматривать оптимизационные задачи с различным числом ограничений не менее двух, не переходя к эквивалентной системе равенств-неравенств. Во-вторых, используются в общем случае ограничения типа неравенств, что не позволяет получить простого аналитического выражения для решения, на чем основаны метода типа Качмажа. В третьих, при па-раллелизации алгоритма рассматриваются представления допустимого множества задачи (1) в наиболее общем представлении: в виде пересечения подмножеств, соответствующих нескольким ограничениям исходной задачи, каждое из которых есть пересечение линейных полупространств, соответствующих отдельным ограничениям.
Параллельный алгоритм проекций для задачи (1) в рамках первого подхода к представлению ограничений имеет вид
т
х5+1 = X] Р.(х5 — Л5д5), в = 0,1,..., где Р. - оператор про-.7 = 1
екции вектора х5 — на множество С., j = 1..., т, д5 -градиент функции Л в точке х5. Проекции Р. вектора х5 — на множества С. осуществляются на параллельных процессорах.
Выбор весовых коэффициентов может быть тривиальным
(" = 1/т) или более сложным [2, 9, 11].
Недостатком первого подхода к представлению допустимого множества задачи (4) является последовательный характер алгоритма по вектору связывающих переменных у, что затрудняет распараллеливание процесса вычислений.
Для параллельного решения задачи условной оптимизации (1) можно использовать второй подход с введением искусственного множества Ст+1. В этом случае вычисление проекций на отдельные множества Сг, т.е. решение задач (3), осуществляется параллельно на т процессорах. Процесс т + 1, реализующий решение задачи проекции на множество Ст+1, по сути осуществляет лишь пересчет искусственно введенных векторов переменных гг, уже полученных ранее в результате решения задач проекции на множества Сг,г € /т. Таким образом, каждая из т задач поиска проекции содержит только один блок ограничений, определяемых индексом г, что позволяет надеяться на уменьшение объема вычислений и количества используемой памяти для хранения данных по каждому блоку ограничений по сравнению с решением задачи (4). С другой стороны, платой за уменьшение количества ограничений при решении каждой из задач проекции является необходимость решения квадратичной, а не линейной задачи.
Используя второй подход к декомпозиции ограничений блочной задачи условной оптимизации (4) параллельный алгоритм имеет вид х5+1 = Р.(х5 — Л5с),в = 0,1,..., где Р. - оператор проекции вектора х5 — Л5с на множество С., j € 1т+1. На (в + 1)-м шаге определяются оптимальные значения в ре-
зультате решения задачи проекции вектора х5 — Л5с на множество Сг. Тогда для г € 1т получаем ¿г5+1 = одя* + (1 — — Л5с),
и с учетом задачи проекции на множество Ст+1 имеем =
т
+ £ "и + ^тт1 £ и*.
.=г .=1
В вычислительных экспериментах на задачах линейного программирования с ограничениями типа неравенств (2) использовались случайно генерируемые данные разных размерно-
стей для тестовой задачи, описанной в работе [6, 13]. В этих вычислительно-трудоемких тестах система неравенств строилась как набор m случайных опорных плоскостей к n-мерной сфере с центром в случайной точке x° с нормально распределенными координатами и радиусом r = О||ж°||, О G (0,1).
Алгоритм реализован [7] с использованием свободно распространяемого программного обеспечения Octave [21]. Параллельная реализация алгоритма для многопроцессорной вычислительной архитектуры с распределенной памятью реализована с использованием свободно распространяемого программного обеспечения: библиотеки функций интерфейса передачи сообщений OpenMPI [23] и пакета MPI для Octave [24]. Расчеты проводились на вычислительных мощностях центра коллективного пользования ДВО РАН «Дальневосточный вычислительный ресурс» [17]. Для решения задач проекции на множество использовался алгоритм решения задачи с квадратичной целевой функцией на линейных ограничениях (функция qp в Octave).
На рис. 1 показана вычислительная сложность алгоритма, которая близка к полиномиальной.
Рис. 1. Вычислительная сложность алгоритма
Характер сходимости параллельного алгоритма при простом
выборе шаговых множителей вида Л5 = и весовых коэффициентов щ = 1/т показан на рис. 2. Он демонстрирует важный эффект ускорения сходимости на последних шагах алгоритма по мере приближения к оптимуму, чего лишены классические последовательные алгоритмы проекции градиента. По оси ординат в логарифмической шкале расположены значения = |Л,(х5) — Л,(х*)| абсолютного отклонения точного решения от приближенного, получаемого на итерациях алгоритма.
Рис. 2. Сходимость алгоритма. в - номер итерации, - абсолютная погрешность приближенного решения
Величина эффективности Ега параллельного алгоритма определяется как отношение времени выполнения алгоритма на одном процессоре к времени выполнения алгоритма на п параллельных процессорах.
На рис. 3 показана фактическая величина показателя эффективности Дга параллельного алгоритма в зависимости от числа используемых процессоров п.
Рис. 3. Эффективность параллельного алгоритма: п - число процессов, Ега - отношение времени работы алгоритма с одним процессом к времени выполнения параллельного алгоритма
4. Выводы
В работе рассматривались подходы к построению параллельных алгоритмов на основе совмещения двух идей: декомпозиции ограничений задач условной оптимизации и организации вычислений на основе метода проекций градиента, применяемых для решения независимых задач оптимизации на подмножествах множества ограничений исходной задачи.
Разработан и реализован параллельный алгоритм для многопроцессорных вычислительных систем с распределенной памятью с использованием функций интерфейса передачи сообщений ОрепМР1. Хорошие показатели качества алгоритма показаны на специальном вычислительно-трудоемком наборе тестовых исходных данных для задачи линейного программирования большой размерности по числу ограничений.
Наиболее эффективной оказывается стратегия декомпозиции допустимого множества V в задаче (1) на небольшое число подмножеств. Декомпозиция по отдельным ограничениям задачи, когда подмножества содержат по одному ограничению и число про-
цессов принимает большое значение, не является эффективной с точки зрения показателя Дп, как это может показаться на первый взгляд.
Литература
1. БЕРДНИКОВА Е.А., ЕРЕМИН И.И., ПОПОВ Л.Д. Распределенные фейеровские процессы для систем линейных неравенств и задач линейного программирования // Автоматика и телемеханика. - 2004. - №2. - С. 16-32.
2. БРЭГМАН Л.М. Релаксационный метод нахождения общей точки выпуклых множеств и его применение для задач выпуклого программирования // Журн. вычисл. матем. и матем. физики. - 1967. - Т. 7, №3. - С. 620-631.
3. ВАСИН В.В., ЕРЕМИН И.И. Операторы и итерационные процессы фейеровского типа. Теория и приложения. -Екатеринбург: УрО РАН, 2005. - 200 с.
4. ВЕЛИЧКО А.С. Двойственный алгоритм для задач регуляризации с недифференцируемыми стабилизаторами // Вычислительные технологии. - 2014. - №2. - С. 14-20.
5. ВЕЛИЧКО А.С. Об алгоритме двойственных отсечений для задачи двуэтапного стохастического программирования // Известия ВУЗов. Математика. - 2006. - №4. -С. 78-81.
6. ВЕЛИЧКО А.С. О выборе шага в проективных алгоритмах для задач линейного программирования большой размерности // Дальневосточный математический журнал. -2012. - №2.-С. 160-170.
7. ВЕЛИЧКО А.С. Решение задач оптимизации методом последовательных проекций с различными способами выбора начального приближения и шаговых множителей : св. о гос. рег. прог. для ЭВМ Росс. Фед. 2011614018 от 24.05.11; заявл. 06.04.11; опубл. 20.09.2011. В бюлл.: ЯИ ОБПБТ, №3. С. 319.
8. ГОЛЬШТЕЙН Е.Г., ТРЕТЬЯКОВ Н.В. Модифицированные функции Лагранжа. Теория и методы оптимизации. -М.: Наука, 1989. - 400 с.
9. ГУРИН Л.Г., ПОЛЯК Б.Т., РАЙК Э.В. Методы проекций для отыскания общей точки выпуклых множеств // Журн. вычисл. матем. и матем. физики. - 1967. - Т. 7, №6. - С. 1211-1228.
10. ДОЛГИЙ Д.В., НУРМИНСКИЙ Е.А. Ускоренный параллельный метод проекций для решения задачи о наименьшем расстоянии // Вычислительные методы и программирование. - 2006. - Т. 7. - С. 273-277.
11. КАРМАНОВ В.Г. Оценки сходимости итерационных методов минимизации // Журн. вычисл. матем. и мат. физики. - 1974. - Т. 14, №1. - С. 3-14.
12. НЕСТЕРОВ Ю.Е. Введение в выпуклую оптимизацию. -М.: Издательство МЦНМО, 2010. - 280 с.
13. НУРМИНСКИЙ Е.А. Фейеровские процессы c малыми возмущениями // Доклады АН. - 2008. - Т. 422, №5. -С. 601-605.
14. ПОЛЯК Б.Т. Введение в оптимизацию. - М.: Наука, 1983. -384 с.
15. ПОЛЯК Б.Т. Минимизация негладких функционалов // Журн. вычисл. матем. и матем. физики. - 1969. - Т. 9, №3. - С. 509-521.
16. ПОЛЯК Б.Т. Рандомизированные алгоритмы решения выпуклых неравенств // Стохастическая оптимизация в информатике. - 2005. - Т. 1, №1. - С. 150-169.
17. Центр коллективного пользования ДВО РАН «Дальневосточный вычислительный ресурс». - [Электронный ресурс] - URL: http://www.cc.dvo.ru (дата обращения: 15.01.2016).
18. BAUSCHKE H.H., BORWEIN J.M. On projection algorithms for solving convex feasibility problems // SIAM Review. - 1996. - Vol. 38, №3. - P. 367-426.
19. BERTSEKAS D.P., TSITSIKLIS J.N. Parallel and
Distributed Computation: Numerical Methods. - Athena Scientific, 1997. - 735 p.
20. CENSOR Y., ZENIOS S.A. Parallel optimization: theory, algorithms, and applications. - Oxford University Press, 1997. - 541 p.
21. GNU Octave. - [Электронный ресурс] - URL: http://www.octave.org (дата обращения: 15.01.2016).
22. KACZMARZ S. Approximate solution of systems of linear equations // International Journal of Control. - 1993. -Vol. 57, №6. - P. 1269-1271.
23. OpenMPI. - [Электронный ресурс] - URL: http://www.open-mpi.org (дата обращения: 15.01.2016).
24. MPI package for Octave. - [Электронный ресурс] -URL: http://octave.sourceforge.net/mpi/ (дата обращения: 15.01.2016).
PARALLEL ALGORITHMS WITH FEASIBLE SET PARTITIONING FOR LARGE-SCALE OPTIMIZATION PROBLEMS
Andrei Velichko, Institute for Automation and Control Processes of FEB RAS, Vladivostok, Cand.Sc., researcher (vandre@dvo.ru).
Abstract: We consider parallel algorithms for the class of constrained optimization problems. The algorithms are based on the two ideas: gradient projection method and decomposition of the set of constraints into disjoint subsets with different number of constraints. The proposed approach is used for a class of large-scale linear programming problems. A number of simulations was made to analyze computational complexity, convergence speed and the parallelization efficiency. The algorithms show a good performance for the special special computationally difficult benchmark problem. The most efficient decomposition strategy is to partition the feasible set into a relatively small number of subsets. The decomposition by separate constraints demonstrates less efficiency of parallelization.
Keywords: parallel algorithm, gradient projection method, decomposition, large-scale problem.
Статья представлена к публикации членом редакционной коллегии В.И. Зоркальцевым.
Поступила в редакцию 01.02.2016. Дата опубликования 31.07.2016.
