Параллельный двухуровневый индексный алгоритм глобальной оптимизации Текст научной статьи по специальности «Математика»

Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 2012, № 5 (2), с. 208-213

УДК 541.186

ПАРАЛЛЕЛЬНЫЙ ДВУХУРОВНЕВЫЙ ИНДЕКСНЫЙ АЛГОРИТМ ГЛОБАЛЬНОЙ ОПТИМИЗАЦИИ

© 2012 г. С.В. Сидоров

Нижегородский госуниверситет им. Н.И. Лобачевского

81^гоу. sergey@gmail.com

Поступила в редакцию 25.09.2012

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

Ключевые слова: глобальная оптимизация, индексный метод, параллельное программирование, развертки типа Пеано.

1. Введение

Данная работа продолжает развитие известного подхода к минимизации многоэкстремальных функций при невыпуклых ограничениях, описанного в работах [1-7, 9, 10, 13, 14] и получившего название индексного метода глобальной оптимизации. Подход основан на раздельном учете каждого ограничения задачи и не связан с использованием штрафных функций. При этом решение многомерных задач сводится к решению эквивалентных им одномерных. Соответствующая редукция основана на использовании кривых Пеано (называемых также развертками Пеано), однозначно отображающих единичный отрезок вещественной оси на гиперкуб. Используется схема построения множества кривых Пеано («вращаемые развертки»), которую можно применять при решении задачи на кластерных системах с десятками и сотнями процессоров. Данная схема распараллеливания дополняется одновременной обработкой нескольких точек испытания на каждом узле кластера с использованием многоядерности (двухуровневая параллельная схема). Приведены результаты экспериментов по оценке масштабируемости двухуровнего параллельного индексного метода. Эксперименты выполнены на вычислительном кластере ННГУ им. Н.И. Лобачевского, установленном в ходе выполнения национального проекта «Образование».

2. Постановка задачи

Рассмотрим задачу глобальной оптимизации вида

ер* = ф(у*) = min (ф(у): у & D, gj(y) < 0, 1 <j < т}, D = {у е lf: а{<у{< bu 1 <i<N), (1)

где целевая функция ср(у) (в дальнейшем обозначаемая также gm+\(y)) и левые части ограничений gj(y), 1 <j<m, удовлетворяют условию Липшица с соответствующими константами Lj, \ <j <т + \,& именно

\ё}(У1)-ё}(У2)\^Ц\У\-У2\, l<j<rn+ 1 ,y\j>2^D.

Используя кривые типа развертки Пеано y(x), однозначно отображающие отрезок [0,1] на TV-мерный гиперкуб Р

P={yeBN: -2~х <Уг< 2~\ 1 <i<N} =

= {у(х): 0<х<1}, исходную задачу можно редуцировать к следующей одномерной задаче:

ф(у(х*)) = тт{ф(у(х)): х е [0,1],

gj(y(x))<0, 1 <j<m}. (2)

Рассматриваемая схема редукции размерности сопоставляет многомерной задаче с липши-цевой минимизируемой функцией и липшице-выми ограничениями одномерную задачу, в которой соответствующие функции удовлетворяют равномерному условию Гельдера (см. [2]), то есть

х',х" є [0,1], 1 <] < т + 1, где N есть размерность исходной многомерной задачи, а коэффициенты Кі связаны с константами Липшица Ц исходной задачи соотношениями К<4Ь] У^.

Различные варианты индексного алгоритма для решения одномерных задач и соответствующая теория сходимости представлены в работах [1, 4, 7].

3. Использование множественных отображений

Редукция многомерных задач к одномерным с помощью разверток имеет такие важные свойства, как непрерывность и сохранение равномерной ограниченности разностей функций при ограниченности вариации аргумента. Однако при этом происходит потеря части информации

о близости точек в многомерном пространстве, так как точка х є [0, 1] имеет лишь левых и правых соседей, а соответствующая ей точка у(х) є

гМ

