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

МЕНЮ

Искусственный интеллект
Поиск
Регистрация на сайте
Помощь проекту

ТЕМЫ

Новости ИИРазработка ИИВнедрение ИИРабота разума и сознаниеМодель мозгаРобототехника, БПЛАТрансгуманизмОбработка текстаТеория эволюцииДополненная реальностьЖелезоКиберугрозыНаучный мирИТ индустрияРазработка ПОТеория информацииМатематикаЦифровая экономика

Авторизация

RSS

RSS новости
  1. Хаято Гото *,
  2. Косукэ Тацумура и
  3. Александр Р. Диксон

Абстрактный

Комбинаторные задачи оптимизации встречаются повсеместно, но их трудно решить. Аппаратные средства для решения этих задач в последнее время разрабатываются различными подходами, в том числе и квантовыми компьютерами. Вдохновленный недавно предложенной квантовой адиабатической оптимизацией с использованием нелинейной осцилляторной сети, мы предлагаем новый алгоритм оптимизации, имитирующий адиабатические эволюции классических нелинейных гамильтоновых систем, проявляющих бифуркационные явления, которые мы называем моделируемой бифуркацией (SB). СБ основан на адиабатической и хаотической (эргодической) эволюциях нелинейных гамильтоновых систем. SB также подходит для параллельных вычислений из-за его одновременного обновления. Реализуя SB с помощью программируемой по полю стробоскопической матрицы, мы демонстрируем, что машина SB может получить хорошие приближенные решения все-ко-всем подключенной 2000-узловой задачи MAX-CUT за 0,5 мс, что примерно в 10 раз быстрее, чем современная лазерная машина, называемая когерентной машиной Изинга. СБ позволит ускорить крупномасштабную комбинаторную оптимизацию, использующую цифровые компьютерные технологии, а также предложить новое применение вычислительной и математической физики.

ВВЕДЕНИЕ

В социальной жизни и промышленности часто возникают проблемы с поиском наилучшего сочетания среди огромного числа кандидатов. Эти комбинаторные оптимизационные задачи математически формулируются как минимизация или максимизация некоторых функций дискретных переменных, которые называются целевыми функциями или функциями затрат. Эти проблемы заведомо сложны из-за комбинаторного взрыва ( 1 , 2 ). Таким образом, специальные аппаратные средства для этих проблем, как ожидается, будут полезны. В частности, “Изинг-машины” предназначены для нахождения основного состояния модели спина Изинга ( 3) в последнее время привлекают большое внимание, потому что многие комбинаторные оптимизационные задачи могут быть сопоставлены с проблемой Изинга ( 4), включая очень крупномасштабную конструкцию интегральных схем ( 5), дизайн лекарств ( 6) и управление финансовым портфелем ( 7 ). Эти машины были разработаны различные подходы: квантовые компьютеры на квантовых отжига (8, 9) или квантовой адиабатической оптимизации (10–12) выполнены с сверхпроводящих проводов (13, 14), согласованный Изинга машины (Кисм) выполнены с лазерными импульсами (15–20), а Изинга машин на основе алгоритма имитации отжига (SA) и (1, 21, 22) реализовано с помощью цифровых схем, таких как полевые программируемые вентильные матрицы (FPGA) ( 23-29 ).

Вдохновленный недавно предложенной квантовой адиабатической оптимизацией с использованием нелинейной осцилляторной сети ( 30-34), здесь мы предлагаем новый эвристический алгоритм для задачи Изинга. В этом алгоритме, который мы называем моделируемой бифуркацией (SB ), мы численно моделируем адиабатические эволюции классических нелинейных гамильтоновых систем , проявляющих бифуркационные явления ( 35), где две ветви бифуркации в каждом нелинейном осцилляторе соответствуют двум состояниям каждого Изинговского спина ( 30, 34). [CIMs также используют две ветви бифуркации для двух состояний Изинг-спина ( 15-20), но они не являются Гамильтоновыми системами.] В отличие от моделирования квантовых систем, мы можем эффективно моделировать Гамильтоновы системы ( 36 ) с использованием современных цифровых компьютеров. Это позволяет рассматривать крупномасштабные проблемы с плотной связью, которые являются сложными для ближайших квантовых компьютеров. Кроме того, SB позволяет одновременно обновлять переменные на каждом временном шаге. Это в отличие от SA, поскольку он обычно требует индивидуального обновления, чтобы гарантировать конвергенцию. [Одновременное обновление разрешено только для изолированных спинов (22, 24).] Это различие между SB и SA означает, что ускорение SB при массивно-параллельной обработке легче, чем у SA. Используя эти преимущества, мы реализуем сверхбыструю все-ко-всем подключенную 2000-спиновую Ising машину на основе SB с одним FPGA и демонстрируем, что наша машина примерно в 10 раз быстрее, чем современный CIM ( 17, 19). Кроме того, мы также решаем полностью связанную проблему 100 000-спинового Изинга с непрерывными параметрами SB с помощью кластера графических процессоров (GPU), что примерно в 10 и 1000 раз быстрее, чем наши лучшие SA и программное обеспечение SA, предоставляемые ( 22), соответственно. Это говорит о том, что СБ может ускорить крупномасштабную комбинаторную оптимизацию без специально разработанных аппаратных устройств или пользовательских схем.

