Представлены результаты исследования параллельного информационностатистического алгоритма глобальной оптимизации, разработанного в ННГУ им. Н.И. Лобачевского. Данный метод комбинируется с блочной схемой редукции размерности, что позволяет решать многомерные задачи путем их сведения к параллельному решению информационно-независимых подзадач меньшей размерности. Новым элементом, реализованным в рамках проводимого исследования, является использование нескольких графических ускорителей, задействованных на разных вычислительных узлах. Приводятся результаты решения тестовых задач на суперкомпьютере «Лобачевский» с использованием десятков тысяч GPU ядер.
Решение задач глобальной оптимизации на гетерогенных кластерных системах* на основе конечного числа k вычислений значений оптимизируемой функции. Относительно класса рассматриваемых задач предполагается выполнение двух важных условий.
Во-первых, предполагается, что оптимизируемая функция (y) может быть задана не
аналитически, а некоторым алгоритмом вычисления ее значений в точках области D; при этом
испытание (вычисление одного значения) является вычислительно-трудоемкой операцией.
Во-вторых, будем предполагать, что (y) удовлетворяет условию Липшица
(
Рассмотрим теперь основные способы распараллеливания вычислений, которые могут быть применены при решении задач глобальной оптимизации.
Во-первых, можно организовать разделение области решения между процессорами и параллельно решать подзадачи в этих подобластях. Однако такой подход обладает низкой эффективностью, поскольку при разделении области поиска только небольшая часть процессоров (в худшем случае – только один из них) будет решать задачу в подобласти с искомым глобальным минимумом; остальные процессоры будут работать в подобластях, в которых отсутствует решение исходной задачи.
Во-вторых, можно распараллеливать вычисление целевой функции, описывающей оптимизируемый объект. Данный путь может давать ускорение, но является специфичным для каждой конкретной решаемой задачи.
В-третьих, можно распараллелить реализацию вычислительных правил алгоритма, обеспечивающих выбор точки проведения очередного испытания. В этом случае способ распараллеливания будет зависеть от конкретного класса алгоритмов и, кроме того, часто эти правила достаточно просты и распараллеливать их нецелесообразно (накладные расходы на организацию параллелизма могут свести к нулю возможное ускорение).
Наконец, можно изменить схему алгоритма с целью параллельного выполнения нескольких испытаний (именно этот подход и будет рассматриваться в данной работе). Он является наиболее перспективным, т.к. характеризуется эффективностью (распараллеливается именно та часть вычислительного процесса, в котором выполняется основной объем вычислений) и общностью (применим для широкого класса характеристических алгоритмов многоэкстремальной оптимизации). Одновременно данный подход позволяет эффективно задействовать гетерогенные вычислительные ресурсы современных суперкомпьютеров (ядра на центральном процессоре, графические ускорители, математические сопроцессоры).
Данная статья продолжает серию исследований, начальные результаты которых были отражены в [1, 2]. 2. Базовый параллельный алгоритм глобального поиска
Для снижения сложности алгоритмов глобальной оптимизации, формирующих неравномерное покрытие области поиска, широко используются различные схемы редукции размерности, которые позволяют свести решение многомерных оптимизационных задач к семейству задач одномерной оптимизации. Поэтому в качестве базовой задачи мы будет рассматривать одномерную задачу многоэкстремальной оптимизации
* (x*) min{ (x) : x [0,1]}, в которой целевая функция (x) удовлетворяет условию Липшица. Дадим детальное описание параллельного алгоритма глобального поиска (ПАГП), применяемого к ее решению.
Пусть в нашем распоряжении имеется p 1 вычислительных элементов. Тогда на данной итерации можно провести одновременно p испытаний. Тогда общее число испытаний, выполненных после n параллельных итераций, составит k pn .
Предположим, что выполнено n 1 итераций метода (в качестве точек x1,...,x p первой итерации выбираются произвольные различные точки отрезка [0,1]). Тогда точки xk 1,...,xk p текущей (n+1)-ой итерации определяются по следующим правилам.
Правило 1. Перенумеровать точки множества
Xk {x1,...,xk} {0} {1} которое включает в себя граничные точки интервала [0,1], а также точки предшествующих испытаний, нижними индексами в порядке увеличения значений координаты, т.е.
0 x0 x1 ... xk1 1 Правило 2. Полагая zi (xi ),1 i k , вычислить величины r , 0, max | zi zi1 | , M 1ik i 1, 0, , где r 1 является заданным параметром метода (параметр надежности), а i xi xi1.
Правило 3. Для каждого интервала (xi1, xi ),1 i k 1, вычислить характеристику в соответствии с формулами
z z
R(
M M
2 R(i) i (zi zi1) 2 zi zi1 , 1 i k 1
M 2i M R(t1) R(t2) R(tk1) R(tk1) xk j xt j xt j 1 , t j 1 , t j k 1
2 xk j xt j xt j 1 zt j zt j 1 , 1 t j k 1 .
2 2M
(
Правило 5. Провести новые испытания в точках xk j , 1 j p , вычисленных по формулам
Алгоритм прекращает работу, если выполняется условие t j хотя бы для одного
номера t j , 1 j p ; здесь 0 есть заданная точность. В качестве оценки глобально-оптимального
решения задачи (
Одним из подходов к решению многомерных задач глобальной оптимизации является сведение их к одномерным и использование эффективных одномерных алгоритмов глобального поиска к редуцированной задаче. Данном разделе будут кратко изложены два известных способа редукции размерности, а также их обобщение. 3.1 Редукция размерности с использованием кривых Пеано
Первым из рассматриваемых способов редукции размерности является использование кривой Пеано y(x), однозначно отображающей отрезок вещественной оси [0,1] на n-мерный куб y RN : 21 yi 21,1 i N y(x) : 0 x 1, Вопросы численного построения отображений типа кривой Пеано и соответствующая теория подробно рассмотрены в [3]. Здесь же отметим, что численно построенная развертка является приближением к теоретической кривой Пеано с точностью порядка 2m , где m – параметр построения развертки.
Использование подобного рода отображений позволяет свести многомерную задачу (
( y*) ( y(x*)) min{ ( y(x)) : x [0,1]}.
Важным свойством является сохранение ограниченности относительных разностей
функции: если функция (y) в области D удовлетворяла условию Липшица (
H 4Ld N , d max{bi ai :1 i N}.
Соотношение (8) позволяет модифицировать приведенный в разделе 2 алгоритм решения
одномерных задач для решения многомерных задач, редуцированных к одномерным. Для этого
длины интервалов i , участвующие в правилах (
i (xi xi1)1 N , а вместо формулы (7) вводится выражение 1 | zt j zt j 1 |N xk j xt j xt j 1 sign (zt j zt j 1) 2r , 1 t j k 1 .
2 3.2 Рекурсивная схема редукции размерности Схема рекурсивной оптимизации основана на известном (см. [5]) соотношении min{ ( y) : y D} min min ... min ( y) ,
a1 y1b1 a2 y2 b2 aN yN bN
которое позволяет заменить решение многомерной задачи (
Введем в рассмотрение множество функций
N ( y1,..., yN ) ( y1,..., yN ) , i ( y1,..., yi )
min ai1 yi1bi1
i1( y1,..., yi , yi1) , 1 i N 1. 1( y1*) min{1( y1) : y1 [a1,b1]}. Однако при этом каждое вычисление значения одномерной функции 1( y1) в некоторой фиксированной точке предполагает решение одномерной задачи минимизации
2( y1, y2*) min{2( y1, y2) : y2 [a2,b2]}, и так далее до вычисления N согласно (10). 3.3 Блочная рекурсивная схема редукции размерности
Для изложенной выше рекурсивной схемы предложено обобщение (блочная рекурсивная схема), которое комбинирует использование разверток и рекурсивной схемы с целью эффективного распараллеливания вычислений.
Рассмотрим вектор y как вектор блочных переменных
y ( y1, y2,...,yN ) (u1,u2,...,uM ) , где i-я блочная переменная ui представляет собой вектор размерности Ni из последовательно взятых компонент вектора y, т.е. u1 ( y1, y2,...,yN1 ) , u2 ( yN1 1, yN1 2,..., yN1 N2 ) , … , u2 ( yN NM 1, yN NM 2,..., yN ) , причем N1 N2 ... NM N .
С использованием новых переменных основное соотношение многошаговой схемы (9) может быть переписано в виде min ( y) min min ... min ( y) , yD u1D1 u2D2 uM DM (13) где подобласти Di ,1 i M , являются проекциями исходной области поиска D на подпространства, соответствующие переменным ui ,1 i M .
Формулы, определяющие способ решения задачи (
При этом принципиальным отличием от исходной схемы является тот факт, что в блочной схеме вложенные подзадачи являются многомерными, и для их решения может быть применен способ редукции размерности на основе кривых Пеано.
Число векторов и количество компонент в каждом векторе являются параметрами блочной
многошаговой схемы и могут быть использованы для формирования подзадач с нужными
свойствами. Например, если M N , т.е. ui yi ,1 i N , то блочная схема идентична
исходной; каждая из вложенных подзадач является одномерной. А если M 1 , т.е. u u1 y , то
решение задачи эквивалентно ее решению с использованием единственной развертки,
отображающей [0,1] в D; вложенные подзадачи отсутствуют.
4. Организация параллельных вычислений на гетерогенном кластере
Для организации параллельных вычислений будем использовать небольшое (
M i 1 j ,
i1 j1 параллельно работающих процессоров/ядер.
Используя различные параметры вектора распараллеливания можно адаптировать алгоритм для работы на гетерогенной вычислительной системе. При этом операцией, которую можно эффективно реализовать на ускорителе, является параллельное вычисление сразу многих значений целевой функции. Пересылки данных от CPU к GPU будут минимальные: требуется лишь передать на GPU координаты точек испытаний, и получить обратно значения функции в этих точках. Обработка результатов испытаний в соответствии с алгоритмом, требующая работы с большим объемом накопленной поисковой информацией, может быть эффективно реализована на CPU.
Общая схема организации вычислений с использованием нескольких узлов кластера и нескольких GPU приведена на рис. 1; процессы параллельной программы будут образовывать дерево, соответствующее уровням вложенных подзадач. В соответствии с данной схемой вложенные подзадачи при i1,…,M–2 решаются только с использованием CPU. Непосредственно в данных подзадачах вычислений значений оптимизируемой функции не происходит: вычисление значения функции i (u1,...,ui ) это решение задачи минимизации следующего уровня. Каждая подзадача решается в отдельном процессе; обмен вычисленными значениями организован с помощью MPI.
Подзадача последнего (M–1)-го уровня M 1(u1,...,uM 1) min M (u1,...,uM )
uM DM отличается от всех предыдущих подзадач – в ней происходит вычисление значений оптимизируемой функции, т.к. M (u1,...,uM ) ( y1,..., yN ) . Данная подзадача также решается на CPU, но вычисление значений функции производится на ускорителе. При этом на CPU выполняются правила 1 – 4 параллельного алгоритма глобального поиска. Координаты p точек испытаний, вычисленные на шаге 4 алгоритма, накапливаются в промежуточном буфере, а затем передаются на графический процессор. На GPU происходит вычисление значений функции в этих точках, после чего результаты испытаний (снова через промежуточный буфер) передаются на CPU.
Самый простой вариант использования данной схемы будет соответствовать двухкомпонентному вектору распараллеливания ( 1, 2 ) . Здесь 1 1 будет соответствовать числу MPI-процессов, а 2 количеству используемый ядер на GPU; тем самым общее количество задействованных ядер (CPU и GPU) будет определяться как 1 1 1 2 .
Отметим также, что в случае невозможности эффективно реализовать процесс вычисления значения оптимизируемой функции на GPU на последнем уровне распараллеливания можно использовать ядра центрального процессора или ускорителя Xeon Phi в многопоточном режиме.
1(u1) … … …………………... 3(u1,u2,u3) … 3(u1,u2,u3) 3(u1,u2,u3) … 3(u1,u2,u3) M1(u1,..,uM1) M(u1,…,uM)
…………………...
Рис. 1. Схема организации параллельных вычислений на кластере 5. Результаты вычислительных экспериментов
Вычислительные эксперименты проводились на суперкомпьютере «Лобачевский» (операционная система – CentOS 6.4, система управления – SLURM). Один узел суперкомпьютера располагает 2-я процессорами Intel Sandy Bridge E5-2660 2.2 GHz, 64 Gb RAM, и 2-я GPU NVIDIA Kepler K20Х. Центральный процессор является 8-и ядерным (т.е. всего на узле доступно 16 ядер CPU), а на графическом процессоре доступны 14 потоковых мультипроцессоров (2688 CUDA-ядер). Использовался компилятор Intel C++ 14.0.2 и CUDA Toolkit 6.0.
В качестве тестовых задач были выбраны задачи, порождаемые с помощью GKLSгенератора [6]. Данный генератор позволяет получать задачи многоэкстремальной оптимизации с заранее известными свойствами: количеством локальных минимумов, размерами их областей притяжения, точкой глобального минимума, значением функции в ней и т.п. С целью имитации вычислительной трудоемкости, присущей прикладным задачам оптимизации, расчет целевой функции во всех проводимых экспериментах был усложнен дополнительными вычислениями, не меняющими вид функции и расположение ее минимумов (суммированием ряда из 20 тыс. элементов).
Для примера на рис. 2 приведены линии уровня двумерной функции, порождаемой GKLSгенератором. Темные точки соответствуют 464 испытаниям, потребовавшимся для решения данной задачи с точностью 0.01 по координате; при этом число итераций составило 58 (на каждой итерации проводилось 8 испытаний). Линии уровня наглядно демонстрируют многоэкстремальность задачи, а расположение точек – неравномерное покрытие, сгущающеюся только в районе глобального минимума. Рис. 2. Решение двумерной задачи Эксперимент с использованием одного графического ускорителя проведем для решения серии из 100 шестимерных задач. Глобальный минимум y* в тестовой задаче будем считать найденным, если алгоритм сгенерировал точку испытания yk в -окрестности глобального минимума, т.е. yk y* . Размер окрестности 0.01 b a , где a и b – границы области поиска D. Параметр надежности метода был выбран r 4.5; параметр построения развертки – m 10 . Максимально допустимое число параллельных итераций составляло Kmax = 10 000. Так как задача решалась с использованием единственного ускорителя, то разбиения размерностей не проводилось.
Среднее время решения задачи с использованием GPU составило 10.78 с., тогда как среднее время решения задачи с использованием всех 16 ядер CPU на узле – 53.8 с.; наблюдается почти пятикратное ускорение.
Более масштабный вычислительный эксперимент был проведен для задач размерности N = 8 и N = 10: было использовано 12 узлов кластера с 36 графическими ускорителями (по три ускорителя на узел), тем самым было задействовано 96 768 CUDA-ядер.
В соответствии с блочной рекурсивной схемой (13) было использовано два уровня подзадач с размерностями N1 N2 4 для восьмимерной задачи и N1 N2 5 для десятимерной. Вектор распараллеливания (15) был выбран как (36,2688) в соответствии с общим числом используемых GPU (компонента 1 ) и CUDA-ядер на одном ускорителе (компонента 2 ). Время решения восьмимерной задачи составило 405.6 с., ускорение по сравнению с CPU версией алгоритма составляет 5.9 раз. Время решения десятимерной задачи составило 2055.8 с. Десятимерная задача ввиду большой сложности вычислений на CPU не решалась, поэтому оценить ускорение не представляется возможным. 6. Заключение
Результаты проведенных экспериментов на серии тестовых задач разной размерности показывают, что предложенная блочная многошаговая схема редукции размерности в сочетании с параллельным алгоритмом глобального поиска эффективно реализуется на современных вычислительных системах. Данная схема обладает: высоким запасом параллелизма (было задействовано порядка 105 ядер); малой информационной зависимостью параллельно выполняемых вычислений (между параллельными процессами нужно было передавать лишь точки и результаты испытаний, полученные на текущей итерации метода).
Полученные результаты позволяют сделать предположение, что предложенная схема организации параллельных вычислений будет эффективна и при использовании большего (порядка 106) числа ядер. Дальнейшее развитие также будет заключаться в использовании локальной информации о поведении оптимизируемой функции [7], учете ограничений в задачах условной оптимизации [8], адаптации алгоритмов для решения многокритериальных задач [9]. Литература Solving the global optimization problems on heterogeneous cluster systems Konstantin Barkalov, Victor Gergel, Ilya Lebedev and Alexander Sysoyev Keywords: global optimization, parallel algorithms, heterogeneous computing The paper presents the results of investigation of parallel information-statistical algorithm of global optimization that developed in University of Nizhny Novgorod. This method is combined with multilevel scheme of dimension reduction. It allows us to solve multidimensional problems by reducing it to parallel solving a series of data-independent subtasks of lower dimension. New suggestion developed in the framework of current research is the method of using several graphics accelerators on different nodes of cluster. Results of solving a series of test problems on supercomputer Lobachevsky using tens of thousands GPU cores are given.