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

УДК 51:004+53:004

АЛГОРИТМ ВЫДЕЛЕНИЯ ТРЕНДА ЗАШУМЛЕННЫХ БОЛЬШИХ ВРЕМЕННЫХ РЯДОВ

М.Л. Бахмутский, к.ф.-м..н. (НИИСИРАН, г. Москва, mbakhmut@mail.ru, mbakhmut@ftiisi.ras.ru)

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

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

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

Сингулярный спектральный анализ

При изучении временных рядов, помимо анализа Фурье и вейвлет-анализа, широко применяется метод сингулярно-спектрального анализа (ССА) [1]. Данный метод основан на построении по исходному временному ряду, по так называемой траекторной матрице, ее аппроксимации при помощи сингулярного разложения матриц. Поэтому его можно рассматривать как исследование временного ряда хорошо известным в статистике методом главных компонент (другое название -метод естественных ортогональных функций). Для построения траекторной матрицы временной ряд {х 1(.. .,Хм) последовательно, со сдвигом на один отсчет, проходится окном размером т отсчетов и строится множество векторов задержки размерности т о), = (х ¡, х,+1( . . ., х ,+т _ 1 ) Т [2], которые составляют траекторную матрицу размерности т*п

A=

V

Х1 Х2 х3 х„ \

Х2 х3 х4 • хп+1

Х3 Х4 Х5 Хп+2 (1)

ХП1 + 1 хт+2 • xN /

где х1,х2,. . .,Хм - значения ряда в моменты времени 1=1, 2 ,..., К, п=К-т+1, т < р^-1] .

Для траекторной матрицы находится ее сингулярное разложение, то есть матрица А представляется в виде А=и£ УТ, где и - матрица размерности т*п с ортонормальными столбцами, и={ир); У={ур) - матрица размерности пхп с ортонормальными столбцами; £=diag (стр) - матрица размерности т*п с неотрицательными диагональными элементами (сингулярными числами). Матрица аппроксимируется при помощи г троек { стр, ир,ур], отвечающих г наибольшим сингулярным числам. Иными словами, матрица А заменяется матрицей А :

А=ст1 * Е1 + ст2 * Е2 + —I- стг * Ег, (2)

где Е,=и, хут - внешнее произведение Ьго столбца матрицы и и соответствующего столбца матрицы V; Е , - простая матрица ранга 1 с единичной нормой. Временной ряд восстанавливается путем диагонального усреднения матрицы , то есть матрицы, в которой оставлены только г более гладких компонент.

=

1 "i.j-i+l.lsjsm J

-S™iai,j-i+i, m < j < n,

(3)

N-j+l .

VN-j + 1

я:

*i+j-ii,n-i+l,n<j<N.

Если параметр m выбран таким, что m « n, можно решить спектральную задачу для матрицы C=A*AT=U*(£*£T)*U т.

Матрица U - это матрица собственных векторов матрицы C, а сингулярные числа матрицы A являются положительными корнями собственных чисел матрицы C. Умножив скалярно матрицу U т на A, получим матрицу главных компонент:

ГА

