ПОИСК ВЕСОВ ДЛЯ ЗАДАЧИ ОЦЕНИВАНИЯ СИГНАЛА КОНЕЧНОГО РАНГА В ПРИСУТСТВИИ СЛУЧАЙНОГО ШУМА

Звонарев Никита Константинович

Вестник СПбГУ. Математика. Механика. Астрономия. 2021. Т. 8 (66). Вып. 3 УДК 519.246.8+519.853.4 МЯС 62М10, 90С30

Поиск весов для задачи оценивания сигнала конечного ранга в присутствии случайного шума*

Н. К. Звонарев

Санкт-Петербургский государственный университет,

Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9

Для цитирования: Звонарев Н. К. Поиск весов для задачи оценивания сигнала конечного ранга в присутствии случайного шума // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2021. Т. 8(66). Вып. 3. С. 417-429. https://doi.org/10.21638/spbu01.2021.304

Рассматривается задача взвешенной аппроксимации временного ряда рядом конечного ранга с целью оценивания сигнала в модели «сигнал плюс шум», где обратная ковариационная матрица шума является (2р + 1)-диагональной. Решается задача поиска весов, которые приводят к улучшению точности оценок сигнала. Построен и теоретически обоснован эффективный численный метод поиска весов. Проведено численное моделирование для различных примеров шума, показывающее улучшение точности метода оценивания сигнала.

Ключевые слова: временные ряды, метод Кэдзоу, ряды конечного ранга, взвешенный метод наименьших квадратов, дивергенция Кульбака — Лейблера, методы оптимизации.

1. Введение. 1.1. Сигналы конечного ранга. Рассмотрим задачу извлечения сигнала S = (si,..., sn)t из зашумленного временного ряда X = S + R, где S управляется линейной рекуррентной формулой (ЛРФ) порядка r: sn = aisn-i,

n = r+1,..., N, ar =0, R — шум с нулевым матожиданием и ковариационной матрицей X = Sn G RnxN. Хорошо известен результат, который задает явный вид управляемых ЛРФ рядов в параметрическом виде: sn = i Pi(n) exp(ajn) cos(2nwjn +1i), где Pi(n) — многочлены от n (см. [1, теорема 5.3] и [2]).

Известно, что для выделения такого рода сигналов хорошо работают методы, основанные на оценке сигнального подпространства (subspace-based methods), и, в частности, анализ сингулярного спектра (см. [2-4]). Идея методов состоит в следующем: зафиксируем длину окна L, 1 < L < N, положим K = N — L +1 и определим траекторную матрицу для ряда S G XN (XN — множество вещественных временных рядов длины N):

S = Tl( S) =

/si

S2

52

53

sk \

sk+1

\sl sl+1 ... sn j

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант №20-01-00067).

