CC BY
23
УДК 519.246.8+519.853.32
Вестник СПбГУ. Сер. 1. Т. 3(61). 2016. Вып. 4
ПОИСК ВЕСОВ В ЗАДАЧЕ ВЗВЕШЕННОЙ АППРОКСИМАЦИИ ВРЕМЕННЫМ РЯДОМ КОНЕЧНОГО РАНГА*
Н. К. Звонарев
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Рассматривается задача взвешенной аппроксимации временного ряда рядом конечного ранга с целью оценивания сигнала. Решается проблема поиска весов, которые приводят к улучшению точности оценок. Строится и теоретически обосновывается эффективный метод численного расчета поиска весов с помощью теории квадратичной оптимизации. Для того чтобы получить эффективный алгоритм, задача квадратичной оптимизации с большим числом ограничений сводится к последовательности задач с меньшим числом ограничений и критерию остановки. Для обоснования алгоритма доказывается эквивалентность различных формулировок исходной задачи оптимизации. Проводится численное моделирование, подтверждающее эффективность алгоритма и улучшение точности метода оценивания сигнала. Библиогр. 10 назв. Ил. 2. Табл. 1.
Ключевые слова: временные ряды, метод Cadzow, ряды конечного ранга, взвешенный метод наименьших квадратов, сингулярное разложение, квадратичное программирование.
Введение. Рассмотрим задачу извлечения сигнала S = (si,..., sn) из зашум-ленного временного ряда X = S + N, где S управляется линейной рекуррентной формулой (ЛРФ) порядка r: sn = Y1 ¿=1 aiSn-i, n = r + 1,...,N, ar = 0, N —белый шум. Ряды, управляемые ЛРФ, могут быть записаны в параметрическом виде sn = i Pi(n) exp(ain) cos(2nwin + фi), где Pi(n) —многочлены от n [1].
Известно, что для выделения такого рода сигналов хорошо работают методы, основанные на оценке сигнального подпространства (subspace-based methods) и, в частности, анализ сингулярного спектра (см. [1-4] и [5, глава 4]). Идея методов состоит в следующем: зафиксируем длину окна L, 1 < L < N, положим K = N — L +1 и определим траекторную матрицу для ряда S G XN (xn — множество вещественных временных рядов длины N):
где Tl обозначает биекцию между Xn и H, и H — множество ганкелевых матриц порядка L х K с одинаковыми значениями на побочных диагоналях i + j = const.
Скажем, что ряд S, траекторная матрица S которого имеет неполный ранг r, называется рядом L-ранга r. Если L-ранг не зависит от L при достаточно больших N и L, ряд называется рядом конечного ранга. Рассмотрим задачу аппроксимации временного ряда рядом L-ранга по взвешенной норме [6, 7]:
/ Si S2 S2 S3
sk+i
S = Tl(S)
\sl sl+1
sn
N
(1)
* Работа выполнена при финансовой поддержке РФФИ (грант №16-04-00821). (¡5 Санкт-Петербургский государственный университет, 2016
где X = (xixN) € XN — исходный временной ряд, Y = (yi,..., yN) — требуемая аппроксимация, Y = TL (Y) = T(Y) — соответствующая траекторная матрица, qi > 0, i = 1,...,N, — веса (в дальнейшем называемые весами ряда).
Полученное решение Y* задачи (1) можно использовать как оценку сигнала S рядом конечного ранга. Более того, полученная оценка позволяет решить множество связанных задач: оценка коэффициентов ЛРФ и параметров ряда, прогнозирование ряда, заполнение пропусков и так далее [1]. Заметим, что в случае белого шума естественно выбирать веса qi равными.
Обозначим соответствующее скалярное произведение в Xn как (Ж, Y)q — qixiyi. Для решения задачи (1) в [6] рассматривается эквивалентная задача структурной аппроксимации матрицей неполного ранга (Structured Low-Rank Approximation, SLRA):
||X - Y||c ^ min (2)
rank y <r
Yen
где X = T(X), а матричная норма || • ||c (C —диагональная, положительно определенная матрица) порождена следующим скалярным произведением:
l к
(Y, Z)c = tr(YCZT) = ]T Y. Ck,
l=1 k=1
где C = diag(ci,...,cK), Y = (yl}k), Z = (zl}k). Между весами qi в задаче (1) аппроксимации ряда и весами ci в задаче (2) (в дальнейшем называемые матричными весами) существует соотношение [6, предл. 4] для эквивалентности норм.
Основная проблема в [6] состояла в нахождении весов Ci, соответствующих равным весам qi, естественных для модели с белым шумом. Тем самым решение задачи (2) соответствовало оценке ряда, являющейся решением задачи невзвешенных наименьших квадратов в пространстве рядов. Однако опять же согласно [7] и [8], если это и возможно, то только при вырожденной матрице C, что недопустимо при решении задачи методами типа Oblique Cadzow [6, замечание 4]. Следовательно, равномерностью весов ряда qi приходится жертвовать и использовать матричные веса ci такие, чтобы C была невырожденной, а веса qi лишь близкими к равномерным.
Цель статьи — разработка эффективного алгоритма поиска весов при зафиксированной мере невырожденности матрицы C. В качестве такой меры ограничим число обусловленности C снизу. Так же как ив [6], введем параметр 0 < а < 1 и потребуем, чтобы все ci удовлетворяли неравенствам ci > 0 и
mini c. , . -— > 3
maxi c
что эквивалентно неравенству ^nd(C) < 1/а. После этого поставим задачу аппроксимации единичных весов ряда по наиболее естественной евклидовой норме. В этом и состоит ключевое отличие от предыдущего подхода, описанного в [6, предл. 5] и неприменимого в большинстве практических ситуаций. Например, в [6] требуется, чтобы длина ряда N была кратной длине окна L. Данная работа посвящена решению задачи в такой постановке, что требует применение теории экстремальных задач и квадратичного программирования (КП). Поэтому получение эффективного по скорости работы алгоритма поиска весов является нетривиальной задачей.
В [6, раздел 5] показано, что чем равномернее веса ряда в задаче (1), тем точнее получается решение. Более того, как было показано в [6], параметр а влияет на скорость сходимости метода Oblique Cadzow: чем больше Cond(C) (что равнозначно меньшему значению а), тем медленнее сходится метод. В разделе 6 показано, что эти соотношения выполняются и для весов, полученных новым методом.
1. Постановка задачи аппроксимации весов. Для удобства переформулируем предложение [6, предл. 4], связывающее ci и qi, в матричном виде. Для этого рассмотрим матрицу B = (bitj) порядка N х K, имеющую следующий вид:
il, для г = j,...,j + L - 1,
bi,j = 1 п (4) I U, в противном случае.
Заметим, что B — матрица полного ранга. Предложение выглядит следующим образом.
Предложение 1. Пусть Y = T(Y), Z = T(Z), Q = (qx,...,qN)T, C = (ci,...,ck )T, C = diag(C ). Тогда (Y, Z)q = (Y, Z)c для любых Y, Z G xn тогда и только тогда, когда Q = BC.
В общем виде задачу аппроксимации весов сформулируем следующим образом.
Задача 1.
р* = min p(C), где p(C) = ||BC - 1n||2, (5a)
=C
а>0, г = 1 ,...,К, (5Ь)
тахг вг )
где = (1,..1)Т € — требуемые веса ряда (вектор из N единиц), 0 < а < 1 — параметр, регулирующий степень вырожденности матрицы С, || • || —евклидова норма
Замечание 1. В дальнейшем случай а = 1 опустим. При а =1 все вг = в равны между собой, и задача 1 сводится к поиску множителя в такого, что ||вВ1к — ||2 достигает своего минимума. Он достигается в точке в = ВТ1^/(1^ВТВ1к).
Решение задачи 1 можно использовать для решения задачи (2) и тем самым для решения задачи (1).
2. Эквивалентные формулировки. Рассмотрим выражение ц>(С)/2, где ц>(С) определено в (5а). Преобразуем его, не учитывая константное слагаемое, а множество С из (5Ь) запишем в виде набора линейных ограничений. Получим следующую эквивалентную формулировку.
Задача 2.
/* = тт/(С), где /(С) = ^СТ8С-ЫТКС, (6а)
С = {С I вг - авэ > 0, г = у, 1 < г,] < К} , (6Ь)
в RN
где S = (si,j) = BTB — положительно определенная матрица размера K х K с элементами
jL — \i — j\, \i — j\<L,
si ,j = (7) I U, в противном случае,
а вектор L1k = BT1W = (L,..., L)T G rk.
Теперь задача (2) представлена в том виде, для которого разработана теория квадратичного программирования [9, с. 111]. Справедливо следующее утверждение.
Предложение 2.
1. Задача 2 имеет единственное решение C*.
2. Ее решение C* = (с*,... ,c*K)T является палиндромом, то есть для любого индекса 1 < i < K справедливо с* = c*K_i+1.
Доказательство.
1. Задача 2 — задача КП с набором из K2 линейных ограничений (6b) и целевой функцией (6a), квадратичная форма в которой положительно определена, поэтому решение такой задачи единственно [10, с. 452].
2. Достаточно рассмотреть вектор C** = (cK,...,c*) и заметить, что f (C*) = f (C**) и C** gC в (6b); значит, C* = C**, что и требовалось доказать. ■
Для решения задачи 2 можно использовать методы квадратичного программирования, но, к сожалению, K2 линейных ограничений являются серьезной проблемой при решении задачи на практике. Рассмотрим еще две эквивалентные формулировки, свойства которых используются при решении задачи.
Сделаем дополнительное предположение: пусть в точке j веса достигают своего максимума (тем самым снизим число линейных ограничений до линейного по K размера), а само j будем перебирать в цикле от 1 до \K/2]. При этом воспользуемся тем фактом, что решение является палиндромом. Формально задача выглядит следующим образом.
Задача 3.
f ** = min f*, (8a)
fj* = min f (C), (8b)
Cj = {C \ а = ск-i+i, i =1,..., K/2]; (8c)
ci — acj > U, i = 1,...,j — 1,j + 1,..., [K/21; (8d)
cj — ci > U,i =1,...,j — 1,j + 1,..., [K/2] }. (8e)
Чтобы не перебирать все возможные j, удобно уметь проверять, дает ли полученное на очередной итерации решение глобальный минимум. Для такой проверки рассмотрим еще одну эквивалентную формулировку задачи 1. Введем дополнительную переменную cmax, хранящую максимальный вес, и перейдем к новым переменным Ci = cmax — ci, i = 1,...,K, определяющим разницу между максимальным среди всех и текущим весами. При этом выполняются неравенства cmax > U, ci > U. Между векторами C = (c1,..., cK)T и C = (c1 ,..., cK, cmax)T существует простое линейное соответствие: C = HC7, где H G rkх(к+1) — матрица следующего вида:
( —1 о 1 \
H = .. . I . (9)
V 0 —1 1 /
Условие (3), устанавливающее границу снизу для весов, в новых обозначениях записывается как (1 — a)cmax — с\ > 0, г = 1,...,K.
Получаем следующую задачу квадратичного программирования с линейными ограничениями, но с нестрого выпуклой целевой функцией.
Задача 4.
/* = min f(Ô), где f(C) = f(UC) = l-CTUTSUC-LtTKUC\ (10а)
сей 2
C ={ C\èj > 0, (1 — a)cmax — Cj > 0, j = 1,...,k}. (10b)
Теорема 1. Задачи 1, 2, 3 и 4 эквивалентны.
Доказательство. Эквивалентность задач 1 и 2 очевидна, так как задача 2 является переформулировкой задачи 1. Для оставшихся задач приведем идею доказательства.
Требуется доказать f * > f ** > f* > f * и тот факт, что по решению одной из задач можно построить решение любой другой.
Докажем неравенство f * > f **. Пусть C* — решение задачи 2, то есть C* = argminCeC f (C ). Положим j* = argmaxici. Имеем C* G C и C* —палиндром, следовательно, C* G Cj*. Цепочка неравенств f * = f (C*) > f* = minc eCj* f (C ) > min j=i..^K/2~\ minCeCj f (C ) = f ** завершает доказательство неравенства.
Докажем неравенство f** > f*. Пусть j * = argminj=1,_jK/2\ f*, C** = argminc * f (C ) —решение задачи 3. Положим cmax = Cj*j*, ci = cmax — Cj* ,i,
г = 1,...,K, где Cj** = (cj*i,...,Cj*,K)T. Нетрудно заметить, что C G C, из чего
получим f* < f (C) = f (C** ) = f **.
Докажем неравенство f * > f *. Пусть C* = argmince ¿f (C) — решение задачи 4. Положим ci = cmax — с и заметим, что C G C, из чего получим f * < f (C ) = f (C) = f*.
3. Общий алгоритм решения. Используя эквивалентность, доказанную в теореме 1, можем сформулировать следующий алгоритм решения задачи 1. Алгоритм 1. Вход: Параметры Ь, К, а. Результат: Вектор оптимальных весов С*.
1. Положить у = 1.
2. Решить подзадачу КП С* = (в^,..., )Т = а^ттс еС. / (С).
3. Задать эквивалентный вектор (О, взяв втах = вjtj, сг = втах — вj .
4. Проверить, является ли вектор С решением задачи 4. Если да, положить С* = С* решением задачи, иначе взять у ^ у + 1 и перейти к пункту 2.
Таким образом, если встретился нужный индекс ], алгоритму не понадобится перебирать оставшиеся индексы и решать подзадачу (8Ь) лишний раз. Практические эксперименты показывают, что максимальный вес всегда находится на краях: алгоритм 1 останавливается уже при у = 1, то есть фактически сводится к задаче КП с положительно определенной квадратичной формой в целевой функции.
Для реализации алгоритма 1 необходимо разработать алгоритмы решения задач из пункта 2 и 4, что и будет сделано в следующих разделах работы. Сначала в разделе 4 рассмотрим пункт 4, а в разделе 5 — пункт 2.
4. Проверка вектора на решение задачи 4. Предложим быстрый алгоритм, проверяющий, является ли заданный вектор С точкой, в которой достигается глобальный минимум в задаче 4, то есть реализующий пункт 4 из алгоритма 1. Для этого применим теорему о необходимом и достаточном условии минимума в задаче квадратичного программирования для задачи 4.
Теорема 2.
1. Рассмотрим вектор Я = (т\,... ,тк+1)т = НТ8ИС — ЬИТ1к. Тогда С является решением задачи 4 в том и только в том случае, если выполняются условия:
(1а) ЯТС = 0,
(1Ь) существует вектор V = (и1,...,ик)т € Рк такой, что имеют место неравенства и^ > 0, г = 1,...,К, щ > —т^, г = 1,...,К, (1 — а) ^ к=1 щ < г к+1.
2. Положим щ = тах(0, —т^), г = 1,...,К. Тогда С является решением задачи 4 в том и только в том случае, если выполняются условия:
(2a) RTC = 0,
(2b) (1 - Ki ui < tk+1.
Доказательство. Первый пункт данной теоремы является переформулировкой теоремы [9, теорема 9.2] для задачи 4, поэтому перейдем сразу к доказательству второго пункта. Первое условие в обоих пунктах одинаковое. Докажем эквивалентность второго условия.
Достаточность. Коэффициенты ui выбраны так, чтобы условия ui > 0 и ui > —ri выполнялись. Если удовлетворено и второе условие, выполняются все требования первого пункта.
Необходимость. Очевидно, что щ выбраны наименьшими из всех тех, которые удовлетворяют условиям ui > 0 и ui > -ri. Пусть имеем неравенство (1 — а) ^^ ui > tk+i. Рассмотрим любой вектор U = (ui), удовлетворяющий условию (2b), то есть (1 — а) ^k=1 ui < tk+1 . Тогда существует индекс i такой, что ui < ui. Следовательно, либо ui < 0, либо ui < —Ti, из чего следует, что не существует такого вектора U, который бы удовлетворял условиям первого пункта. ■
Рассмотрим следующий алгоритм, базирующийся на теореме 2.
Алгоритм 2. Вход: предполагаемое решение С, параметр а.
Результат: Булево значение: является ли С решением задачи 4.
1. Вычислить R = HTSHC — LHT1K.
2. Если RTC = 0, вернуть FALSE, иначе перейти к следующему пункту.
3. Выбрать в качестве вектора U = (ui,...,uk)T следующие значения: ui = max(0, —Tj), i = 1,...,K.
4. Если (1 — a)Y,K1 ui < tk+1, вернуть TRUE, иначе FALSE.
5. Задача квадратичного программирования специального вида. Теперь переходим к пункту 2 алгоритма 1. Для решения подзадачи (8b) при зафиксированном индексе j воспользуемся так называемым «Primary Active Seto-методом решения задачи КП, описанным в [10]. Специфика задачи позволяет эффективно реализовать алгоритм.
5.1. Общая схема «Active Set»-Memoda. В общем случае задача квадратичного программирования выглядит следующим образом.
Задача 5.
min -XTGX — VTX, (IIa)
xex 2
x = {X | AjX = pi, г e y; (11b)
AjX > pi, г e Z}, (11c)
где G e rk xK —произвольная положительно определенная матрица, V e rk —произвольный вектор, y и z — множества индексов, векторы Ai e rk вместе с pi e r задают ограничения.
Суть любого «Active Seto-метода состоит в последовательном переборе подмножества ограничений, которые выполняются как равенство для промежуточной точки — кандидата в решение задачи квадратичного программирования. Такое множество называется рабочим множеством и обозначается w С y U z, при этом всегда y С w, а ограничения, лежащие в рабочем множестве, называются активными. Ниже представлена схема этого метода для решения задачи 5.
Алгоритм 3 [10, с. 472]. Вход: параметры задачи квадратичного программирования: матрица G, вектор V, множества y, z и коэффициенты условий Ai, pi.
Результат: Решение X* задачи 5.
1. Найти начальную точку Xo, удовлетворяющую условиям задачи (11b), (11c), и положить wo — множество активных ограничений в этой точке.
2. Положить к = 0.
3. Положить Gk = V — GXk.
4. Решить подзадачу квадратичного программирования и найти множители Ла-гранжа ui, i e Wk:
Pfc* = argminPfcePfcipfcTGPfc - GjPk, (12a)
pk = {Pk | AjPk =0, г e wk}, (12b)
после чего найти ui из системы уравнений ieWk uiAi = GP* — Gk .
5. Если P* = 0, и все ui > 0, i e wk П z, положить X* = Xk и STOP.
6. Если P* = 0, но существует i такое, что ui < 0, взять j = argminjeWfcnZuj, положить wk+i = wk \ {j}, увеличить к на единицу и перейти к п. 3.
7. Положить ак = min (l, min^Wfci АтР*<0 ) •
8. Положить Хк+1 = Хк + а.кР*.
9. Если ак < 1, положить ] = а^тш^^ а^р*<оЕ£аТр^± и = ^ и Ш> иначе Wk+1 = Wk.
10. Положить к ^ к +1 и перейти к пункту 3.
5.2. Специфика задачи. Подзадача (8Ь) при фиксированном у переписывается в терминах задачи 5 следующим образом: X* = С*, X = С, С = Я, V = Ь1К, (11Ь) состоит из (8с), (11с) состоит из (8^ и (8е).
Таким образом, необходимо объяснить, как находить начальную точку (п. 1 алгоритма 3), решение подзадачи квадратичного программирования и множители Лагранжа (п. 4 алгоритма 3) применительно к частному случаю задачи (8Ь).
Обозначим вкратце те особенности, которые помогают получить быстрое решение.
Для выполнения пункта 4 алгоритма 3 требуется уметь решать задачу (12а) с ограничениями (12Ь). Положив А = [Аг : г € Wk], А € ^Ккт, ограничения (12Ь) можно записать как АТРк = 0, где т — количество активных ограничений.
Для решения поставленной вспомогательной задачи (12а) есть явная формула обобщенного метода наименьших квадратов:
Р? = А (АТС А! 1 АТС
к = А ^А СА^ А ^к,
где матрица А € ^Кх(к-т) СОстоит из столбцов, составляющих базис ортогонального дополнения к базису столбцов матрицы А. Обычно т велико, поэтому (К — т) мало, что позволяет быстро искать решение подзадачи.
В случае задачи 3 матрица С имеет простой вид, и ее можно умножать на вектор за время О (К). Матрица А — разреженная, в ней содержится максимум 2 ненулевых коэффициента в каждом столбце. За счет этого матрицу А можно также быстро найти и хранить, используя О(К) памяти. Из-за быстрого умножения матрицы А для вычисления формулы обобщенного метода наименьших квадратов удобно и эффективно использовать метод сопряженных градиентов, вычислительная сложность которого равна О(К(К — т)).
Аналогично за счет разреженности матрицы А можно быстро находить множители Лагранжа за время О (К). Таким образом, вычислительная сложность одного шага алгоритма 3 при решении подзадачи (8Ь) равна О(К(К — т)).
Есть две стратегии выбора начальной точки: первая применяется при малом размере задачи и заключается в том, чтобы назначить первым £ точкам, где £ перебирается по целочисленной сетке до |_К/2_|, максимальный вес (то есть (8е)), а оставшимся — минимальный (то есть (8^), после чего выбрать то где достигается наименьшее значение целевой функции ](С).
Вторая стратегия используется при ] = 1 и основана на следующем наблюдении: решения задач при примерно одинаковом отношении N к Ь и одинаковом а схожи. Например, на рис. 1, а, б изображены два отнормированных (умноженных на N) решения задачи 2 при N = 200, Ь = 60, а = 0.1 и N = 1000, Ь = 300, а = 0.1 соответственно.
Схема эвристики следующая: зафиксируем параметр масштаба 0 < ^ < 1 (на практике хорошими значениями 7 являются 0.5-0.7), найдем решение задачи при
2015105-
~Г
"Г
"Г
"Г
"Г
"Г
20 40 60 80 100 120 140
20-
100 200 300 400 500 600 700
Index
Рис. 1. Два решения задачи 2 при одинаковом отношении N к Ь.
N « YN, « а7 = а, после чего используем решение задачи меньшего порядка для выбора начального рабочего множества и начальной точки. Формально, алгоритм выглядит следующим образом.
Алгоритм 4■ Вход: параметр 7.
Результат: Начальное рабочее множество ограничений Шо и стартовая точка Х0.
1. Положить Ky = [7K], Ly = [7L], NY = LY + KY — 1.
2. Определить функцию s(i) = T^rf—— + 1] перехода к уменьшенной размерности.
3. Найти масштабированное решение CY = (с7д,.cy,ky)T задачи (8b) при K = К7, L = Lj, a = a, j = jj = 1.
4. Добавить к рабочему множеству w индексы всех ограничений вида (8c).
5. Для i = 2, 3,..., \К/2\
1) вычислить iY = s(i);
2) если cY,iY — acYjY = 0, добавить в w индекс ограничения вида [ci — acj > 0} из (8d) и перейти к следующей итерации;
3) если cY,iY = cYjY, добавить в w индекс ограничения вида [cj — ci > 0} из (8e) и перейти к следующей итерации.
6. Для к =1, 2,...
1) решить следующую задачу квадратичного программирования, вычислить значение точки минимума C*:
C* = argminCeCfcf (C), Ck = [C | ATC = 0, i e W},
где f (C) определено в (6a);
2) проверить все условия вида (8d) и (8e) для точки C*; если все они выполняются, ответ Xo = C*, wo = w и STOP, иначе добавить к рабочему множеству W индекс любого ограничения, которое не выполняется.
Заметим, что шаг 6 алгоритма 4, необходимый для того, чтобы получить корректную начальную точку, совершит не более \K/2~\ итераций вследствие того, что каждая итерация уменьшает размерность подпространства Ck, в которой лежит начальная точка C*, на единицу.
Таким образом, при использовании данной эвристики получаем рекурсивный алгоритм, который уменьшает размерность задачи в y раз до тех пор, пока K не станет достаточно малым для использования первой стратегии.
6. Численный эксперимент. Ниже приведем численное сравнение эффективности полученных с помощью алгоритма 1 весов в задаче оценки сигнала рядом конечного ранга методом Oblique Cadzow [6]. Сравнение было проведено на примере синусоидального сигнала.
Был взят сигнал S = (si,...,sn) длины N = 40 и ранга r = 2, имеющий вид
Sfc = 5 sin ■
2-кк 6~'
Рассматривался ряд вида X = S + N, где N — гауссовский белый шум с нулевым средним и единичной дисперсией. Точность оценки сигнала S измерялась с помощью корня из среднего по точкам ряда и по 1000 реализациям ряда квадрата отклонения от сигнала S. Эту меру будем называть RMSE (root mean-square error) оценки сигнала. В качестве меры скорости сходимости алгоритма Oblique Cadzow используется среднее число итераций до остановки. Длина окна зафиксирована равной L = 20. Был использован следующий критерий остановки метода Oblique Cadzow: ||T-1(Yk) -T-1(Yk+1)\\2/N < 10-8. Рассматривались веса, полученные алгоритмом 1 (тип I) и прежним методом Cadzow(a) [6, предл. 5] (тип II). Одинаково от-нормированные веса в пространстве рядов при а = 2-2 и а = 2-4 изображены на рис. 2, а, б соответственно.
1.4-
- New approach (I)
- - - Cadzow (а) (П) Unit weights
Index
Рис. 2. Полученные веса в пространстве рядов при различных а: а = 2 2 (а), а = 2 4 (б).
Результаты сравнения приведены в таблице. Заметим, что результаты при а = 1 и а = 2~5, 2~6 совпадают, так как численно равны матричные веса, полученные двумя различными методами. В целом методы дают сравнимые результаты, а при а = 0.5 и а = 0.25 новый подход дает оценку лучше, чем метод Cadzow(a). Также заметим, что заявленные соотношения между а и точностью оценивания, а также а и скоростью сходимости метода Oblique Cadzow, выполняются и для нового метода. При малых а точность оценки существенно лучше, чем при а =1, при этом точность убывает с ростом а.
Сравнение методов по ИМБЕ и среднему числу итераций при различных а и алгоритмах получения весов (I, II)
а RMSE, I RMSE, II Iterations, I Iterations, II
1 0.3811 0.3811 8.44 8.44
2-1 0.3615 0.3682 8.67 8.71
2~'2 0.3475 0.3529 9.76 9.55
0.3371 0.3388 11.97 11.61
2-4 0.3300 0.3293 14.73 15.95
2-ь 0.3247 0.3247 24.37 24.37
2"ь 0.3231 0.3231 39.95 39.95
Заключение. Для решения задачи поиска весов было применено несколько нетривиальных переходов. Сформулированная в терминах минимизации квадратичного функционала задача 2 содержит K2 линейных ограничений, что не позволяет эффективно искать решение на практике. Поэтому была выведена эквивалентная задача 3, решение которой — минимум из решений K/2 задач квадратичного программирования с количеством ограничений порядка K. Преимущество формулировки 3 состоит в том, что перебор задач КП можно остановить заранее, используя разработанный метод остановки, реализованный в алгоритме 2. На практике перебор ограничивается уже на первой итерации, вследствие чего решение находится эффективно.
Качество полученных алгоритмом 1 весов было проверено численным экспериментом. Для длины ряда N, кратной длине окна L, результаты сравнимы с прежним алгоритмом Cadzow(a), но новый алгоритм позволяет находить решение для произвольных N и L.
Литература
1. Golyandina N., Nekrutkin V., Zhigljavsky A. Analysis of Time Series Structure: SSA and Related Techniques. Chapman & Hall/CRC, 2001.
2. Broomhead D., King G. Extracting qualitative dynamics from experimental data // Physica D. 1986. Vol. 20. P. 217-236.
3. Va,uta,rd, R., Yiou P., Ghil M. Singular-Spectrum Analysis: a Toolkit for short, noisy chaotic signals // Physica D. 1992. Vol.58. P. 95-126.
4. Eisner J.B., Tsonis A. A. Singular Spectrum Analysis: a New Tool in Time Series Analysis. Plenum, 1996.
5. Stoica P., Moses R.L. Spectral analysis of signals. Pearson Prentice Hall Upper Saddle River, New Jersey, 2005. Vol. 452.
6. Zvonarev N., Golyandina N. Iterative algorithms for weighted and unweighted finite-rank time-series approximations // Statistics and Its Interface. 2017. Vol. 10. P. 5-18.
7. Gillard J., Zhigljavsky A. A. Stochastic algorithms for solving structured low-rank matrix approximation problems // Communication in Nonlinear Science and Numerical Simulation. 2015. Vol. 21, N 1. P. 70-88.
8. Zhigljavsky A., Golyandina N., Gryaznov S. Deconvolution of a discrete uniform distribution // Statistics & Probability Letters. 2016. Vol.118. P. 37-44.
9. Гавурин М. К., Малоземов В. Н. Экстремальные задачи с линейными ограничениями. Учебное пособие. Л.: Изд-во Ленингр. ун-та, 1984.
10. Nocedal J., Wright S. Numerical optimization. Springer Science & Business Media, 2006.
Статья поступила в редакцию 11 марта 2016 г.
Сведения об авторе
Звонарев Никита Константинович —аспирант; nikitazvonarev@gmail.com
SEARCH OF WEIGHTS IN PROBLEM OF
WEIGHTED FINITE-RANK TIME-SERIES APPROXIMATION
Nikita K. Zvonarev
St. Petersburg State University, Universitetskaya nab., 7—9, St. Petersburg, 199034, Russian Federation; nikitazvonarev@gmail.com
The problem of weighted finite-rank time-series approximation is considered for signal estimation. Finding of weights, which lead to improvement of estimate accuracy, is examined. An effective method for numerical search of weights is constructed and proved. The theory of quadratic optimization is used. For effective algorithm, the problem of quadratic optimization with large number of linear constraints reduced to a sequence of problems with smaller number constraints, which is completed by a stopping criterion. For the algorithm deduction, equivalence of different statements of the original optimization problem is proved. A numerical simulation is performed to approve the algorithm efficiency and improvement of the estimate accuracy. Refs 10. Figs 2. Table 1.
Keywords: time series, time series of finite rank, Cadzow's iterative method, weighted low-rank approximation, oblique singular value decomposition, quadratic programming.
References
1. Golyandina N., Nekrutkin V., Zhigljavsky A., Analysis of Time Series Structure: SSA and Related, Techniques (Chapman & Hall/CRC, 2001).
2. Broomhead D., King G., "Extracting qualitative dynamics from experimental data", Physica D. 20, 217-236 (1986).
3. Vautard R., Yiou P., Ghil M., "Singular-Spectrum Analysis: a Toolkit for short, noisy chaotic signals", Physica D. 58, 95-126 (1992).
4. Elsner J.B., Tsonis A. A., Singular Spectrum Analysis: a New Tool in Time Series Analysis (Plenum, 1996).
5. Stoica P., Moses R. L., Spectral analysis of signals 452 (Pearson Prentice Hall Upper Saddle River, New Jersey, 2005).
6. Zvonarev N., Golyandina N., "Iterative algorithms for weighted and unweighted finite-rank time-series approximations", Statistics and Its Interface 10, 5-18 (2017).
7. Gillard J., Zhigljavsky A. A., "Stochastic algorithms for solving structured low-rank matrix approximation problems", Communication in Nonlinear Science and Numerical Simulation 21(1), 7088 (2015).
8. Zhigljavsky A., Golyandina N., Gryaznov S., "Deconvolution of a discrete uniform distribution", Statistics & Probability Letters 118, 37-44 (2016).
9. Gavurin M.K., Malozemov V. N., Extremal problems with linear constraints (Izdat. Leningrad State Univ., Leningrad, 1984) [in Russian].
10. Nocedal J., Wright S., Numerical optimization (Springer Science & Business Media, 2006).
Для цитирования: Звонарев Н. К. Поиск весов в задаче взвешенной аппроксимации временным рядом конечного ранга // Вестник Санкт-Петербургского университета. Серия 1. Математика. Механика. Астрономия. 2016. Т.3(61). Вып. 4. С. 570-581. DOI: 10.21638/11701/spbu01.2016.406
For citation: Zvonarev N. K. Search of weights in problem of weighted finite-rank time-series approximation. Vestnik of Saint Petersburg University. Series 1. Mathematics. Mechanics. Astronomy, 2016, vol. 3(61), issue 4, pp. 570-581. DOI: 10.21638/11701/spbu01.2016.406