Y=UT*A=( ; I , где Y, , i=1, 2 ,..., m - строки дли-

W

ны n.

Тогда А=Л х uj xY, и временной ряд восстанавливается по формуле (3). Этот подход довольно широко используется в задачах нахождения периодических зависимостей во временных рядах наблюдений, фильтрации шумов и сглаживания временных рядов, так как использование сингулярного (спектрального) разложения матрицы позволяет выделить наиболее значимые составляющие ряда и отсеять случайные возмущения. Однако при непосредственном применении этого подхода к обработке очень больших рядов наблюдений возникают трудности. В изложенном методе ССА есть два существенных параметра - размерность окна, то есть размерность траекторной матрицы, и количество сингулярных троек (главных компонент), используемых для аппроксимации траекторной матрицы. Количество используемых главных компонент в данной задаче можно определять интерактивно.

Выбор длины скользящего окна определяется основной идеей метода ССА, то есть соображени-

ем о том, что, если исходный временной ряд имел какую-то структуру, то и его отрезки (столбцы траекторной матрицы) имеют такую же структуру. Если стоит задача выделения составляющих с известным периодом, размер окна равен или кратен этому периоду. Поэтому в задаче о выделении глобального тренда и сглаживании ряда рекомендуемый размер окна примерно равен половине длины ряда. При этом выбор размера движущегося окна оказывается очень существенным. На рисунке 1 изображен временной ряд, полученный в ходе реального эксперимента по гидродинамическому прослушиванию скважин. На рисунках 2 и 3, иллюстрирующих влияние выбора размера окна, приведены результаты сглаживания (белая линия) исходного временного ряда (черная линия) с использованием окон в 100 и 1 000 отсчетов соответственно. Были проведены численные эксперименты по последовательному применению метода ССА. Исходный ряд сглаживался методом ССА с использованием окна размером 100, 500, 5 000 отсчетов. На рисунке 4 приведено сопоставление результатов последовательного сглаживания для различных окон (100, 500 и 5 000 отсчетов). Из этих расчетов следует, что для надежного и корректного нахождения тренда временного ряда необходимо использовать траекторную матрицу с длиной столбца, равной примерно половине длины ряда. Заметим, что для вычисления сингулярного разложения матрицы размерности порядка 5 000*5 000 на компьютере c CPU Intel Core i5-750 (2,67 GHz) потребовалось несколько суток непрерывной работы. Это вполне понятно, так как вычисление SVD-разложения матрицы размерности m*n требует порядка 4*п 2 * ш — ^ * п 3 +

+О (п2) операций [3]. Отсюда следует, что непосредственно применить метод ССА к выделению тренда и сглаживанию больших временных рядов не удается, так как невозможно вычислить сингулярное разложение матриц размерности порядка 65 000*65 000 элементов, во всяком случае, для компьютеров традиционной последовательной архитектуры. Поэтому предлагается сочетать метод ССА и применение вейвлет-пакетов.

Обобщение метода ССА на случай больших временных рядов

Допустим имеется большой, порядка десятков тысяч отсчетов, временной ряд длины N. Будем полагать без потери общности, что N=2 k* N ^ Пусть N=2*M. Вместо временного ряда X= = рассмотрим

ряд Z, состоящий из двух временных рядов S и D. Z={S | D}, где S={s1( s2,. . ., sM} и D={ di , d2 ,. . ., dM}.

Переход от исходного ряда X к ряду Z происходит следующим образом: si=S±2llS*i; di=^™ i=i, 2, M. (4)

Рис. 1. Исходный временной ряд 130 816 значений

давления, измеренных с интервалом 1 минута

Рис. 2. Исходный ряд (черная линия) и его сглаживание (светлая линия) методом ССА с движущимся окном в 100 отсчетов

Рис. 3. Исходный ряд и его сглаживание методом ССА с движущимся окном в 1 000 отсчетов

.,.1 , . „

Г

iWr' ™ '

k*

sr-.rs::.™——- !

™

Рис. 4. Вторичное сглаживание исходного ряда методом ССА (сплошная черная линия - окно 100 отсчетов, сплошная светлая - окно 500 отсчетов, светлые точки - 5 000 отсчетов)

Обратное преобразование от S и D к X имеет следующий вид:

х 2 * i - !=lf; X 2 * i = s2f , i=1, 2, . M; (5)

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

Б - уточняющие коэффициенты для элементов ряда X.

Фактически на основе (4) раскладывается ряд X по вейвлетам Хаара [4].

Преобразования (4), (5) являются взаимно-однозначными, как легко убедиться непосредственно, и сохраняющими нормы. То есть, если I I X | | 2 = £ 2=5* (х,) 2, | | 5 I I 2 = £ Г=1 (5 ,) 2, | | О | | 2 = , то .

Таким образом, пренебрежение какой-либо величиной в преобразованном временном ряде означает (после обратного преобразования) пренебрежение такой же малой величиной в исходном. На следующем этапе применяем введенное преобразование к рядам 8 и Б. То есть от ряда Z перейдем к ряду Ъ ( 2 ) = (5 5 О 5 | 5 О О О ) по тем же правилам, что и при переходе от ряда X к ряду Z

ss, = sd, =

S2,i-1 + S2!

V2

d2,i-i + d2,i V2

1 ; dSj =

; ddj =

s2 + i—1—S2+j

H ' d2,i-i-d2)i

V2

Обратный переход:

_ ssj + dsj

S2ti"1 - ~JT

_ sdj+ddj

S2»i =

ssj—dsj

V2 '

Mi =

, в (1: — с! с!; . . _ __

-—; d2* , = ' ' , здесь 1=1, 2, ..., М1,

V 2 V 2

М_ N 2 " 4 .

Все то, что справедливо при переходе от ряда X к ряду Z, справедливо и для перехода от ряда Z к ряду Ъ ( 2 -1. Ряд 88 является временным рядом локально сглаженных, усредненных значений элементов ряда 8, а ряд Б8 состоит из уточняющих коэффициентов. Аналогичное справедливо и для рядов 8Б и ББ. При этом ряд 88 описывает самую гладкую часть исходного временного ряда. Сохраняя взаимно-однозначное соответствие, удалось перейти от исходного временного ряда к ряду, состоящему из двух временных рядов (гладкой части и поправок), а затем к ряду, состоящему из 4 временных рядов (гладкой части и поправок различной гладкости). Последовательно применяя этот прием (называемый вейвлет-пакетным разложением [4]) к раз получаем временной ряд, состоящий из временных рядов, имеющих длину N. Для каждого из этих 2 к временных рядов построим свою траекторную матрицу (1) размерно-

[(N1+1)1 .. , .

стью т*п, где т = I—-—], п = N 1 - т + 1 . Вычислим сингулярные разложения этих матриц. Заметим, что эта часть алгоритма может быть эффективно выполнена как на традиционных последовательных, так и на параллельных системах. Анализируя сингулярное разложение траекторной матрицы наиболее гладкой части, определяем количество сингулярных троек (или главных компонент), описывающих тренд, и тем самым определяем тот порог, ниже которого пренебрегаем сингулярными числами всех траекторных матриц. После этого формируем аппроксимации (2) траек-

торных матриц, по формулам (3) восстанавливаем временные ряды и, последовательно применяя (5), получаем сглаженную версию исходного временного ряда.

Обратимся к примеру. Исходный ряд состоит из 130 816 отсчетов, то есть N=130 816, k=8, N j =511. На 8-м уровне вейвлет-пакетного разложения имеется 256 временных рядов длиной 511 значений, каждый из которых соответствует (как гладкая часть и поправочные коэффициенты) всему исходному временному ряду. На рисунке 5 приведен график гладкой части вейвлет-разложе-ния, на рисунках 6 и 7 следующие 2 временных ряда поправочных коэффициентов. Заметим, что график на рисунке 5 внешне похож на график рисунка 4, соответствующий последовательному двойному сглаживанию методом ССА с окном в 100 отсчетов, поскольку последовательное применение первой формулы преобразования (4) можно рассматривать как своеобразное усреднение в окне 256 отсчетов.

Рис. 5. Самая гладкая часть 8-го уровня

вейвлет-преобразования

Рис. 6. Следующий после самой гладкой части массив поправочных коэффициентов 8-го уровня

вейвлет-преобразования

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

На рисунках 6 и 7 приведены следующие после самой гладкой части (рис. 5) массивы вейвлет-ко-эффициентов. Легко заметить, что их нерегулярность увеличивается, то есть растет влияние шумовой компоненты. Однако просматривается также закономерность. Влияние шума на графике на рисунке 7 больше, чем на рисунке 6, и тем более, чем на рисунке 5. На рисунке 8 приведен график 20 наибольших значений сингулярных чисел для 1-й и 2-й траекторных матриц. Предположим, что в разложении (2) решено ограничиться только первым членом. Тогда, отбрасывая все члены в разложении (2) для всех рядов коэффициентов, сингулярные числа которых меньше использованного, и к раз проделав обратное вейвлет-преобраз-ование, получим тренд (светлая линия) на рисунке 9. Последовательно добавляя в разложение матрицы (2) новые члены, повышаем детальность воспроизведения временного ряда. На рисунке 10 приведен график первого ряда вейвлет-коэффици-ентов, восстановленного с использованием 10 первых членов в (2). А на рисунке 11, где приве-

Рис. 11. Тренд первого ряда вейвлет-коэффициентов. Построен при использовании 11 членов разложения (2). Видно появление шумовой компоненты

Рис. 12. Исходный ряд (черная линия) и тренд, найденный с использованием 10 сингулярных троек

(10 членов аппроксимации (2))

Рис. 13. Исходный ряд (черная линия) и тренд, найденный с использованием 50 сингулярных троек. Шумовая компонента усиливается

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

Рис. 8. Первые 20 сингулярных чисел из графика рисунка 10

Рис. 9. Исходный ряд (черная линия) и тренд, найденный с использованием 1-й сингулярной тройки

Рис. 10. Тренд первого ряда вейвлет-коэффициентов. Построен при использовании 10 членов разложения (2)

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

7 4

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

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

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

N1

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

пользовании предлагаемого алгоритма уменьшается в 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. Матрица А симметричная, с нулевой диагональю, а ее элементы - случайные независимые величины. Количество минимумов функционала E( Я) велико и экспоненциально растет с ростом размерности задачи. Так, уже при размерности вектора Я, равной 200, число минимумов ~100 000, и поэтому нахождение достаточно глубокого минимума сложная задача.

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

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

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

Й = -B + Ax s , (2)

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

S = signH • (3)

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

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

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

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