(¡5 Санкт-Петербургский государственный университет, 2021

где Tl обозначает биекцию между Xn и H, и H — множество ганкелевых матриц порядка L х K с одинаковыми значениями на побочных диагоналях i + j = const.

Скажем, что ряд S, у которого траекторная матрица S имеет неполный ранг r для любого r < L < N — r + 1, называется рядом ранга r. Известно, что определение рядов ранга r эквивалентно равенству ранга r у траекторной матрицы с L = r + 1, если N > 2r [5, следствие 5.1]. Множество временных рядов ранга r длины N обозначим как D>n = Dr.

1.2. Сведение к задаче аппроксимации. Рассмотрим задачу аппроксимации временного ряда рядом ранга r по косоугольной евклидовой норме:

Y* = argminYe—||Х — Y||w, (1)

где X, Y G RN — временные ряды длины N, ||Z|| w — косоугольная евклидова норма, порожденная скалярным произведением

(X, Y)w = XTWY, (2)

W G RN xN — симметричная положительно определенная матрица весов (назовем ее «весовой матрицей»), Т>г — замыкание множества Т>г.

Полученное решение Y* задачи (1) можно использовать как оценку сигнала S рядом конечного ранга. Более того, полученная оценка позволяет решить множество связанных задач: оценка коэффициентов ЛРФ и параметров ряда, прогнозирование ряда, заполнение пропусков и так далее (см. [1]). Заметим, что при W = Х-1 полученная оценка сигнала Y* в предположении о гауссовости шума R является оценкой максимального правдоподобия: логарифм функции правдоподобия в многомерной нормальной модели выглядит как — у log(27r) + ^ det(W) — — (например,

см. [6, формула (8.19)]), где ряд Y удовлетворяет параметрической модели Y € Т>г; выражение с точностью до константных слагаемых, множителя i и знака минус совпадает с (1).

1.3. Проблема подбора весов. Для решения задачи (1) рассмотрим связанную с ней задачу ганкелевой структурной аппроксимации матрицей неполного ранга (Structured Low-Rank Approximation, SLRA):

Y* = argminYeMrпн/ь,н(П /l,r(Y) = ||X — Y||l,r, (3)

где ||X||l,r — порожденная скалярным произведением

(X, Y)l,r = tr(LXRYT) (4)

матричная норма [7], Mr С RLxK — множество матриц ранга, не превосходящего r, L G RLxL, R G RKxK — положительно определенные матрицы весов, а X и Y — матрицы, связанные с задачей (1) соотношениями X = 7L(X) и Y = TL(Y).

Между весовой матрицей W в задаче (1) и матрицами L, R в (3) существует соотношение, при котором нормы || • ||w и || • ||l,r являются эквивалентными [8].

Основная проблема состоит в том, чтобы найти матрицы L, R (назовем их «матричными весами»), которые бы соответствовали матрице Wo = Е-1, что естественно для модели с такой ковариационной матрицей шума и приводит к задаче оценивания сигнала по методу максимального правдоподобия. Для случая авторегрессионного шума (и для белого шума в частности) в [8] показано, что если это

и возможно, то только при вырожденности либо матрицы L, либо R, что недопустимо при решении задачи методами типа Oblique Cadzow (см. [9, замечание 4]). Следовательно, приходится жертвовать точным равенством W и Х-1 и использовать невырожденные матрицы L, R, которые дают матрицу W, приближенную к Е-1. Частный случай этой задачи для модели с белым шумом был разобран в [10].

Цель статьи — рассмотреть подходы к поиску матриц L, R в случае, когда матрица Е-1 является ленточной и содержит (2p + 1) ненулевые диагонали, p > 0. Такие матрицы возникают при нестационарном обобщении авторегрессионной дискретной модели [11]. Важные частные случаи такой модели — классический стационарный процесс авторегрессии порядка p (AR(p)) [12] и модель нестационарного белого шума.

1.4. Подход к подбору весов. Известно, что если взять (2p + 1)-диагональную матрицу L и диагональную матрицу R, то соответствующая им матрица W будет (2p + 1)-диагональной [8]. Этим фактом можно воспользоваться, чтобы упростить подход. В предположении, что матрица L зафиксирована, положительно определена и имеет сходную с матрицей Е-1 структуру, задача поиска весов сводится только к подбору матрицы R. Например, если Е = XN(b), где Ем (b) G RMxM — матрица ковариации стационарного процесса AR(p) длины M с коэффициентами b G то зафиксируем L = E-1(b). Для матрицы R ограничим число обусловленности сверху, чтобы гарантировать ее невырожденность. Как и в работах [9, 10], мы рассмотрим диагональную матрицу R = diag(r1,..., г к ), введем параметр 0 < а < 1 и потребуем, чтобы выполнялось неравенство г, > 0 и

> а, (5)

max,- г,-

что эквивалентно неравенству со^ И. < 1/а. После этого сформулируем задачу поиска нужной матрицы И при помощи минимизации дивергенции Кульбака — Лейблера (являющейся естественной функцией различия между распределениями, задающими требуемую и приближенную модели шума) между многомерными нормальными распределениями с нулевыми средними и ковариационными матрицами Е и W-1.

Полученные в результате матричные веса Ь, И применяются в методе итераций Кэдзоу (Cadzow) [13], который дает приближенное численное решение задачи БЬЯЛ (3). Определим следующую последовательность матриц, начинающуюся с заданной матрицы X:

Хо = X, Хк+1 = ПнПМгХк, к > 0, (6)

Пн и Пмг — проекторы на соответствующие множества по матричной норме У • ||ь,н, порожденной скалярным произведением (4). В [9] показано, что параметр а влияет на скорость сходимости метода Кэдзоу (6): чем больше со^(И) (что равнозначно меньшему значению а), тем медленнее сходится метод.

Ключевое отличие от предыдущего подхода, описанного в [10], состоит в обобщении на случай неединичной матрицы Е и замене целевой функции на выражение, обоснованное с точки зрения математической статистики и определенное для произвольных положительно определенных матриц. Решение задачи сводится к оптимизации гладкой функции в многомерном кубе. Для подобных задач созданы эффективные численные алгоритмы, например алгоритм Ь-БЕСВ-Б [14].

1.5. Структура работы,. Раздел 2 содержит известную информацию о связи между матричными весами и матрицей W, необходимую для изучения задачи поиска весов. В разделе 3 формулируется задача поиска весов как оптимизационная задача и рассматривается выбор матрицы Ь на основе обобщения модели авторегрессии. Раздел 4 посвящен вычислительному аспекту решения задачи поиска весов: задача сводится к оптимизации гладкой функции в многомерном кубе, выведены выражения для вычисления значения и градиента функции в точке. Приведена временная трудоемкость одной итерации для метода оптимизации при различных матрицах Я. Раздел 5 содержит численное сравнение полученных весов с результатом работы [10] по ошибке оценки сигнала с помощью метода Кэдзоу на модельных примерах. В [10] для случая единичных весов показано, что чем ближе веса к требуемым, тем точнее решение задачи. В разделе 5 показано, что эти соотношения выполняются и для неединичных весов, полученных предложенным методом.

2. Соотношение между весами рядов и матричными весами. Рассмотрим соотношение, связывающее (полу)норму ||7||-^, порожденную скалярным произведением (2) и (полу)норму ||Ъ||ь,н, которая порождена скалярным произведением (4) (полунормы получаются, если не требовать невырожденность матриц W, Ь, И.). Для эквивалентности задач (1) и (3) нужно, чтобы отображение 7ъ сохраняло норму, то есть чтобы выполнялись равенства ||7||-^ = ||Ъ||ь,н для любого 7, Ъ = 71(7), где Ь = ), И = (г.,.) — матричные веса, W = (ш.,.) — весовая матрица.

Определим свертку матриц С = А * В для А е К1хЬ, В е ККхК, С е х№:

^^ а»1>Л Ь»2.2,

¿1 +¿2-1=«л +¿2-1=.;

где А = («¿1 ,¿1 ^ В = (Ь.2,.2 ^ С = К.).

Известно соотношение, связывающее матричные веса с весовой матрицей.

Теорема 1 [8]. ||7||-^ = ЦЪЦь.и. для любого ряда 7 е тогда и только тогда, когда W = W(L, И) = Ь * И.

Рассмотрим интересующий нас случай диагональной матрицы И. При этом пусть р обозначает количество ненулевых главных диагоналей матрицы Ь в ее верхнем или нижнем треугольнике, то есть Ь является (2р + 1)-диагональной.

Определим матрицы М. = Ь * Е.,., где Е.,. е ККхК — матрица, содержащая единицу по индексу (г, г) и ноль по всем остальным индексам. Заметим, что М. е х:м содержит матрицу Ь как подматрицу в строчках и столбцах с г-го по (г + Ь — 1)-й, при этом все остальные элементы нулевые. Представим матрицу И как

и = 52 К=1 г.Ем.

Следствие. Рассмотрим (2р + 1)-диагональную матрицу Ь и диагональную И = diag(r), г = (г1,..., гк)т. Тогда весовая матрица W является (2р + 1)-диагональной, W = 52¿=1 г.М., где М. = Ь * Е.,..

3. Эквивалентные формулировки задачи поиска весов. 3.1. Постановка задачи. Сформулируем задачу поиска матричных весов И при фиксированной матрице Ь и ограничении на вырожденность матрицы И, таких что W(L, И) равна заданной обратной матрице ковариаций Wo = Е-1 или приближена к ней.

Во введении было упомянуто, что уже для авторегрессионного случая Wо = Е^(Ь) не существует таких невырожденных Ь, И, что W(L, И) = Ь * И = Wо [8]. Поэтому необходимо поставить задачу поиска Ь и И, которые бы дали матрицу W(L, И), схожую с Wо по некоей мере близости этих матриц.

Рассмотрим разумную функцию различия между обратной ковариационной матрицей Wо = Е-1 и весовой матрицей W. В качестве такой меры выберем дивергенцию Кульбака — Лейблера [15] между многомерными нормальными распределениями с нулевыми средними и ковариационными матрицами Е = W—1 (соответствующей модели реального шума в предположении гауссовости) и W-1 = (W(L, И)) (соответствующей модели с ковариационной матрицей, обратной весовой матрице). По определению дивергенция выглядит следующим образом [16]:

АаХ^Ол^о-1)!!!^,1^-1)) = 1 ИЛУ1^1) - N - 1<^(сИ\¥\¥0-1))) .