В то время как квантовая адиабатическая оптимизация основана на квантовой адиабатической теореме ( 37 , 38 ), СБ основан на адиабатических и хаотических (эргодических) эволюциях классических нелинейных гамильтоновых систем. Мы обсудим этот механизм СБ, используя результаты моделирования простого примера вместо того, чтобы дать математически строгое доказательство принципа действия СБ, который остается интересной задачей для дальнейшей работы.

РЕЗУЛЬТАТЫ

SB алгоритм

Задача N-спинового Изинга без внешних магнитных полей формулируется следующим образом. Безразмерная энергия Изинга задается с помощью

E Ising (s) = - 1 2 ? i = 1 N j j = 1 N J i, j s i s j

(1)где ыя обозначает яче Изинга спин, который принимает значение 1 или -1, х = (х1 х2 ... хП) является векторное представление спиновой конфигурации, и Джейя,J в(= дждж,я) это соединение коэффициент между Я- й и J Вй спинов (J нея,Я = 0). Задача состоит в том, чтобы найти спиновую конфигурацию, минимизирующую энергию Изинга. Эта задача математически эквивалентна известной комбинаторной оптимизационной задаче с именем MAX-CUT ( 5 , 17 , 18 , 39 ): разделите узлы взвешенного графа на две группы, максимизирующие общий вес (называемый “значением среза”) ребер, вырезанных при делении. Установив коэффициенты связи как J i, j = - w i, j ( w i, j = w j, i-вес между i-м и j-м узлами), мы можем решить MAX-CUT с помощью Ising-машин ( см. Также дополнительные материалы) (5 , 17 , 18 ). Поскольку общее число спиновых конфигураций равно 2 N, чрезвычайно трудно найти точные решения для больших N . Если топология графа не имеет некоторой специальной структуры, такой как двумерные решетки, эта задача, как известно, является недетерминированным полиномиальным временем (NP)–hard ( 2, 3). Следовательно, существует твердое убеждение, что для этой задачи не существует эффективного точного алгоритма. На практике часто используют эвристические алгоритмы, такие как SA, чтобы найти приближенные, но практически полезные решения как можно быстрее. Здесь мы предлагаем SB в качестве нового варианта для этих эвристических алгоритмов. [Для получения точных решений малоразмерных задач машина SA, называемая "цифровым Отжигателем" в ( 29), может быть самой быстрой до сих пор.]

Для задачи Изинга недавно был предложен новый вид квантовой адиабатической оптимизации с использованием сети Керр-нелинейных параметрических осцилляторов (KPO) (30 – 34 ). Квантово-механический гамильтониан в этом подходе задается выражением ( 30, 34)

Чщ(Т)=??я=1Н[К2в†2- яа2- я?П(Т)2(а†2я+а2я)+?яв†яэтоя]???0?я=1n и?с J=1НJ ия,Джейв†явДж

(2) где ?-приведенная постоянная Планка, a†i

и этоя - создание и уничтожение операторов, соответственно, Для я- го осциллятора, к - положительное Керр коэффициент, Р(Т) является зависимым от времени параметрических двухфотонной накачки амплитуда, ?- я положительная расстройка частот между резонансной частотой яче осциллятора и частотой накачки, а ?0 - это положительная константа с размерностью частоты. Начальное состояние каждого КПО устанавливается в вакуумное состояние, и амплитуда накачки p ( t ) постепенно увеличивается от нуля до достаточно большого значения. Постоянная ? 0 устанавливается на достаточно малое значение, такое , что состояние вакуума является основным состоянием исходного гамильтониана ( 30, 34 ). Затем каждый КПО окончательно переходит в когерентное состояние с положительной или отрицательной амплитудой через квантовую адиабатическую бифуркацию ( 30, 34), и знак конечной амплитуды для i-го КПО обеспечивает i-й Изинг-спин основного состояния модели Изинга , что гарантируется квантовой адиабатической теоремой ( 30, 34).