є к имеет соседей по 2 направлениям, а при использовании отображений типа кривой Пеано близким в М-мерном пространстве образам у ’, у" могут соответствовать достаточно далекие прообразы X, х" на отрезке [0, 1]. Как результат, единственной точке глобального минимума в многомерной задаче соответствует несколько (не более 2м) локальных экстремумов в одномерной задаче, что, естественно, ухудшает свойства одномерной задачи.

Сохранить часть информации о близости точек позволяет использование множества отображений

ВД = (у^х),..., /(х)} (3)

вместо применения единственной кривой Пеано у(х) (см. [3, 5]). Каждая кривая Пеано уг(х) из УЬ(х) может быть получена в результате некоторого сдвига вдоль главной диагонали гиперинтервала В. Таким образом сконструированное множество кривых Пеано позволяет получить для любых близких образов у' , у", отличающихся только по одной координате, близкие прообразы х', х" для некоторого отображения уг(х).

3.1. Вращаемые развертки. К числу недостатков ставшей уже классической схемы построения множественных разверток (далее будем называть их сдвиговыми развертками, или С-развертками) можно отнести, во-первых, наличие дополнительного ограничения, порождающего сложную допустимую область на одномерных отрезках. А во-вторых, при построении С-разверток число разверток Ь (а следова-

тельно, и число параллельно решаемых задач) зависит от требуемой точности е поиска решения задачи. С позиции экономии вычислительных ресурсов невыгодно использовать число разверток большее, чем Г^2 (е-1)~|. Например, при решении задачи с точностью 10 3 по координате целесообразно выбирать число разверток не больше 10.

Преодолеть эти недостатки, сохранив информацию о близости точек в М-мерном пространстве, позволяет схема построения множественных отображений, предложенная в [8]. Отличительной чертой этой схемы является построение множества кривых Пеано не с помощью сдвига вдоль главной диагонали гиперкуба, а поворотом развертки вокруг начала координат. При этом найдется отображение у(х), которое точкам многомерного пространства у’, у", которым при исходном отображении соответствовали достаточно далекие прообразы на отрезке [0, 1], будет сопоставлять более близкие прообразы х’, х".

Развертки, порождаемые в соответствии с новой схемой, будем называть вращаемыми развертками, или В-развертками.

Максимальное число различных поворотов развертки, отображающей М-мерный гиперкуб на одномерный отрезок, составляет 2м. Использование всех из них является избыточным, требуется выбрать лишь часть из всех возможных вариантов. В предложенной схеме преобразование развертки осуществляется в виде поворота на угол ±п/2 в каждой из координатных плоскостей. Число подобных пар поворотов определяется числом координатных плоскостей пространства, которое равно СМ = N(N -1) / 2 , а общее число преобразований будет равно N(N-1). Учитывая исходное отображение, приходим к заключению, что данный способ позволяет строить до М(М-1)+1 развертки для отображения М-мерной области на соответствующие одномерные отрезки. При этом дополнительное ограничение, которое возникает при построении С-разверток [3], отсутствует. В случае необходимости данный способ построения множества отображений может быть легко «от-масштабирован» для получения большего (вплоть до 2м) числа разверток.

3.2. Параллельный двухуровневый индексный алгоритм решения задач глобальной оптимизации. Ввиду того, что современные высокопроизводительные системы имеют иерархическую архитектуру - набор многоядерных процессоров, соединенных локальной сетью, -

в работе предлагается параллельный двухуровневый индексный алгоритм, учитывающий эту особенность.

На уровне процессоров будет организовано параллельное вычисление очередных точек испытаний на отдельных вычислительных ядрах.

На уровне межпроцессорного взаимодействия будет организована асинхронная схема обмена поисковой информацией.

3.2.1. Описание алгоритма параллельного индексного метода на межпроцессорном уровне. Использование множественных отображений (3) позволяет решать задачу (1) путем параллельного решения Ь задач вида (2) на наборе отрезков [0, 1]. Каждая одномерная задача решается на отдельном процессоре с использованием развертки у'5, 1 < 5 < Ь. Результаты иск

