Южно-Уральский государственный университет В работе описывается параллельный алгоритм решения нестационарных задач линейного программирования большой размерности, ориентированный на кластерные вычислительные системы. В основе алгоритма, получившего название «следящий», лежат фейеровские отображения. Алгоритм отслеживает изменения исходных данных и вносит корректировки в вычислительный процесс. При этом задача разбивается на большое количество подзадач, которые могут решаться независимо без обменов данными. Приводятся диаграммы деятельности UML, описывающие реализацию следящего алгоритма. Ключевые слова: нестационарная задача линейного программирования, фейеровские отображения, следящий алгоритм, диаграммы деятельности UML, массовый параллелизм, кластерные вычислительные системы.
которая будет использоваться в следящем алгоритме.
Обозначим ϕ (x ) = x − λ m max { ai, x − bi, 0}
∑ m i =1 a 2 i
⋅ ai , ϕs(x) = ϕ …ϕ(x). Под ϕ -проектированием (псевдопроектированием) точки x ∈ ℝn на многогранник M поs Под фейеровским процессом, порождаемым отображением ϕ при произвольном начальном приближении x0 ∈ ℝn , будем понимать последовательность {ϕs(x0)} s =0 занный фейеровский процесс сходится к точке, принадлежащей множеству M: +∞ {ϕs(x0 )}s =0 → x ∈ M . Будем кратко обозначать это следующим образом: lim ϕs(x0) = x . s→∞ +∞ . Известно, что ука2. Постановка задачи Пусть задана задача линейного программирования
max { c, x | Ax ≤ b,x ≥ 0} Определим фейеровское отображение ϕ : ℝn → ℝn следующим образом:
m ϕ (x ) = x − ∑ αiλi i =1 max { ai, x − bi, 0}
⋅ ai гая в формуле (2) λi = λ и αi = 1 m (i = 1,…,m) , получаем формулу
n−1 kα = ∑ αini .
i=0
(
s→∞ 3. Построение следящей области
Без ограничения общности мы можем считать, что все процессы происходят в положительной области координат. Пусть n – размерность пространства решений, K – количество ячеек в следящей области по одному измерению. Пусть P – количество MPI-процессов, используемых для распараллеливания вычислений. Будем предполагать, что всегда выполняется равенство:
K n = P , (5)
то есть, количество ячеек следящей области равно количеству MPI-процессов. Зададим в
пространстве целочисленных координат u0,…,un линейную нумерацию ячеек следящей области
следующим образом. Пусть ячейка α имеет целочисленные координаты (α0,…, αn ). Тогда ее
номер kα вычисляется по формуле:
u1
6
3
0
7
4
1
8
5
2
На рис. 1 приведен пример такой линейной нумерации при n = 3 . Например, ячейка с
номером 19 на рис. 1 имеет целочисленные координаты (
n α2 = kα − α0 − α1n mod n ;
n2 . . . . . .
α0 = kα − (kα ÷ n) ⋅ n 1. α1 = kα − (kα − (kα ÷ n) ⋅ n) mod n
n (kα ÷ n) ⋅ n = mod n
n
= (kα ÷ n )mod n.
α1 = kα ÷ n − ((kα ÷ n) ÷ n) ⋅ n .
(7)
(8)
(9)
(11)
(
i Пусть g = (g0,…, gn−1) – нулевая вершина куба следящей области. Пусть y = (y0,…,yn−1) – нулевая вершина произвольной ячейки α . Выразим координаты точки y 4. Пересечение многогранника M с ячейкой
Пусть y – нулевая вершина ячейки . Тогда область внутри ячейки задается системой из 2n неравенств: (включая границы) −x0 x0 Эта же система в матричной форме: −x1 x1 ⋯ ⋯ ≤ ≤ ⋯ −xn−1 ≤ ≤ ≤ ⋯ ≤ xn−1 −y0 −y1 ⋯ −yn−1 y0 + s y1 + s
⋯ yn−1 + s где (для n = 3 )
Aαx ≤ bα , (14) (15) (19) (20) Положим
−01 −01 00 −−yy10 Aα = 0 0 −01 , bα = −y2 .
1 0 y0 + s 00 10 10 yy21 ++ ss A′ = AAα , b′ = .
b bα A′x ≤ b′ , (21) (22) (23) Тогда пересечение многогранника M с ячейкой форме задается системой неравенств в матричной где A′ – расширенная матрица размера (m + 2n) × n , b′ – расширенный столбец свободных членов. Расширенный столбец b′ в соответствии с формулой (22) имеет инвариантную часть b, не зависящую от координат нулевой вершины ячейки , и вариативную часть bα , зависящую от координат нулевой вершины ячейки . Элементы расширенной матрицы A′ не зависят от координат нулевой вершины ячейки . 5. Реализация следящего алгоритма
В данном разделе описывается полная реализация следящего алгоритма в виде диаграмм деятельности UML. 5.1 Схема головной подпрограммы
Общая схема головной подпрограммы следящего алгоритма приведена на рис. 2. На шаге 1
выполняется подпрограмма init (см. раздел 5.2), выполняющая инициализацию переменных.
Затем в цикле until с меткой 2 выполняется корректировка следящей области в соответствии с
описанием идеи алгоритма в [11]. Одна итерация соответствует одной корректировке. Головная
подпрограмма следящего алгоритма оформляется в виде независимого процесса, который
выполняется до тех пор, пока переменная stop не примет значение true (истина). Начальную
установку переменной stop в значение false (ложь) осуществляет головной процесс,
соответствующий основной программе. Он же присваивает переменной stop значение true , когда
вычислительный процесс нужно остановить. В качестве текущего приближения решения задачи
(
В теле цикла until выполняются следующие действия. На шаге 3 организуется K параллельных потоков управления (нитей), которые независимо друг от друга вычисляют псевдопроекции из целевой точки z на пересечение i-той ячейки с многогранником M (i = 0,…, P − 1). Напомним, что P равно количеству MPI-процессов, и в соответствии с формулой (5) равно количеству ячеек в кубической следящей области. Схема подпрограммы вычисления псевдопроекции детально описана в разделе 5.3.
В цикле for с меткой 5 для полученных на шаге 3 точек псевдопроекций x0,…, xP−1 вычисляется номер kα ячейки, на которой достигается максимум C целевой функции. Для корректной работы цикла 5 переменным kα и C на шаге 4 присваиваются начальные значения MinInt и MinFloat соответственно. Значение MinInt соответствует минимальному машинному значению целого типа, MinFloat – минимальному машинному значению вещественного типа. Подпрограмма, вычисляющая точку xk = (xk0,…, xk,n−1) псевдопроекции точки z на пересечение многогранника M с ячейкой с номером kα , присваивает координате xk0 значение −1 в том случае, когда точка xk псевдопроекции не принадлежит многограннику M . Эта ситуация возникает в случае, когда пересечение многогранника M с ячейкой с номером kα является пустым.
2
1 0..n-1 z[0..n-1], ⋅ ⋅ ⋅ 0..n-1 z[0..n-1],P 3 5 k =MinInt; C=MinFloat 4 σ = <c[0..n-1];xk[0..n-1]> σ C
C = σ; k =k Рис. 2. Головная подпрограмма следящего алгоритма.
8 0..n-1 k
k =MinInt xk[0]=-1 ‖q-q'‖<¼r z[0..n-1]
0..n-1 13 s= 0..n-1 k
τ=k ÷ n [0] = [0] + (k − τ∗n)∗s [i]= [i]+(τ mod n)∗s
τ=τ ÷ n Рис. 3. Схема подпрограммы zero вычисления нулевой
вершины ячейки с порядковым номером k . Если же xk принадлежит многограннику, то в силу предположения о том, что все процессы находятся в положительной области координат (см. раздел 3), значение xk0 не может быть отрицательным. Указанное условие проверяется на шаге 6. Случаи xk0 = −1 из рассмотрения исключаются. Если при выполнении цикла 5 оказывается, что ни одна из ячеек следящей области не имеет непустого пересечения с многогранником M , то в переменной kα сохраняется значение MinInt . Этот факт проверяется на шаге 7. В этом случае шаг сетки s , длина r ребра следящей области и координаты целевой точки z увеличиваются в w раз, где w – положительная константа, являющаяся параметром алгоритма (шаг 13 на рис. 2).
Если на шаге 7 выясняется, что kα ≠ MinInt , значит найдена ячейка с номером kα , имеющая непустое пересечение с многогранником, на которой достигается максимум целевой функции. В этом случае на шаге 8 вычисляется вектор q′ , представляющий нулевую вершину новой центральной ячейки следящей области. Схема подпрограммы вычисления нулевой вершины ячейки с порядковым номером k приведена на рис. 3. Вычисления осуществляются с использованием формул (11), (15) и (16).
На шаге 9 (рис. 2) анализируется, насколько новая центральная ячейка далеко отстоит от 3 предыдущей. Если расстояние между новой и старой центральными ячейками превышает r 4 (где r – длина ребра кубической следящей области), то длина ребра кубической следящей области r , шаг сетки s и координаты целевой точки z на шаге 10 увеличиваются в 1.5 раза. Если 1 расстояние между новой и старой центральными ячейками меньше r , то длина ребра кубиче4 ской следящей области r , шаг сетки s и координаты целевой точки z на шаге 11 уменьшаются в 2 раза. Величина 1 2 , используемая в шагах 10 и 11, в общем случае является параметром горитма.
На шаге 12 следящая область сдвигается по вектору (q′ − q) , и в качестве текущей нулевой вершины q центральной ячейки следящей области берется точка q′ . input A[0..n-1,0..m-1], b[0..m-1], c[0..n-1], G[0..n-1] 3 input n, m, R, p, K, u, L, T
1 n mod u = 0 2 5 =rank
Kn=P [0..m-1]= [0..m-1] Aα[0..2n-1,0..n-1]=0 7
Aα[i,i]=-1 Aα[i+n,i]=1 6 4
8 [0..m-1,0..n-1]=A[0..m-1,0..n-1]
9 [m..m+2n-1,0..n-1]=Aα[0..2n-1,0..n-1] 10 r=R; s=r/K
11 g[0..n-1]=G[0..n-1]
12 z[0..n-1]=c[0..n-1]∗T
13 0..n-1 [0..n-1]+s∗γ[0..n-1]
14 γ[0..n-1]=(6K/27,...,6K/27)
15 Рис. 4. Схема подпрограммы инициализации переменных init. 5.2 Схема подпрограммы инициализации переменных init
Подпрограмма инициализации переменных init выполняет ввод исходных данных и инициализацию переменных. Схема подпрограммы init приведена на рис. 4. На шаге 1 вводятся значения переменных: n – размерность пространства решений; m – число неравенств в системе ограничений; R – начальное значение длины ребра следящей области, обеспечивающее покрытие многогранника M; p – количество итераций при построении псевдопроекции, выполняемое между обновлениями входных данных (этот параметр используется в подпрограмме dataChange обновления исходных данных); K – количество ячеек в следящей области по одному измерению; u – размерность подвектора; L – число независимых фейеровских итераций на подвекторах в подпрограмме вычисления псевдопроекции (см. рис. 5); T – масштабирующий коэффициент для вычисления координат целевой точки z. На шаге 2 проверяется условие (assert) n mod u = 0 , означающее, что размерность пространства решений n кратно размерности подвектора u . На шаге 3 осуществляется ввод исходных данных задачи линейного программирования: A – матрица коэффициентов неравенств; b – столбец свободных членов; c – вектор коэффициентов целевой функции. Также здесь вводится вектор G, содержащий начальные координаты нулевой вершины кубической следящей области. Шаг 4 присваивает переменной P значение, равное количеству доступных MPI-процессов и вычисляемое с помощью системной функции rank . На шаге 5 проверяется условие (assert) K n = P , означающее, что общее количество ячеек следящей области равно количеству MPI-процессов. На шаге 6 формируется инвариантная часть расширенного столбца b′ свободных членов, определяемого по формуле (22). На шаге 7 происходит инициализация матрицы Aα путем присваивания всем ее элементам нулевых значений. На шаге 8 определяются ненулевые элементы матрицы Aα в соответствии с системой неравенств (19). На шагах 9 и 10 строится расширенная матрица A′ , определяемая формулой (22). На шаге 11 в качестве начального значения длины r ребра следящей области определяется значение R, обеспечивающее покрытие многогранника M, и вычисляется значение s длины ребра ячейки. На шаге 12 в качестве начальной нулевой вершины g кубической следящей области берется точка G, определяющая такое положение следящей области, при котором она полностью покрывает многогранник M. На шаге 13 вычисляются координаты целевой точки z по формуле z = Tc . На шаге 14 вычисляется нулевая вершина q центральной ячейки следящей области по формуле (18). На последнем шаге вычисляется вектор γ начальных целочисленных координат центральной ячейки по формуле (17). 5.3 Схема подпрограммы вычисления псевдопроекции
На рис. 5 приведена схема подпрограммы вычисления псевдопроекции x = π(z,kα) из
целевой точки z на пересечение многогранника M с ячейкой следящей области, имеющей
порядковый номер kα , где kα вычисляется по формуле (
На шаге 3 организуется итерационный процесс вычисления псевдопроекции. Для обеспечения высокой масштабируемости процедуры вычисления псевдопроекции используется метод разбиения вектора x на h подвекторов размерности u , предложенный в работе [13]. Мы здесь предполагаем, что n = h ⋅ u . На каждом v-том подвекторе делается L независимых итераций вида (xvu,…, x(v+1)u−1) := (xvu,…, x(v+1)u−1) − − λ m max { (ai,vu,…,ai,(v+1)u−1),(xvu,…, x(v+1)u−1) − bi, 0}
∑ m i =1 a 2 i
⋅ (ai,vu,…,ai,(v+1)u−1). z[0..n-1],k
2 1 x[0..n-1]=z[0..n-1]; [0..n-1] = zero(k ); [0..n-1];
[0..n-1]+s = Рис. 5. Схема подпрограммы π вычисления псевдопроекции. q[i]-ε<x[i]<q[i]+s+ε 0..n-1
k Рис. 6. Схема подпрограммы in, проверяющей принадлежность точки x ячейке с номером k . Указанная формула получается из формулы (3) путем ограничения действия фейеровского отображения ϕ на соответствующий подвектор. Подпрограмма dataChange вносит изменения в исходные данные с периодом в t секунд (t – положительное число, которое может принимать значения меньше единицы).
Итерационный процесс заканчивается, когда расстояние между двумя последними приближениями x и x ′ будет меньше ε . На четвертом шаге подпрограмма in (см. рис. 6) проверяет принадлежность найденной точки псевдопроекции x ячейке с номером kα . Если x не принадлежит ячейке с номером kα , то x[0] присваивается значение (−1) .
Переменная s в подпрограмме in на рис. 6 задает длину ребра ячейки. Ее значение определяется на шаге 11 подпрограммы init (см. рис. 4.). Константа ε задает малое положительное число, позволяющее корректно обрабатывать приближенные значения. 6. Заключение
В работе описана параллельная реализация следящего алгоритма для решения нестационарных задач большой размерности на кластерных вычислительных системах. Данный алгоритм использует подход, основанный на применении фейеровских отображений для построения псевдопроекции на многогранник. Приведена формальная постановка задачи, и дано описание фейеровского процесса. С помощью математических формул определены следящая область и пересечение многогранника, задаваемого системой ограничений, с произвольной ячейкой следящей области. В плане дальнейших исследований – реализация параллельного алгоритма на языке C++ с использованием технологий параллельного программирования MPI и OpenMP, и проведение вычислительных экспериментов на искусственных и реальных данных. Литература 5. Дышаев М.М., Соколинская И.М. Представление торговых сигналов на основе адаптивной скользящей средней Кауфмана в виде системы линейных неравенств // Вестник ЮжноУральского государственного университета. Серия: Вычислительная математика и информатика. 2013. Т. 2, № 4. С. 103-108. 7. Sokolinskaya I.M., Sokolinskii L.B. Parallel algorithm for solving linear programming problem under conditions of incomplete data // Automation and Remote Control. 2010. Vol. 71, No. 7. P. 1452-1460. 11. Sokolinskaya I., Sokolinsky L. Solving unstable linear programming problems of high dimension on cluster computing systems // 1st Russian Conference on Supercomputing Days 2015, RuSCDays 2015; Moscow; Russian Federation; 28 September 2015 through 29 September 2015.
CEUR Workshop Proceedings. V. 1482, CEUR-WS.org 2015. P. 420-427. Implementation of Parallel Pursuit Algorithm for Solving Unstable Linear Programming Problems
I.M. Sokolinskaya, L.B. Sokolinsky
South Ural State University The paper describes an implementation of the parallel pursuit algorithm for solving unstable linear programming problems of high dimension on cluster computing systems. This algorithm uses Fejer,s mappings for building pseudo-projection on polyhedron. The algorithm tracks changes in input data and corrects the calculation process. This task is divided into set of independent subtasks, which can be processed in parallel. The UML activity diagrams describing the algorithm implementation are presented.
Keywords: unstable linear programming problem, Fejer’s mappings, pursuit algorithm, UML activity diagrams, massive parallelism, cluster computing systems. 7. Sokolinskaya I.M., Sokolinskii L.B. Parallel algorithm for solving linear programming problem under conditions of incomplete data // Automation and Remote Control. 2010. Vol. 71, No. 7. P. 1452-1460. 8. Rechkalov T.V., Zymbler M.L. Accelerating Medoids-based Clustering with the Intel Many Integrated Core Architecture // Proceedings of the 9th International Conference on Application of Information and Communication Technologies (AICT'2015), October 14–16, 2015, Rostov-on-Don, Russia. IEEE, 2015. P. 413–417. 9. Zymbler M.L. Best-match Time Series Subsequence Search on the Intel Many Integrated Core Architecture // Proceedings of the 19th East-European Conference on Advances in Databases and Information Systems, ADBIS 2015 (Poitiers, France, September 8–11, 2015). Lecture Notes in Computer Science. Vol. 9282. Springer, 2015. P. 275–286.