Соответствующая классическая Гамильтонова система получена из классической аппроксимации, где математическое ожидание a i аппроксимируется как комплексная амплитуда x i + iy i ( где i обозначает мнимую единицу) (30 , 34 ). Здесь действительная и мнимая части, x i и y i, являются парой канонических сопряженных переменных, таких как положение и импульс, для i-го КПО. Классический механический гамильтониан и уравнения движения для этой классической системы задаются формулами ( 30, 34)

Hc(x,y,t)=?i=1N[K4(x2i+y2i)2?p(t)2(x2i?y2i)+?i2(x2i+y2i)]??02?i=1N?j=1NJi,j(xixj+yiyj)

(3)

x.i=?Hc?yi=[K(x2i+y2i)+p(t)+?i]yi??0?j=1NJi,jyj

(4)

y.i=??Hc?xi=?[K(x2i+y2i)?p(t)+?i]xi+?0?j=1NJi,jxj

(5 )где точки обозначают дифференцирование по времени t. Эмпирически было установлено, что эта классическая система может также находить хорошие приближенные решения со значительной вероятностью , где знак конечного x i обеспечивает i-й спин решений ( 30, 34 ). Этот результат предполагает новый подход к проблеме Изинга. В отличие от квантовых компьютеров, нам не нужно строить реальную машину, описанную выше Гамильтонианом, потому что вместо этого мы можем эффективно моделировать такую машину, используя классические компьютеры.

Однако, Эквалайзер. 4 и 5 не подходят для быстрого численного моделирования. Игнорируя некоторые термины, пропорциональные моменту y , который изменяется около нуля ( 30, 34), мы упрощаем эквалайзер. 4 и 5 следующим образом

HSB(x,y,t)=?i=1N?2y2i+V(x,t)=?i=1N?2y2i+?i=1N[K4x4i+??p(t)2x2i]??02?i=1N?j=1NJi,jxixj

(6)

x.i=?HSB?yi=?yi

(7)

Y.i = - ? H SB ? x i = - [ K x 2 i-p ( t ) + ? ] x i + ? 0 j j = 1 N J i, j x j

(8 ) где V ( x , t) обозначает потенциальную энергию, и все расстройки были приняты за одно и то же значение ?. Результирующая классическая система представляет собой сеть генераторов Дуффинга (35, 40) с массой ? -1 и коэффициентами связи {? 0 J i , j }. [Аналогичное упрощение для моделирования CIM недавно было сообщено в ( 41 ).]

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

Алгоритм SB основан на Эквалайзерах. 7 и 8 . Отметим, что новый гамильтониан в уравнении Eq. 6 является сепарабельным по отношению к позициям и моментам, в отличие от предыдущего в эквалайзере. 3. Следовательно, Эквалайзер. 7 и 8 могут быть решены явным симплектическим методом Эйлера (см. Также методы) (36) вместо неявных методов или, например, явного метода Рунге-Кутты четвертого порядка. Простота явного симплектического метода Эйлера имеет решающее значение для монтажа алгоритма SB с пользовательскими схемами, такими как FPGA. Его стабильность также важна, потому что это позволяет установить шаг времени на большое значение, что приводит к высокоскоростному моделированию. [Симплектический метод Эйлера является точным первого порядка. Явные симплектические методы более высокого порядка, такие как метод Стермера-Верле ( 36), можно также использовать, но это не эффективно не только для скорости и стабильности настоящего моделирования, но и для точности решения задачи Изинга. Следовательно, мы принимаем самый простой метод первого порядка. Дополнительные сведения см. В разделе методы.]

Алгоритм SB выглядит следующим образом. Все переменные, x и y , изначально устанавливаются вокруг нуля. Постепенно увеличивая p (t ) от нуля, мы решаем эквалайзеры. 7 и 8 численно с помощью явного симплектического метода Эйлера. (Для дальнейшего ускорения мы модифицируем метод из его стандартной формы. Дополнительные сведения см. В разделе методы.) Знак конечного x i обеспечивает i-й спин приближенного решения задачи Изинга.

SB имеет дополнительные алгоритмические преимущества: эквалайзер. 7 и 8 имеют только одну операцию суммирования продукта относительно матрицы связи J, которая является наиболее вычислительно-интенсивной частью в этом алгоритме, в отличие от эквалайзеров. 4 и 5 включая две такие операции Product-sum; мы можем обновить переменные одновременно, и поэтому мы можем полностью использовать массовую параллельную обработку для ускорения SB; и SB включает только сложение и умножение, а не вероятностные процессы, и поэтому может быть легко запрограммирован, например, с FPGA.

Механизм SB

Механизм СБ качественно объясняется следующим образом. Когда p ( t ) постепенно увеличивается от нуля до достаточно большого значения p f, каждый осциллятор проявляет бифуркацию с двумя устойчивыми ветвями, положения которых приблизительно задаются путем s i ( p f-? ) / K? ? ? ? ? ? ? ? ? ? ?