С учетом условия на невырожденность матрицы И, ограничим ее число обусловленности: сопё И < 1/а, 0 < а < 1. Получим оптимизационную задачу в следующем виде:

И* = а^шшсопан^/«^^^, W0-1)||N(0N, W-1)). (7)

w=l*r

r=diag(r)

Введем функцию, которая представляет собой упрощение дивергенции Кульба-ка — Лейблера:

^0 = tr(WW01) - ^(ёе^)), (8)

после чего сформулируем следующую задачу:

И* = а^шт^ r<l/aDwo (9)

w=l*r

r=diag(r)

Предложение 1. Задачи (7) и (9) эквивалентны.

Доказательство. Целевые функции задач отличаются только на константные множитель и слагаемое. □

3.2. Выбор левой весовой матрицы Ь. В общем случае выбор подходящей матрицы Ь по заданной матрице Е (или Е-1) — нетривиальная задача. Ниже мы приведем подход, основанный на статье [11] и позволяющий выбрать Ь в некоторых частных случаях.

Используя разложение Холецкого, любую положительно определенную ковариационную матрицу Е, чья обратная является (2р +1)-диагональной, можно записать как

Е = (I - Л)-1Б ((I - Л)-1)Т ,

или через ее обращение:

Е-1 = (1 - Л)ТБ-1(1 - Л), (10)

