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

На правах рукописи

Карандашев Яков Михайлович

ИССЛЕДОВАНИЕ МЕТОДОВ ТРАНСФОРМАЦИИ ЭНЕРГЕТИЧЕСКОЙ ПОВЕРХНОСТИ В ЗАДАЧАХ БИНАРНОЙ МИНИМИЗАЦИИ КВАДРАТИЧНОГО ФУНКЦИОНАЛА

05.13.01 - «Системный анализ, управление и обработка информации» (информационно-вычислительное обеспечение)

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Москва 2013

21 НОЯ 2013

005539407

Работа выполнена в Федеральном государственном бюджетном учреждении науки Научно-исследовательском институт системных исследований Российской академии наук в Центре оптико-нейронных технологий

Научный руководитель: Крыжановский Борис Владимирович,

чл.-корр. РАН, доктор физико-математических наук, руководитель Центра оптико-нейронных технологий НИИСИ РАН

Официальные оппоненты: Литвинов Олег Станиславович,

доктор физико-математических наук, профессор кафедры физики Московского государственного технического университета им. Н.Э. Баумана

Доленко Сергей Анатольевич,

кандидат физико-математических наук, старший научный сотрудник Научно-исследовательского института ядерной физики им. Д.В. Скобельцына Московского государственного университета им. М.В. Ломоносова (НИИЯФ МГУ)

Ведущая организация: Федеральное государственное автономное

образовательное учреждение высшего профессионального образования "Национальный исследовательский ядерный университет «МИФИ»"

Защита состоится «24» декабря 2013 года в 11:00 часов на заседании диссертационного совета Д-002.086.02 Федерального государственного бюджетного учреждения науки Институт системного анализа Российской академии наук (ИСА РАН) по адресу: 117312, Москва, проспект 60-летия Октября, 9.

С диссертацией можно ознакомиться в библиотеке Федерального государственного бюджетного учреждения науки Институт системного анализа Российской академии наук.

Автореферат разослан «11» ноября 2013 года

Ученый секретарь

диссертационного совета Д-002.086.02 доктор технических наук, профессор

Пропой А.И.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы

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

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

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

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

Цели и задачи диссертации

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

Для достижения данной цели ставятся следующие задачи:

1) Исследовать энергетическую поверхность функционала, в частности форму локальных минимумов.

2) Исследовать зависимость вероятности отыскания минимума при случайном поиске от глубины минимума

3) Исследовать влияние трансформации матрицы функционала на изменение глубины минимумов, сдвиг минимумов в пространстве и смещение спектра находимых минимумов.

4) На основе проведенных исследований построить эффективный алгоритм минимизации квадратичного функционала.

Методы исследования

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

эффективности разработанных методов осуществлялась путём численных экспериментов с использованием тестовых задач. Вычислительные эксперименты проводились на персональном компьютере (Intel Core i3 CPU 3.08ГГц, 2 Гб ОЗУ), и обычно занимали от нескольких минут до нескольких дней счёта (в зависимости от размерности задачи и объёма эксперимента).

Основные положения выносимые на защиту

1) Получены выражения для формы минимума, её зависимость от глубины минимума.

2) Показано, что вероятность отыскания минимума при случайном поиске экспоненциально растёт с глубиной минимума.

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

4) На основании предложенной трансформации предложен двухступенчатый алгоритм минимизации, использование которого приводит к экспоненциальному по N росту эффективности случайного поиска по сравнению со стандартным случайным поиском SRS, при этом количество операций, затрачиваемых на один старт, увеличивается вдвое.

Научная новизна

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

Практическая ценность

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

Апробация работы и публикации

По материалам диссертации опубликовано 19 работ, из них 10 - в российских и зарубежных журналах (в том числе 6 из перечня ВАК), 9 - в трудах конференций. Основные результаты работы докладывались на конференциях:

1. XI Всероссийская научно-техническая конференция Нейроинформатнка 2009. Москва, МИФИ (2009)

2. ICONS 2009: The 2nd IFAC International Conference on Intelligent Control Systems and Signal Processing. Turkey, Istanbul (2009)

3. ИКТМР-2009: Международная конференция Актуальные проблемы информационно-компьютерных технологий мехатроники и робототехники. Дивноморск (2009)

4. МСО - 2009: Всероссийская конференция Методы и средства обработки информации. Москва, МГУ (2009)

5. ММРО-14: 14-ая Всероссийская конференция Математические Методы Распознавания Образов. Суздаль (2009)

6. 52-я научная конференция МФТИ. Долгопрудный, МФТИ (2009)

7. X Всероссийская научно-техническая конференция Нейроппформатика-2010. Москва, МИФИ (2010)

8. IJCNN 2010: The International Joint Conference on Neural Networks. Spain, Barcelona (2010)

9. ICANN 2010: The 20 International Conference on Artificial Neural Networks, Greece, Thessaloniki (2010)

10. ICANN 2011: The 21st International Conference on Artificial Neural Networks. Finland, Espoo (2011).

11.1CCG1 2012: The Seventh International Multi-Conference on Computing in the

Global Information Technology. Italy, Venice (2012) 12. ICANN 2012: The 22 International Conference on Artificial Neural Networks. Switzerland, Lausanne (2012)

Объём и структура диссертации.

Диссертация состоит из Введения, 4 глав, Заключения и Списка литературы. Список литературы содержит 85 наименований. Текст диссертации содержит 119 страницу машинописного текста, включая 31 рисунок и 9 таблиц.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

В первой главе даётся введение в задачу минимизации квадратичного бинарного функционала £(5):

построенного на заданной N х N -матрице Т в ТУ-мерном пространстве состояний $ = с бинарными переменными ^ =±1 ; су -

стандартное отклонение матричных элементов Ту.

В качестве основной процедуры локальной минимизации выбрана асинхронная динамика модели Хопфилда. Стандартная (асинхронная) динамика минимизации состоит в следующем: вычисляется величина локального поля А, •—8£(5) / дя,, действующего на произвольно выбранный /-й спин:

Если спин направлен вдоль поля (hlsl >0), то его состояние не изменяется. Если спин находится в неустойчивом состоянии (Л,.^ < 0), то он разворачивается вдоль поля, принимая состояние si = sign А,. Эта процедура последовательно применяется ко всем спинам до тех пор, пока система не конвергирует в устойчивое состояние. Уменьшение энергии при каждом перевороте неустойчивого спина гарантирует нам, что процесс, описываемый асинхронной динамикой, за конечное число шагов приводит систему в устойчивое положение с минимумом энергии. Конечно, найденный минимум скорее всего окажется локальным, в то время как нам бы хотелось найти самый глубокий (глобальный) минимум функционала. Стандартным приёмом здесь является алгоритм случайного поиска (Standard Random Search (SRS)), суть которого состоит в многократном повторении описанных выше спусков из различных случайных стартовых конфигураций.

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