(s i = ±1) (34). Соответственно, потенциальная энергия в конечном итоге имеет 2 n минимумов, соответствующих 2 n спиновых конфигураций, которые приблизительно задаются с помощью

В(С(пф??)/к???????????,т)=?п(пф??)24к??02пф??K и?я=1н?с J=1пJ ия,J ВыясДж

(9)

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

Вместо того, чтобы показать общую справедливость приведенного выше качественного описания, здесь мы показываем ситуацию простого примера: модель двухспинового Изинга с ферромагнитной связью ( J 1,2 = J 2,1 = 1), основными состояниями которой являются (вверх, вверх) и (вниз, вниз). Это даст нам интуитивное понимание SB, которого достаточно для использования SB в качестве эвристического алгоритма.

На рис .1 (А и Б) приведены результаты моделирования для этой задачи с линейно увеличенной амплитудой накачки p ( t ) = ? p t. Как показано на этом рисунке, после двух бифуркаций потенциальная энергия в конечном итоге имеет четыре минимума, соответствующие четырем конфигурациям двухспиновой модели Изинга.

Инжир. 1 динамика в СБ с двумя осцилляторами для двухспиновой задачи Изинга с ферромагнитной связью.

А) зависящая от времени амплитуда накачки, p (t ) = ? p t (? p = 0,01), при моделировании (сверху) и временных эволюциях x и y (осциллирующие тонкие линии в средней и нижней панелях). Жирные линии представляют собой x стабильных неподвижных точек (минимумов потенциальной энергии): x 1 = x 2 = 0 ( p ? ?-? 0 = 0.4), x 1 = x 2 = ± (p-? + ? 0 ) / K? ? ? ? ? ? ? ? ? ? ? ? ? ?

(p > ?-?> 0 = 0,4), и x 1 = - x 2 = ± (p-?-? 0 ) / K? ? ? ? ? ? ? ? ? ? ? ? ? ?

(p > ? + ?0 = 0.6). (Б) траектории (циркуляционных линий в белом) в трех временных интервалах, (0, 40) (Топ), (40, 60) (в центре), и (60, 200) (внизу) и потенциальной энергии в(Х, Т) измеряется от полной энергии е(т) = н(х(т), г(т), т) в заключительные времена интервалов. Границы временных интервалов обозначены в (А) пунктирными ( t = 40) и пунктирными ( t = 60) вертикальными линиями. Петли (пурпурные) показывают границы энергетически допустимых областей, определяемых V ( x , t ) ? E (t ) ? 0 в конечные моменты времени. Другие параметры задаются как K = 1, ? = 0.5, ? 0 = 0.1, x 1 (0) = x 2 (0) = y 1(0) = 0, а y 2 (0) = 0,1. C и D) аналогичные результаты для p (t) приведены на верхней панели в подпункте C). Границы временных интервалов обозначены в (C) пунктирными ( t = 100) и пунктирными ( t = 131.2) вертикальными линиями.

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

Эта адиабатическая эволюция объясняется следующим образом. Во-первых, полная энергия E (t) и потенциальные энергии на минимумах, V ±,± ( t) (+и ? соответствуют вверх и вниз соответственно), формулируются из их временных производных как

