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

Вычислительные технологии

Том 9, № 1, 2004

ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ ЦЕЛОЧИСЛЕННОГО КВАДРАТИЧНОГО ПРОГРАММИРОВАНИЯ*

Г. И. Злвиняко, Е.А. Котельников Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: zabin@rav.sscc.ru

The parallel algorithm of integer and mixed-integer quadratic programming, based on the branch and bound method. The algorithm was realized in FORTRAN using the MPI system of parallel programming. The efficiency of the parallel and the sequential algorithms are compared for test problems.

Введение

Задача целочисленного и частично целочисленного квадратичного программирования (ЦКП) представляется в следующем виде:

минимизировать f (ж) = ^(х,Ох) + (с,х) (1)

при условиях

Ах = Ь, (2)

а < х < в, (3)

х^ — целые числа для ] £ 3. (4)

Здесь А — матрица размера т х п; О — симметричная матрица размера I х I (I < п); О > 0; с,х,а, в £ Яп; Ь € Кт; 3 — список целочисленных переменных.

Алгоритм решения задач ЦКП основывается на методе ветвей и границ с односторонним ветвлением, который представляет собой обобщение алгоритма целочисленного линейного программирования (ЦЛП). Решение исходной задачи сводится к решению последовательности оценочных задач квадратичного программирования.

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 01-07-90367).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

некоторое сочетание алгоритмов с односторонним и одновременным ветвлением [1]. Использование MPI существенно упрощает программирование такого сочетания алгоритмов. Во время исполнения процессы обмениваются короткими сообщениями на фоне выполнения большого числа арифметических операций, что обеспечивает высокую эффективность распараллеливания как на вычислительных системах с общей памятью, так и на кластерах.

1. Решение оценочных задач

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

Для решения оценочных задач используется метод приведенного градиента [2], в котором переменные х подразделяются на базисные хв, супербазисные х^ и небазисные хN. В матрице А, соответственно, выделяются В-базисная т х т, Б-супербазисная т х в и N-небазисная т х (п — (т + в)) матрицы. Предполагается, что В — невырожденная матрица, а переменные XN принимают свои граничные значения.

В текущем подпространстве супербазисных переменных минимизация проводится методом сопряженных градиентов. На к-м шаге метода последовательно определяются:

1) вектор двойственных переменных ук = (В-1)ТV/(хкв);

2) вектор приведенного градиента Ьк = V/(х|) — БТук;

3) вектор направления в подпространстве супербазисных переменных

кк , 1И1 2 „к-1.

— hk I 11 112

ps — -h + Tnjk—ûvïPs

Il hk ||2 l|hk-1 Il?'

4) вектор направления в подпространстве базисных переменных рВ = —В 1Бр|;

г\ ^ к \ (к ,ря) к к к

5) оптимальный шаг вдоль р А к = — , к ~ к, , где вектор р составлен из рВ, р£ и

(р , )

pN = 0. Полагается х к+1 = х к + Ар к, где А = шт{А к, Атах}, а Атах = argmax{А : аг < хк+

Л

Арк < вг}.

Если А к > Атах, то по методу приведенного градиента меняется состав супербазисных и других переменных, процесс по методу сопряженных градиентов начинается в новом подпространстве. В противном случае процесс по методу сопряженных градиентов продолжается до достижения с заданной точностью равенства нулю приведенного градиента.

Для уменьшения влияния ошибок вычислений на процесс по методу сопряженных градиентов осуществляется контроль точности решения систем уравнений ВТу к = V/в (хкв), ВрВ = Бр| и, если необходимо, их решения уточняются.

Против зацикливания метода сопряженных градиентов в плохо обусловленных задачах используется прием из [3], который состоит в проверке через определенное число итераций выполнения неравенства ||кк||2 < с(к)||кк — кк||2, где величина с(к) возрастает с увеличением к, кк — приведенный градиент с непосредственным вычислением V/(х|), а в кк используется рекуррентно вычисленный V/(х|). При выполнении неравенства процесс по методу сопряженных градиентов в данном супербазисном подпространстве завершается.