где Л — нижнетреугольная (р+1)-диагональная матрица с нулевой главной диагональю, Б — положительно определенная диагональная матрица. Пусть X — случайный вектор с нулевым средним и ковариационной матрицей Е. Тогда его можно представить в виде X = (I - Л)-1 V = ЛХ + V, где V — случайный вектор с нулевым средним

и ковариационной матрицей О, независимый с X [11]. Добавление в представление ненулевой матрицы А дает обобщение процесса авторегрессии на «динамический» случай, V — вектор некоррелированных ошибок.

В случае, когда Я соответствует классическому стационарному процессу авторегрессии ЛИ(Ъ), матрица А = А(Ь) является теплицевой, значения на диагоналях соответствуют значениям вектора коэффициентов Ь = (61,..., 6р)т:

а. = (6.,...Л)т е -\ г =1,...,р, (11)

где а. обозначает р-ю главную диагональ матрицы А. Это следует из определения =52Р=1 + е^, е^ — независимые одинаково распределенные величины с

нулевым средним и дисперсией а2, X = (ж )т, при этом добавляется тре-

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

Исходя из этого, мы можем подбирать Ь в представлении, аналогичном (10):

Ь =(1Ь — Аь)тБ-1(1ь — Аь), (12)

где А_1, — матрицы порядка Ь х Ь, имеющие ту же структуру, что и соответствующие матрицы А и О.

В частности, если матрица А теплицева, то матрицу А^ можно тоже выбрать (р + 1)-диагональной нижнетреугольной с теми же значениями на диагоналях. Заметим, что случай теплицевой матрицы А может соответствовать и нестационарному процессу авторегрессии, если хотя бы один корень характеристического полинома с коэффициентами Ь лежит внутри единичного круга. Случай (возможно, нестационарного) белого шума соответствует нулевой матрице А. Если матрица О содержит все одинаковые элементы, кроме, возможно, первых р, то в качестве можно взять диагональную матрицу, у которой на диагонали находятся первые Ь элементов диагонали матрицы О.