пытания в точке х , полученные конкретным процессором для решаемой им задачи, интерпретируются как результаты испытаний во всех остальных задачах (в соответствующих точках

к1 кЬ\

х , ., х ) и рассылаются другим процессорам. При таком подходе испытание в точке хк є є [0, 1], осуществляемое в 5-й задаче, состоит в последовательности действий:

1. Определить образ ук = у5(хк) при соответствии у5(х).

2. Проинформировать остальные процессоры

о начале проведения испытания в точке ук (бло-

кирование точки у ).

3. Вычислить ве значения индекса V < т определяются условия-

(у1(хк1)) = ... = gv(yL(xkL)) = /.

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

3.2.2. Описание алгоритма параллельного индексного метода на уровне отдельного процессора. Начальная итерация осуществляется в произвольной точке (начальные точки для всех процессоров задаются различными). Выбор точки хц+1, ц > 1, любого последующего испытания определяется следующими правилами.

Правило 1. Изъять из очереди, закрепленной за данным процессором, записанные для него результаты, включающие множество Yq = {у: 1 < г < sq} точек итераций в области Б и вычисленные в них значения индекса и величин из (3.1). Определить множество Хц = {хг: 1 < г < ¿ц} прообразов точек множества Уц при соответствии у6'(х). Мощность множеств Уц и Хц есть

> q,

3. Вычислить величины gl(yk),...,gV(yk), где

g}(yk) < 0.1 <7 < V, gv(yk) > 0, V < т. Выявление первого нарушенного ограничения прерывает испытание в точке у к. В случае,

к

когда точка у допустима, испытание включает вычисление значений всех функционалов задачи, при этом значение индекса принимается равным величине V = т + 1. Тройка

у*(хк), V = v(xk), / = gv(ys(xk)) (4)

является результатом испытания в точке хк.

4. Определить прообразы хк1 є [0, 1], 1 < I < Ь,

к

точки у и интерпретировать испытание, проведенное в точке у є В, как проведение испытаний в Ь точках

к1 к1 /С\

х , ., х (5)

с одинаковыми результатами

v(xk1) = . = v(xkL) = v(xk),

Проинформировать остальные процессоры о

к

результатах испытания в точке у , разослав им

кк тройки (у , V, 2 ).

Каждый процессор имеет свою копию программных средств, реализующих вычисление

так как эти множества содержат точки, вычисленные на данном процессоре и полученные с других процессоров.

Правило 2. Точки множества итераций {х1} и Хц

перенумеровать нижними индексами в порядке увеличения значений координаты

0 = х0 < х1 < ... < Хг < ... < хк < Хк+1 = 1, (6)

где к = + 1, и сопоставить им значения гг =

= gv(xг), V = v(xг■), 1 < г < к, вычисленные в этих точках. При этом индекс блокированной точки хг (т.е. точки, в которой начато проведение испытания другим процессором) полагается равным -1, т.е. v(xг■) = -1, значение гг является неопределенным. Точки Х0, Хк+1 введены дополнительно для удобства последующего изложения, индекс данных точек полагается равным -2, т.е. v(x0) = v(xk+l) = -2, а значения г0, гк+1 являются неопределенными.

Правило 3. Провести классификацию номеров г, 1 < г < к, точек из ряда (4) по числу ограничений задачи, выполняющихся в этих точках, путем построения множеств

1-2 = {0, к + 1},

1-1 = {г: 1 < г < к, v(xг■) = -1},

!,= {/: 1 < г < к, v(xг) = V}, 1 < V < т + 1, содержащих номера всех точек хг, 1 < г < к, имеющих индексы, равные одному и тому же значению v. Определить максимальное значе-

ми

ние индекса

М = шах^ = v(xг■), 1 < г < к}.

Правило 4. Вычислить текущие нижние границы

