Рассмотрен параллельный алгоритм решения задач многоэкстремальной оптимизации. Описывается реализация алгоритма на современных вычислительных системах с использованием сопроцессора Xeon Phi. Рассматриваются два подхода к распараллеливанию алгоритма, учитывающие информацию о трудоемкости вычисления значений оптимизируемой функции. Приводятся результаты вычислительных экспериментов, полученные на суперкомпьютере «Лобачевский». Демонстрируется, что реализация для Xeon Phi опережает версию для CPU. Результаты подтверждают ускорение алгоритма с использованием Xeon Phi по сравнению с алгоритмом, реализованным только на CPU. Ключевые слова: глобальная оптимизация, многоэкстремальные функции, редукция размерности, параллельные алгоритмы, Intel Xeon Phi.
программных средств для современных вычислительных систем является актуальной задачей. Особый интерес представляет разработка схем распараллеливания, позволяющих эффективно использовать ускорители вычислений, такие, как сопроцессор Intel Xeon Phi.
В рамках проводимого исследования мы будем предполагать, что время проведения одного испытания (вычисления значения функции в области поиска) может значительно отличаться в различных решаемых задачах. Это определяет разработку разных подходов к распараллеливанию в задачах с «легкими» и «сложными» критериями. В статье приведено описание универсального подхода к распараллеливанию алгоритма глобального поиска, который охватывает оба этих случая. Указанный подход реализован в разработанной в ННГУ им. Н.И. Лобачевского параллельной программной системе решения задач глобальной оптимизации. 2. Параллельный алгоритм глобального поиска
Рассмотрим задачу поиска глобального минимума N-мерной функции ϕ(y) в гиперинтервале D = {y ∈ RN : ai ≤ yi ≤ bi ,1 ≤ i ≤ N}. Будем предполагать, что функция удовлетворяет условию Липшица с априори неизвестной константой L.
ϕ ( y*) = min{ϕ ( y) : y ∈ D} , ϕ ( y1) −ϕ ( y2 ) ≤ L y1 − y2 , y1, y2 ∈ D, 0 < L < ∞ .
Существует ряд способов адаптации эффективных одномерных алгоритмов для решения многомерных задач, см., например, методы диагонального [6] или симплексного [7] разбиения области поиска. В данной работе мы будем использовать подход, основанный на идее редукции размерности с помощью кривой Пеано y(x), непрерывно и однозначно отображающей отрезок вещественной оси [0,1] на n-мерный куб
{y ∈ RN : −2−1 ≤ yi ≤ 2−1,1 ≤ i ≤ N}= {y(x) : 0 ≤ x ≤ 1}.
Вопросы численного построения отображений типа кривой Пеано и соответствующая теория подробно рассмотрены в [4]. Здесь же отметим, что численно построенная развертка является приближением к теоретической кривой Пеано с точностью порядка 2−m , где m – параметр построения развертки. Использование подобного рода отображений позволяет свести многомерную задачу к одномерной задаче
ϕ ( y*) =ϕ ( y(x*)) = min{ϕ ( y(x)) : x ∈[0,1]}.
Важным свойством является сохранение ограниченности относительных разностей функции: если функция ϕ(y) в области D удовлетворяла условию Липшица, то функция ϕ(y(x)) на интервале [0,1] будет удовлетворять равномерному условию Гельдера
ϕ ( y(x1)) −ϕ ( y(x2 )) ≤ H x1 − x2 1 N , x1, x2 ∈[0,1] , где константа Гельдера H связана с константой Липшица L соотношением
H = 4Ld N , d = max{bi − ai :1≤ i ≤ N} .
Поэтому, не ограничивая общности, можно рассматривать минимизацию одномерной функции f (x) =ϕ(y(x)) , x∈[0,1] , удовлетворяющей условию Гельдера.
Рассматриваемый алгоритм решения данной задачи предполагает построение последовательности точек xk, в которых вычисляются значения минимизируемой функции zk = f(xk). Процесс вычисления значения функции (включающий в себя построение образа yk=y(xk)) будем называть испытанием, а пару (xk,zk) – результатом испытания. Множество пар {(xk,zk)}, 1≤k≤n составляют поисковую информацию, накопленную методом после проведения n шагов. В нашем распоряжении имеется p ≥1 вычислительных элементов и в рамках одной итерации метода мы будем проводить p испытаний одновременно. Обозначим k(n) общее число испытаний, выполненных после n параллельных итераций.
На первой итерации метода испытание проводится в произвольной внутренней точке x1
интервала [0,1]. Пусть выполнено n ≥1 итераций метода, в процессе которых были проведены
(
Шаг 1. Перенумеровать точки множества Xk = {x1,..., xk }∪{0}∪{1}, которое включает в себя граничные точки интервала [0,1], а также точки предшествующих испытаний, нижними индексами в порядке увеличения значений координаты, т.е.
Шаг 2. Полагая zi = f (xi ),1 ≤ i ≤ k , вычислить величины
,
где r > 1 является заданным параметром метода, а ∆i = (xi − xi−1)1 N .
ствии с формулами
Шаг 3. Для каждого интервала (xi−1, xi ),1 ≤ i ≤ k +1, вычислить характеристику в
соответz z
R(
M M R(i) = ∆i + (zi − zi−1)2 − 2 zi + zi−1 . , 1< i < k +1.
M 2∆i M
Шаг 4. Характеристики R(i),1≤ i ≤ k +1, упорядочить в порядке убывания
(
2 xk+ j = xtj +2xtj −1 − sign(zt j − ztj −1) 21r | ztj −µztj −1 | N , 1 < t j < k +1 (7) Алгоритм прекращает работу, если выполняется условие ∆t j ≤ ε хотя бы для одного номера t j , 1 ≤ j ≤ P ; здесь ε > 0 есть заданная точность. В качестве оценки глобально-оптимального решения задачи выбираются значения fk* = min f (xi ) , xk* = arg min f (xi ) 1≤i≤k 1≤i≤k (8) Теоретическое обоснование данного способа организации параллельных вычислений изложено в [5]. 3. Обобщенная схема редукция размерности
Одним из подходов к решению многомерных задач глобальной оптимизации является
сведение их к одномерным и использование эффективных одномерных алгоритмов глобального
поиска к редуцированной задаче. В предыдущем разделе была изложена идея редукции
размерности с использованием кривых Пеано. Ниже излагается обобщенный способ редукции
размерности, комбинирующий использование разверток и схему вложенной (рекурсивной)
оптимизации.
3.1 Рекурсивная схема редукции размерности
Схема рекурсивной оптимизации основана на известном [11] соотношении
min{ϕ ( y) : y ∈ D} = min min ... min ϕ ( y) ,
a1≤y1≤b1 a2≤y2≤b2 aN ≤yN ≤bN
(
Введем в рассмотрение множество функций
ϕ N ( y1,..., y N ) = ϕ ( y1,..., yN ) ,
ϕ i ( y1,..., yi ) = min ϕ i+1( y1,..., yi , yi+1) ,1 ≤ i ≤ N − 1 . (
ai+1≤yi+1≤bi+1
Тогда, в соответствии с соотношением (
ϕ1( y1*) = min{ϕ1( y1) : y1 ∈[a1,b1]}.
Однако при этом каждое вычисление значения одномерной функции ϕ1( y1) в некоторой фиксированной точке предполагает решение одномерной задачи минимизации ϕ2( y1, y2*) = min{ϕ2( y1, y2) : y2 ∈[a2,b2 ]},и так далее до вычисления ϕ N согласно (10). 3.2 Блочная рекурсивная схема редукции размерности
Для изложенной выше рекурсивной схемы предложено обобщение (блочная рекурсивная схема), которое комбинирует использование разверток и рекурсивной схемы с целью эффективного распараллеливания вычислений.
Рассмотрим вектор y как вектор блочных переменных
y = ( y1, y2,..., yN ) = (u1,u2,...,uM ) , где i-я блочная переменная ui представляет собой вектор размерности Ni из последовательно взятых компонент вектора y, т.е.
u1 = ( y1, y2,..., yN1 ) , u2 = ( yN1+1, yN1+2,..., yN1+N2 ) ,…, uM = ( yN−NM +1, yN−NM +2,..., yN ) , причем N1 + N2 + ... + NM = N .
С использованием новых переменных основное соотношение многошаговой схемы (
Формулы, определяющие способ решения задачи на основе соотношений (13) в целом
совпадают с рекурсивной схемой (10)−(
При этом принципиальным отличием от исходной схемы является тот факт, что в блочной схеме вложенные подзадачи ϕ i (u1,..., ui ) = min ϕ i+1(u1,..., ui , ui+1) ,1 ≤ i ≤ M − 1 , (14)
ui+1∈Di+1 являются многомерными, и для их решения может быть применен способ редукции размерности на основе кривых Пеано.
Число векторов и количество компонент в каждом векторе являются параметрами блочной
многошаговой схемы и могут быть использованы для формирования подзадач с нужными
свойствами. Например, если M = N , т.е. ui = yi ,1 ≤ i ≤ N , то блочная схема идентична
исходной; каждая из вложенных подзадач является одномерной. А если M = 1 , т.е. u = u1 = y , то
решение задачи эквивалентно ее решению с использованием единственной развертки,
отображающей [0,1] в D; вложенные подзадачи отсутствуют.
(10)
(
В2012 году компания Intel представила первый сопроцессор с архитектурой Intel MIC (Intel® Many Integrated Core Architecture). Архитектура MIC позволяет использовать большое количество вычислительных ядер архитектуры x86 в одном процессоре. В результате для параллельного программирования могут быть использованы стандартные технологии, такие как OpenMP и MPI.
Архитектура Intel Xeon Phi поддерживает несколько режимов использования сопроцессора, которые можно комбинировать для достижения максимальной производительности в зависимости от характеристик решаемой задачи.
В режиме Offload процессы MPI выполняются только на CPU, а на сопроцессоре происходит запуск отдельных функций, аналогично использованию графических ускорителей.
В режиме MPI базовая система и каждый сопроцессор Intel Xeon Phi рассматриваются как отдельные равноправные узлы, и процессы MPI могут выполняться на центральных процессорах и сопроцессорах Xeon Phi в произвольных сочетаниях. 4.1 Режим Offload
Вначале рассмотрим ситуацию, когда проведение одного испытания является трудоемкой операцией. В этом случае ускоритель Xeon Phi может быть использован в режиме Offload для параллельного проведения сразу многих испытаний на одной итерации метода. Пересылки данных от CPU к Xeon Phi будут минимальны – требуется лишь передать на сопроцессор координаты точек испытаний, и получить обратно значения функции в этих точках. Функции, определяющие обработку результатов испытаний в соответствии с алгоритмом и требующие работы с большим объемом накопленной поисковой информацией, могут быть эффективно реализованы на CPU. Общая схема организации вычислений с использованием Xeon Phi будет следующий.
На CPU выполняются шаги 1 – 4 параллельного алгоритма глобального поиска из п. 2. При этом на каждой итерации происходит накопление координат точек испытания в буфере, и этот буфер передается на сопроцессор. На Xeon Phi выполняется параллельное вычисление значений функции в этих точках (шаг 5 алгоритма). Используется OpenMP распараллеливание цикла, на каждой итерации которого происходит вычисление значений функции. По завершению испытаний происходит передача вычисленных значений функции на CPU. 4.2 Режим MPI
В случае, если проведение одного поискового испытания является относительно простой операцией, параллельное проведение многих итераций в режиме Offload не дает большого ускорения (сказывается влияние накладных расходов на передачу данных). Однако здесь можно увеличить вычислительную нагрузку на Xeon Phi, если применить блочную схему редукции размерности из п. 3.2, а для решения возникающих подзадач использовать сопроцессор в режиме MPI.
Для организации параллельных вычислений будем использовать небольшое (
Общая схема организации вычислений с использованием нескольких узлов кластера и нескольких сопроцессоров состоит в следующем. Процессы параллельной программы образуют дерево, соответствующее уровням вложенных подзадач, при этом вложенные подзадачи ϕ i (u1,..., ui ) = min ϕ i+1(u1,..., ui , ui+1)
ui+1∈Di+1 при i=1,…,M–2 решаются только с использованием CPU. Непосредственно в данных подзадачах вычислений значений оптимизируемой функции не происходит: вычисление значения функции ϕi (u1,...,ui ) − это решение задачи минимизации следующего уровня. Каждая подзадача решается в отдельном процессе; обмен данными осуществляется лишь между процессамипредками и процессами-потомками. Подзадача последнего (M–1)-го уровня процессов будет определяться как 1+π1 +π1π 2 .
Отметим, что синхронный параллельный алгоритм глобальной оптимизации, описанный в п.3, в сочетании с блочной схемой редукции размерности обладает существенным недостатком, связанным с возможными простоями всех процессов, кроме корневого. Простои могут возникать в том случае, если часть потомков некоторого процесса закончили решение своих подзадач и отправили данные родителю раньше остальных. Для преодоления данной проблемы был применен предложенный в [12] асинхронный вариант параллельного алгоритма. Конкретные детали реализация блочной схемы в сочетании с асинхронным алгоритмом описаны в [12]. 5. Результаты численных экспериментов
Вычислительные эксперименты проводились на узле суперкомпьютера «Лобачевский», установленного в ННГУ им. Н.И. Лобачевского, на сегменте под управлением операционной системы семейства Windows (HPC Server 2008). Узел располагает двумя процессорами Intel Sandy Bridge E5-2660 2.2 GHz, 64 Gb RAM, и двумя сопроцессорами Intel Xeon Phi 5110P.
В работе [6] описан GKLS-генератор, позволяющий порождать задачи многоэкстремальной оптимизации с заранее известными свойствами: количеством локальных минимумов, размерами их областей притяжения, точкой глобального минимума, значением функции в ней и т.п. Тестовые задачи, порождаемые данным генератором, характеризуются малым временем вычисления значений целевой функции. Поэтому с целью имитации вычислительной трудоемкости, присущей прикладным задачам оптимизации, расчет критерия был усложнен дополнительными вычислениями, не меняющими вид функции и расположение ее минимумов. Рис. 1. Решение двумерной задачи Для примера на рис. 1 приведены линии уровня двумерной функции, порождаемой GKLSгенератором. Темные точки соответствуют 464 испытаниям, потребовавшимся для решения данной задачи с точностью 0.01 по координате; при этом число итераций составило 58 (на каждой итерации проводилось 8 испытаний). Линии уровня наглядно демонстрируют многоэкстремальность задачи, а расположение точек – неравномерное покрытие, сгущающееся только в районе глобального минимума.
Ниже приведены результаты численного сравнения трех последовательных алгоритмов – DIRECT [9], DIRECTl [10] и алгоритма глобального поиска (АГП) (результаты работы первых двух алгоритмов приводятся по работе [6]). Численное сравнение проводилось на классах функций Simple и Hard размерности 4 и 5 из [6]. Глобальный минимум y* считался найденным, если алгоритм генерировал точку испытания yk в δ-окрестности глобального минимума, т.е.
y k − y* ≤ δ . При этом размер окрестности выбирался (в соответствии с [6]) как δ = b − a N ∆ , здесь N – размерность решаемой задачи, a и b – границы области поиска D, параметр ∆=10−6 при N=4 и ∆=10−7 при N=5. При использовании метода АГП для класса Simple выбирался параметр r=4.5, для класса Hard − r=5.6; параметр построения кривой Пеано был фиксированный m=10. Максимально допустимое число итераций составляло Kmax = 1 000 000.
В таблице 1 отражено среднее число итераций kav, которые выполнил метод при решении серии задач из данных классов. Символ «>» отражает ситуацию, когда не все задачи класса были решены каким-либо методом. Это означает, что алгоритм был остановлен по причине достижения максимально допустимого числа итераций Kmax. В этом случае значение Kmax=1 000 000 использовалось при вычислении среднего значения числа итераций kav, что соответствует нижней оценке этого среднего значения. Количество нерешенных задач указано в скобках.
Таблица 1. Среднее число итераций N 4 5 Класс функции
DIRECT
>47282 (4)
>95708 (7)
>16057 (
DIRECTl 18983 68754 16758 >269064 (4) АГП 11953 25263 15920 >148342 (4) Как видно из таблицы 1, АГП превосходит методы DIRECT и DIRECTl на всех классах задач по среднему числу итераций. При этом в классе 5-Hard каждый из методов решил не все задачи: DIRECT не решил 16 задач, DIRECTl и АГП – по 4 задачи. 5.1 Задачи с трудоемким критерием
Вначале приведем результаты экспериментов, которые соответствуют задачам с большим временем проведения одного поискового испытания. Оценим ускорение параллельного алгоритм глобального поиска (ПАГП), реализованного на CPU, в зависимости от использованного количества ядер p. В таблицах 2 и 3 приведено ускорение по времени S(p) и по итерациям s(p). Ускорение вычислено относительно однопоточного запуска.
p
1.15 2.82 3.47
Hard 1.32 2.59 5.34 Таблица 3. Ускорение по итерациям s(p) на CPU p Результаты экспериментов показывают значительное ускорение ПАГП при использовании CPU. Лучшие результаты получены при использовании всех ядер процессора.
Далее изложим результаты экспериментов, проведенных с использованием Intel Xeon Phi. Вначале рассмотрим эксперименты на одном Xeon Phi в режиме Offload. Варьировалось количество потоков на сопроцессоре, все остальные параметры совпадают с предыдущими запусками. Приведено (таблицы 4 и 5) ускорение относительно восьмипоточного запуска на центральном процессоре.
Таблица 4. Ускорение по времени S(p) на Phi Таблица 5. Ускорение по итерациям s(p) на Phi Результаты экспериментов показывают, что только на классе Simple при N=4 реализация на Xeon Phi медленнее, чем CPU реализация, на классе Hard при N=4 и на классе Simple при N=5 версии для Xeon Phi и CPU демонстрирует примерно равную эффективность, а на классе Hard при N=5 версия для Xeon Phi существенно превосходит CPU-реализацию. Наибольшее преимущество реализация для Xeon Phi показывает при 120 потоках для четырехмерных задач и при 240 потоках для пятимерных. При этом на всех классах наблюдается значительное ускорение по итерациям, а также практически линейная масштабируемость от числа потоков.
Далее рассмотрим эксперименты с запуском сопроцессора в режиме MPI. Число уровней вложенности разбиваемой задачи равно двум. Число процессов, запущенных на сопроцессоре, – 30 (таблицы 6 и 7) и 60(таблицы 8 и 9), что соответствует ситуации, когда у корня дерева имеется соответственно 30 и 60 потомков. Процессы, работающие на сопроцессоре, использовали OpenMP для параллельного вычисления значений функции. Ускорение приведено относительно восьмипоточного запуска на CPU.
Таблица 6. Ускорение по времени S(p) на Phi, 30 MPI процессов Таблица 7. Ускорение по итерациям s(p) на Phi, 30 MPI процессов p
5,15 3,35 4,61
9,27 8,24 8,28 N=5 N=5
14,65 10,49 14,36 Hard 27,38 23,85 22,12 Таблица 8. Ускорение по времени S(p) на Phi, 60 MPI процессов Таблица 9. Ускорение по итерациям s(p) на Phi, 60 MPI процессов 5.2 Задачи с легко вычисляемым критерием
Рассмотрим теперь решение задач, в которых поисковые испытания проводятся быстро. Вначале приведем результаты OpenMP версии. В таблице 10 приведено ускорение относительно однопоточного запуска.
Таблица 10. Ускорение по времени S(p) на CPU
N=4
N=5 Эксперименты показывают незначительное ускорение, а во многих случаях – и замедление при вычислениях только на центральном процессоре.
Далее рассмотрим ускорение при использовании Xeon Phi в режиме MPI. Результаты приведены относительно однопоточного запуска на CPU. В таблице 11 приведено ускорение по времени S(p) для запуска 30 процессов на Xeon Phi, в таблице 12 для 60 процессов. Таблица 11. Ускорение по времени S(p) на Phi, 30 MPI процессов p
0,56 0,74 0,45 0,76 0,22
0,37 0,39 0,54 0,40 N=5
0,57 0,66 1,77 0,66 0,55 Hard 1,92 1,23 2,63 1,53 Таблица 12. Ускорение по времени S(p) на Phi, 60 MPI процессов p Результаты экспериментов показывают ускорение на Xeon Phi, превосходящее ускорение на центральном процессоре. Лучшие результаты получены при использовании 30 процессов на Xeon Phi только для простых четырехмерных задач, все остальные задачи лучше всего решаются при использовании 60 процессов.
Теперь перейдем к эксперименту с шестимерными задачами простого класса. При решении шестимерной задачи с использованием одного процесса (вложенные подзадачи отсутствуют) и сохранением точности построения развертки m=10 требуется использовать числа с плавающей запятой расширенной точности, что ведет к значительному увеличению времени работы алгоритма. При этом в случае использования хотя бы одного уровня вложенных подзадач, задачи на каждом уровне будут решаться без подключения расширенной точности, что обеспечивает дополнительное преимущество блочной многошаговой схемы.
Рассмотрим результаты эксперимента с использованием Xeon Phi в режиме MPI, таблица 13. Число уровней вложенности разбиваемой задачи равно двум. Приведены результаты экспериментов для разбиений 3:3 (три переменных на первом уровне, три – на втором) и 4:2 (4 переменных на первом уровне, две – на втором). Число процессов, запущенных на сопроцессоре равно 30 и 60, что соответствует тому, что у корня дерева имеется соответственно 30 и 60 потомков. Процессы, работающие на Xeon Phi, использовали OpenMP для параллельного вычисления значений функции. Ускорение приведено относительно однопоточного запуска на CPU.
Таблица 13. Ускорение по времени S(p) на Phi p Результаты экспериментов показывают значительное ускорение запуска на сопроцессоре по сравнению с запуском только на центральном процессоре. Лучшее ускорение наблюдается при шестидесяти MPI процессах на Xeon Phi по 8 потоков на каждый процесс, время работы в 25 раз меньше, чем решение задачи только на центральном процессоре. 6. Заключение
В работе рассмотрен параллельный алгоритм глобального поиска, разработанный в рамках информационно статистического подхода. Предложен подход к распараллеливанию данного алгоритма, который может быть эффективен как в задачах с трудоемким критерием оптимизации, так и в задачах с относительно простым критерием. Проведены эксперименты с реализацией предложенного параллельного алгоритма на суперкомпьютере «Лобачевский» с использованием сопроцессоров Intel Xeon Phi. Результаты вычислительных экспериментов подтверждают высокую эффективность распараллеливания для Xeon Phi во всех классах рассмотренных задач. Литература 7. Žilinskas J. Branch and bound with simplicial partitions for global optimization. Mathematical
Modelling and Analysis. 2008. vol. 13(
Global Optimization. 2001. vol. 21(
Implementation of a parallel algorithm for searching the global extremum of a function on Intel Xeon Phi * K.A. Barkalov, I.G. Lebedev, V.V. Sovrasov, A.V. Sysoyev
Lobachevsky State University of Nizhni Novgorod The paper considers parallel algorithm for solving multiextremal optimization problems. The issues of implementation of the algorithm on state-of-the-art computing systems using Intel Xeon Phi coprocessor are examined. Two approaches for algorithm parallelization, which take into account information about laboriousness of the objective function computing, are considered. Speed up of the algorithm using Xeon Phi compared to the algorithm using CPU only is experimentally confirmed. Computational experiments are carried out on Lobachevsky supercomputer. 4. Strongin R.G., Sergeyev Ya.D. Global Optimization with Non-convex Constraints. Sequential and Parallel Algorithms. Kluwer Academic Publishers. 2000. 704 p. 7. Žilinskas J. Branch and bound with simplicial partitions for global optimization. Mathematical
Modelling and Analysis. 2008. vol. 13(
Gaviano M., Lera D., Kvasov D.E., Sergeyev Ya.D. Software for generation of classes of test functions with known local and global minima for global optimization // ACM Transactions on Mathematical Software. 2003. Vol. 29. P. 469−480. 10. Gablonsky J.M., Kelley C.T. A Locally-Biased Form of the DIRECT Algorithm // Journal of
Global Optimization. 2001. vol. 21(