Замечание 1. Такой выбор матриц А^ и для случая стационарного процесса авторегрессии эквивалентен выбору Ь = Е-1(Ь) — обратной автоковариационной матрицы процесса авторегрессии длины Ь с теми же коэффициентами авторегрессии Ь, что и у Е, где Е-1(Ь) можно напрямую вычислить по формулам, приведенным

в [12].

Таким образом, мы обобщили выбор матрицы Ь на класс ковариационных матриц Е, которым соответствует теплицева матрица А и для которых диагональ матрицы О содержит равные элементы после первых р элементов. Помимо стационарных, такой класс включает и нестационарные процессы авторегрессии, когда нестационарность достигнута за счет выбора коэффициентов полинома Ь.

Остается непокрытым случай, когда матрица А не является теплицевой или диагональ матрицы О содержит разные элементы и вне первых р (например, случай нестационарного белого шума). Без дополнительной информации можно лишь предложить эвристики, которые точно дадут положительно определенную матрицу Ь. К примеру, можно выбрать А^ теплицевой, где элементы на диагоналях равны усредненным значениям элементов диагонали А, а матрицу выбрать соответствующей таковой при стационарном варианте, если матрица А теплицева и коэффициенты Ь отвечают стационарному случаю, либо единичной, если подобные условия не выполнены.

4. Задача поиска весов как задача минимизации гладкой функции. 4-1. Подход к оптимизации. Основная трудность при решении задачи (7) и эквивалентной ей задачи (9) — набор ограничений на элементы матрицы R. Сформулируем эквивалентную форму задачи (9), в которой вместо вектора г, R = diag(r), используется вектор Г, оптимизация по которому происходит в многомерном кубе. Рассмотрим следующую задачу:

Г* = arg min )T D* (Г), где D*(r) = - log(det(L * diag(r))) + N log (r, q(W-1)),

fi <1,'i='l,...,K fi>a, i=1,...,K

(13)

q(Z) = (q1 (Z),..., qK(Z))T : qi(Z) = tr(MiZ) = (Mi; Z)F — вектор, определенный для симметричной положительно определенной матрицы Z G RNxN, (•, -)f и (•, •) обозначают соответственно фробениусово и обычное евклидово скалярные произведения.

Теорема 2. Задачи (9) и (13) эквивалентны.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Доказательство. Параметризуем вектор г = сГ с учетом ограничения на число обусловленности матрицы R, где c > 0 — свободный множитель, при этом ri < 1, ri > a, i = 1,..., K, Г = (r1,..., гк)T. Тогда задачу (9) с учетом следствия теоремы 1 можно переписать следующим образом:

Г* = argminfi<1, i=1,...,K V cfitr(MiW-1) - log ( det ( V cTiMj J . (14) fi >a,i=i,...Kt1 \ \t1 ))

Преобразуем целевую функцию из задачи (14):

V cfitr(MiW-1) - log |det cfiM^ =

= с(Г, q(W-1)) - N log с - 1og(det(L * diag(Г))).

Заметим, что легко выполнить оптимизацию по множителю с, если найти выражение, при котором целевая функция достигает минимума. Достаточно взять у выражения производную по с и приравнять ее к нулю. Получим оптимальное зна-тт-. Подстановкой с = с* в (14) и сокращением константного

N

(г.ч^-1))'

слагаемого получим выражение (13). □

Предложение 2. Функция ^|(г) из (13) является гладкой, ее градиент определяется 'равенством