^= шах ] 1 ^ г'1/ м : ] < г, ] е Iv, г е Iv\ (7)

I (Х - Х1 ) ]

для относительных разностей функций gv, 1 < < v < т + 1. Если множество Д, содержит менее двух элементов или если из (7) оказывается равным нулю, то принять = 1. Из (7) следует,

что оценки являются неубывающими начиная с момента, когда (7) порождает первое положительное значение ц„.

Правило 5. Для всех непустых множеств Д,,

1 < v < т + 1, вычислить оценки

-sv, v< M,

min{gv(): i єіv}, v= M,

>v\ г/ ' V)

где вектор

Єк = (єь ..., Єт), (8)

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

Правило 6. Для каждого интервала (х-1, х),

1 < г < к + 1, вычислить характеристику К(г), где

К(г) = Дг + (ВДг-і£ - 2-^ + ^ - 2^)2

і 2 2 *

ГКЛ-v = v( х-і) = v( х),

rvKv

z _ Z

R(i) = 2Д -4^-v, v(х-і) < v(х) = v,

rvKv

*

z ___ z

R(i) = 2Д _4-i_l—^, v = v(х_і) > v(х), rvKv

Д = (х _ х_і)17N.

Величины rv > 1, 1 < v < m + 1, являются параметрами алгоритма.

Правило 7. Определить 5 (s _ число вычислительных ядер процессора) интервалов (х, х+1)ь

1 < k < s, которым соответствует s наибольших характеристик

Rk(t) = max {R(i): 1 < і < k}, 1 < k < s. (9)

Правило 8. Провести параллельно на s вычислительных ядрах очередное испытание в серединной точке интервала (хь х1+1)ь 1 < k < s, если индексы его концевых точек не совпадают, то есть

хЧ+1= (X + Х_1)k , v(Xt_1) ^v(xt ),1 < k < s.

В противном случае провести испытание в точке

q+1 = (хі + хі _1)k

Xі =

_sign( zt _ zt _1)k-1-2 2rv

I zt _ zt _1 Ik

Результаты испытания занести в очередь, закрепленную за данным процессором. Увеличив ц на единицу, перейти к новой итерации.

Описанные правила можно дополнить условием остановки, прекращающим испытания, если Д < е, где t из (9) и е > 0 имеет порядок желаемой покоординатной точности в задаче.

Ввиду того, что поисковая последовательность строится по принципу, описанному в [1], сходимость метода на каждом вычислительном узле обеспечивается теоремой из [11].

4. Реализация в программной системе С1оЬа1ЕхреГ+

На кафедре математического обеспечения ЭВМ ННГУ им. Н.И. Лобачевского разрабатывается программный комплекс параллельной глобальной оптимизации О1оЬа1Бхрег1;, одной из особенностей которого является быстрая работа с большими объёмами поисковой информации (т. е. совокупностью всех точек, в которых вычислялись функции задачи), включающая собственную подсистему подкачки, учитывающую особенности алгоритмов глобального поиска.

Автором за основу был взят комплекс ОЬЬаШхрей и разработана новая программная система О1оЬа1Бхре1!;+, предназначенная для решения трудоемких прикладных задач глобальной оптимизации с использованием высокопроизводительных вычислений.

5. Результаты экспериментов

Исследование эффективности и масштабируемости параллельного двухуровневого индексного алгоритма с вращаемыми развёртками (ЛО-Я) проведём при решении многомерных задач, построенных на основе функции Растри-гина [12].

Рассмотрим задачу условной многоэкстремальной оптимизации на основе одной из функций Растригина, описанную в [12]:

f(у) = X( у- _10cos (2пУ )) ^ min

i=1

§1(у) = X Уі _ °.5 < 0,

i =1

g-(у) = j^y-2 _ 1.0 < 0,

i =1 n

g3(У) = X sin (іпУі ) _0.3 < 0,

(10)

i =1

v (х_1) = v (х),1 < k < s.

у є [_1.0,1.0]N, где N - размерность поискового пространства.