После завершения оптимизации методом сопряженных градиентов в текущем супербазисном подпространстве на основе двойственных оценок небазисных столбцов проверяется

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

2. Алгоритм целочисленного квадратичного программирования

Оценки границ целевой функции (1) с учетом (4) получаются из решений оценочных задач квадратичного программирования. Выбор переменной ветвления в ЦЛП часто осуществляют с использованием штрафов [1, 2]. В рассматриваемом алгоритме ЦКП также делается попытка, когда это возможно, для выбора переменной ветвления использовать штрафные оценки.

Вначале рассмотрим линейный случай. Значение базисной переменной Xj для ] £ 3 в оптимальном базисе очередной оценочной задачи представим в виде Xj = [х,] + Vj, где [xj] — целая часть х,. Штраф за увеличение х, на величину 1 — V, обозначим через Р+, а за уменьшение на величину Vj — через Р-.

На примере Р+ покажем алгоритм определения штрафных оценок. Номера небазисных переменных зададим в двух списках: в I- поместим номера переменных, фиксированных на нижней границе, а в 1+ — на верхней границе. Если переменная х, занимает в базисе р-ю позицию, то вычислим г = е^В-1, ер — р-й орт в Кт. Тогда

P+ = min < min |(1 — Vj)(—min |(1 — Vj)(— q

I q:q€/-, apq<0 ^ q:q£l+, apq>0 ^ ^ a;

■pq

где apq = (z,Aq), Aq — q-й небазисный столбец матрицы A; dq — его приведенная оценка. Если q-я небазисная переменная целочисленная, то величины (1 — Vj)(dq/ — apq) или (1 — Vj)(—dq/apq), которые меньше |dq |, заменяются на |dq |. Начальные значения Pj- и Pj+ полагаются равными

Штрафы для квадратичной функции f заменим штрафами для линейной функции g(x) = (Vf(ж), ж — ж*), где x* — решение очередной оценочной задачи квадратичного программирования. Поскольку исходная функция выпуклая, штрафные оценки для g(x) будут не больше штрафов для f (ж). В нелинейном случае при вычислении штрафов необходимо дополнительно учесть супербазисные переменные (если они присутствуют в оптимальном решении оценочной задачи).

Величины dq для супербазисных переменных равны нулю (dq — q-я компонента приведенного градиента), и штрафы при рассмотрении очередной j-й переменной не обнуляются, если apq = 0. Для любой j-й супербазисной переменной Pj+ = Pj- = 0.

Пусть Pmax = max max {Pj-, P+}. В задачах ЦЛП Pmax = 0 получается в оценочных

задачах, которые имеют не единственный оптимальный базис. В нелинейном случае возможность получить Pmax = 0 увеличивается за счет супербазисных переменных. Однако численные эксперименты показали, что использование штрафов в отдельных задачах ЦКП позволяет сократить время решения в несколько раз.

В случае Pmax = 0 для ветвления выбирается базисная или супербазисная переменная Xj, которая имеет значение, наиболее уклоняющееся от целочисленного. При этом только для определения переменной ветвления временно полагаются P+ = 1 — Vj и Pj- = Vj.

Схемы метода ветвей и границ с односторонним ветвлением позволяют использовать компактную форму списков оценочных задач. Пусть на некотором уровне к для ветвления выбрана переменная х , со штрафами Р- и Р+. Если Р- < Р+, то в качестве очередной

^ и и и и

оценочной задачи выбирается задача, соответствующая Р-, а в списке оценочных задач (во вспомогательном массиве к) необходимо запомнить информацию об альтернативной задаче. Для этого производятся присвоения: к(1, к) = к(2, к) = в,, где в, — верхняя граница переменной х, в г-й оценочной задаче на предшествующем уровне (к—1); к(3, к) = 1, если /г + Р+ > гг и к(3,к) = 0 в противном случае, здесь /г — оптимальное значение / в г-й оценочной задаче, а гг — текущее значение рекорда; к(4,к) = /г; к(5, к) = Р+. При Р+ < Р, выбирается оценочная задача, отвечающая штрафу Р+, а в массиве к запоминаются величины: к(1, к) = —к(2, к) = а,, к(3, к), равная 0 либо 1 в зависимости от выполнения неравенства /г + Р- > гг, к(4, к) = /г, к(5, к) = Р-.

В дальнейшем при рассмотрении параллельного алгоритма ветвь, которой соответствует к(3,к) = 1, будем называть помеченной и не помеченной при к(3, к) = 0.

Алгоритм.

Шаг 0. Положить г = 0, к = 0, значение рекорда г0 = а0 = а и в0 = в.

Шаг 1. Решить текущую задачу квадратичного программирования. Пусть хг — ее оптимальное решение и /г = / (хг).

А. Если /г > гг или система ограничений несовместна, то присвоить гг+1 = гг, г = г + 1 и перейти на шаг 3.

Б. Пусть /г < гг. Если вектор хг нецелочисленный, то перейти на шаг 2. Если решение хг целочисленное, то присвоить гг+1 = /г, г = г + 1. При /г = /0 перейти на шаг 5, иначе на шаг 3.

Шаг 2. Для тех базисных и супербазисных переменных х,, ] £ 3, у которых х, — нецелое, вычислить штрафы Р+ , Р,-. Среди них найти минимальный штраф Рт;п:

А. Если /г + Рт;п > гг, то перейти к шагу 3.

Б. При /г + Рт;п < гг осуществить ветвление по базисной или супербазисной переменной х,, которой соответствует максимальный штраф Р+ или Р,-. Из двух величин Р+, Р- выбрать меньшую. Пусть Р- < Р+, тогда выбрать очередную оценочную задачу, соответствующую Р - . Присвоить к = к + 1, записать в список задачу, отвечающую штрафу Р+. Изменить верхнюю границу переменной х,, положив ее равной [х, ]. Перейти на шаг 1.

(Если Р+ < Р-, то в списке к запоминается задача, отвечающая Р-, а в формируемой задаче нижняя граница переменной х, полагается равной [х,] + 1.)

Шаг 3. Если к = 0, то перейти к шагу 5.

Если к > 0 и к(3, к) = 1 или к(4, к) + к(5, к) > гг, то перейти к шагу 4.

Иначе, используя массив к, сформировать задачу, альтернативную образованной в шаге 2 Б. Присвоить к(3, к) = 1 и перейти на шаг 1.

Шаг 4. Установить двусторонние ограничения на переменную х,, используя значения к(1, к), к(2, к). Присвоить к = к — 1 и перейти на шаг 3.

Шаг 5. Остановиться.

Схемы ветвей и границ с односторонним ветвлением позволяют с небольшими затратами осуществить переход к очередной оценочной задаче. Рассмотренная реализация схемы во многом сходна с реализацией для ЦЛП в [4].

3. Параллельный алгоритм

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

На начальном этапе данные о задаче считываются с дисковой памяти нулевым процессором и пересылаются всем остальным. Далее на нулевом процессоре решается нулевая оценочная задача (исходная задача без учета условий целочисленности), а остальные процессоры находятся в ожидании.

Для загрузки любого процессора в его оперативную память пересылаются следующие данные: целочисленные массивы 1в,1и — соответственно списки базисных, супербазисных и небазисных переменных; часть массива Н (часть списка оценочных задач); х^ — массив значений супербазисных переменных. Значения элементов массивов 1в ,1я ,1и, х^ и какая часть массива Н пересылается, определяется уровнем к, с которого начинается решение первой после загрузки оценочной задачи. На основе списка 1в строится матрица, обратная к базисной, и начинается исполнение алгоритма ветвей и границ с уровня к.

Для упорядочения загрузки процессорных элементов на нулевом процессоре хранится массив пар (г, кг), где г — номер процессора, а кг — уровень непомеченной ветви на г-м процессоре. В списке пар для каждого процессора г сохраняется кг с наименьшим значением (что соответствует наиболее высокой непомеченной ветви на данном процессоре). Если некоторый г-й процессор простаивает, то это помечается присвоением кг = —1.

Рассмотрим случай, когда г-й процессор завершил выполнение алгоритма (достигнут нулевой уровень) или процессор изначально еще не загружен. Если номер процессорного элемента г > 0, то нулевому процессору г-й процессор посылает запрос на загрузку. После получения запроса на нулевом процессоре выбирается пара (], kj) с минимальным значением kj > 0 и на ]-й процессор отправляется сообщение о необходимости загрузки г-го процессора. В результате на ]-м процессоре ветвь уровня kj помечается и отправляются данные для загрузки на г-й процессор. Если не существует kj > 0, то г-й процесс находится в ожидании. При завершении исполнения алгоритма на нулевом процессоре аналогично определяется номер ]-го процессора для загрузки.

Если на некотором процессоре получено новое значение рекорда г, то это значение пересылается всем остальным процессорам. При получении значения рекорда на каждом из процессоров проверяется выполнение неравенства Н(4, к) — Н(5, к) > г, где к — уровень, с которого начиналось исполнение алгоритма на данном процессоре. Если неравенство выполняется на каком-либо из процессоров, то на этом процессоре завершается исполнение алгоритма ветвей и границ и формируется запрос на загрузку.

Полученное значение г используется для обновления информации о помеченных и непомеченных ветвях в массиве Н. Может оказаться, что на некотором из процессоров ] помечается не помеченная ветвь с минимальным значением kj. В этом случае на нулевой процессор посылается информация об изменении kj.

После решения очередной оценочной задачи на каждом из процессоров производится обращение к функции МР1_1РИ,ОВЕ [5] для проверки, имеются ли сообщения для данного процесса. Кроме того, обращение к этой функции производится после выполнения определенного числа итераций решения оценочной задачи. Частота обращений к функции за-

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

Задача решена, если на всех процессорных элементах достигнут нулевой уровень.

Передача исходных данных о задаче производится с помощью блокированной функции MPI "один — всем". В дальнейшем все процессы выполняются в асинхронном режиме и обмены осуществляются с помощью неблокированных функций.

4. Решение задач

Тестовые задачи взяты из [6]. В табл. 1 указаны входные параметры задач.

В табл. 2 приведены результаты расчетов в однопроцессорном режиме на системе МВС 1000М с процессорами альфа-21264, с тактовой частотой 830 МГц.

Таблица 1

Входные параметры задач

ь Название

задачи задачи m me n ni nb nz

1 ibell3a 106 — 122 60 31 304

2 ibienstl 576 128 505 28 28 2184

3 ield76 75 75 1898 1898 1898 19111

4 imas284 68 — 151 150 150 9631

5 imisc07 212 35 260 259 259 8619

6 imodOll 4480 4367 10957 97 96 22253

7 iqiu 1192 132 840 48 48 3432

8 iran13 195 26 338 169 169 676

9 iran8 296 40 512 256 256 1024

Примечание. т — общее число ограничений, среди которых те ограничений являются равенствами; п — число переменных, среди которых Пг являются целочисленными и щ переменных — булевыми; пг — число ненулей в матрице А.

Таблица 2

Результаты расчетов в однопроцессорном режиме на системе МВС 1000М

ь задачи it cp t Pi P2 Po

1 354065 35765 66.76 2680 9648 18

2 32317905 71883 45315. 26 2500 30294

3 42003240 91423 51143. 150 4385 10577

4 1054772 53052 380.8 409 16672 127

5 3848924 54826 3057. 202 6797 2086

6 5500000 5411 83258. 15 193 —

7 32000000 150246 85163. 368 16556 4166

8 18496522 1324478 4480. 8847 537387 6604

9 16239520 985215 5751. 4385 400440 15115

Примечание. И — суммарное число итераций в оценочных задачах; ср — количество оценочных задач; Ь — время решения, Р1 указывает, сколько раз в процессе решения задачи выполнилось условие fг + Рт1п > гг; Р2 — число случаев выполнения условия fг + Р, > гг, где Р, = Р+ или

Pj

Р- ; Р0 — указывает, сколько раз в процессе решения задачи было получено Ртах = 0.

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

f - f *

ции, через f * — оптимальное значение и определим относительную ошибку f

If *l

Значения величины 8f для задач 6 и 7 соответственно равны 0.119 и 0.002.

Характеристики процесса решения задач параллельным алгоритмом содержит табл. 3. В параллельном варианте расчетов использовались восемь процессоров. Во всех вариантах расчетов значения управляющих парметров выбирались одни и те же, кроме параметра, определяющего частоту обращения к функции MPI_IPROBE. Этот параметр для разных задач принимался в пределах от 10 до 30 в зависимости от числа строк в матрице ограничений.

Для задач 6 и 7, в которых в однопроцессорном режиме не удалось получить точное решение, были проведены дополнительно расчеты на 12 процессорах. В результате суммарное число итераций сократилось по сравнению с расчетами на восьми процессорах, для задачи 6 — на 4.1 млн итераций и для задачи 7 — на 9.6 млн итераций. Отношение te/t12, где te — время исполнения на восьми процессорах и t12 — на двенадцати, для задач 6 и 7 соответственно составило 1.9 и 1.7.

Анализ результатов численного эксперимента показывает, что параллельный алгоритм обеспечивает достаточно равномерную загрузку процессоров. С увеличением размерности, например, при числе строк, равном нескольким десяткам тысяч, можно ожидать возрастания простоев на начальном этапе у процессорных элементов с номерами i > 0. Здесь был бы полезен вариант MPI с динамическим определением числа используемых процессоров [7].

Параллельный вариант программы включен в пакет программ ЦКП. В него кроме программ решения задач входят программы обработки данных, которые осуществляют логический контроль данных, их перевод из входного MPS-формата [2] во внутренний разреженный формат и другие преобразования. Пакет организован аналогично пакету ЦЛП [4].

Таблица 3

Характеристики процесса решения задач параллельным алгоритмом

Номер задачи

1 2 3 4 5 6 7 8 9

it1 48243 3650064 7246301 123501 339065 2319167 15065957 2018111 1466794

it2 47246 3770439 4379823 116647 324842 2500000 14285510 2005200 1490414

i GÛ 48071 3504662 4669789 121271 333464 2348354 14172405 2040903 1486702

it4 48566 3947817 4523608 115781 315766 2361807 13030163 2039631 1423969

it5 47621 3745081 3431145 113520 348888 2320268 14535137 1993103 1460954

ite 46654 3707780 3473066 115535 309939 2326606 14784738 2019802 1449547

it7 47678 3413254 3294001 115900 307952 2412781 13322902 2010772 1467747

ite 45102 3532041 3281568 119760 302719 2317997 14299376 1963402 1529396

Sit 379181 29271138 34299301 941915 2582635 18906980 113496188 16090924 11775523

tp 9.18 5187. 8851 45.49 293.4 40439. 38292 502.7 528.6

t/tp 7.27 8.74 5.78 8.37 10.42 — — 8.91 10.88

Примечание: г^ — суммарное число итераций в оценочных задачах, выполненных на г-м процессоре; Бц — суммарное число итераций, выполненное на всех процессорах; Ьр — время решения задач параллельным алгоритмом, Ь/Ьр — отношение времени решения задачи последовательным алгоритмом ко времени решения параллельным алгоритмом.

Заключение

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

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

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

[1] Ковалев М.М. Дискретная оптимизация (целочисленное программирование). Минск: Изд-во Белорус. ун-та, 1977.

[2] Муртаф Б. Современное линейное программирование. Теория и практика. М.: Мир, 1984.

[3] Уилкинсон Дж.Х., РАйнш К. Справочник алгоритмов на языке Алгол. Линейная алгебра. М.: Машиностроение, 1976.

[4] Zabinyako G.I., Kotel'nikov E.A. Linear optimization programs // NCC Bulletin, Series Num. Anal. Novosibirsk: NCC Publisher, 2002. Issue 11. P. 103-112.

[5] Корнеев В.Д. Параллельное программирование в MPI. Новосибирск: Изд-во СО РАН, 2000.

[6] FTP://PLATO.LA.ASU.EDU

[7] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. СПб.: БХВ Санкт-Петербург, 2002.

Поступила в редакцию 1 октября 2003 г.