{ПШ' = * с!1аё(г))-1) +

(г, q(Wo ))

Доказательство. Достаточно продифференцировать выражение (13) по г, используя известную формулу logdetZ = trlogZ, где Z = Ь * diag(r) € КДГхЛГ. □

Обоснуем необходимость введения задачи поиска весов в форме (13). Согласно утверждению теоремы 2 и предложений 1 и 2, исходная задача минимизации (7)

может быть записана в эквивалентной форме оптимизации гладкой функции в многомерном кубе. Поэтому для численного решения задачи (13) можно использовать, например, достаточно эффективный численный метод Ь-БЕСВ-Б [14], разработанный для функций такого вида.

4-2. Исследование трудоемкости. Для работы метода требуется уметь вычислять как значение целевой функции Д|(г) в заданной точке, так и ее градиент.

Пусть Wo и Ь являются (2р +1)-диагональными матрицами, где р фиксировано. Рассмотрим асимптотическую сложность вычисления функции Б? (Г) и ее градиента по N. Применим следующую вычислительную схему.

Предварительно вычислим q(W—1), что сводится к обращению матрицы Wo и вычислению К фробениусовских скалярных произведений. Это потребует временной сложности О^2), так как (2р + 1)-диагональную матрицу можно обратить за время О^2) с использованием разложения Холецкого [17]. Дальнейшее вычисление К скалярных произведений с М1,..., Мк потребует О^2) времени, так как каждая из матриц М. содержит О^) ненулевых элементов, и К < N.

Вычисление свертки Ь * diag(f) = 52¿=1 г.М., которая является (2р + 1)-диагональной матрицей, тоже занимает О^2) времени (сложение К матриц, каждая из которых содержит О^) ненулевых элементов). Взятие логарифма ее определителя делается через разложение Холецкого и занимает О^) времени. Трудоемкость вычисления градиента, требующего нахождения —q((Ь * diag(f))-1), эквивалентна трудоемкости вычисления q(W-1).

Таким образом, временная сложность вычисления функции и ее гради-

ента в случае (2р + 1)-диагональных матриц весов составляет О^2), что допустимо для многих практических приложений.

Следующее замечание устанавливает порядок трудоемкости для частного случая р = 0.

Замечание 2. В случае диагональных матриц Wo и L асимптотическая сложность вычисления функции (г) и ее градиента может быть снижена до O(N log N) за счет применения преобразования Фурье для вычисления L * diag(r). В случае, когда L = II, асимптотика может быть снижена до O(N) за счет вычисления сверток с использованием накопленных сумм.

Таким образом, был получен метод поиска весов, эффективный с вычислительной точки зрения и учитывающий вид ковариационной матрицы шума.

5. Численный эксперимент. Ниже приведем численное сравнение эффективности полученных с помощью решения задачи (13) весов и метода нахождения весов, упомянутых в работе [10]. Сравнение было проведено на примере синусоидального сигнала.

Предложенный в [10, алгоритм 3] метод (в дальнейшем назовем его «квадратичный метод») решал иную задачу. Во-первых, там использовалась минимизация квадратичного функционала вместо дивергенции Кульбака — Лейблера. Во-вторых, метод никак не учитывал вид матриц Ь и Wo — можно считать, что веса искались в предположении, что Ь = и Wo = являются единичными матрицами соответствующих размеров. Новый подход учитывает форму данных матриц.

Для того чтобы изучить разницу между полученными весами, был взят сигнал S = (si,..., sn) длины N = 40 и ранга r = 2, имеющий вид

sfc = 5sm-, k=l,...,N,

6

и рассматривались три ряда Xj = S + Rj, Rj — случайный гауссовский шум с нулевым средним, i = 1, 2, 3. В качестве первой модели шума Ri была взята стационарная авторегрессия с нулевым средним, p =1, b = (0.9)T, а = 0.1. Варианты R2 и R3 используют нестационарную авторегрессию. Согласно представлению ковариационной матрицы (10), при теплицевости матрицы A нестационарности можно добиться как выбором коэффициентов b характеристического полинома b, где A и b связаны соотношением (11), так и выбором матрицы D с неравными коэффициентами. Для варианта R2 возьмем b = (0.9)T, D = 0.01 • diag(12, 22,..., 402), нестационарность достигается за счет матрицы D. Для шума R3 возьмем b = (1)T, D = diag(0.01,..., 0.01) в представлении (10) матрицы Я, что соответствует модели броуновского шума, который представляет из себя нестационарную авторегрессию за счет выбора коэффициентов b.

