Высокопроизводительные парал. вычисления на кластерных системах. Воеводин (2005) (Высокопроизводительные парал. вычисления на кластерных системах. Воеводин (2005).pdf), страница 16

А именно: если в одномерной задаче для достижения точности решения ε требуется p вычислений функции, то в задаче с размерностью N для решения с той же точностью необходимо осуществитьαpN испытаний, где α зависит от целевой функции, допустимой области и используемого метода. В свою очередь многоэкстремальностьограничений может порождать области сложной структуры (невыпуклые, неодносвязные и даже несвязные), существенным образомвлияющие на процесс поиска решения.Для анализа задач рассматриваемого класса в данной статье предлагается подход, основанный на идеях редукции сложности, когда решение исходной задачи заменяется решением одной или несколькихболее простых задач. Такая редукция осуществляется с помощью схемы вложенной оптимизации, позволяющей заменить решение многомерной задачи решением семейства рекурсивно вложенных одномерных подзадач, в сочетании с индексным методом учета ограничений [1,2], обобщенным на многомерный случай [1, 3, 4].Индексная функцияВведем обозначение Q0 = D и определим набор вложенных множеств75Q j = { y ∈ Q j −1 : g j ( y ) ≤ 0}, 1 ≤ j ≤ m,(4)для которых справедливо следующее включение:D = Q0 ⊇ Q1 ⊇ K ⊇ Qm = Q ,(5)причем все непустые множества Qj являются замкнутыми вследствиенепрерывности функций gj(y), обеспечиваемой их липшицевостью.Введем понятие индекса v(y) точки y ∈ D либо как номера первогонарушенного в этой точке ограничения (в соответствии с порядкомнумерации, установленным в (2)), т.е.g j ( y ) ≤ 0, 1 ≤ j ≤ ν( y ) − 1, g ν ( y ) ( y ) > 0 ,(6)либо принимающего значение v(y) = m + 1, если в точке y выполненывсе ограничения.Определим максимальное значение M среди индексов точек гиперпараллелепипеда D, т.е.M = max ν( y ) .y∈ DТогда из (4), (5) следует, что если допустимая область Q не пуста, тоM = m + 1, а в противном случаеQM −1 ≠ ∅, Q j = ∅, M ≤ j ≤ m.Переобозначим целевую функцию как gm+1(y) = f (y) и введем вспомогательную задачуg M ( y ) → min, y ∈ QM −1 ,(7)которая всегда имеет решение, поскольку gM(y) непрерывна, а QM–1непусто, замкнуто и ограничено.Сконструируем индексную функцию⎧0, ν( y ) < m + 1,(8)Φ( y) = g ν( y ) ( y) − ⎨ *⎩ f , ν( y ) = m + 1,где f * – минимальное значение функции f (y) в области Q.

Очевидно,что индексная функция обладает свойствами(i) Φ ( y ) > 0 , когда v(y) < m + 1, т.е. точка y является недопустимой;(ii) Φ ( y ) = 0 , когда v(y) = m + 1 и gm+1(y) = f *;(iii) Φ ( y ) > 0 , когда v(y) = m + 1 и gm+1(y) = f *.76Рассмотрим новую задачуΦ ( y ) → min, y ∈ D .(9)Вследствие указанных свойств индексной функции исходная задача (1)–(3) и задача (9) являются эквивалентными в том смысле, чтомножества их глобально-оптимальных решений совпадают, когда область Q не пуста. Таким образом, мы можем заменить отыскание глобального минимума функции f (y) в сложной области Q решением задачи оптимизации индексной функции в гиперпараллелепипеде D.Данный подход, который носит название индексного метода [1, 2],может рассматриваться как вариант метода штрафных функций, однако, в отличие от классической схемы индексный метод не требует подбора каких-либо параметров типа константы штрафа или выравнивающих коэффициентов при ограничениях.