N

Как видно из (10), глобальный минимум в этой задаче достигается в точке (0,...,0)N со значением f = 0.0.

Результаты решения задачи при N = 8 последовательным и параллельным индексным методом (56 процессоров, 1 развертка на процессор, 4 ядра) приведены в таблице 1. В обоих случаях метод завершал работу по достижении заданной относительной точности р = 0.0035.

Из табл. 1 видно, что параллельный алгоритм AG-R имеет ускорение 150.37 по числу итераций, 149.95 - по числу испытаний и 128.8 - по количеству вычислений критерия.

раметры надёжности (г=3.0 для ЛО, г=2.5 - для ЛО-Я) и разную относительную точность в условиях остановки.

Рассмотрим задачу безусловной многоэкстремальной оптимизации на основе одной из функций Растригина:

f (У) = Х (УІ - C0S(18Уг ))

г=1

(11)

у е[-1.0,1.0]Д,

где Д - размерность поискового пространства.

Минимум этой функции достигается в точке (0,.,0)Д с оптимальным значением / = —Д.

Таблица 1

Сравнение последовательного и параллельного алгоритма AG-R на задаче с ограничениями

Последовательный, 56 развёрток Параллельный, 56 процессоров, 1 развертка на процессор, 4 ядра

всего (максимум на один процесс) ускорение с поправкой на разное число испытаний*

Количество испытаний 609640 444356 (2967)

Количество вычислений критерия 63874 42722 (495)

Ускорение по испытаниям 205.43 149.95

Ускорение по вычислениям критерия 179.30 128.8

Время вычислений, с 93900.6 626.8

Ускорение по времени 150.37 150.37

Найденное минимальное значение целевой функции 0.016412 0.022182

Замечание. В задачах большей размерности сравнение последовательного и параллельного алгоритма ЛО-Я усложняется, поскольку в последовательном случае вся поисковая информация сосредоточена на одном вычислительном узле. Последовательному алгоритму, выполняемому на одном узле, не хватает ресурсов оперативной памяти, в результате чего начинается использование медленной внешней памяти. В параллельном алгоритме все данные распределяются по процессам, что позволяет не использовать внешнюю память.

Сравним последовательный алгоритм ЛО (с одной развёрткой) и параллельный ЛО-Я (с р развёртками, гдер - число процессоров, 4 ядра). В этом случае необходимо брать различные па-

Результаты минимизации мно

В данном эксперименте измерение ускорения по числу испытаний не представляется возможным в силу различия в параметрах методов. Однако можно измерить ускорение по времени и эффективность нахождения глобального оптимума при одинаковом условии остановки по числу испытаний.

В табл. 2 приведены результаты безусловной минимизации задачи (11) при Д = 11, 15 и 20 для последовательного алгоритма и параллельного для 50 и 55 процессоров с 4 ядрами. Индексный метод с вращаемой разверткой, выполняемый на более чем 50 четырехъядерных процессорах, показывает близкое к линейному ускорение по времени выполнения вычислений, что доказывает его высокую эффективность.

Таблица 2

мерных функций Растригина

AG AG-R (50, 55 процессоров), 4 ядра

Кол-во испытаний Время, с Оценка f* Кол-во испытаний Время процессов, с Оценка f* Ускорение по времени

N = 11 5000000 544068.8 -10.993 4993599 2720. (55) -10.99 200.01

N = 15 3000000 605788.4 -14.838 3000016 3346.6(50) -14.92 181.7

N = 20 3000000 605956.7 -19.98 3000056 3172. (50) -19.04 191.06

* Поправка в данном случае - это умножение на величину p/s, где p - число испытаний, сделанных параллельным методом, а s - последовательным.

6. Заключение

В рамках данной работы реализована двухуровневая схема параллельных вычислений, позволяющая на уровне отдельных вычислительных узлов с общей разделяемой памятью исключить передачу данных, а на уровне всей вычислительной системы обеспечивающая асинхронное выполнение параллельных вычислений и балансировку вычислительной нагрузки. Проведено исследование эффективности предложенного алгоритма. Результаты наглядно подтверждают близкое к линейному ускорение.

