Применение чисел малой разрядности в задаче бинарной оптимизации Текст научной статьи по специальности «Математика»

Сравним стоимость 5УО-разложения траектор-ной матрицы исходного ряда в предположении равенства длины движущегося окна половине длины ряда и траекторных матриц предлагаемого алгоритма. Стоимость £КО-разложения матрицы

7 4

размерности ш*п примерно равна 4 * п 2 * щ — - *

* п 3 + О (п2) операций [3]. Если длина движущегося окна примерно равна половине ряда, то п « щ « 2 к_ 1 * М и стоимость 5УО-разложения

~ 2з_ * м 13 + О (2 2 к"2 * N 1 2 ) операций. Для данного алгоритма стоимость £КО-разложения 2 к

N1

траекторных матриц с п « щ « примерно равна 2 к * + О (М-) | операций. Таким образом, стоимость вычисления 5УО-разложения при ис-

пользовании предлагаемого алгоритма уменьшается в 2 2 k раза.

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

Автор благодарен С.Г. Вольпину, предоставившему экспериментальный материал.

Литература

1. Vautard R., Yiou P., Ghil M. / Physica D. 1992. Vol. 58.

2. Broomhead D.S., King G.P. / Physica D. 1986. Vol. 20.

3. Деммель Дж. Вычислительная линейная алгебра. М.: Мир, 2001.

4. Столниц Э., ДеРоуз Т., Салезин Д. Вейвлеты в компьютерной графике. Ижевск: НИЦ «Регулярная и хаотическая динамика», 2002.

УДК 004.8:9

ПРИМЕНЕНИЕ ЧИСЕЛ МАЛОЙ РАЗРЯДНОСТИ В ЗАДАЧЕ БИНАРНОЙ ОПТИМИЗАЦИИ

М.В. Крыжановский, к.ф.-м..н.; М.Ю. Мальсагов

(НИИСИ РАН, г. Москва, iont.niisi@gmail.com, magomed.malsagov@gmail.com)

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

Ключевые слова: нейронные сети, модель Хопфилда, дискретизация, минимизация, квадратичный функционал.

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

Е= -(я, А« )+2-(В , я), (1)

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

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

минимума [2, 3]. Это означает, что нейросеть с подавляющей вероятностью находит если не оптимальное решение (глобальный минимум), то одно из субоптимальных (локальный минимум).

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

Й = -B + Ax s , (2)

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

S = signH • (3)

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

Для уменьшения объема вычислений в [4]

предложен метод дискретизации матрицы A , то есть выделение из матричных элементов среднего

значения A0 и замена остатка матрицей C , нормированные элементы которой имеют целочисленные значения в диапазоне [-m; +m], где m -число градаций. При таком подходе поиск минимума E( S) заменяется поиском минимума дискре-

тизованного функционала е = -(s,Cs) + 2-(b,s), если A0=0.

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

h = -b+C • s, (4)

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

S = signh . (5)

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

Минимум функционала. В данной работе применяется двухэтапный алгоритм [5]. Процесс минимизации функционала (1) начинается с некоторой случайной точки пространства S. Подчиняясь решающему правилу (5), сеть конвергирует в некоторое устойчивое состояние s"0*, являющееся

минимумом функционала е. Если из этой точки продолжить спуск с решающим правилом (3), то сеть конвергирует в состояние S o, соответствующее минимуму функционала (1). На всем пути S ^ s* ^ S o величина ошибки (вероятность ошибки P в направлениях градиентов) только уменьшается. С ростом числа градаций m величина ошибки убывает как P = ——1-г-.

я( 2m +1)

Эффективность минимизации определяется величиной 5E = (E* - E„) / E„, где E* = E( so*) -

энергия минимума s"0* дискретизированного функционала; E0 = E(so) - энергия ближайшего к

3

1 \ s°

_в

Mill Г S

\

Рис. 1. Распределение плотности вероятности состояний по энергии для разного числа градаций

Рис. 2. Функции распределения состояний

нему минимума ¥ ю в который можно было спуститься при использовании исходного решающего правила (3). Показано [4], что эффективность минимизации 8E«P, (6) из которой следует, что разница в энергиях с ростом m убывает как 8E~m-1.

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

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

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

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

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

чисел матрицы А будем называть исходным форматом, а p обозначим количество чисел малой разрядности, которое упаковывается в исходный 4-байтный формат.

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

Ускорение алгоритма определяется величиной 0=10Я, где Х0 - исходное время работы алгоритма; I - время работы алгоритма при упакованных элементах матрицы.

0,4 06 0,8

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

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

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

0 а.

\

2000 4000 6000 8000 10000 12000 N 10 20 30 40

л Рис. 4. Количество Рис. 3. Ускорение переворотов спинов алгоритма при на каждой итерации использовании упаковки 7 (размерность нейросети чисел малой разрядности

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

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

Т 100 80 60 40 20 0

0 2000 4000 6000 8000 10000 12000

,-/100 80 60 40 20 0

—-

0 2000 4000 6000 8000 10000 1200 N

б

Рис. 5. Вычислительная сложность алгоритма

при изменении состояния нейрона вектор Н модифицируется по правилу Н = Н +2( А )ь если направление спина положительно/отрицательно;

(А ) - вектор строки с номером текущего нейрона. Объем вычислений компонент Н в исходном состоянии равен При перевороте каждого спина производится N операций. Число переворотов спинов не превышает N (рис. 5б), поэтому вычислительная сложность такого алгоритма ~(1+а)№2, а<1. Отношение объемов вычислений такого метода к исходному методу ~Т, (Т>>1). Далее в качестве исходного алгоритма используется метод с обновлением компонент Н.

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

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

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

I оН

чае 0 =10 - 0

I о+о2

На первом этапе число итераций при минимизации е( Я) равно числу итераций исходным методом [5], поэтому объем вычислений на первом этапе составляет 01=0н/5. На втором этапе число итераций Т мало [5], соответственно мало и число изменений в состояниях нейронов, поэтому основной объем вычислений приходится на вычисление градиента Н (в промежуточных состояниях и составляет 02=0н/2. В этом случае ускорение алгоритма составит 0ю1,4.

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

оН

ускорение составит 0 я-ю 2,2 .

О+°»5о2

Конечно, при таком подходе велика вероятность получить несколько худшее решение, оценить которую можно экспериментально (рис. 7).

N

а

1

—

4

ч ГЁ2оь

0 .00 ,000 ,000 .00 ,.0 ,0. ,000 8000 ,000 N ... ... ,.. .. .... ... N

Рис. 6. Ускорение, Рис. 7. Вероятность получаемое попадания в самый при упаковке чисел (р=8) глубокий минимум из по сравнению левой и правой областей с обычнът методом распределения I*

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

Чтобы определить границы области, состояния

—*

¥ которых используются на втором этапе алгоритма, и тем самым уменьшить объем вычислений, определена интегральная эффективность минимизации. То есть определялся минимум величины 5Е = (Е^ - е) / Е^ на отрезке, начало

которого соответствует наиболее глубокому минимуму дискретизированной сети. С помощью 8E можно оценить хеммингово расстояние D между глобальным минимумом Е* и самым глубоким минимумом E, найденным при стартах из отрезка. Действительно, при изменении состояния одного нейрона энергия меняется в среднем на величину e. Поскольку все минимумы при N>100 находятся в области Ес+3а и отношение ае/Ес<0,01, то

Е* ~ Е Ес~е^ Тогда 5Е ~ ^ = ^ .

Таким образом, достаточно подобрать область

^ * ^

состояний , для которой выполняется условие

D=8E•N<<1. (7)

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

г SE ■ 0,045 ■ 0,04 ■ 0,035 ■ 0,03 ■ 0,025 ■ 0,02 ■ 0,015 ■ 0,01 ■ 0,005

\N = 200 /

СГ -3,5 -3 -2,5 -2 -1,5 -1 -0 5

Рис. 8. Интегральная эффективность минимизации от размера стартовой области (размер стартовой области задается в единицах стандартного

—* \

отклонения энергии е в промежуточных точках )

распределения точек е( s0*), а для N=200 - глубже 2,6сте; сте - полуширина функции распределения дискретизированных минимумов е(?0*) . С ростом

размерности N наиболее глубокая величина е( so*)

уменьшается, поэтому граница отодвигается влево. Значение границы, равное 2,0сте, можно принять для всех размерностей нейронной сети, больших 100.

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

OH

работы алгоритма составит 8 »-» 5,0.

о + 0,0502

На основании изложенного можно сделать вывод о том, что предложенная процедура дискретизации позволяет использовать числа малой разрядности для ускорения минимизации квадратичного функционала. Авторами разработан двух-этапный алгоритм на основе нейросетевой модели Хопфилда, дающий возможность быстро находить решение поставленной задачи. Ускорение работы алгоритма составляет 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. N° 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.