Следует также отметить, чтодля построения индексной функции не требуется, чтобы функции gj(y),1 ≤ j ≤ m + 1 существовали во всех точках области D; достаточно, чтобы каждая функция gj(y) была определена только в области Qj–1. Задачи такого рода называются задачами с частично-вычислимыми ограничениями, и для их решения классический метод штрафных функцийприменить невозможно.С другой стороны заметим, что, несмотря на липшицевость функций исходной задачи, функция Φ(y) является, вообще говоря, разрывной, что требует применения специальных методов оптимизации, учитывающих данное свойство.Вторая особенность индексной функции состоит в том, что в ееформировании участвует величина f *, которая, разумеется, неизвестна.Поэтому в численном эксперименте вместо f * можно использовать ееоценку f k* по результатам вычисления целевой функции, т.е.f k* = min f ( y i ) ,(10)y i ∈Qгде y1, y2, …, yk – точки области D, в которых проведены вычислениязначений функции Φ(y).

Это означает, что при решении задачи (9), мыиспользуем адаптивное семейство функций⎧⎪0, ν( y ) < m + 1,Φ k ( y) = g ν( y) ( y) − ⎨ *⎪⎩ f k , ν( y ) = m + 1,(11)Укажем, что для функции Φk(y) сохраняется свойство (i), а для то77чек yi ∈ Q справедливы свойства (ii)-(iii) при замене f * на f k* .Схема вложенной оптимизацииВведем обозначенияui = ( y1 ,..., yi ), vi = ( yi +1 ,..., y N ),позволяющие при 1 ≤ i ≤ N − 1 записать вектор y y в виде парыy = (ui , vi ) , и примем, что y = v0 при i = 0 и y = uN при i = N.Построим семейство функцийΦ N ( y ) ≡ Φ ( y ),Φ i (ui ) =minyi +1∈[ ai +1 ,bi +1 ]Φ (ui +1 ), 0 ≤ i ≤ N − 1,(12)определенных на множествахDi = {ui ∈ R i : as ≤ ys ≤ bs , 1 ≤ s ≤ i} .Тогда имеет место соотношение [3,4]min Φ ( y ) = miny∈Dmin .

. .y1∈[ a1 ,b1 ] y2∈[ a2 ,b2 ]miny N ∈[ a N ,bN ]Φ ( y ).(13)Как следует из (13), для решения задачи (9) достаточно решить одномерную задачуΦ1 ( y1 ) → min,y1 ∈ [a1 , b1 ] ⊆ R1 .(14)1При этом каждое вычисление функции Φ (y1) в некоторой фиксированной точке y1 ∈ [a1, b1] представляет собой согласно (12) решениеодномерной задачиΦ 2 ( y1 , y 2 ) → min,y2 ∈ [a2 , b2 ] ⊆ R1.Эта задача является одномерной задачей минимизации по y2, т.к.

y1фиксировано.В свою очередь, каждое вычисление значения функции Φ2(y1, y2)при фиксированных y1, y2 требует решения одномерной задачиΦ 3 (u 2 , y3 ) → min, , y3 ∈ [a3 , b3 ] ⊆ R1 ,и т.д. вплоть до решения задачиΦ N (u N −1 , y N ) = f (u N −1 , y N ) → min,78y N ∈ [ a N , b N ] ⊆ R1при фиксированном uN–1.Окончательно решение задачи (9) сводится к решению семейства«вложенных» одномерных подзадачΦ i (ui −1 , yi ) → min,yi ∈ [ai , bi ] ⊆ R1 ,(15)где фиксированный вектор ui–1 ∈ Di–1.Решение исходной многомерной задачи (9) через решение системывзаимосвязанных одномерных подзадач (15) называется многошаговойсхемой редукции размерности.Для решения одномерных подзадач (15) в следующем разделепредлагается параллельный одномерный алгоритм, обобщающий последовательный индексный метод, предложенный Д.Л.Маркиным иР.Г.Стронгиным [1, 2].Параллельный индексный метододномерной оптимизацииРассмотрим одномерную задачу (1)–(3), когда функции gi(y),1 ≤ j ≤ m + 1, зависят от единственной скалярной переменной y и область D является отрезком, т.е.