Работа выполнена при поддержке Российского фонда фундаментальных исследований (номер проекта 11-07-97017-р_поволжье_а), а также гранта № 14.B37.21.0878.

Список литературы

1. Стронгин Р.Г. Численные методы в многоэкстремальных задачах (Информационно-статистические алгоритмы). М.: Наука, 1978.

2. Стронгин Р.Г. Поиск глобального оптимума. М.: Знание, 1990.

3. Стронгин Р.Г. Параллельная многоэкстремальная оптимизация с использованием множества разверток // ЖВМиМФ. 1991. Т. 31. № 8. С. 1173-1185.

4. Стронгин Р.Г., Баркалов К.А. О сходимости индексного алгоритма в задачах условной оптимизации с s-резервированными решениями // Математические вопросы кибернетики. М.: Наука, 1999. С. 273-288.

5. Strongin R.G., Sergeyev Ya.D. Global optimization with non-convex constraints. Sequential and parallel algorithms. Kluwer Academic Publishers, Dord-

recht, 2000.

6. Баркалов К.А., Стронгин Р.Г. Метод глобальной оптимизации с адаптивным порядком проверки ограничений // ЖВМиМФ. 2002. Т. 42. № 9. C. 1338-1350.

7. Баркалов К.А. Ускорение сходимости в задачах условной глобальной оптимизации. Н. Новгород: Изд-во Нижегородского гос. ун-та, 2005.

8. Баркалов К.А., Рябов В.В., Сидоров С.В. Использование кривых Пеано в параллельной глобальной оптимизации // Матер. IX Междунар. конф.-семинара «Высокопроизводительные параллельные вычисления на кластерных системах». Владимир, 2009. С. 44-47.

9. Стронгин Р.Г., Гергель В.П., Баркалов К.А. Параллельные методы решения задач глобальной оптимизации // Изв. вузов. Приборостроение. 2009. Т. 52. № 10. С. 25-32.

10. Баркалов К.А., Рябов В.В., Сидоров С.В. Масштабируемые параллельные алгоритмы глобальной оптимизации со смешанной локально-глобальной стратегией // Матер. Междунар. науч. конф. «Параллельные вычислительные технологии 2010». Уфа, 2010. С. 402-409.

11. Стронгин Р.Г., Маркин Д.Л. Минимизация многоэкстремальных функций при невыпуклых ограничениях // Кибернетика. 1986. № 4.

12. Törn A. and Zilinskas A. «Global Optimization». Lecture Notes in Computer Science, № 350, SpringerVerlag, Berlin, 1989.

13. Gergel V.P., Sergeyev Ya.D. Sequential and parallel algorithms for global minimizing functions with Lipschitzian derivatives // Computers & Mathematics with Applications. 1999. Vol. 37. № 4-5. P. 163-179.

14. Gergel V.P., Strongin R.G. Parallel computing for globally optimal decision making on cluster systems // Future Generation Computer Systems. 2005. Vol. 21. № 5. P. 673-678.

TWO-LEVEL PARALLEL INDEX ALGORITHM FOR GLOBAL OPTIMIZATION

S.V. Sidorov

We continue to develop information-statistical approach to minimization of multiextremal functions in the case of non-convex constraints. The proposed approach is called the index method of global optimization. The solution of multi-dimensional problems is reduced to their one-dimensional equivalents. This reduction is based on singlevalued Peano curve mapping a hypercube onto the unit segment on the real axis. The scheme of constructing a set of Peano curves («rotated evolvents») is also used which can be effectively applied to the problem solution on a cluster with tens and hundreds of processors. Special attention is paid to the two-level parallel index algorithm for global optimization that takes into account the hierarchical structure of modern cluster systems.

Keywords: global optimization, index method, parallel programming, Peano scans.