Пусть S0 =(s\°),s2)- конфигурация, соответствующая некоторому минимуму £0 = E(S0). Совокупность точек, находящихся на расстоянии п (по Хеммингу) от S„, назовём п-окрестностью точки S0. Тогда, объём области притяжения минимума (число всех точек пространства, стартуя с которых спускаемся в S0) равен:

О)

N

(2)

где Рп - вероятность попадания из «-окрестности в минимум SQ, которую мы оценим ниже.

Вычислим вероятность перехода из «-окрестности в (и-1)-окрестность, которая совпадает с вероятностью события А( (Sn)s-0) > 0. Как показано в диссертации, величина АД5„)л'<|) представляет собой нормально распределенную величину со средним {N -2п)\Ей\ и дисперсией NOj. Соответственно, для искомой вероятности перехода из п- в (п -1) -окрестность получим Рг{/г( {S, )sf11 > 0} = Ф(у„), где:

<D(jg = _L \e^dz, r„=y(l-2n/N), (4)

yl2n i

где у = I^In/T/ - глубина минимума S0. Тогда выражение для вероятности конвергенции из 5„ в минимум S0 примет вид: Р„ =Ф(У„)- ®()/„_i) • — • Ф^)-Подставляя это выражение в (3), используя соотношение Стирлинга и переходя от суммирования к интегрированию по переменной х = п/ N для объема области притяжения получим:

0.5 NF( I)

V = j2N Г ,е dx, (5)

QJ s] 7ix(l-x)

где

X

F(x) = -xlnx-(l-x) In(l -x)+ Jin [Ф(гг ))dz, yT = у (1 - 2z). (6)

о

Интеграл (6) оценим методом перевала. Приравняв производную F(x) к нулю, получим уравнение для седловой точки хд:

-^ = Ф[К1-2х0)] (7)

1 хп

и выражение (5) преобразуется к виду:

V = 2eNFW (8)

Величина х0 фактически представляет собой радиус области притяжения минимума, поскольку объем области притяжения с точностью до несущественного коэффициента совпадает с величиной х". Анализ выражений (7-8) показывает, что с ростом глубины минимума у радиус области притяжения возрастает. Одновременно экспоненциально возрастает и объем области притяжения. Чтобы убедиться в этом, получим из (7-8) более простые оценочные выражения. Эксперимент показывает, что глубина среднестатистического минимума у отличается от глубины глобального минимума у0 всего лишь на 10-15%. Это позволяет для относительно глубоких минимумов (у < у <у0) представить функцию ха=х0(у) в виде разложения в окрестности значения у, которое наиболее легко определяется в эксперименте. Ограничиваясь линейным членом разложения, получим оценочное выражение:

хъ*х+кх(у-у) (9)

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

Р{у) = У/2К получим:

Р(у) = РекЩг-7) (10)

где Р-Р(у) есть вероятность отыскания минимума среднестатистической глубины у.

Для некоторого среднестатистического минимума, характерного для матриц модели Шеррингтона-Киркпатрика (8К-модель), глубина равна у = 1.36. Из (7) и (9) получаем для него * « 0.384 иЬ 0.076. Это означает, что при случайном поиске большая часть точек, из которых система конвергирует в данный минимум, расположена в его и-окрестности, где п и 0.3 847^. Соответственно, вероятность отыскания минимума глубины у равна

РяР . Как видим, чем глубже минимум, тем больше вероятность его

отыскания, причём зависимость экспоненциальная. При небольших значениях /V данную зависимость можно подтвердить экспериментально.

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

<">

где Еа(п) и <т0(п) - соответственно среднее значение и стандартное отклонение энергии в п-окрестности минимума. Энергетическая поверхность (II) очень гладкая, поскольку дисперсия энергии в п-окрестности с ростом размерности стремится к нулю как а ~1/ N. Более того, все минимумы имеют одинаковую форму, задаваемую выражением (11), причём ширина минимума (на полувысоте) прямо пропорциональна его глубине Е0. Это и означает, что число точек внутри области притяжения экспоненциально возрастает с ростом глубины минимума.

Известно (КгугЬапоувку, 2008), что любую матрицу можно представить в виде взвешенного разложения по внешним произведениям случайных конфигурационных векторов:

Т^тЬЛ8т , (12)

лмО

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

где ас » 0.138 - критическое значение параметра загрузки в модели Хопфилда. Это утверждение прежде всего относится к точке 50, которая по определению является минимумом функционала (1), и для нее справедливы соотношения:

1 >га>гс, -1 <Е0<ЕС, Ес ~ —гс (14)

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

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

, с.—^ (15)

где стм - стандартное отклонение элементов матрицы М. Все рассуждения проведем на примере трансформации глобального минимума. Очевидно, что при трансформации поверхности глобальный минимум сдвигается в пространстве, одновременно изменяются его глубина и ширина области притяжения. Ниже мы покажем, что при относительно небольших значениях параметра к {2<к<5) указанная трансформация приводит к значительному увеличению глубины минимума при относительно небольшом сдвиге.