Точность оценки сигнала S измерялась с помощью среднего по точкам ряда и по I = 10000 реализациям ряда квадрата отклонения от сигнала S:

11 1 ~

j=i

где Sj получена с помощью применения алгоритма Кэдзоу (6) к j-й реализации ряда Xj. Эту меру оценки сигнала будем называть MSE (mean-squared error). Длина окна зафиксирована равной L = 20. В каждом из вариантов у метода Кэдзоу было произведено 40 итераций.

Рассматривались следующие три типа весов. Тип I соответствует методу Oblique Cadzow из работы [10], в нем L = Il, R получена квадратичным методом, то есть информация о виде матрицы Я никак не используется. Тип II — промежуточный вариант, в котором R также получена квадратичным методом без использования матрицы Я, но для метода Кэдзоу матрица L выбрана уже не единичной, а дву-диагональной согласно подразделу 3.2. Тип III соответствует этой работе, в нем матрица L выбрана согласно подразделу 3.2, для выбора R используется решение задачи (13).

Результаты сравнения различных типов весов для рядов Xi и X2 представлены на рис. 1, для ряда X3 — на рис. 2.

Для всех трех моделей шума тип весов I провалился и показал большую ошибку (в случае рядов X2 и X3 даже возрастающую при уменьшении а). Это связано с тем, что в нем никак не учитывается авторегрессионная природа шума. Тип II показал очень близкий к типу III результат в случае стационарной авторегрессии из-за схожести полученных в результате поиска матричных весов, что видно на рис. 1, a. Однако в нестационарных случаях X2 и X3 тип II отличается от типа III. На рис. 1, б видно, что полученные им оценки хуже, чем у типа III при всех а < 1, а на рис. 2, б — что ошибка начинает возрастать при уменьшении а на промежутке а < 0.05 и становится значительно хуже, чем у метода III. Это объясняется тем, что квадратичный метод поиска весов никак не учитывает информацию о виде ковариационной матрицы Е.

Рис. 1. Сравнение методов по MSE при различных а и различных алгоритмах получения весов (I, II, III) для рядов Xi (а) и X2 (б).

а

0.05 0.10 0.20

--о-.о-о-о^в"0

I I I I

0.01 0.02 0.05 0.10 0.20 0.50 1.00

Рис. 2. Сравнение методов по МЯЕ при различных а и различных алгоритмах получения весов для ряда Хз.

a

а

6. Заключение. Таким образом, в работе удалось улучшить метод, предлагаемый в [10], тем самым еще больше увеличив преимущество взвешенных вариантов метода Кэдзоу по сравнению со стандартно используемым невзвешенным вариантом [13, 18]. А именно, с помощью нового подхода к поиску весов, состоящего в решении задачи (7) по заданной обратной ковариационной матрице шума Х-1, удалось получить более точные, чем в работе [10], оценки сигнала конечного ранга по зашумленному ряду, что подтверждено численным экспериментом. Также удалось расширить область применения метода Кэдзоу на случай модели авторегрессионного шума, в том числе и на нестационарный случай.

Для решения задачи поиска весов был применен обобщенный подход, сводящийся к поиску весовой матрицы, близкой к обратной от требуемой ковариационной матрицы шума в смысле дивергенции Кульбака — Лейблера между соответствующими распределениями. Принципиальное отличие от предыдущего подхода, изложенного в [10], состоит в том, что вид ковариационной матрицы учитывается при поиске весов. Для исходной задачи (7) построена эквивалентная задача (13), которая является задачей оптимизации гладкой функции в многомерном кубе, а для решения подобных задач существуют эффективные методы оптимизации. Одна итерация оптимизационного метода имеет временную трудоемкость 0(Ж2) в слу-

чае (2p + 1)-диагональной матрицы Е 1, но может быть сокращена до O(N log N) и O(N) в случае диагональной и единичной матриц соответственно.

Литература

1. Golyandina N., Nekrutkin V., Zhigljavsky A. Analysis of Time Series Structure: SSA and Related Techniques. Chapman & Hall/CRC (2001).