D = [a, b]. В этом случае функция Φ(y)представляет собой в общем случае разрывную функцию, дуги котороймежду двумя соседними точками разрыва или крайними точками a, bудовлетворяют условию Липшица с константой, которая зависит оттого, какая из функций gi(y) образовала данную дугу.Опишем вычислительную схему параллельного алгоритма, предназначенного для поиска глобального минимума задачи (1)–(3).Назовем испытанием операцию вычисления в некоторой точкеy ∈ [a, b] индекса ν( y ) и значения gv(y)(y), и будем использовать термин итерация для одновременного (параллельного) выполнения нескольких испытаний (каждое на своем процессоре), обозначая числомp(n) количество испытаний n-й итерации, а величиной k(n) – суммарное количество испытаний, выполненных в процессе всех n итераций.На начальной итерации поиска проведем первые испытания вp = p(1) ≥ 2 точках y1, y2, …, yk отрезка [a, b], среди которых в обязательном порядке должны присутствовать концы отрезка a и b, причемэти испытания могут быть проведены параллельно (каждое на отдельном процессоре).

По результатам испытаний будут вычислены индексы vi = v(yi) и значения g νi ( y i ), 1 ≤ i ≤ p. После этого установим номеритерации n = 1 и положим k = k(1) = p.79Общее правило получения координат новой (n + 1)-й итерации(n ≥ 1) состоит в следующем.Шаг 1. Координаты проведенных испытаний y1, y2, …, yk, k = k(n),перенумеровываются нижним индексом в порядке возрастания, т.е.a = y0 < y1 < K < y τ−1 < y τ = b,(16)где τ = k(n) – 1.Шаг 2. Каждой точке yi множества (16) ставится в соответствиеиндекс vi = v(yi) и значение zi функции (11), т.е. zi = Φk(yi).Шаг 3.

Для функций gj(y), 1 ≤ j ≤ m + 1 вычисляются величины⎧⎪ z q − z s⎫⎪λ j = max ⎨: 0 ≤ s < q ≤ k, ν s = νq = j⎬ .(17)⎪⎩ y q − y s⎪⎭Если значение λj не может быть вычислено, оно полагается равнымнулю.Шаг 4. На базе величин (17) определяются оценки µj для константЛипшица Lj функций gj(y), 1 ≤ j ≤ m + 1:⎧⎪rλ j , λ j > 0 ,µj =⎨(18)⎪⎩ 1, λ j = 0где r > 1 – параметр метода.Шаг 5. Каждому подынтервалу (yi–1, yi), 1 ≤ i ≤ τ, ставится в соответствие число⎧( zi − zi −1 ) 22( zi + zi −1 )−, ν i −1 = ν i = j⎪ yi − yi −1 + 2µjµ j ( yi − yi −1 )⎪⎪⎪4z2( yi − yi −1 ) − i , ν i −1 < ν i = j ,(19)R(i ) = ⎨µj⎪⎪4z2( yi − yi −1 ) − i −1 , ν i < ν i −1 = j ,⎪µj⎪⎩называемое характеристикой подынтервала.Шаг 6.

Характеристики R(i), 1 ≤ i ≤ τ, упорядочиваются по убыванию:R(t1 ) ≥ R(t 2 ) ≥ K ≥ R(t τ−1 ) ≥ R(t τ )80(20)Шаг 7. В последовательности (20) выбираются p = p(n + 1) ≤ τнаибольших характеристик с номерами ti, 1 ≤ i ≤ p, и в интервалах,соответствующих этим характеристикам, формируются точки испытаний yk+1, …, yk+p новой (n + 1)-й итерации согласно выражениямy k +szts − zts −1⎧1, ν ts −1 = ν ts = j ,⎪ ( yts + yts −1 ) −2µ j⎪2=⎨⎪1⎪⎩ 2 ( yts + yts −1 ), ν ts −1 ≠ ν ts ,(21)для 1 ≤ s ≤ p.Шаг 8.

Еслиmin( yt s − yt s −1 ) ≤ ε(22)1≤ s ≤ pгде ε > 0 – заданная точность поиска, то поиск прекращается и в качестве решения принимается величина f k* из (10), а также ее координата.Если же (22) не выполняется, осуществляется параллельная итерация метода, т.е. вычисление в точках yk+s, 1 ≤ s ≤ p, индексов этих точек v(yk+s) и значений функции Φk(yk+s) – каждое испытание на своемпроцессоре.