Чтобы нивелировать сдвиг, предлагается следующий двухэтапный алгоритм минимизации. На первом этапе осуществляется спуск по поверхности Ек(Б) и находится конфигурация являющаяся минимумом функционала £¿.(5). На втором этапе проводится коррекция: из точки по поверхности £(5) спускаемся в ближайший минимум 5т функционала £'(5'). Алгоритм спуска по поверхности Ек(Б) обычный: рассчитывается локальное поле й/*' = -дЕк (8) / дз, на ¡-м. спине:

Т=2>л об)

и при й^'л, < 0 состояние спина обновляется по правилу 5,. = /?,<1>.

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

DDK (Double Descent Algorithm с параметром к = 2,3,...) и сравнивать его эффективность со стандартным алгоритмом случайного поиска SRS. Во избежание недоразумений отметим, что при к = 1 трансформированный функционал Ек (S) совпадает с исходным £(5), а предлагаемый алгоритм DDK в этом случае есть не что иное, как стандартный алгоритм случайного поиска SRS {DD\=SRS).

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

Углубление минимумов. Рассмотрим энергию E2Q =£2(50) в точке Sa. Как показано в диссертации, в пределе N»1 энергию Е10 можно рассматривать как нормально распределенную величину со средним Е10 и относительно малым стандартным отклонением аЕ:

£20=-V^£02, сге =(l-r02)/ N (17)

Поскольку Ё201Е0 = r0-jN > rc-jN «1.35, то следует ожидать, что с ростом N вероятность углубления Рг{£'20 <£<,} асимптотически стремится к единице. Иными словами, трансформация поверхности с подавляющей вероятностью приводит к существенному увеличению глубины минимума (примерно в 1.35 раз), что в соответствии со сказанным в главе 2 должно приводить к экспоненциальному по N возрастанию вероятности отыскания глобального минимума.

Сдвиг минимума. Среднюю величину сдвига можно представить в виде dk = NP, где Р = Рф,'0^"0 <01 Л,5,<0) >0} - вероятность того, что направления спина j((0) и локального поля h't] в конфигурации минимума Sa не совпадают. Для рассматриваемого здесь случая к = 2 эта вероятность выражается через функцию ошибок:

Р = \e(l-r)1 ег/с{4гух)с1х (18)

Pj* о

где Pa=\ + erfy и y = -J~Nr0.

Рис. la показывает наблюдаемое углубление глобальных минимумов. Как видим, глобальные минимумы становятся глубже для к = 2,3,4 и даже для к = 5, однако с увеличивающейся дисперсией. Наибольшее среднее углубление наблюдается для к = 3. Вызываемый трансформацией сдвиг минимума в пространстве относительно невелик (рис. 1А). Видно, что наименьший сдвиг также наблюдается при к = 3.

Рис. 1 Углубление а) и сдвиг Ь) глобальных минимумов при трансформации поверхности при возведении матрицы в степень к = 2,3,...,7. Матрицы 20 ЕА размера N = 100,196,484. Планки слева показывают ст. отклонение.

Эффективность алгоритма DDK. Для иллюстрации работоспособности предложенных идей мы выбрали матрицы модели 3-мерной модели Эдвардса-Андерсона размера N = Z? при L = 8,10,20. Эксперименты включали в себя

106 стартов для динамики Хопфилда и 104 стартов для динамики Хоудайера-Мартина (НМ) (см. дальше). Результаты были усреднены по 100 матрицам каждого размера. На рис. 2а показана типичная картина распределения плотности вероятности для алгоритмов DD3 и SRS. Как видно из рисунка, кривые, полученные для достаточно больших значений N (больше 500), хорошо аппроксимируются гауссовым распределением с параметрами: 1) для DD3 среднее значение Е} ~ -1.29 / -J~N и стандартное отклонение сг3 я 0.41 / N; 2) для SRS среднее а -1.16 / \[n и стандартное отклонение сг, » 0.66 / N.

Алгоритм DD3 имеет существенное преимущество перед SRS. При JV >1000, ни в одном из экспериментов алгоритм SRS не нашёл минимум с энергией Е<Ег. Этот результат вполне ожидаем поскольку вероятность такого события оказывается <1.3-10"'° для уже при N = 1000, т.е. выигрыш составляет 10 порядков!

Ещё раз подчеркнём, что предлагаемый в диссертации метод основан не на динамике, а на изменении статики, поскольку трансформирует энергетическую поверхность. Однако можно было бы попробовать также модифицировать динамику спуска по поверхности. С этой целью мы попробовали динамику Хоудайера-Мартина (НМ) вместо динамики Хопфилда.

Р{Е)

so

70

so 50 40

зо 20 10 о

Р(£)

W-/WM

\

N-JJ2

SRS

г. s -жю

-1.4

-1.1

Рис. 2. а) Плотность вероятности Р(Е) нахождения минимумов с помощью алгоритма SRS (точечная кривая) и с помощью алгоритма DD3 (сплошная кривая); Ь) Сплошные кривые показывают плотность вероятности Р{Е), получаемую модифицированным алгоритмом случайного поиска DD3-HM; для сравнения, точечная кривая показывает плотность вероятности, получаемая стандартным алгоритмом SRS. Данные для 3D ЕА модели с L - 8,10,20. По оси абсцисс отложена энергия помноженная на 4Ñ.

Алгоритм DD3-HM позволяет найти глобальный минимум для 2D ЕА матриц размера N = 100 с вероятностью примерно 6%, которая в 30 раз выше чем вероятность (порядка 0.2%) даваемая алгоритмом DD3 и в 25000 раз больше чем вероятность (~ 2.4 • 10"6) даваемая алгоритмом SRS. Результаты для DD3-HM для больших размерностей показаны на рис. 2Ь. Пики описываются гауссовым распределением со средним и стандартным отклонением:

Е t_

где £ =1.331,1.327,1.315 и t = 0.37, 0.40,0.46для iV = 512,1000,8000 соответственно. Ясно, что алгоритм DD3-HM показывает гораздо лучшие результаты чем DD3. Однако с ростом размерности это преимущество исчезает: пик кривой плотности вероятности сдвигается в сторону пунктирной линии, которая соответствует центру распределения DD3 на рис. 2а (£ —>1.29).

Применения алгоритма к задаче разбиения графа Суть задачи о разбиении графа в следующем: пусть имеется некоторый граф G, у которого N вершин F = (v,,v2,...,vJV), а веса рёбер описываются матрицей Т размера NxN. Требуется разбить вершины графа на две равные непересекающиеся части У, и V2 так, чтобы сумма весов рёбер, связывающих разные части, была минимальна. Легко показать, что эта задача сводится к минимизации

N

квадратичного бинарного функционала (1) при ограничении ^s¡ =0.

i=i

Алгоритм DDK для этой задачи на первом шаге ищет разбиение графа с матрицей связи Тк, где Т - матрица связи исходного графа, к - некоторая степень (к = 2,3,...,5). Полученное на первом шаге разбиение используется как начальное на втором шаге. Как показывают эксперименты, такой двухступенчатый алгоритм существенно повышает эффективность алгоритма разбиения графа А именно, такая трансформация связей графа понижает

-12-

среднюю энергию находимых минимумов и разница между энергией глобального минимума и средней энергией находимых минимумов уменьшается почти в 3 раза. Как следствие, вероятность отыскания близких к глобальному минимумов повышается на порядки (для матриц Изинга уже при Лг = 100, а также для матриц с равномерным распределением при больших размерностях начиная с N = 200).

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

M = —(]-z)T + —zT2¡ (19)

ат а1Т

где агт - среднеквадратичное отклонение элементов матриц Т1. При изменении параметра г от 0 до 1 происходит переход от матрицы М = Т к матрице М = Т2. Соответственно, поверхность, описываемая функционалом £(5), трансформируется в поверхность Е, (£):

где аи = <yT\j~2 +(1 - г)2 - стандартное отклонение элементов матрицы М.

Как и в главе 3, на основе предложенной трансформации (19), был разработан аналогичный DDK двухэтапный алгоритм, который мы назвали ММ (Mix-Matrix algorithm). Минимизация функционала (20) на первом шаге производится как обычно. Чтобы обосновать работоспособность алгоритма ММ, мы показали, что трансформация приводит к углублению минимумов при небольшом сдвиге.

Углубление минимума. Рассмотрим энергию Ег11 = Ем(S„) в точке S0. После ряда преобразований получаем:

Е,а=Е0 , А^^^Щ^Т), (21)

V(1 -z)2+z2 Г

где q - число отличных от нуля элементов в строке матрицы Т (q = N -\ в SK-модели и q - 4 в модели 2D ЕА), а\ - дисперсия величины Я, = hisj / в минимуме.

Из (21) следует, что имеется оптимальное значение параметра z

v-lfy <22)

при котором углубление достигает максимального значения:

(£г0/£„)„,, =^/í77 (23)

Для матриц 2D ЕА экспериментально найдено, что глубина глобальных минимумов j'»1.31, о, и0.69 и из (21)-(23) находим ~ор, =0.51,

(Ег0/£■„) , »1.45. Аналогично, для SK-модели / ~ 1.51, ov, я0.76, значит

0.55, (Ег0 / Ед)ар1 « 1.59. Видим, что оптимальное значение г находится

вблизи 0.5. Эксперимент подтверждает полученные соотношения (рис. За). Сдвиг минимума. Аналогично (18) получаем вероятность сдвига:

Р=-

1

■\dxi

(г-Г)1

1-Ф

ах + а] -1

—^

где Ф(0 - интеграл вероятностей, а = £(1 - z)\j{q -1) / q + -У ]/г •

(24)

E,JE,

//

а)

SK 500

/.' 2DEA48Î4 Ц

Ч \

Ч1

Рис. 3 Углубление а) и сдвиг Ь) глобальных минимумов в зависимости от параметра г. Маркеры с планками - экспериментальные кривые: вК. 500 - для N = 500 модели Шеррингтона-Киркпатрика; 2Г) ЕА 484 - для матриц 20 модели Эдвардса-Андерсона размерности N = 22 х 22 . Штриховые лиши (слева) - теоретические, построенные по формуле (21) с параметрами у ~ 1.51, ак а0.76 для ЯК 500, и у к 1.31, <у1 »0.69 для 20ЕА484.

Как и следовало ожидать, формула (24) описывает монотонное увеличение величины с1 с ростом 2. Это подтверждается экспериментом (см. рис. 36). Более того, с ростом г возрастают и флуктуации величины сдвига <1. Это ожидаемый результат, поскольку аналогичные вычисления показывают, что дисперсия сдвига описывается выражением а2(Ы) = с1. Значения сдвига, достигаемые при г = ~ор,, приблизительно в 3 раза меньше чем при г = 1: = 0.015ЛГ для ЭК-модели и йор, = 0.0 Ш для модели 20 ЕА. Результаты для алгоритма ММ. Для проверки эффективности предложенного алгоритма в процессе численных экспериментов строилась микс-матрица при значениях г от 0 до 1 с интервалом Дт = 0.05. Результаты экспериментов усреднялись по ансамблю из 100 матриц. Каждый эксперимент включал в себя по Ы^, = 10е стартов. В экспериментах помимо описанной выше микс-матрицы (19) мы также попробовали использовать микс-матрицу типа М = аТ + ЬТЪ. _

На рис. 4 показано как меняется средняя энергия Е найденных минимумов при различных г. Результаты для ЗО Изинга сливаются с результатами для 2Б Изинга и потому не приведены на графике. Самым важным является то, что

независимо от размерности задачи алгоритм ММ позволяет находить минимумы, всего лишь на несколько процентов отличающиеся от энергии глобального минимума. Так, для матриц двумерной и трёхмерной модели ЕА алгоритм SRS (г = 0) даёт минимумы с энергией на 17% хуже глобального, в то время как алгоритм ММ (при г = 0.7) даёт минимумы отличающиеся от энергии глобального минимума лишь на 6.5 - 7%. Аналогично, для матриц SK модели, алгоритм SRS не доходит до глобального минимума примерно на 10%, в то же время средняя глубина минимумов, получаемых алгоритмом ММ при z = 0.7, оказывается лишь на 4% меньше глубины глобального минимума.

Пусть PSRS и Рш - вероятности отыскания минимумов с энергией £є[-1;-0.99] с помощью алгоритмов SRS и ММ, соответственно. Разница между Рш и P,:RS даже при относительно небольшой размерности задачи составляет примерно 3 порядка. При больших размерностях распределение плотности вероятности аналогично представленным на рис. 2.

Сравнение с алгоритмами для задачи MAXCUT. Суть задачи о максимальном разрезе в следующем: для заданного взвешенного неориентированного графа G(V,IV) требуется найти такое разбиение его вершин V на два непересекающихся подмножества Vi и К2, чтобы сумма весов ребёр соединяющих вершины из разных подмножеств была максимальна, т.е. max У V w„„. Легко показать, что эта задача сводится к минимизации

17 ГГ t I I J UV

квадратичного функционала (1), построенного на матрице TtJ = -wtj / 4.

Мы сравнили эффективность алгоритма ММ с алгоритмами, наиболее хорошо

зарекомендовавшими себя для решения задачи MAX-CUT. В частности, это алгоритмы: Scatter Search (Marti et al, 2009) и CirCut (Burer et al, 2000). На сайте (http://www.optsicom.es') выложены наборы задач MAX-CUT и результаты для этих алгоритмов, протестированных на данных наборах на компьютере примерно той же вычислительной мощности что и наш (процессор Intel 3.06 GHz).

Сравнение результатов показано в Таблице 1. В первом столбце показаны названия матриц. Это матрицы того же сорта, что и матрицы 3D ЕА с периодическими

Рис. 4. Среднее значение энергии локальных минимумов ЕI Еа, найденных с помощью алгоритма ММ для матриц моделей ЙК и 20 ЕА Нормировка на глубину глобального минимума Еа введена для исключения зависимости от N, Пунктирные линии - шк-матрица с Т}, сплошные - гшх-матрица с Т1.

(тороидальными) граничными условиями. Первые две задачи с гауссовыми связями Т, , следующие две со связями Ту = ±1. Размерность задач N = Lx Lx L

при L = 15 и 8. Второй и третий столбцы показывают результаты (значение разреза и время в секундах), приведенные на сайте для алгоритмов Scatter Search и CirCut. В последнем столбце показаны усреднённые значения наилучших разрезов, найденные нашим алгоритмом за 10 стартов.

Таблица 1.. Сравнительная таблица для 4-х примеров задачи MAXCUT.

Scatter Search CirCut Mix-Matrix z=0.7, 104 стартов

Problem Obj Val time Obj Val time Obj Val time

toursg3-15 281 029 888 1023 268 519 648 788 (2.73 ± 0.004)*108 8.8

toursg3-8 40 314 704 65 41 684 814 53 (4.10 ± 0.010)*107 2.3

tourspm3-15-50 2 964 333 2 895 427 2 889.7 +4.6 8.9

tourspm3-8-50 442 48 454 38 451.9 ± 1.3 2.2

Как видно из таблицы, наш алгоритм находит решение, сравнимое с одним из алгоритмов (СлгСЩ или 55) и не более чем на 3% хуже другого, при том что СнШ и Э5 затрачивают в 20-100 раз больше времени. Наилучший результат оказался для 4-ой задачи, где наш алгоритм также нашёл в некоторых запусках решение 454, соответствующее лучшему из ББ и С1гСи1. Более того, при увеличении числа стартов, удалось найти решение с величиной разреза 456, которое не было до сих пор получено другими алгоритмами.

В пятой главе подробно рассматривается проблема ограничения на глубину минимумов, вкратце рассмотренной в главе 2. Известно (КгугЬапоувку, 2008), что любая матрица может быть представлена в виде разложения по конфигурациям минимумов (12). При этом конфигурация может быть минимумом, только если её вес гт в разложении не меньше некоторого критического значения гс. К сожалению, теоретически точно оценить значение гс не удаётся. Мы пошли обратным путём и рассмотрели квазихеббовскую матрицу, построенную на случайных паттернах. В этом случае, пользуясь методами статистической физики, удаётся для самого общего вида распределения весов методом седловой точки получить уравнение, позволяющее исследовать свойства энергетической поверхности вблизи паттернов, в том числе получить оценки для ёмкости памяти модели Хопфилда и связанной с ней критическим весом гс.

В качестве модели рассмотрена ситуация, когда только один вес г, = г отличается от всех остальных весов, тождественно равных 1: г^ = 1, ц > 2 и г Ф1. Также исследованы три частных случая различных весов:

1) Для весов в виде гармонического ряда г,.=1/к максимальное число запомненных паттернов невелико но сравнению с размерностью N, всего *м«>/ЛГ/21 пЛГ.

2) Для весов в виде геометрической прогрессии г, = цк было показано, что при некотором оптимальном значении д, достаточно близком к 1 (дт ~ ф - <5^, где <50 «1/(0.165-ЛО), число распознаваемых паттернов будет пропорционально N, кт »0.05Л'.

3) Для весов в виде арифметической прогрессии, при оптимальном значении параметра шага прогрессии, получается максимальное значение загрузки кт «0,067/.

На основании числа кт, можно сделать оценку критического значения веса гс, ниже которого паттерны не распознаются, т.е. не являются минимумами. Так, для весов в виде гармонического ряда гс равно 1 / кт, что в точности совпадает с результатом, полученным методом «сигнал/шум»: Однако максимальное значение гс будет в случае одинаковых весов паттернов (стандартная модель Хопфилда). Таким образом, если статвес паттерна гт больше критического величины: гс—\1 , где ас & 0.14 - критическое значение параметра загрузки стандартной модели Хопфилда, то он точно будет минимумом.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

1) Исследованы свойства энергетической поверхности, описываемой квадратичным функционалом в Л^-мерном пространстве состояний с бинарными переменными. Получены выражения для формы минимумов. Показано, что средняя энергия в л-окрестности минимума описывается формулой Е(и)» Е0(1-2и/ АО2, где Е0 - энергия минимума, при этом разброс относительно среднего стремится к нулю как 1/М. Таким образом, форма минимума определяется лишь глубиной минимума.

2) Показано, что радиус области притяжения минимума х0 полностью

определяется глубиной минимума у = из уравнения:

-А- = ф[у(1-2*0)],

1 ,х0

где Ф(-) - интеграл вероятностей. Поскольку спектр энергий минимумов имеет узкое распределение и глубина минимумов отличается от среднего значения глубины у не более чем на 10-20%, то в интересующей нас области глубоких минимумов можно считать зависимость радиуса от глубины линейной:

*0 в1х+кх(у —у),

где Зс и к - константы (х а 0.384 и к « 0.076 для у = 1.35). Из этого следует, что вероятность отыскания минимума при случайном поиске, равная относительному объёму области притяжения минимума во всём пространстве конфигураций, экспоненциально зависит от глубины минимума Р{у) = х" ! 2" * Ре1"{>-7\

3) Показано, что возведение матрицы Т в степень к (к = 2,3,4....) меняет энергетическую поверхность функционала, однако наиболее глубокие минимумы исходного функционала Е остаются минимумами в новом функционале Ек, либо присутствуют минимумы в их малой ^-окрестности (при к =3 сдвиг с! « 0.026Аг ) и при этом оказываются глубже (для матриц 2Б ЕА при к- 3 происходит углубление минимума, соответствующего глобальному на исходном функционале, примерно на 15%).

4) При сложении матрицы Г и матрицы Т2 с весами (1-г) и г, соответственно, получается новый функционал Е,, который обладает теми же свойствами, что и Ек, только в этом случае можно контролировать сдвиг и углубление минимумов, варьируя параметр г от 0 до 1. Получены выражения для величины углубления минимума и его сдвига в пространстве в зависимости от г. Показано, что существует некоторое оптимальное значение ^«0.51 (для матриц 2В ЕА), при котором происходит максимальное углубление (порядка 45%), при этом сдвиг остается относительно малым с/ч„ = 0.0 ЬУ.

5) На основании двух предложенных выше трансформаций функционала, в диссертации предложены два алгоритма минимизации (ББК и ММ), состоящие

- 18-

из двух шагов. На первом шаге происходит минимизация трансформированного функционала (спуск по поверхности Ек в алгоритме DDK и спуск по поверхности Ег в алгоритме ММ). Далее из найденной точки процедура минимизации продолжается на исходном функционале (продолжение спуска по поверхности Е). Такая двухэтапная процедура алгоритмов позволяет находить минимумы исходного функционала, нивелируя сдвиги, введённые трансформацией поверхности.

6) Показано, что использование предложенных алгоритмов (DDK и ММ) приводит к существенному сдвигу спектра находимых минимумов в сторону глобального минимума и к экспоненциальному по N росту эффективности алгоритма случайного поиска, при этом количество операций, затрачиваемых на один старт, увеличивается вдвое. Дано строгое обоснование эффективности предложенных алгоритмов для полносвязных матриц модели Шеррингтона-Киркпатрика; применение алгоритмов для другого типа матриц является эвристическим. Эффективность предложенных двухэтапных алгоритмов проверена на нескольких типах случайных матриц и продемонстрирована сравнением с известными алгоритмами на задачах MaxCut и Graph В ¡partitioning.

РАБОТЫ АВТОРА ПО ТЕМЕ ДИССЕРТАЦИИ

1. Я.М. Карандашев, Б.В. Крыжановский. Увеличение глубины глобального минимума квадратичного функционала путём возведения матрицы в степень. // XI Всероссийская научно-техническая конференция Нейроинформатика-2009. Москва, МИФИ (2009)

2. J.M. Karandashev, B.V. Kryzhanovsky. An Effective Increase of the Basin of a Quadratic Functional Global Minimum. // ICONS 2009: The 2nd IFAC International Conference on Intelligent Control Systems and Signal Processing. Turkey, Istanbul (2009)

3. Карандашев Я.М., Крыжановский Б.В. Эффективное увеличение области притяжения глобального минимума квадратичного бинарного функционала при нейросетевом поиске. // Международная конференция Актуальные проблемы информационно-компьютерных технологий мехатроники и робототехники (ИКТМР-2009). Дивноморск. (2009)

4. Я.М. Карандашев, Б.В. Крыжановский. Эффективное увеличение области притяжения глобального минимума квадратичного бинарного функционала. // Всероссийская конференция Методы и средства обработки информации (МСО - 2009). Москва, МГУ (2009)

5. Я.М. Карандашев, Б.В. Крыжановский. Эффективное увеличение области притяжения глобального минимума квадратичного бинарного функционала при нейросетевом поиске. //14-ая Всероссийская конференция Математические Методы Распознавания Образов (ММРО-14), Суздаль (2009)

6. Я.М. Карандашев, Б.В. Крыжановский. Улучшение нейросетсвой оптимизации бинарного квадратичного функционала. // 52-я научная конференция МФТИ. Долгопрудный, МФТИ (2009)

7. Я.М. Карандашев, Б.В. Крыжановский. О трансформации энергетической поверхности в задаче бинарной оптимизации. //ДАН. Т. 429, №4, с. 465-469

(2009)

8. Я.М. Карандашев, Б.В. Крыжановский. Эффективное увеличение области притяжения глобального минимума квадратичного бинарного функционала при нейросетевом поиске. // Искусственный интеллект, № 4, с. 37-44 (2009)

9. Я.М. Карандашев, Б.В. Крыжановский. Повышение эффективности алгоритма Кернигана-Лина путём модификации матрицы межсвязей. // Всероссийская конференция Нейроинформатика-2010. Москва, МИФИ

(2010)

10. Ya. М. Karandashev, В. V. Kryzhanovsky. Binary Optimization: Efficient Increasing of Global Minimum Basin of Attraction. // Optical Memory & Neural Networks, Vol. 19, No. 2, pp. 110-125 (2010)

11. Ya. M. Karandashev and В. V. Kryzhanovsky. Efficient Energy Landscape Transformation in the Problem of Binary Minimization. // The 2010 International Joint Conference on Neural Networks (1JCNN), Spain, Barcelona, pp. 1-6 (2010)

12. Yakov Karandashev, Boris Kryzhanovsky. Binary Minimization: Increasing the Attraction Area of the Global Minimum in the Binary Optimization Problem. // ICANN (2) 2010, LNCS, v. 6353/2010, pp. 525-530 (2010).

13. Iakov M. Karandashev, Boris V. Kryzhanovsky. Transformation of Edge Weights in a Graph Bipartitioning Problem. //Proceedings of the 21st international conference on Artificial neural networks ICANN (Part 2), LNCS, v. 6792/2011, pp. 25-31 (2011).

14. Я.М. Карапдашев, Б.В. Крыжановский, JI.Б. Литинский. Сильная неустойчивость спектра минимумов бинарного квадратичного функционала. ДАН. - Т. 436, N 6. С. 733-737 (2011)

15. Iakov Karandashev, Boris Kryzhanovsky and Leonid Litinskii. Weighted Patterns as a Tool to Improve the Hopfield Model. // Physical Review E 85, 041925 (2012).

16. Я.М. Карандашев. Трансформация матрицы связей в задаче разбиения графа.// Вестник компьютерных и информационных технологий. №1, 33-37 (2012)

17.1. Karandashev and В. Kiyzhanovsky . Mix-matrix Method in Problem of Discrete Optimization. // ICCGI 2012 : The Seventh International MultiConference on Computing in the Global Information Technology, Italy, Venice (2012)

18.1. Karandashev and B. Kryzhanovsky. The Mix-Matrix Method in the Problem of Binary Quadratic Optimization. // Proc. of ICANN 2012, Part I, LNCS 7552, pp. 41^48 (2012)

19. Karandashev and B. Kryzhanovsky. Increasing the attraction area of the global minimum in the binary optimization problem. // J. Glob. Optim., V. 56, 3 , pp. 1167-1185 (2013), DOI: 10.1007/s 10898-012-9947-7

Подписано в печать:

08.11.2013

Заказ № 9052 Тираж -100 экз. Печать трафаретная. Типография «11-й ФОРМАТ» ИНН 7726330900 115230, Москва, Варшавское ш., 36 (499) 788-78-56 www.autoreferat.ru

Федеральное государственное бюджетное учреждение науки Научно-Исследовательский Институт Системных Исследований РАН

04201454513 На правах рукописи

Карандашев Яков Михайлович

ИССЛЕДОВАНИЕ МЕТОДОВ ТРАНСФОРМАЦИИ ЭНЕРГЕТИЧЕСКОЙ ПОВЕРХНОСТИ В ЗАДАЧАХ БИНАРНОЙ МИНИМИЗАЦИИ КВАДРАТИЧНОГО ФУНКЦИОНАЛА

05.13.01 — «Системный анализ, управление и обработка информации» (информационно-вычислительное обеспечение)

Диссертация на соискание ученой степени . кандидата физико-математических наук

Научный руководитель: член-корреспондент РАН доктор физико-математических наук Крыжановский Борис Владимирович

Москва —2013

СОДЕРЖАНИЕ

СОДЕРЖАНИЕ........................................................................................................................................2

1. ВВЕДЕНИЕ...........................................................................................................................................3

1.1 Актуальность темы.........................................................................................................................3

1.2 Постановка задачи..........................................................................................................................3

1.3 Обзор литературы.........................-..................................................................................................4

1.4 Цели и задачи диссертации..........................................................................................................11

1.5 Основные положения выносимые на защиту...............................................................................12

1.6 Научная новизна...........................................................................................................................12

1.7 Практическая ценность.................................................................................................................13

1.8 Методы исследования...................................................................................................................13

1.9 Апробация работы и публикации.................................................................................................14

1.10 Структура диссертации...............................................................................................................17

2. ОБЛАСТЬ ПРИТЯЖЕНИЯ МИНИМУМА....................................................................................18

2.1 Область притяжения минимума...................................................................................................18

2.2 Форма минимума..........................................................................................................................23

2.3 Ограничение на глубину минимума.............................................................................................24

2.4 Выводы по главе 2........................................................................................................................25

2.5 Приложение А. Форма энергетической поверхности в п-окрестности.......................................26

3. АЛГОРИТМ С ВОЗВЕДЕНИЕМ МАТРИЦЫ В СТЕПЕНЬ (DDK).............................................29

3.1 Стандартный алгоритм случайного поиска (SRS).......................................................................29

3.2 Базовые соотношения...................................................................................................................32

3.3 Трансформация поверхности........................................................................................................33

3.4 Обоснование алгоритма DDK.......................................................................................................35

3.4а) Углубление минимумов.........................................................................................................35

3.46) Сдвиг минимума...................................................................................................................36

3.4в) Выводы и их проверка...........................................................................................................37

3.5 Эффективность алгоритма DDK..................................................................................................40

3.5а) Результаты для матриц двумерной модели Эдвардса-Андерсона (2D ЕА).......................40

3.56) Результаты для матриц модели Шеррингтона-Киркпатрика (SK)..................................43

3.5в) Улучшение динамики. Алгоритм Хоудайера-Мартина (НМ)..............................................45

3.5г) Результаты для матриц 3-мерной модели Эдвардса-Андерсона (3D ЕА)..........................48

3.6 Применения алгоритма к задаче разбиения графа (graph bipartitioning).....................................49

3.7 Выводы по главе 3........................................................................................................................56

4. АЛГОРИТМ MIX-MATRIX (ММ)...................................................................................................59

4.1 Углубление минимума..................................................................................................................60

4.2 Сдвиг минимума...........................................................................................................................62

4.3 Описание алгоритма ММ.............................................................................................................63

4.4 Результаты для алгоритма ММ.....................................................................................................63

4.5 Сравнение с алгоритмами для задачи MAX CUT.........................................................................70

4.6 Выводы по главе 4........................................................................................................................71

4.7 Приложение А...............................................................................................................................74

5. ОГРАНИЧЕНИЯ НА ВЕСА ПАТТЕРНОВ В КВАЗИХЕББОВСКОЙ МАТРИЦЕ (ДОПОЛНИТЕЛЬНАЯ ГЛАВА).....................................................................................................................75

5.1 Квазихеббовская модель...............................................................................................................75

5.2 Основное уравнение.....................................................................................................................76

5.3 Простейший случай......................................................................................................................79

5.4 Произвольное распределение весов.............................................................................................87

5.5 Выводы по главе 5......................................................................................................................105

6. ЗАКЛЮЧЕНИЕ................................................................................................................................108

СПИСОК ЛИТЕРАТУРЫ..................................................................................................................114

1. ВВЕДЕНИЕ

1.1 Актуальность темы

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

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

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

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

1.2 Постановка задачи

Настоящая работа направлена на повышение эффективности процедуры случайного поиска, применяемой при решении задач бинарной минимизации. Это класс задач, решение которых сводится к минимизации квадратичного функционала Е(Б):

(1-1)

■/V С7Г ,=1 У=1

построенного на заданной Л'" х N -матрице Т в N -мерном конфигурационном пространстве состояний £ = (5,, я2 с дискретными

переменными ^=±1, / = 1, А^; сгт - стандартное отклонение матричных

элементов Ти. Матрицу Т без ограничения общности будем полагать

симметричной (Т =Т.) с нулевой диагональю (Тц =0). Нормировка на от

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

1.3 Обзор литературы

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

называется спиновым стеклом и характеризуется беспорядком и фрустрацией, приводящей к очень сложной энергетической поверхности [25]. Эти физические модели также формулируются как задачи оптимизации, которые применимы во многих областях научной деятельности [3,5]. К сожалению, проблема поиска основного состояния является ЫР-сложной [6]. Поэтому разнообразные алгоритмы адаптировались на протяжении многих лет специально для поиска основных состояний и для исследования спиновых стёкол при конечных температурах и при переходе к пределу нулевой температуры [7,8]. Ввиду большой сложности задач разрабатывались специализированные компьютеры для их решения [9].

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

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

В случае когда система является ферромагнетиком, т.е. взаимодействия Ти имеют все одинаковый положительный знак, задача отыскания основного

состояния является тривиальной, поскольку решением является состояние, в котором все спины сонаправлены (все +1 или все -1). Даже в случае добавления к функционалу (1.1) линейного члена (внешнее случайное поле), задача для ферромагнетика остаётся полиномиально разрешимой при помощи алгоритмов, основанных на методах максимального потока в сети [10,11].

Обычно система становится фрустрирована если существует замкнутый цикл, в котором нечётное число отрицательных связей. В этом случае, спиновое стекло не может быть преобразовано в ферромагнетик [11]. По этой причине быстрые алгоритмы поиска основного состояния не могут быть применены, а проблема становится ЫР-сложной, т.е. для её решения известны только алгоритмы с экспоненциальным временем работы (более

подробно о проблеме Р=ЫР, о классе ЫР-сложных задач, и о полиномиальной сводимости задач друг к другу см. [12-14]).

Как исключение, в случае двумерного спинового стекла без внешнего поля и с периодическим граничным условием не более чем в одном измерении существуют эффективные полиномиальные алгоритмы отыскания основного состояния [15], многие из которых основаны на паросочетаниях в графах [16-18]. Полиномиально разрешим также случай двумерной решётки с периодическими граничными условиями в обоих направлениях, если связи имеют вид ±3 [19]. Для гауссовых связей вопрос до сих пор открыт. Как только мы добавляем в систему внешнее поле, т.е. линейный член в (1.1), задача становится ЫР-сложной для всех типов взаимодействий. Для

трёхмерной кубической решётки задача является NP-сложной даже в отсутствии внешнего поля.

Известно, что задача квадратичной бинарной оптимизации напрямую связана с задачей разбиения графа. В своей базовой версии, эта задача состоит в том, чтобы разбить вершины графа на две непересекающиеся части так, чтобы вес рёбер соединяющих этих части был минимален или максимален. В первом случае задача обозначается как MIN-CUT, а во втором MAX-CUT. Задачи о разрезах имеют множество применений, например в разметке электронных схем, в повышении надёжности сетей, распознавании зрительных сцен и др. Множество важных комбинаторных оптимизационных задач может быть естественным образом сформулировано как задачи квадратичной бинарной оптимизации с ограничениями, причём, применение знаний, полученных из исследования случая без ограничений, часто позволяет существенно повысить скорость конечных алгоритмов [20].

Для неотрицательных весов рёбер, задача MIN-CUT может быть решена при помощи метода потоков в сети благодаря дуальности максимального потока и минимального разреза в сети (алгоритм Форда-Фалкерсона [10]), либо при помощи алгоритма Штора-Вагнера [21].

Для общего случая весов рёбер, задача MAX-CUT (как и задача MINCUT, получающаяся инверсией знаков весов рёбер) является NP-сложной. Для некоторых классов графов, например, как сказано выше, для планарных графов [15], для графов не содержащих К5 (полный граф из пяти вершин)[23], для слабо двудольных графов, в случае целочисленных весов когда род графа ограничен константой и значения весов рёбер ограничены по абсолютному значению полиномом от размера графа [22], а также для некоторых других случаев существуют полиномиальные алгоритмы решения этой задачи. Детальные исследования различных случаев граничных условий и весов рёбер, делающих задачу легко решаемой или NP-сложной, приведены в [24-25].

Строгие алгоритмы. Для решения задачи в общем случае предложено несколько алгоритмов (обзор дан в [26-27]). Простейший метод состоит в перенумерации всех 2N возможных состояний и очевидно имеет экспоненциальное время работы. Даже кубическая решётка состоящая из 4J элементов уже оказывается слишком большой. Базовая идея метода ветвей и границ [28] состоит в том, чтобы исключить из просмотра те области пространства состояний, в которых состояния с низкой энергией не могут быть найдены, и таким образом чтобы все низкоэнергетические состояния системы размера 4J уже можно было посчитать [29]. Также к спиновым стёклам размера 43 применялся метод трансфер-матриц [30].

Для решения NP-сложных задач предложены алгоритмы s-аппроксимации, т.е. полиномиальные алгоритмы, строго дающие решение, отличающееся не более чем на е от энергии глобального минимума. Выдающимся результатом здесь является работа [31], в которой был предложен алгоритм, который за полиномиальное время даёт решение задачи MAX-CUT со значением, составляющим не менее 0.878 от величины максимального разреза, т.е. отличающимся менее чем на 13% от оптимума. С другой стороны, как следует из работы [32], в предположении P^NP, не существует полиномиального алгоритма, который доказуемо даёт решение со значением разреза, составляющим хотя бы 98% от оптимального.

Исследование алгоритмов е -аппроксимации привело к открытию

методов сведения задачи бинарной оптимизации к линейному и

полуопределённому программированию (так называемые LP- и SDP-

релаксации). Такой подход активно используется для поиска верхних и

нижних границ в методе ветвей и границ и позволяет существенно сузить

полный перебор и найти точное решение, пусть не за полиномиальное, но за

умеренное время, для более широкого круга задач. В работе [33] описан

алгоритм, основанный на LP-релаксации, который переписывает

квадратичную функцию энергии в виде линейной функции с

дополнительными набором неравенств, которые должны быть выполнены

7

для всех допустимых решений. Поскольку не все неравенства известны заранее, метод итерационно решает линейную задачу, ищет неравенства, которые нарушены, и добавляет их к набору до тех пор пока решение не будет найдено. Поскольку число неравенств растёт экспоненциально с размером задачи, то же наблюдается и со временем работы алгоритма. Однако данный метод за несколько минут рабочего времени обычного компьютера позволяет находить основное состояние для задач 2D и 3D моделей Эдвардса-Андерсона [34] с периодическими граничными условиями и гауссовыми связями для размерности задачи вплоть до нескольких тысяч спинов. В работах [35-36] описаны методы SDP-релаксации задачи MAXCUT, позволяющие решать задачи с более плотными матрицами связей, чем у 2D и 3D решёток. Однако размер задач, с которыми данный алгоритм справляется, существенно уменьшается с увеличением плотности ненулевых матричных элементов.

Эвристические методы, локальный поиск. Для решения задач больших размерностей, а также с более плотными матрицами (например, модель твёрдого тела Шеррингтона-Киркпатрика [37]) применяются эвристические методы.

Широко распространённым в физике подходом является эволюция системы при некоторой температуре используя метод Монте-Карло (МС) [38]. Для динамики Метрополиса обычно просто выбирают на каждом шаге случайный спин и переворачивают его с вероятностью min(l,exp(-/?AE)), где (5 - обратная температура, АЕ - изменение энергии (1.1) к которому приведёт переворот данного спина. Теоретически можно достичь равновесия при низкой температуре, следовательно основное состояние будет найдено. На этом основаны методы симуляции отжига [39-42]. Тем не менее, для таких систем достижение равновесия в действительности занимает очень долгое время, поскольку динамика застревает в метастабильных состояниях из-за существования многих низкоэнергетических локальных минимумов.

Для борьбы с этим можно использовать алгоритм МС3 (Metropolis-coupled

8

Markov-chain Monte-Carlo) [43], также известный в физике как репличный обмен или параллельная закалка (Parallel Tempering) [44].

Предлагались подходы основанные на нейронных сетях [45]. С разным успехом применялись варианты генетических алгоритмов [46-49]. В работе [1] была предложена комбинация генетического алгоритма с рекурсивной процедурой ренормализации. Основная идея состояла в том, чтобы делить задачу на подзадачи меньшего размера и решать их тем же путём. Данная техника применялась к конечномерным системам с гауссовым распределением весов, с размерами вплоть до 103, а также к другим комбинаторным задачам типа коммивояжёра [50].

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