2. Golyandina N., Zhigljavsky A. Singular spectrum analysis for time series. 2nd ed. Berlin, Heidelberg, Springer (2020).

3. Van Der Veen A.-J., Deprettere E. F., Swindlehurst A.L. Subspace-based signal analysis using singular value decomposition. Proceedings of the IEEE 81 (9), 1277-1308 (1993).

4. Gonen E., Mendel J.M. Subspace-based direction finding methods. In: Madisetti V. K., Williams D. B. (eds.) The Digital Signal Processing Handbook, vol. 62. Boca Raton, CRC Press LLC (1999).

5. Heinig G., Rost K. Algebraic Methods for Toeplitz-like Matrices and Operators. In: Operator Theory: Advances and Applications. Birkhauser, Verlag (1985).

6. Adachi K. Matrix-based introduction to multivariate data analysis. Springer (2016).

7. Allen G.I., Grosenick L., Taylor J. A generalized least-square matrix decomposition. Journal of the American Statistical Association 109 (505), 145-159 (2014).

8. Golyandina N., Zhigljavsky A. Blind deconvolution of covariance matrix inverses for autoregressive processes. Linear Algebra and its Applications 593, 188-211 (2020).

9. Zvonarev N., Golyandina N. Iterative algorithms for weighted and unweighted finite-rank time-series approximations. Statistics and Its Interface 10 (1), 5-18 (2017).

10. Звонарев Н. К. Поиск весов в задаче взвешенной аппроксимации временным рядом конечного ранга. Вестник Санкт-Петербургского университета. Серия 1. Математика. Механика. Астрономия 3 (61), вып. 4, 570-581 (2016). https://doi.org/10.21638/11701/spbu01.2016.406

11. Wiesel A., Globerson A. Covariance estimation in time varying ARMA processes. IEEE 7th Sensor Array and Multichannel Signal Processing Workshop (SAM) / IEEE, 357-360 (2012).

12. Verbyla A. A note on the inverse covariance matrix of the autoregressive process. Australian & New Zealand Journal of Statistics 27 (2), 221-224 (1985).

13. Cadzow J. A. Signal enhancement: a composite property mapping algorithm. IEEE Trans. Acoust. 36 (1), 49-62 (1988).

14. Byrd R. H., Lu P., Nocedal J., Zhu C. A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing 16 (5), 1190-1208 (1995).

15. Kullback S., Leibler R. A. On information and sufficiency. The annals of mathematical statistics 22 (1), 79-86 (1951).

16. Duchi J. Derivations for linear algebra and optimization (2007). Доступно на: https://web.stanford.edu/~jduchi/projects/general_notes.pdf (дата обращения: 01.04.2021).

17. Golub G. H., Van Loan C. F. Matrix computations. 3rd ed. Baltimore, MD, USA, Johns Hopkins University Press (1996).

18. Gillard J. Cadzow's basic algorithm, alternating projections and singular spectrum analysis. Stat. Interface 3 (3), 335-343 (2010).

Статья поступила в редакцию 11 декабря 2020 г.;

после доработки 11 января 2021 г.; рекомендована в печать 19 марта 2021 г.

Контактная информация:

Звонарев Никита Константинович — канд. физ.-мат. наук; nikitazvonarev@gmail.com

Search of weights in the problem of finite-rank signal estimation in presence of noise*

N. K. Zvonarev

St. Petersburg State University, 7—9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation

For citation: Zvonarev N. K. Search of weights in the problem of finite-rank signal estimation in presence of noise. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2021, vol. 8(66), issue 3, pp. 417-429. https://doi.org/10.21638/spbu01.2021.304 (In Russian)

The problem of weighted finite-rank time-series approximation is considered for signal estimation in "signal plus noise" model, where the inverse covariance matrix of noise is (2p + 1)-diagonal. Finding of weights, which improve the estimation accuracy, is examined. An effective method for the numerical search of the weights is constructed and proved. Numerical simulations are performed to study the improvement of the estimation accuracy for several noise models.

Keywords: time series, Cadzow's iterative method, time series of finite rank, weighted low-rank approximation, Kullback — Leibler divergence, optimization methods.

References