E (t ) = E 0-? p 2 ? t 0 / x (t ' ) | 2 d t '

(10)

V ±, ± (t) = - ? p 2 ? t 0 | X ±, ± (t ' ) | 2 d t '

(11 ) где E 0-начальная полная энергия, которая считается малой, и X ±,± ( t) обозначают положения минимумов. Перед соответствующей бифуркацией X ±, ± (t) определяются как ноль. Заметим, что полная энергия монотонно уменьшается из-за второго члена в уравнении Eq. 10 Описание отрицательной работы, выполненной зависимым от времени термином накачки. После первой бифуркации X +, + (t ) и X?, ? (t ) покидают начало координат, и x ( t ) начинает циркулировать вокруг них, вместо начала координат. Это приводит к увеличению | x ( t)| и поэтому дальнейшее уменьшение полной энергии. В результате энергетически допустимая область, определяемая V ( x , t ) ? E (t ) ? 0 (внутри петель на фиг. 1Б) разделяется на две части, как показано на средней панели фиг. 1B. В настоящем моделировании x ( t ) ограничивается допустимой областью вокруг X?, ? (t ). После этого x ( t ) обращается вокруг X ?,? ( t ), что приводит к ?| x ( t )| 2 dt ? ? | x ?,? ( T )| 2 dt . Отсюда и эквалайзер. 10 и 11, разница между E ( t ) и V?, ? (t ) сохраняется постоянной при небольшом значении. Таким образом, x (t ) сохраняется в допустимой области вокруг X?, ? (t), как показано на нижней панели рис. 1B. Таким образом, SB находит основное состояние, адиабатически следуя одной из стабильных неподвижных точек (минимумов потенциальной энергии), соответствующих основным состояниям, которые появляются в первой точке бифуркации.

Вышеизложенный механизм SB предполагает, что решения, полученные с помощью SB, определяются в первой точке бифуркации. Поскольку x (t ) в первой точке бифуркации задается максимумом-собственным значением собственного вектора матрицы связи J (см. Дополнительные материалы), SB может обеспечить признаки элементов собственного вектора в виде приближенного решения, которое фактически является приближенным решением, полученным методом непрерывной редукции ( 30, 34). Однако, как показано в следующем подразделе, SB может обеспечить гораздо лучшие решения, чем это. Этот факт подразумевает, что может существовать еще один механизм СБ; то есть СБ может продолжать поиск лучших решений после первой бифуркации.

Чтобы исследовать этот новый механизм SB, мы также провели еще одно моделирование, где p (t ) приведено в верхней панели фиг. 1С . Результаты представлены на Рис.1. 1 (C и D). В этом случае начальное значение p (t ) больше двух точек бифуркации, и поэтому даже в начальный момент времени существует четыре минимума потенциальной энергии. В то время как p ( t ) является постоянной величиной, x ( t ) движется в энергетически допустимой области, включая четыре минимума, как показано на верхней панели фиг. 1D. После этого допустимая область разделяется на четыре части, как показано на средней панели фиг. 1D. Тогда, x ( t) ограничена в одном из четырех регионов. В конечном счете, SB находит основное состояние (вверх, вверх) успешно.

Так как x (t ) также движется вокруг других минимумов потенциальной энергии, как показано на верхней и средней панелях Рис.1. 1D, вероятность отказа может быть конечной. Чтобы оценить вероятность успеха, мы повторяем это моделирование 10 4 раза с различными начальными условиями, где x 1 (0) = x 2 (0) = 0 и y 1 (0) и y 2(0) задаются случайным образом из интервала (-0.1, 0.1). В результате мы получаем основные состояния с вероятностью около 0,83. Этот результат указывает на то, что более низкие энергетические состояния (лучшие решения) получены с более высокими вероятностями даже в начальном состоянии, когда энергетически допустимая область включает несколько минимумов потенциальной энергии.

Этот интересный результат можно объяснить, предположив эргодичность этой нелинейной Гамильтоновой системы, подразумевая, что траектория в конечном итоге посетит все состояния в допустимой области фазового пространства ( 36 ). Об этом свидетельствует верхняя панель на фиг. 1D, где x (t ) движется вокруг двух низких минимумов чаще, чем два высоких минимума. Из эргодичности следует, что вероятность успеха может быть пропорциональна объему фазового пространства допустимой области, удовлетворяющей x 1 x 2 > 0. Вероятность успеха, численно оцененная по объему фазового пространства для средней панели фиг. 1D это примерно 0.841, что удивительно близко к фактической вероятности успеха 0.83. Этот численный факт поддерживает механизм, основанный на эргодичности. В общем случае эргодичность предполагает, что x ( t) движется вокруг минимумов низкой потенциальной энергии за более длительное время, чем высокие минимумы, что приводит к более низким энергиям Изинга. В следующем подразделе мы приводим дополнительные численные доказательства этого механизма; то есть SB может найти гораздо лучшие решения для задачи спинового Изинга 2000, чем те, которые были получены при первой бифуркации. Таким образом, мы приходим к выводу, что SB находит более низкие минимумы среди многих минимумов, использующих эргодичность. [Совсем недавно был предложен новый подход к задаче Изинга с использованием классической Гамильтоновой системы с” классическими спинами", имитирующими квантовый отжиг ( 42), где фиксированная точка отслеживается методом, называемым сокращением до адиабатичности. Хотя этот подход также использует классическую Гамильтонову динамику, он сильно отличается от SB, потому что он не будет использовать эргодический поиск, а просто отслеживать неподвижную точку. Кроме того, до сих пор было неясно, может ли он достичь высокоскоростных вычислений со своими сложными уравнениями.]

Решение проблемы максимального сокращения всех подключенных 2000 узлов с помощью одной FPGA SB-машины

Исходя из преимуществ и свойств SB, мы ожидаем, что SB будет полезен для приближенного решения крупномасштабных проблем Ising или MAX-CUT с плотной связью как можно быстрее. Чтобы проверить это ожидание, мы сравниваем нашу машину SB с современным CIM ( 17 , 19 ), потому что эта машина достигла выдающейся производительности для этих проблем. Для этого сравнения мы решили все-ко-всем подключенную 2000-узловую задачу MAX-CUT с именем K 2000 , которая была решена CIM ( 17 , 19 ).

Мы внедрили SB-based 2000-spin Ising machine с одним FPGA. Операция product-sum, ? N j = 1 J i, j x j

, в эквалайзере. 8, которая является наиболее вычислительно-интенсивной частью в SB, эффективно выполняется путем массивной параллельной обработки. J i , j x j слагаемые до 8192 обрабатываются за один такт, что примерно в четыре раза больше суммарных спинов. Смотрите дополнительные материалы для получения дополнительной информации.

Одним из главных результатов в (17) является то, что ММК достиг цели энергия дается Goemans-Уильямсон полуопределенном программировании (ГВт-СДП) алгоритм (17, 38) около 70 МКС в лучшем случае из 100 проб, что составляет около 50 раз короче, чем на СА очень настроены в (17). [SA в ( 17, 19 ) очень быстро, поставив все элементы матрицы сопряжения J на кэш одного ядра. Таким образом, параллельные вычисления с несколькими ядрами не могут ускорить SA из-за накладных расходов на связь.] По сравнению с этим, наша машина SB одиночн-FPGA достигает значение GW-SDP в только 58 µs даже в самом плохом случае среди 100 проб, как показано в Fig. 2А .

Инжир. 2 производительность машины SB, реализованной с помощью FPGA для K 2000 .

А) временные эволюции энергий Изинга. (Время вычислений не включает загрузку данных и вывод результатов; см. Дополнительные материалы.) Константы в SB задаются как K = ? = 1 и ?0=0.7??N?

(? = 1-это SD элементов из J .) p (t ) линейно увеличивается от 0 до 1, в то время как число временных шагов увеличивается от 1 до N шагов . Жирные линии: 100-пробные средние энергии SB с шагом N = 40 (слева), CIM из ( 17, 19) (середина) и SA из ( 19) (справа). [Данные CIM являются 26-пробным средним, исключающим 74 испытания из-за нестабильности оптических параметрических осцилляторов ( 17 , 19 ).] Пунктирные тонкие и пунктирные тонкие линии, расположенные между жирными линиями: нижняя (в лучшем случае) и верхняя (в худшем случае) огибающие 100 трасс. Пунктирная тонкая горизонтальная линия: GW-SDP ( 17). Пунктирная тонкая горизонтальная линия: 100-пробная средняя энергия HNN после сближения при локальном минимуме. Таблица: расчетные времена для достижения энергий HNN и GW-SDP. N / A, не применимо. (B ) гистограммы 100 значений среза, полученных с помощью SB с шагом N = 186 через 0,5 мс ( сверху), CIM через 5 мс из (17) (середина) и SA через 50 мс из (17) (снизу). Цифры в панелях представляют собой средние значения сокращений по 100 испытаниям.

Другой основной результат в ( 17) заключается в том, что CIM потребовалось всего 5 мс для получения хороших решений (больших значений разреза), как сильно настроенный SA, полученный за 50 мс, где значение разреза задачи MAX-CUT задается с помощью ? E Ising 2-1 4 ? N i = 1 ? N j = 1 J i, j

с параметрами Ising (см. Дополнительные материалы). Рисунок 2B показывает, что наша машина SB занимает всего 0,5 мс для получения лучших решений в среднем, чем CIM в 5 мс и SA в 50 мс.

На Фиг. 2A, мы также показываем энергию, полученную нейронной сетью Хопфилда ( HNN) двух состояний нейронов (19 , 43 ). HNN является самым наивным алгоритмом локального поиска, и поэтому более сложные алгоритмы должны превосходить HNN. (HNN эквивалентен SA при нулевой температуре; см. Дополнительные материалы.) Обратите внимание, что энергия HNN ниже, чем энергия GW-SDP. Наша машина SB достигает энергию HNN около 13 времен более быстро чем CIM на среднем. Из этих результатов мы делаем вывод, что наша одно-FPGA SB машина примерно в 10 раз быстрее, чем CIM, разработанная в ( 17 ).

Как уже упоминалось в предыдущем подразделе, в первой точке бифуркации SB может найти приближенное решение, определяемое максимумом собственного значения собственного вектора матрицы связи J ( 30, 34). В случае K 2000 это приближенное решение равно E Ising = -55724, или отрезанное значение 27 342. Это намного хуже, чем те, которые фактически получены SB, как показано на фиг. 2. Этот факт поддерживает механизм СБ, основанный на эргодичности, рассмотренной в предыдущем подразделе.

Решение все-ко-все подключенные 100,000-узел MAX-CUT проблемы с помощью GPU-кластера SB машины

Далее мы решили неразрывно связанную задачу максимального сокращения на 100 000 узлов с непрерывными весами, заданными случайным образом из интервала [-1, 1], который называется M 100000 . [См. дополнительные материалы для его определения с помощью генератора псевдослучайных чисел в ( 44).] Обратите внимание, что в M 100000, существует около 5 миллиардов ребер, размер данных которых слишком велик для ПЛИС. Вместо этого мы используем две большие вычислительные системы и сравниваем производительность каждой из них. Система 1-это кластер графических процессоров, использующий восемь графических процессоров, взаимодействующих друг с другом по быстрой ссылке. Система 2 представляет собой кластер ПК, состоящий из нескольких узлов с двумя 14-ядерными процессорами, каждый из которых взаимодействует друг с другом через быструю связь. Мы также решили M 100000 с помощью hnn и SA, разработанных нами на основе ( 22). Наш SA является удивительно быстрым, например, примерно в 100 раз быстрее, чем программное обеспечение, предоставляемое ( 22). Смотрите дополнительные материалы для получения дополнительной информации.

Сравнение времен вычислений для SB и SA показано на фиг. 3А . Самая быстрая реализация-это кластер GPU для SB и 50-ядерный кластер ПК для SA. Это отличие связано с тем, что массовая параллельная обработка кластера GPU эффективна для SB, но не для SA, что обусловлено отсутствием распараллеливаемости спинового обновления и вытекающих из этого больших коммуникационных накладных расходов. Рисунок 3B показывает время вычислений SB и SA в их самых быстрых реализациях. Пунктирные линии примерно дают огибающие кривых эволюции времени. Это означает, что пунктирные линии указывают наилучшие случаи относительно N шага и N развертки . Из этого сравнения при тех же энергиях Изинга ниже значения HNN получается, что машина GPU-кластера SB примерно в 10 раз быстрее нашей самой быстрой SA или в 1000 раз быстрее программного обеспечения, предусмотренного ( 22). Это показывает, что SB позволяет ускорить крупномасштабную комбинаторную оптимизацию без специально разработанных аппаратных устройств или пользовательских схем.

Инжир. 3 производительность машин SB реализована с кластером GPU и кластером ПК для M 100000 .

А) время вычисления SB (N шаг = 100) и SA ( N шаг = 100) для M 100000 . (Время вычислений не включает загрузку данных и вывод результатов; см. Дополнительные материалы.) Константы в SB задаются как K = ? = 1 и ?0=0.7??N?

(?=1/3–?

является SD элементов из J .) p (t ) линейно увеличивается от 0 до 1, в то время как число временных шагов увеличивается от 1 до N шагов . В СА скорость отжига регулируется общим числом "стреловидностей" N стреловидности , а обратная температура линейно увеличивается от 0 до 0,05. Полный и пустой круги: SB и SA, соответственно, с использованием кластера ПК, число ядер которого задается горизонтальной осью (каждый узел кластера ПК использует 1 ядро или 25 ядер). Пунктирные и пунктирные линии: SB и SA, соответственно, используя кластер GPU с восемью графическими процессорами. (B) Временные эволюции энергий Изинга (тонкие линии). Полные и пустые окружности, Соединенные пунктирными линиями: конечные значения SB и SA, соответственно, с различным шагом N и N разверткой (оба являются 10, 20, 50, 100, 200, 500, и 1000 слева). Пунктирная горизонтальная линия: 10-пробное среднее значение HNN после сведения к локальному минимуму.

ОБСУЖДЕНИЕ

Мы предложили новый эвристический алгоритм для задачи Изинга, которую мы называем SB, вдохновленный квантовой адиабатической оптимизацией с использованием нелинейной осцилляторной сети. Используя преимущества SB, такие как одновременное обновление и простые вычисления, мы реализовали сверхбыстрый все-ко-всем подключенный 2000-спиновый Изинг-автомат на основе SB с использованием одного FPGA, что примерно в 10 раз быстрее, чем современный CIM. Мы также решили все-ко-все подключенные 100,000-узловые задачи максимального сокращения с непрерывными весами, используя машину SB GPU-кластера, которая примерно в 10 раз быстрее, чем наша самая быстрая SA. Мы обнаружили, что скорость вычислений машины GPU-cluster SB ограничена полосой пропускания памяти каждого GPU (см. Дополнительные материалы). Это указывает на то, что аппаратные устройства с более широкой полосой пропускания памяти, специально разработанные для SB, такие как специализированные для приложений интегральные схемы или мульти-FPGA системы, значительно улучшат производительность. Механизм СБ основан на адиабатической и хаотической (эргодической) эволюциях нелинейной Гамильтоновой системы. Мы обсудили механизм, используя результаты моделирования простой модели. Механизм СБ, по-видимому, связан с адиабатическими инвариантами классических интегрируемых систем ( 36, 40, 45, 46 ) , а также с эргодическими адиабатическими инвариантами классических неинтегрируемых систем (47-49). Математически строгая обработка механизма с использованием этих понятий оставлена для дальнейшей работы. Другим интересным направлением исследований является распространение SB на термические состояния при конечной температуре, например, с использованием метода носе-Гувера ( 36 ), который будет связывать SB с неравновесной статистической механикой ( 50 ).

МЕТОДЫ

Модифицированный явный симплектический метод Эйлера

Симплектический метод ( 36)является стабильным численным методом для решения гамильтоновых уравнений движения. В дальнейшем отображение потока по Гамильтониану H с шагом по времени ? t обозначается через ? (? t, H ). Поскольку карта течения только по кинетической энергии ?(? t , H ? V ) или только по потенциальной энергии ?(? t , V), очевидно , симплектична , их состав ?(? t, V ) ? ?(? t, H ? V ) также симплектичен ( 36 ). Поскольку настоящий гамильтониан сепарабелен относительно x и y , мы получили явный симплектический метод Эйлера для настоящего гамильтониана

x i (t n + 1 ) = x i (t n ) + ? y i (t n ) ?t

(12)

yi(tn+1)=yi(tn)?[Kx3i(tn+1)+(??p)xi(tn+1)??0?j=1NJi,jxj(tn+1)]?t

(13), где t n = ? t n-n-е дискретизированное время и p = p ( t n + 1). Суть этого метода заключается в том, что моменты обновляются с использованием обновленных позиций. Этот метод не только прост, но и стабилен для значительно больших ?t . Таким образом, мы можем быстро решить уравнения Гамильтона, используя большой ?t .

Для реализации дальнейшего ускорения мы модифицировали этот метод следующим образом. В алгоритме SB операция product-sum, ?Nj=1Ji,jxj

, является наиболее вычислительно-интенсивной частью. Чтобы смягчить это вычисление, мы разделяем гамильтониан как HSB=MHSB?HJM+HJ, где M-целое число, большее 2 и HJ=??02?Ni=1?Nj=1Ji,jxixj. Поскольку карты течения ? (? t, (H SB ? H J )/ M ) = ?(? t / M , H SB ? H J ) и ? (? t, H J), очевидно, симплектические, их состав также симплектичен. Таким образом, мы получили модифицированный явный симплектический метод Эйлера для SB

x (m + 1) i = x (m ) i + ? y (m ) i ? t

(14)

y ( m + 1) i = y (m ) i - [ K x (m + 1) 3 i + (?-p) x ( m + 1) i] ?t

(15)

xi(tn+1)=x(M)i

(16)

yi(tn+1)=y(M)i+?0?j=1NJi,jx(M)j?t

(17) где ? t = ? t / M, m = 0,..., M ? 1, x(0)i=xi(tn), y(0)i=yi(tn), и p = p ( t n +1 ), игнорируя его детальную временную зависимость, потому что p (t ) изменяется медленно. Заметим, что карта течения ?(? t , H SB ? H J ), которая включает нелинейные члены и поэтому индуцирует нестабильность, заменяется M-временным повторением ? (? t , H SB ? H J ) с меньшим временным шагом ? t = ? t / M . Таким образом, данный метод является стабильным по сравнению с исходным. Мы обнаружили, что исходный и модифицированный методы являются стабильными при ? t ? 0,5 и ? t ? 1 (для больших M ), соответственно, где K = ? = 1 и ?0=0.7??N?

, и p (t ) увеличивается линейно от 0 до 1, как и в настоящем моделировании. Таким образом, используя модифицированный метод, мы можем установить ? t на величину, большую в 2 раза, чем исходный метод, и, следовательно, ускорение в 2 раза может быть достигнуто. В настоящем моделировании мы установили ? t = 0.9 и M = 2 на фиг. 2 и ? t = 0.9 и M = 5 на фиг. 3.

Симплектический метод второго порядка, называемый методом Штермера-Верле, задается картой потока ? (? t /2 , H SB ? V ) ? ?(? t, V ) ? ?(? t /2, H SB ? V ) ( 36 ). Из-за ?(? t , V) устойчивость этого метода эквивалентна явному симплектическому методу Эйлера. Мы обнаружили, что этот метод стабилен для M 100000, когда ?t ? 0.5. Таким образом, этот метод не способствует ускорению. Мы также обнаружили, что точность второго порядка для гамильтоновых уравнений не способствует точности решения задачи Изинга. Это может быть связано с тем, что нам нужны только знаки позиций, которые будут устойчивы к численным ошибкам. Поэтому мы приняли приведенный выше модифицированный явный симплектический метод Эйлера.


Источник: advances.sciencemag.org

Комментарии: