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

E= –(, )+2×(, ),                                      (1)

в котором компоненты вектора  принимают значения ±1. Матрица  симметричная, с нулевой диагональю, а ее элементы – случайные независимые величины. Количество минимумов функционала E() велико и экспоненциально растет с ростом размерности задачи. Так, уже при размерности вектора , равной 200, число минимумов ~100 000, и поэтому нахождение достаточно глубокого минимума сложная задача.

Алгоритм минимизации. За основу процедуры минимизации взята нейросетевая модель Хопфилда [1]. Это полносвязная рекуррентная нейронная сеть из N нейронов, имеющих два состояния si=±1, . Энергия сети задана выражением (1). Такую сеть можно рассматривать как систему, решающую задачу бинарной минимизации: конвергируя в устойчивое состояние, сеть находит конфигурацию, соответствующую минимуму энергии E. Показано, что при случайном поиске вероятность отыскания какого-либо минимума экспоненциально растет с увеличением глубины этого минимума [2, 3]. Это означает, что нейросеть с подавляющей вероятностью находит если не оптимальное решение (глобальный минимум), то одно из субоптимальных (локальный минимум).

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

,                                                    (2)

и компоненте si присваивается значение

.                                                           (3)

Эта процедура последовательно применяется ко всем компонентам  до тех пор, пока сеть не конвергирует в устойчивое состояние, соответствующее минимуму энергии.

Для уменьшения объема вычислений в [4] предложен метод дискретизации матрицы , то есть выделение из матричных элементов среднего значения A0 и замена остатка матрицей , нормированные элементы которой имеют целочисленные значения в диапазоне [–m; +m], где m – число градаций. При таком подходе поиск минимума E() заменяется поиском минимума дискретизованного функционала , если A0=0.

Оптимизация функционала (4) проводится аналогично оптимизации функционала (1). На каждом такте работы вычисляется одна из компонент локального поля

= –+×,                                                        (4)

и ей присваивается значение

.                                                            (5)

Качество такой аппроксимации будет зависеть от выбранного числа градаций и определяется функцией распределения. С увеличением числа градаций максимум функции распределения смещается влево (рис. 1), в сторону более глубоких минимумов, что увеличивает вероятность их нахождения. Так, для энергии dE=0,05 имеем: при переходе от m=1 к m=8 плотность состояний увеличится в 1,5 раза. На рисунке 1 представлены функции распределения плотности вероятности по энергии функционала (1), построенные на минимумах e() для разного числа градаций. По оси X отложена энергия, определяемая величиной dE=(E0–E)/E0, где E0 – энергия глобального минимума.

Минимум функционала. В данной работе применяется двухэтапный алгоритм [5]. Процесс минимизации функционала (1) начинается с некоторой случайной точки пространства . Подчиняясь решающему правилу (5), сеть конвергирует в некоторое устойчивое состояние , являющееся минимумом функционала e. Если из этой точки продолжить спуск с решающим правилом (3), то сеть конвергирует в состояние o, соответствующее минимуму функционала (1). На всем пути ®®o величина ошибки (вероятность ошибки P в направлениях градиентов) только уменьшается. С ростом числа градаций m величина ошибки убывает как .

Эффективность минимизации определяется величиной , где  – энергия минимума  дискретизированного функционала;  – энергия ближайшего к нему минимума o, в который можно было спуститься при использовании исходного решающего правила (3). Показано [4], что эффективность минимизации dE»P,                                                         (6)

из которой следует, что разница в энергиях с ростом m убывает как dE~m–1.

При этом хеммингово расстояние между минимумами функционала E и его аналога e определяется простым соотношением d=NP.

Хеммингово расстояние не превышает значения d=0,11N, когда m=1, и снижается до значения d=0,02N при m=16.

Сказанное подкрепляется результатами эксперимента (см. рис. 2). На рисунке кривые 1, 2, 4 – функции распределения начального состояния , дискретизированного  и конечного o. Кривая 3 – распределение состояний исходного функционала. На рисунке они же отмечают точки максимумов функций плотности. Расстояние между этими максимумами на графике соответствует в точности формуле (6). Видно, что первый участок пути ® значительно больше второго участка ®o. Преимущество двухэтапного алгоритма в том, что первая часть пути ® проходится с большей скоростью, вследствие чего достигается общее увеличение скорости работы алгоритма. Вероятность нахождения глубоких минимумов двухэтапным алгоритмом выше, чем простым методом.

Работа с укороченной арифметикой

Полученные результаты дают основу для применения метода дискретизации, что означает возможность применения чисел малой разрядности. При этом для хранения матричного элемента при выборе параметра m=1 достаточно 2 битов, а для m=15 – половины байта. В дальнейшем размер чисел матрицы  будем называть исходным форматом, а p обозначим количество чисел малой разрядности, которое упаковывается в исходный 4-байтный формат.

Основной объем вычислений по алгоритму Хопфилда приходится на вычисление градиента, то есть на вычисления скалярных произведений двух векторов. Поэтому возможность такого представления элементов  позволяет увеличить быстродействие алгоритма. Например, сложение пары вещественных чисел в 4-байтном формате требует одного такта процессорного времени. За это же время можно сложить 4 пары целых чисел в байтном формате.

Ускорение алгоритма определяется величиной θ=I0/I, где I0 – исходное время работы алгоритма; I – время работы алгоритма при упакованных элементах матрицы.

На рисунке 3 приведен пример увеличения скорости работы алгоритма (m=1, параметр упаковки p=8). По оси ординат отложено ускорение q, а по оси абсцисс – размерность нейросети N. С увеличением размерности нейросети ускорение алгоритма увеличивается и достигает своего предельного значения q»7,3.

На скорость работы сети также влияет исходный формат чисел. Например, если элементы исходной матрицы  вещественные (4 байта), то время работы возрастет до 1,4, а для вещественных 10-байтных чисел – до 6,3.

Алгоритм расчета нейронной сети Хопфилда основан на последовательном вычислении компонент градиента  при каждом изменении состояния нейрона.

На рисунке 4 показано уменьшение числа нейронов с каждой итерацией (одна итерация соответствует проверке всех N нейронов). На вертикальной оси отложена доля измененных нейронов  n в процентах. Как видно из рисунка, с увеличением числа шагов число нейронов, изменяющих свое направление, становится малым. Это означает, что направление компонент также редко изменяется. Поэтому вычисление на каждом шаге компоненты вектора  неэффективно. Вычислительная сложность алгоритма ~2TN2, где T – число итераций. С увеличением размерности сети число итераций T растет (рис. 5а). Например, при размерности N=100 число шагов T=5, а при N=10 000 число шагов T=80.

В настоящее время используется другой метод расчета. В исходном состоянии  вычисляются все компоненты . На каждом шаге процедуры при изменении состояния нейрона вектор  модифицируется по правилу =±2()i, если направление спина положительно/отрицательно; ()i – вектор строки с номером текущего нейрона. Объем вычислений компонент  в исходном состоянии равен N2. При перевороте каждого спина производится N операций. Число переворотов спинов не превышает N (рис. 5б), поэтому вычислительная сложность такого алгоритма ~(1+a)N2, a<1. Отношение объемов вычислений такого метода к исходному методу ~T, (T>>1). Далее в качестве исходного алгоритма используется метод с обновлением компонент .

Для модификации алгоритма обновления при использовании чисел малой разрядности необходима их упаковка в исходный формат. В результате параллельно выполняются операции сложения и умножения, что приводит к увеличению скорости работы. На рисунке 6 показано ускорение алгоритма при упаковке (p=8, m=1) для нейронных сетей различной размерности. Установлено, что достигаемое ускорение q»5. При упаковке p=4 в исходный формат ускорение составляет q»3,8 (m<15), а при p=16 (m=1) – q»8.

Алгоритм случайного поиска

Проведем анализ быстродействия двухэтапного алгоритма при использовании упакованных данных (m=1, p=8). на первом этапе. В этом случае .

На первом этапе число итераций при минимизации e() равно числу итераций исходным методом [5], поэтому объем вычислений на первом этапе составляет O1=OH/5. На втором этапе число итераций T мало [5], соответственно мало и число изменений в состояниях нейронов, поэтому основной объем вычислений приходится на вычисление градиента  в промежуточных состояниях и составляет O2=OH/2. В этом случае ускорение алгоритма составит q»1,4.

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

Конечно, при таком подходе велика вероятность получить несколько худшее решение, оценить которую можно экспериментально (рис. 7). Для этого распределение минимумов  разбивалось на левую и правую области относительно пика распределения и определялась вероятность найти самый глубокий минимум при стартах из каждой области. На рисунке видно, что самый глубокий минимум с вероятностью 96 % находится из состояний левой области, в то время как правая область дает вероятность, близкую к нулю. Поэтому состояния правой области могут не участвовать на втором этапе поиска.

Чтобы определить границы области, состояния  которых используются на втором этапе ал- горитма, и тем самым уменьшить объем вычис- лений, определена интегральная эффективность минимизации. То есть определялся минимум величины  на отрезке, начало которого соответствует наиболее глубокому минимуму дискретизированной сети. С помощью dE можно оценить хеммингово расстояние D между глобальным минимумом  и самым глубоким минимумом E, найденным при стартах из отрезка. Действительно, при изменении состояния одного нейрона энергия меняется в среднем на величину e. Поскольку все минимумы при N>100 находятся в области Ec±3s и отношение se/Ec<0,01, то Ec~eN. Тогда .

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

D=dE×N<<1.                                                            (7)

Поведение интегральной эффективности с увеличением длины отрезка показано на рисунке 8. На горизонтальной оси отложено расстояние от вершины функции распределения состояний  (начало координат) до границы области (состояния, лежащие ниже, используются на втором этапе). На рисунке представлены данные для размерностей сети N=100 и N=200. Из рисунка 8 и с учетом условия (7) можно заключить, что граница области находится на расстоянии, превышающем 2se, от максимума функции распределения. Так, для нейронной сети размерности N=100 можно использовать минимумы глубже 2,1se от пика распределения точек , а для N=200 – глубже 2,6se; se – полуширина функции распределения дискретизированных минимумов . С ростом размерности N наиболее глубокая величина  уменьшается, поэтому граница отодвигается влево. Значение границы, равное 2,0se, можно принять для всех размерностей нейронной сети, больших 100.

Распределение состояний  подчиняется нормальному закону. Доля минимумов состояний , лежащих глубже 2,0se, к общему их числу составляет 0,05. В этом случае увеличение скорости работы алгоритма составит .

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

Литература

1.     Hopfield J.J. Neural networks and physical systems with emergent collective computational abilities // Proc. Nat. Acad. Sci. USA. 1982. Vol. 79, pp. 2554–2558.

2.     Крыжановский Б.В., Магомедов Б.М., Микаэлян А.Л. Взаимосвязь глубины локального минимума и вероятности его нахождения в обобщенной модели Хопфилда // ДАН. 2005. Т. 405. № 3. С. 320–324.

3.     Kryzhanovsky B.V. The shape of a local minimum and the probability of its detection in random search // Lecture Notes in Electrical Engineering. 2009. Vol. 24, pp. 51–61.

4.     Крыжановский Б.В., Крыжановский М.В., Мальса- гов М.Ю. Дискретизация матрицы в задаче бинарной минимизации квадратичного функционала // ДАН. 2011. Т. 438. № 3. С. 312–317.

5.     Kryzhanovsky M.V., Malsagov M.Yu. Clipping procedure in optimization problems and its generalization // Optical Memory & Neural Networks. 2009. Vol. 18. № 3, pp. 181–187.