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

УДК 519.853.4

Р. Г. Стронгин, В. П. Гергель, К. А. Баркалов

ПАРАЛЛЕЛЬНЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ ГЛОБАЛЬНОЙ ОПТИМИЗАЦИИ

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

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

Введение. В настоящей работе продолжено исследование в области нового подхода к минимизации многоэкстремальных функций при невыпуклых ограничениях, описанного в [1—6] и получившего название индексного метода глобальной оптимизации. Метод основан на раздельном учете каждого ограничения задачи, он не требует использования штрафных функций. При этом решение многомерных задач сводится к решению эквивалентных им одномерных задач. Соответствующая редукция основана на использовании разверток единичного отрезка вещественной оси на гиперкуб. Роль таких разверток играют непрерывные однозначные отображения типа кривой Пеано, называемые также кривыми, заполняющими пространство. Предложена новая схема построения множества кривых Пеано („вращаемые развертки"), достоинством которой (в отличие от схемы „сдвиговых разверток") является простота в построении и использовании. Описан параллельный индексный алгоритм, основанный на одновременном использовании множества вращаемых разверток. Приведены результаты экспериментов, подтверждающие достоинство новой схемы построения множественных отображений.

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

Ф(/) = тт{ф(>0: у е ° gJ■(у) < 0, 1 < у < (1)

где область поиска О задана как единичный гиперкуб 1-мерного евклидова пространства

О = {у е Я1 : -2-1 < у < 2-1, 1 < I < 1}.

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

этом решение задачи (1) сводится к построению оценки у* е Q, отвечающей некоторому по-

нятию близости к точке у * (например, чтобы у - у* <в или ф(у )-ф(у*)

<в , где в > 0

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

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

обозначаемая также gm+l) и левые части ограничений gJ■ (у) (1 < у < т ) удовлетворяют условию Липшица с соответствующими константами Ьу (1 < у < т +1), а именно

(У1) - gj (У2)\ < Ьу||у1 - У2||, 1 < у < т +1, У1, У2 еО. (2)

В общем случае все эти функции могут быть многоэкстремальными.

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

ф(У(х*)) = min{ф(У(х)) : х е [0,1], gJ■ (у(х)) < 0, 1 < у < т}. (3)

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

gj (У(х')) - gj (У(х"))| < Ку \х' - х"\ 1/^, х', х" е [0,1], 1 < у < т +1, где N есть размерность исходной многомерной задачи, а коэффициенты К у связаны с константами Липшица Ьу исходной задачи соотношением К у < . Общая теория и вопросы численного построения отображений типа кривой Пеано подробно рассмотрены в работах [1—6]. Отметим, что численно построенная кривая является приближением к теоретической

кривой Пеано с точностью не менее 2-к по каждой координате (параметр к называется плотностью развертки).

Для целей дальнейшего изложения введем классификацию точек У из области поиска О с помощью индекса V = v(y ( х)) , который определяется условиями

gj(у(х)) < 0, 1 < у < V-1, gv(у(х)) > 0, (4)

при этом последнее неравенство несущественно, если V = т +1, и удовлетворяет неравенствам

1 < V = v(y(х)) < т +1.

Данная классификация порождает функцию

I (У( х)) = gv (У( х)), V = К у ( х)), (5)

определенную и вычисляемую всюду в О. Ее значение в точке у есть либо значение левой части ограничения, нарушенного в этой точке (случай, когда V < т ), либо значение минимизируемой функции ( V = т +1). Поэтому определение значения I (у( х)) сводится к последовательному вычислению величин gi(у(х)), 1 < г < V = v(х), т.е. последующее значение gi+l (у(х)) вычисляется лишь в том случае, когда gi (у(х)) < 0 . Процесс вычислений завершается либо в результате установления неравенства gi (у(х)) > 0, либо в результате достижения значения v( х) = т +1.

Описанная процедура, называемая испытанием в точке у, автоматически приводит к определению индекса V этой точки. Пара значений

2 = gv (УХ = v(У), (6)

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

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

многомерном пространстве, так как точка х е [0,1] имеет лишь левых и правых соседей, а соответствующая ей точка у( х) е Я1 имеет соседей по 21 направлениям. Как результат, при использовании отображений типа кривой Пеано близким в 1-мерном пространстве образам у', у" могут соответствовать достаточно далекие прообразы х', х' на отрезке [0,1].

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

Гм (х) = {уЧх),..., ум (х)} (7)

вместо применения единственной кривой Пеано у(х) (см. работы [4, 6]). Каждая кривая

Пеано у1 (х) из Ум (х) может быть получена в результате некоторого сдвига вдоль главной диагонали гиперинтервала О. Таким образом, сконструированное множество кривых Пеано позволяет получить для любых близких образов у , у , отличающихся только по одной

координате, близкие прообразы х ', х ' для некоторого отображения у1 (х).

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

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

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

Для иллюстрации на рис. 1 приведены две кривые, являющиеся приближением к развертке Пеано для случая N = 2. Узлы сетки в 1-мерном пространстве отмечены темными кружками, стрелкой указана одна из пар точек, прообразы которых далеки на одномерной оси при одном отображении и близки при другом.

■ \

Г г

1 1 Л I

Л

<

I - —

(- —»-

<

Л —

Рис. 1

Максимальное число различных поворотов развертки, отображающей 1-мерный гиперкуб на одномерный отрезок, составляет 21. Использование их всех является избыточным, требуется выбрать лишь часть из всех возможных вариантов.

1 0 0 0 0

0 0 0 -1 0

0 0 1 0 0

0 1 0 0 0

0 0 0 0 1

1 0 0 0 0

0 0 0 1 0

0 0 10 0

0 -1 0 0 0

0 0 0 0 1

Предлагается осуществлять преобразование развертки в виде поворота на угол ±п /2 в каждой из координатных плоскостей. В качестве примера на рис. 2 приведены матрицы поворота в плоскости ( у2, у4 ) при N = 5.

Число подобных пар поворотов определяется числом координатных плоскостей пространства:

CN = N (N -1)/2, а общее число преобразований будет равно N^ -1). Учитывая исходное отображение, приходим к заключению, что данный способ позво-Рис.2 ляет строить до N(N -1) +1 разверток для отображе-

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

В случае необходимости данный способ построения множества отображений может

быть легко „отмасштабирован" для получения большего (вплоть до 2N ) числа разверток.

Использование множества отображений YM(х) = {у1(х), ..., ум(х)} приводит к формированию соответствующего множества одномерных многоэкстремальных задач

тт{ф(/ (х)) : х е [0,1], gj (у1 (х)) < 0,1 < у < т}, 1 < l < М . (8)

Каждая задача из данного набора может решаться независимо, при этом любое вычисленное значение 2 = gv(у'), у' = у1 (х') функции gv(у) в г-й задаче может интерпретироваться как вычисление значения г = gv (у'), у' = ( х" ) для любой другой 5-й задачи без повторных трудоемких вычислений функции gv (у). Подобное информационное единство позволяет решать исходную задачу (3) путем параллельного решения индексным методом М задач вида (8) на наборе отрезков [0,1]. Каждая одномерная задача решается на отдельном процессоре. Для организации взаимодействия на каждом процессоре создается М очередей, в которые процессоры помещают информацию о выполненных итерациях. Используемая схема не содержит какого-либо единого управляющего процессора, что повышает надежность выполняемых вычислений.

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

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

2. Проинформировать остальные процессоры о начале проведения испытания в точке ук (блокирование точки ук ).

3. Вычислить величины gl(yk),...,gv(ук), где значения индекса V < т определяются

условиями

gi (ук ) < 0, 1 < у < V, gv (ук ) > 0, V < т .

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

(9)

4. Определить прообразы х е[0,1] (1 < I < М ) точки у и интерпретировать испыта-

,к1

к

ние, проведенное в точке у ёП, как проведение испытаний вМточках

к

(10)

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

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

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

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

Правило 1. Изъять из всех очередей, закрепленных за данным процессором, записанные

для него результаты, включающие множество Yq = {у^ :1 < г < } точек итераций в области Б и вычисленные в них значения индекса и величин из (9). Определить множество Хч = {х^ :1 < г < Sq } прообразов точек множества Ук при соответствии у1 (х) . Правило 2. Точки множества итераций

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

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

{Х1} и Х1 и... и Х(

q

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

0 = Хо < Х1 <... < х1 <... < Хк < Хк +1 = 1,

(11)

I _2 = {0, к + 1}, I_1 = {г: 1 < г < к, v(х ) = -1}, ^ = {г: 1 < г < к, v(xi) = v}, 1 < v < т +1,

30 Р. Г. Стронгин, В. П. Гергель, К. А. Баркалов

содержащих номера всех точек xi (1 < i < k) с индексами, равными одному и тому же значению v. Определить максимальное значение индекса

vmax = тах^ = ЧX X 1 < i < k}.

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

= max

z- — z ■ | - ■ : j <h j е Iv, i е Iv

(X — Xj )VN

(12)

для относительных разностей функций gv (1 < V < т +1). Если множество 1У содержит менее двух элементов или если ^ из (12) оказывается равным нулю, то принять ^ = 1. Из (12) следует, что оценки ^ являются неубывающими начиная с момента, когда (12) порождает первое положительное значение ^.

Правило 4. Для всех непустых множеств ^ (1 < V < т +1) вычислить оценки

z* = I £v, v<v max,

V imin{£v (xi):iеIv } v = vmax,

где вектор

£r = (S1,..., sm ), (13)

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

Правило 5. Для каждого интервала (х-—1, х-) (1 < i < к +1) вычислить характеристику R(i)

R(i) = А — 2(Z + Z'—1 — 2Zv), v = v(x,.—1) = v(х-,

rv2 ^v А-

R(i) =2А- — v(x.—1) <v(x.) = v,

rv Mv

R(i) =2А- — 4( z-—Zv ), v = v(Xi—1) >v(X-),

rv Mv

А- = (х- — X-—1)1/N . Величины rv > 1 (1 < v < m +1) являются параметрами алгоритма.

Правило 6. Определить интервал (xt—1, xt), которому соответствует максимальная характеристика

R(t) = max{R(i) :1 < i < к +1}. (14)

Правило 7. Провести очередное испытание в серединной точке интервала (xt—1, xt), если индексы его концевых точек не совпадают, т.е.

Xq+1 = (Xt + Xt—1) / 2, v(Xt—1) * v(Xt) . В противном случае провести испытание в точке

0+1 = . ^ • . X 1

=(Xt + Xt—1V2—sign( zt— zt—1 )

2r

v

zt — zt—1

v = v( Xt—1) = v( Xt).

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

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

Операционные характеристики и сравнение алгоритмов. Один из известных подходов к оценке эффективности методов безусловной глобальной оптимизации основан на численном решении этими методами всех задач из некоторой случайно генерируемой выборки большого объема (см., например, [7]).

Генераторы таких выборок можно рассматривать как некоторые классы функций / (у) (у ёП) с определенной на них вероятностной мерой. Один из таких генераторов О , предложенный в работе [8], порождает многоэкстремальные функции в соответствии со схемой

Ф( уь у2) = -

( 7 7 V ( 7 7

ИГ

^ г=1 /=1

[А/а/ (у1 ,у2)+В/Ъ/ (л,у2)] + ЕЕ[// (у1 ,у2)-// (у1 ,у2).

г=1 /=1

12

где ау (у1, у2) = ^[п(кгу1)^[п(ц/у2), Ъ/ (у1, у2) = ео8(пгу1) со8(п/у2) определены в области 0 < у1, у2 <1, а параметры -1 < А/, В/, С/, Б/ < 1 являются независимыми случайными величинами, равномерно распределенными в указанном выше интервале. Минимизация подобных функций возникает, например, в задаче оценки максимального напряжения (определяющего прочность) в тонкой пластине при поперечной нагрузке.

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

{ук}. Указанная последовательность прерывается (т.е. процесс решения прекращается)

либо в связи с первым попаданием точки очередного испытания в заданную 5 -окрест*

ность решения х , либо в связи с тем, что в ходе выполнения заданного числа К испытаний такое попадание не имело места. В проведенных численных экспериментах использовалось значение К = 1000.

*

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

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

попадания точек испытаний в заданную 5 - окрестность решения. Такую функцию будем называть операционной характеристикой алгоритма £ .

Эксперименты проводились для описанных выше алгоритмов со сдвиговыми (С) и вращаемыми (В) развертками с плотностью к = 12, число разверток М = 2 . При этом использовались параметры ^ = 2,1, 8 = 0,01 и резервы 8V = 0,05 (0 < v < 1) из (13). Данные параметры

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

Операционные характеристики для обоих методов, полученные на выборках, рождаемых генератором О, представлены на рис. 3. При этом кривая 1 характеризует метод со сдвиговыми развертками, 2 — метод с единственной разверткой, 3 — метод с вращаемыми развертками. Указанное расположение кривых показывает, что алгоритм с В-развертками

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

Рис. 3

С целью иллюстрации использования параллельного алгоритма с множеством В-разверток рассмотрим известную задачу минимизации функции Растригина

N

ф( у) = Е (yi - cos(18 y2)), -1,5 < y < 1,5, 1 < i < N, N = 6, i=1

* *

где минимальное значение ф(у ) = - N достигается в точке y = 0 .

Для решения данной задачи использовались последовательный и параллельный методы с вращаемой разверткой с плотностью к = 10, использовался параметр метода r = 2, и s = 0,05 — в критерии остановки. Число разверток M = 30 соответствовало числу процессоров, задействованных при решении задачи. И последовательный, и параллельный алгоритмы нашли решение с требуемой точностью, при этом последовательный алгоритм выполнил 173 116 итераций, а параллельный — 8535 (максимальное число итераций на одном процессоре). Ускорение по числу итераций составило 20,3, а ускорение по времени — 7,5.

Работа выполнена при финансовой поддержке РФФИ (грант № 07-01-00467-а) и Совета по грантам Президента Российской Федерации по государственной поддержке ведущих научных школ Российской Федерации (грант № НШ-4694.2008.9).

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

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

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

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

4. Стронгин Р. Г. Параллельная многоэкстремальная оптимизация с использованием множества разверток // Журн. вычислительной математики и математической физики. 1991. Т. 31, № 8. С. 1173—1185.

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

6. Strongin R. G., Sergeyev Ya. D. Global optimization with non-convex constraints. Sequential and parallel algorithms. Dordrecht: Kluwer Academic Publishers, 2000.

7. Гергель В. П., Стронгин Р. Г. Абсолют. Программная система для исследований и изучения методов глобальной оптимизации. Н. Новгород: Изд-во Нижегородского ун-та, 1998.

8. Гришагин В. А. Операционные характеристики некоторых алгоритмов глобального поиска // Проблемы случайного поиска. Рига: Зинатне, 1978. № 7. С. 198—206.

Роман Григорьевич Стронгин

Виктор Павлович Гергель

Константин Александрович Баркалов

Сведения об авторах

д-р физ.-мат. наук, профессор; Нижегородский государственный университет им. Н. И. Лобачевского; кафедра математического обеспечения ЭВМ; президент; E-mail: strongin@unn.ac.ru

д-р техн. наук, профессор; Нижегородский государственный университет им. Н. И. Лобачевского; кафедра математического обеспечения ЭВМ; декан; E-mail: gergel@unn.ru

канд. физ.-мат. наук; Нижегородский государственный университет им. Н. И. Лобачевского; кафедра математического обеспечения ЭВМ; ст. препод.; E-mail: barkalov@fup.unn.ru

Рекомендована институтом

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

УДК 681.3

Б. Я. Штейнберг

БЛОЧНО-РЕКУРСИВНОЕ ПАРАЛЛЕЛЬНОЕ ПЕРЕМНОЖЕНИЕ МАТРИЦ

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

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

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

В настоящей работе предложен новый параллельный алгоритм перемножения матриц с количеством пересылок, меньшим, чем у существующих. Для быстрого алгоритма перемножения матриц Штрассена рассмотрена сложность по обращениям к памяти и обоснована нецелесообразность его использования для современных компьютеров; предложена процедура блочно-аффинного размещения данных. Следует отметить, что параллельные алгоритмы перемножения матриц на вычислительных системах с распределенной памятью, описанные в работах [2—5], предполагают размещение матриц полосами или блоками — это частные случаи блочно-аффинных размещений.

Приведенные в [2—5] алгоритмы требуют п / р параллельных умножений. Объем пересылаемых данных каждого процессора для алгоритмов Фокса и Кэннона [4, 6] составляет

П / р1/2. Объем пересылаемых данных каждого процессора в алгоритмах, предложенных в

2 2/3

[3], равен п / р . Объем пересылаемых данных каждого процессора в предлагаемом нами