К анализу метода имитации отжига в многоэкстремальном случае Текст научной статьи по специальности «Математика»

УДК 519.245 МЯО 65С05

Вестник СПбГУ. Математика. Механика. Астрономия. Т. 4 (62). 2017. Вып. 2

К АНАЛИЗУ МЕТОДА ИМИТАЦИИ ОТЖИГА В МНОГОЭКСТРЕМАЛЬНОМ СЛУЧАЕ*

С. М. Ермаков, Д. В. Куликов, С. Н. Леора

Санкт-Петербургский государственный университет,

Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9

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

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

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

Метод имитации отжига [1-3] широко применяется для решения прикладных задач по поиску абсолютного экстремума. Его математическая сторона кратко может быть изложена следующим образом. Пусть /(X), X € В, В € Дя, ограниченная функция, |/(X)| ^ М. Как правило, В также ограничена в Дя. Случай неограниченной В требует отдельного рассмотрения. Пусть также Т — вещественный параметр, Т £ Д1. Можно показать, что функция (плотность распределения)

ЩТ,Х) = С(Т)ехр(-Цр-У (1)

где С — константа нормировки, зависящая от Т, при Т ^ 0 слабо сходится к ¿-функции, сосредоточенной в точке глобального минимума /(X).

Алгоритм состоит в моделировании методом Метрополиса [1] последовательности плотностей (1) с убывающими значениями Т1 > Т2 > ... > Тп. Реализации соответствующих случайных величин, полученные при , служат исходными для алгоритма Метрополиса при Тк+1, к = 1, 2,... ,п — 1.

В [4] предлагалось для нахождения глобального максимума использовать другую функцию типа 9Л(Т, X), а именно

* Работа выполнена при финансовой поддержке РФФИ (грант № 14-01-00271-а). (¡5 Санкт-Петербургский государственный университет, 2017

При п ^ то функция М(Х, п) также сходится к функции глобального максимума. Легко показать, что применение алгоритма Метрополиса приводит к модификации метода имитации отжига, где параметр п играет роль Т. Далее мы будем рассматривать алгоритмы, основанные на моделировании М(Х, п) при п ^ то. При этом в работе специально рассматривается случай двух и более одинаковых максимумов (минимумов) исследуемой функции.

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

Задача 1. Требуется найти приближенные значения вещественных корней системы нелинейных уравнений д;(Х) = 0, X(хх,..., хв),1 = 1, 2,..., в.

Задача будет решена, если найти все точки глобального минимума функции

в

Н (X ) = £ Р1д21(Х) = шш, (3)

1=1

р; — заданные константы, р; > 0. (Очевидно задачу легко переформулировать в задачу нахождения максимума.)

Задача 2. Построение точного Б-оптимального плана регрессионного эксперимента. Самая простая из этих задач в случае двух переменных (построение плоскости по трем точкам (хх,у\), (х2,у2), (хз, уз)) требует нахождения наибольшего значения определителя

||1,х<,у<У3=1. (4)

Если область, где строится оптимальный план, является квадратом [0,1]2, решение известно. Это (0, 0), (1, 0) и любая точка отрезка х € [0,1] при у = 1. Также оптимальными будут планы, полученные из этого плана поворотами на г = 1,2, 3.

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

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

Пусть (А, Б, —пространство с а-алгеброй подмножеств Б и конечной мерой /(X) € Ь(Б,^), / является ^-интегрируемой, ограниченной, неотрицательной функцией в Б, имеющей некоторое множество М глобальных максимумов.

Рассмотрим последовательность распределений, заданных следующим образом:

Гл /"(х)и(^х)

тп{Л) = ! , м , , ^еа, п= 1,2,... (5)

)п /"(х)М«х)

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

I. Пусть ц — дискретная мера, сосредоточенная в точках хх,х2,... Предварительно докажем следующую лемму.

Лемма. Пусть qi — вещественные числа, 0 < qi < 1 и

то

Е qi < то, (6)

i=l

тогда справедливо 'равенство

то

11т £ qk =0. (7)

к^то ' i=1

Доказательство. Для Уе > 0 в силу условия (6) найдется N £ N1, такое что +1 qi < е/3 и, поскольку qi < 1, такое к £ N1, что для любого г < N верно qik < 2— • е/3.

Тогда для выбранных N и к будем иметь

то N то N то то

5>? = 5>? + Е + Т-К- 53 9?< ^

i=1 i=1 i=N +1 i=l i=N +1 i=N +1

Пусть теперь ^ — дискретное мультимодальное распределение с т одинаковыми модами

^к = (Лк к = 1,2,...,

V ск ql' Ск q2 ' ... /

причем ж^, ж^,..., ж^ —одинаковые моды, qil = qi2 = ... = ^, а Ск —константа нормировки, ск = (^qik)-1. Обозначим

1 1/т 1/т ... 1/т Теорема 1. при к ^ то слабо сходится к Qo.

Доказательство. Покажем, что если ^^ £ {qil,..., qim}, то скqj; ^ 0 при к ^ то:

к / ^т + £ (^/q*)k'

i=1 5=1

При к ^ то второе слагаемое знаменателя стремится к нулю по лемме и очевидно вся дробь стремится к нулю также. I

Таким образом, предельным распределением является равномерное распределение на множестве равных экстремумов.

II. Вернемся к общему случаю. Обозначим О = тах(/(X)),Х £ В), и положим

М = {ж £ В : /(X) = О}. (8)

Не умаляя общности, будем считать, что О =1. Справедлива следующая теорема.

Теорема 2. Если ^(М) > 0, то справедливо равенство

¡ЛГ(хЫг1х) = М(АПМ) п^оо ¡в/п(х)ц(<1х) ц(М) ' 1 '

Доказательство. Как и в дискретном случае имеем

JAfn(x)/x(dx) = fAгш fn(xM<b) + 1л\м Id fn(xMdx) /м fn(x)fi(dx) + ¡тм fn(x)n(dx)

«(ЛП M) +SA\M fn(x)«(dx) «(A П M) = lim ---= —---

™ +JDWf"(x)^dx) M(M) '

поскольку интеграл на D \ M стремится к нулю. Напомним, что на множествах A\M и D \ M выполнено f (x) < 1 для всех x(mod«). ■

Случай «(M) = 0 теоретически интересен, но требует отдельных исследований и дополнительных предположений относительно поведения f (X) в окрестности точек экстремума.

Тем не менее, для построения прикладных алгоритмов (решения задач 1 и 2, например) результат теоремы 2 может быть весьма полезным.

Введем в рассмотрение функцию (0 < е < 1)

G Г F(X), F(X) < 1 - е,

G (X )=\ 1 - е, F (X) > 1 - е. (10)

Если «({X : F(X) = 1}) = 0, мы будем предполагать, что при каждом е > 0 справедливо неравенство «({X : Ge(X) = 1 — е}) > 0.

Во всяком случае так будет в случае меры Лебега « для непрерывной F. При этом множество {X : Ge(X) = 1 — е} может состоять из набора изолированных, при достаточно малых е, множеств Mj, i = 1, 2,..., на каждом из которых Ge(X) = const, но /«(Mj) не обязано, вообще говоря, совпадать с «(Mj) при i = j.

В соответствии с теоремой 2 абсолютные экстремумы Ge(X) будут распределены равномерно на множестве Me = {X : Ge (X) = 1 — е}.

Алгоритм, являющийся аналогом метода имитации отжига, предполагает задание последовательности (е&, nk), k =1, 2,... При каждом е^ и nk метод Метрополиса используется для моделирования (Gefc(X))nkCk. Здесь ек —убывающая, а nk —возрастающая последовательности и Ck —константа нормировки. При каждом k заполняется массив m реализаций моделируемых случайных величин.

Здесь, очевидно, возможны различные усовершенствования, зависящие от характера задачи. Более сложная, но во многом аналогичная ситуация возникает при решении задачи 2.

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

Пример 1. Рассмотрим функцию, определенную на промежутке [0; 2п] и имеющую два одинаковых максимума:

f (x) = / Sin(x) 0 < Х < ^

f ( ) [ sin(x - п), n<x < 2п.

Введем функцию согласно (10):

f (x), f (x) < 1 - е,

fe(x) ^ 1 - е, f (x) > 1 - е.

Вестник СПбГУ. Математика. Механика. Астрономия. Т. 4(62). 2017. Вып. 2 223

Рис. 1. Графики функций: /(х) —сплошная линия, /(х) —пунктирная линия.

F

Рис.2. Графики функций Е(га,/(х)), значениям п = 1, 2, 3 соответствуют пунктирные линии, п = 8, 9, 10 — штриховые и п = 18, 19, 20 — сплошные.

Пусть для определенности е = 0.1. Графики функций /(ж) и /е(ж) изображены на рис. 1. Функция имеет два глобальных максимума в точках п/2 и 3п/2.

Рассмотрим класс функций ^(п, /(ж)), где п — параметр, /(ж) —функция, глобальные экстремумы которой подлежат определению

/ "(т)

Р(пЛх))=^Л> . (И)

/ /™(ж)йж

о

Функция вида (11) при п ^ то, как было доказано выше, стремится к равномерному распределению на множестве точек одинаковых глобальных экстремумов. Построенные для ряда значений п графики, приведенные на рис. 2, наглядно это демонстрируют. В соответствии с теоремой 2 абсолютные экстремумы будут распределены на множестве Ме = {ж : /е(ж) = 1 — е}. В данном случае Ме = М1 и М2, причем ^(М1) = ^(М2).

Пример 2. Рассмотрим функцию, определенную на промежутке [0; 5п/4] и имеющую два одинаковых максимума, отличающихся «крутизной» вершин:

, ч _ ( вт(ж), 0 < ж < п,

5(ж) = | вт4(ж — п), п < ж < 5п/4.

В этом случае для поиска экстремумов можно рассматривать класс функций вида

5П(ж)

О(п,^(ж))

5тт/4

о

Построенные для ряда значений п графики функции О(п, /(ж)), приведены на рис. 3.

В данном случае абсолютные экстремумы будут распределены равномерно также на двух изолированных множествах М1 и М2, Ме = М1 и М2, причем ^(М1) = ^(М2). Пример 3. Рассмотрим функцию двух переменных

г(ж, у)

1 + 5(у — ж2)2 + 2(у + ж2 — 1)2'

G

4

Рис.3. Графики функций 0(и, / (х)), значениям п = 1, 2, 3 соответствуют пунктирные линии, п = 8, 9, 10 — штриховые и п = 18, 19, 20 — сплошные.

Рис.4■ График функции Z (1,z(x,y)).

Известно, что данная функция имеет два одинаковых максимума в точках (жтах,

утах) = (±1/^2,1/2).

Пусть В = {(х, у) : —1 < х < 1, —1 < у < 1}. Выберем е = 0.1. Моделирование проводится аналогично двумерному случаю. Рассматривается класс функций %(и,г(х,у)), введенных согласно (10). На рис. 4 приведен график функции %(1, г(х,у)).

Рис. 5. Графики функции Z(n, z(x,y)) при n =5 (слева), n = 15.

На рис. 5 приведены графики функции Z(п, z(x, y)) при различных значениях п. Абсолютные экстремумы распределены равномерно на двух изолированных множествах.

Литература

1. Metropolis V. N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E. Equations of state calculations by fast computing machines //J. Chem. Phys. 1953. Vol.21. P. 1087—1092.

2. Ingber L. Adaptive simulated annealing (ASA): Lessons learned //J. Control and Cybernetics.

1995.

3. Winkler G. Image analysis, random fields and Markov chain Monte Carlo methods: a mathematical introduction. Vol.27. Springer Science and Business Media, 2012.

4. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1975.

Статья поступила в редакцию 26 сентября 2016 г.; рекомендована в печать 22 декабря 2016 г. Сведения об авторах

Ермаков Сергей Михайлович —доктор физико-математических наук, профессор; [email protected]

Куликов Денис Валерьевич — студент; [email protected]

Леора Светлана Николаевна — кандидат физико-математических наук, доцент; [email protected]

TOWARDS THE ANALYSIS OF THE SIMULATED ANNEALING METHOD IN THE MULTIEXTREMAL CASE

Sergey M. Ermakov, Denis V. Kulikov, Svetlana N. Leora

St. Petersburg State University, Universitetskaya nab., 7—9, St. Petersburg, 199034, Russian Federation; [email protected], [email protected], [email protected]

There are many applications, where the number to be calculated global extrema is large, or even infinite. These problems include, for example, some experimental design problems, the problem of solving large systems of equations, and several others. If this is single extreme of a function of several variables, one of the commonly used numerical algorithms is the Simulated Annealing, which is also successfully used in high volume discrete problems (traveling salesman problem). In discrete problems known, that the method of simulated annealing searches equal global extremes with equal probability. Continuous case still has not been investigated. It was assumed that equal extremes to be found consistently, sharing their neighborhood during the computation. It is not always effective method, especially in the case when multiple extrema fills a certain region in Rn.

The results obtained in this article the general outline of the approach to the problem. We give computational examples showing the effectiveness of the approach. It can be used to create programs, algorithms, indicating the localization of the roots of large systems of equations. It can also be noted, that many planning tasks regression experiments are an infinite number of solutions. Refs 4. Figs 5. Keywords: simulated annealing, global extremum, experimental designing.

References

1. Metropolis V. N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E., "Equations of state calculations by fast computing machines", J. Chem. Phys. 21, 1087—1092 (1953).

2. Ingber L., "Adaptive simulated annealing (ASA): Lessons learned", J. Control and Cybernetics (1995).

3. Winkler G., Image analysis, random fields and Markov chain Monte Carlo methods: a mathematical introduction 27 (Springer Science and Business Media, 2012).

4. Ermakov S. M., The Monte Carlo Method and Related Questions (Nauka, Moscow, 1975) [in Russian].

Для цитирования: Ермаков С.М., Куликов Д. В., Леора С. Н. К анализу метода имитации отжига в многоэкстремальном случае // Вестник СПбГУ. Математика. Механика. Астрономия. 2017. Т. 4 (62). Вып. 2. С. 220-226. DOI: 10.21638/11701/spbu01.2017.205

For citation: Ermakov S.M., Kulikov D.V., Leora S. N. Towards the analysis of the simulated annealing method in the multiextremal case. Vestnik SPbSU. Mathematics. Mechanics. Astronomy, 2017, vol. 4(62), issue 2, pp. 220-226. DOI: 10.21638/11701/spbu01.2017.205