CC BY
26Быстрый алгоритм спектрального оценивания с долговременной
численной устойчивостью
В.В.Захаров ([email protected]), Фаусто Каско ([email protected]), Departamento de Ingenieria Eléctrica, Universidad Autonoma Metropolitana-Iztapalapa, Mexico D.F., MEXICO
Синтезирован эффективный рекурсивный алгоритм наименьших квадратов (РНК-алгоритм) для спектрального оценивания временных рядов, позволяющий в отличии от известного алгоритма, повысить его долговременную численную устойчивость за счет использования в алгоритме вычислительной процедуры QR разложения корреляционной матрицы с помощью метода Хаусхолдера. Показано, что использование синтезированного алгоритма расширило диапазон числа обусловленности корреляционной матрицы входной последовательности, при котором сохраняется долговременная численная устойчивость алгоритма, что позволило избежать эффекта искажения спектральных характеристик, возникающего в реальных вычислительных устройствах за счет конечной длины машинных слов.
I. ВВЕДЕНИЕ
Основной задачей спектрального анализа временных рядов является получение достоверных спектральных оценок случайного процесса по его единственной реализации конечной длительности. Из существующих моделей временных рядов особое место занимают авторегрессионные (АР) спектральные оценки. Это объясняется тем, что АР оценки спектральной плотности мощности (СПМ) имеют острые пики, и следовательно более высокое разрешение гармонических составляющих, чем классические методы спектрального оценивания [1]. Рассмотрим АР-модель для оценки линейного предсказания вперед отсчета xn , n = 1,2,...,N+p
Р
xn =-X akxn - к , (1)
к=1
где Xn - оценка отсчета xn , ak - коэффициент линейного предсказания вперед, N -
количество отсчётов временного ряда , p -порядок AP- модели .
Если минимизировать средний квадрат ошибки линейного предсказания
(xn - xn )2 =A2n = тт, (2)
тогда вектор коэффициентов линейного предсказания ak, который минимизирует A^
находится как решение нормальных уравнений. Для этого уравнение (1) с учетом (2) представим в виде
p
A n = xn +Х apxn - к . (3)
к=1
Если учесть ,что n = 1,2,...,N+p и xn =0 при n<1 и n>N, тогда выражение (3) запишем в матричной форме
A = X • A, (4)
где А = А: ,А 2 ,...,А
N+р
Т
А
1, а^, ^2,..., а р
, Т-символ транспонирования,
X - матрица входных отсчетов, размера N+р)х (р+1) со следующей структурой
X
Хл
Хл
о
о
Х
о
о о
ХN ХN-1 ХN-2
о о
Х
о о
N - р
Х
N
(5)
Как показано в [1], уравнение минимизирующее средний квадрат ошибки имеет
вид
ЯЛ = Р.
(6)
где Р = р,о,...,о - вектор размера (р+1)
N+р
Р= 1|А = А~ А,
(7)
п= 1
~ - символ сопряжения и транспонирования, Я = XX - теплицева корреляционная матрица входных отсчетов размера (р+1)х(р+1).
Уравнение (6) по своей структуре является уравнением Юла-Уолкера для автогрессионного процесса и может быть решено относительно коэффициентов а, г=1,2,...,р различными известными алгоритмами, например алгоритмом Левинсона [1]. Однако при обработке отсчетов в реальном времени предпочтительнее являются методы последовательного оценивания АР-параметров, что позволяет после поступления очередного отсчета обновить параметры и график спектральной оценки. При этом спектральная плотность мощности входной последовательности для полученных параметров определяется выражением [1]
Тр
1 +1
п=1
апе
- ]2п/пТ
(8)
е
- ] 2п/пТ
комплексная синусоида частоты /п, Т- интервал
где ап, оценка параметра ап дискритизации.
Такая процедура позволяет отслеживать медленно изменяющиеся параметры нестационарных последовательностей за счет введения некоторого экспоненциального окна, которое движется вдоль записи данных, создавая наименьшее изменение значений текущих ошибок и очень сильно уменьшая значения более старых ошибок, что позволяет уменьшить ошибки за счет нестационарности. Таким образом, выражение (7) представим в виде
N+р
р= X И+р-п|А,
п = 1
где и — положительная действительная скалярная величина, удовлетворяющая условию 0 < и < 1.
Соответственно выражение для матрицы Я запишем в виде
N + р
* =1
^+р—п X и X;, (10)
где X п =
Хп'Хп—1 '■■■'Хп—р
Т
п=1
■ п-ая строка матрицы X, п=1,2,...,К
Одними из наиболее простых процедур оценивания текущих параметров АР-модели являются алгоритмы наименьших средних квадратов (НСК—алгоритмы) такие как градиентные процедуры наискорейшего спуска, например алгоритм Гриффитса [2]. Однако сильная зависимость скорости сходимости этих алгоритмов от характеристик входной последовательности (в частности от обусловленности корреляционной матрицы) не позволяет использовать эти алгоритмы для слежения за параметрами АР—модели, если изменение во времени статистик входной последовательности происходит достаточно быстро, по сравнению со скоростью сходимости этих алгоритмов. В связи с этим, более широкое распространение получили рекурсивные алгоритмы наименьших квадратов (РНК-алгоритмы), позволяющие получить сходимость к оптимальной оценке параметров АР—модели за значительно меньшее число шагов, чем алгоритмы НСК, а также независимость от обусловленности корреляционной матрицы входной последовательности, за счет дополнительных вычислительных затрат. Однако высокая чувствительность к плохой численной обусловленности и вычислениях с конечной длиной машинного слова может привести к тому, что после выполнения определенного числа итераций дисперсия ошибки на выходе фильтра предсказания может резко возрасти (другими словами фильтр не обладает долговременной численной устойчивостью), что может привести к существенным искажениям спектральных оценок.
В [1],[2] для снижения влияния плохой численной обусловленности корреляционной матрицы на получаемые решения к входной последовательности предлагается добавить "белый" шум, или в диагональ матрицы добавить небольшую константу (регуляризовать матрицу). Однако, такие методы приводят к смещению получаемых спектральных оценок, что в некоторых случаях недопустимо. Другим методом повышения численной устойчивости является выбор и<<1, что позволяет резко ограничить величину накопленной ошибки на предыдущих итерациях. Однако это приводит к ухудшению точности вычисляемых значений параметров предсказания (другими словами параметр w обязан отслеживать нестационарность, которая никак не связана с обусловленностью).
По видимому самым перспективным путем повышения численной устойчивости рекурсивных алгоритмов является синтез алгоритмов, использующих более устойчивые к ошибкам округления вычислительные процедуры, чем используемая в известных рекурсивных алгоритмах РНК формула Вудбери для обращения однорангово модифицированных матриц [1] [2].
Фундаментальными трудами по изучению влияния эффектов конечной разрядности в различных вычислительных процедурах являются работы [3] ,[4], из которых следует, что наиболее устойчивыми к вычислительным ошибкам являются методы триангуляризации Хаусхолдера и модифицированный алгоритм Грамма-Шмидта, затем следуют метод Холецкого и метод вращения Гивенса.
Целью настоящей работы является синтез алгоритма спектрального оценивания с линейной вычислительной сложностью имеющего низкую чувствительность к плохой
обусловленности корреляционной матрицы входных отсчетов и обладающего долговременной численной устойчивостью.
Для этого предлагается рекурсивная процедура спектрального оценивания в РНК-алгоритме с использовать наиболее устойчивой к ошибкам вычислений процедуры Хаусхолдера [5].
II. СИНТЕЗ АЛГОРИТМА
Для синтеза алгоритма перепишем уравнение (6) в виде
ЬЬ~Л = Р, (11)
где Ь и Ь соответственно разложение матрицы Я на нижнюю и верхнюю треугольные, такие что Я = ЬЬ .
Из (11) следует, что вектор А может быть найден решением двух систем линейных уравнений
ЬВ = Р, Ь~Л = В. (12)
Известен алгоритм получения матриц Ь и Ь с помощью алгоритма Холецкого [4], однако в этом случае улучшение численной устойчивости незначительное, а матрица Я должна формироваться в явном виде, то есть. под нее должна быть зарезервирована память, после чего может быть произведена ее факторизация. Это требует увеличения вычислительных ресурсов.
Предлагаемый алгоритм спектрального оценивания позволяет произвести рекурсивное обновление коэффициентов АР и спектральной оценки, используя
устойчивую к ошибкам вычислений процедуру получения матриц Ь и Ь , а также
позволяет не формировать матрицу Я в явном виде, а элементы матрицы Ь, Ь получать непосредственно из матрицы входных отсчетов X.
Допустим, что отсчеты входной последовательности поступают в реальном времени. Необходимо после поступления очередного отсчета произвести корректировку СПМ входной последовательности. Представим алгоритм в виде последовательного выполнения следующих шагов.
Шаг 1. Получение входных отсчетов в реальном времени хп, п= 1,2,..
Шаг 2. Формирование и взвешивание матрицы входных отсчетов X, размера пх(р + 1)
X=WX',
(13)
где W =
wn-1 мп - 2,
2 ,w1
X =
Хц
Х-
Х
п -1 Х
0
Х
Х
п-2
Х
п -1
Х
о о
п- р-1
Х
п-р
(1 4)
где п = 1,2, ..., N +р номер текущего елемента. Шаг 3. Вычисление матрицы Ь.
Шаг 4. Решение системы уравнений (12) и вычисление вектора А. Шаг 5. Получение значения СПМ согласно выражению (8).
Рассмотрим рекуррентную процедуру обновления элементов матрицы Ь и решения системы уравнений (12). Предположим, что матрица X сформирована по п
комплексным отсчетам данных. Тогда представим матрицу Ьп в виде
0
= 0Х
(15)
где Q - ортогональная матрица Хаусхолдера размерности п хп [5].
Из (15) следует, что элементы матрицы Ь можно получить с помощью ортогонального преобразования над матрицей X, не формируя в явном виде матрицу Представим матрицу Q в факторизованном виде
Ь
0
Qp+lQp■■■Qlх
I—1
(16)
где Х°=Х, X1 1 = Qi_!. . ОуX = Х{ 1:...:X^ 1 :...:Хр?+11 ,/=1,2,...,р+1,при этом столбцы
X) 1 матрицыX~1J/= 1,2,...,р+1 прошли (/-1) раз ортогонализацию матрицами Ql,Q2,■■■,Qi-l■
Л-1 ■
Матрицы Qi в соответствии с [5] определяются следующим образом
а = 1п
и/С/
/=1,2,...,р+1,
( 1 7)
где 1п - единичная матрица размера п .
В (17) скалярная величина С/ определяется как
С/ = X/ _l~X\—1, /=1,2,...,р+1,
(18)
где Xli 1 = 0,0,... ,0,х/ /1 ,Х/ /+1 ,...,х\ п+/+1 - /-й столбец матрицы X-1. Вектор столбец V в (17) размерности п определяется выражением
V;
X/ ^ +дв;С1
( 1 9)
где
5 =
1, если Яе -1, если Яе
к1} о,
к1 }< о/
Яе{4"'}
- обозначает действительную часть
элемента Х, , е, = 0,...,0,1,0,...,0 - вектор-столбец, размерности п с единицей на /-м месте.
Величина и, в (18) является /-м элементом вектора и,. С учетом (16), 1-й шаг ортогонализации столбцов матрицы X' 1 запишем в виде
Xй = a■x/-l.
(20)
Представим выражение (20) с учетом (17) в виде операций над векторами
( тттт~ ^ ТТ~ Xй-1
/-1
}
X) =
1п
и,и~
и С,
V
1^1
X1- = X)-1
и 1 7
/
и С
) = /, /+1, ... ,р+1. (21)
С учетом (19) представим (21) в виде
Хл Х -1 Лл Хи +С/ /
-- = -- - ---- 0 1{о]1 + СХЦ
X) XГ1 X)-1
где Ол = Xlrl~X¡-1.
Скалярная величина 0, в (22) определяется следующим выражением
1
0 •
и, С,
(23)
(22)
Выражение (22) определяет векторный алгоритм получения элементов Iу и ортогонализации векторов X) с помощью преобразования Хаусхолдера.
Исходя из (22) запишем выражение для элементов ¡у матрицы Ь в виде
¡у = Х}/ = Х)' (х// + С/ )0/ (^л + С/Хл )'
(24)
где )=1,2,...,р+1,)=),...,р+1.
Выражение для ортогонализации векторов запишем в виде
X) = X1- - X/-10 1(в1 + С Х-). (25)
Таким образом, вычислительный процесс рекуррентного получения элементов 1 матрицы Ь с использованием преобразования Хаусхолдера на 1-й итерации (/=1,2,...,р+1) состоит из следующих этапов :
• Формирование скаляра С, в соответствии с выражением (18);
• Формирование вектора и, согласно (19);
• Вычисление ¡у согласно (24);
• Ортогонализация векторов X• в соответствии с (25).
Как следует из выражений (18),(19),(24),(25), выполнение одного шага векторного алгоритма требует приближенно 2п операций умножения и 2п операций сложения. Далее покажем что при получении очередного входного отсчета требуется для обновления матрицы Ь не более, чем 6 операций умножений и 7 операций сложения.
Допустим, что через некоторый промежуток времени Т получен (п+1)-й отсчет. Тогда матрица X на (п+1)-й итерации представится в виде
X =
Хц 0 • 0
х 2 х1! 0
Х п—1 хп-2! хп-Р-1
хп хп-1! хп-р
хп+1 хп! хп- Р+1
(26)
т.е. в матрице добавляется строка (выделена рамкой).
Очевидно, если на п-й итерации по п отсчетам были получены величины Сг, и, ¡'¡,
X., тогда при поступлении (п+1)-го отсчета эти величины могут быть легко
модифицированы на (п+1)-й итерации, исходя из следующих соотношений
Сг
'1(п+1)
С'(п) + хп+2-г хп+2-' , ' = 1-2- - -- р + 1
&р(п+1) = Ср(п) + хП+2-;хп+2-г, . = 1>2>--->р + 1
(27)
(28)
где хП+2- - (п+2-1)-й елемент вектора X] 1.
X;
' -1
и
'( п+1)
х
г-1
п+2-г
+ §е'С'(п+1) ;
(29)
¡'¡(п+1 )= х'п1 (х/'' С1(п+1) ) 0'(п+1) (р¡'(п+1)+ С1(п+1 )х'¡'1); (30)
X1.
где 0'
¡(п+1) 1
г- X!-1
х'-п+2-' х'- п+2-'
г(п+1 )
и'С'(п+1)
0'(п+1) ^¡'(п+1) + С'(п+1 )х¡'
' -1
(31)
Как следует из выражений (27)-(31), для получения элементов ¡¡(п+1) из елементов ¡¡(п) необходимо выполнить 6 операций уменожения и 7 - сложения. Если принять во
внимание, что матрица Ь имеет р /2 элементов, тогда для модификации всей матрицы потребуется 3р2 умножений и 3.5р2 сложений.
По мере получения элементов ¡у матрицы Ь, могут быть также получены элементы векторов В и А из выражения (12).
Ь1 =р //п , (32)
Ь) =
X ¡л
п=1
/)), ) = 2,3,..., р +1
а
= Ь
р+1 и р+1
/
р+1 , р+1
(33)
(34)
а
Ь
л
X1 )пап
/у, ) = р,р -1--1. (35)
п=р+1
Получаемые из выражений (34), (35) коеффициенты используются для обновления спектральных оценок в соответствии с (8). Для решение систем уравнений (32), (33) и (34),(35) требуется дополнительно (р2 + р) операций умножений и сложений. Таким образом, общая вычислительная сложность алгоритма составит
4.5р + р операций
умножения-сложения .
3. ТЕСТИРОВАНИЕ АЛГОРИТМА НА ЧИСЛЕННУЮ УСТОЙЧИВОСТЬ
Для сравнительного известного быстрого алгоритма и синтезированного было проведено их моделирование на ЭВМ при 32-разрядной арифметике и выполнении операций с плавающей запятой. Моделировалось два случая. В первом оценивался средний квадрат ошибки на выходе фильтра (р=10) при подаче на его вход последовательности данных, состоящей из одной "мощной" гармонической составляющей и аддитивного "белого" шума. При этом отношение Ргс/Рш равнялось 60 db, где Ргс -мощность гармонической составляющей, Рш - мощность шума. Этот случай следует расценивать как "худший", поскольку максимальное собственное число матрицы Я в данном случае будет приближенно равно Ргс, а минимальное - Рш, и, следовательно, число обусловленности матрицы равно Ргс/Рш. Результаты зависимости р / р ^ от числа
отсчетов п при числе испытаний равном 100 приведены на рис.1 (кривые 1), где р г
получено из непосредственного решения уравнения (6) с помощью алгоритма Левинсона. Во втором случае последовательность состояла из 8 гармонических составляющих
равной мощности, так что X Рг,
/=1
Рш = 60йЬ. Такой случай следует расценивать, как
"средний" с точки зрения численной устойчивости, поскольку при полностью "разрешенных" спектральных составляющих число обусловленности матрицы будет пропорционально величине Ргс/Рш, что ниже, чем в худшем случае. Результаты приведены на рис. 1 (кривые 2).
Как следует из приведенных результатов моделирования, при фиксированной разрядности машинного слова (что имеет место в реальных вычислительных устройствах) наблюдается повышение дисперсии ошибки на выходе фильтра с увеличением числа отсчетов входной последовательности (штриховая линия). Такой эффект непременно приведет к искажению формы спектра при использовании рекуррентных алгоритмов в
спектральном оценивании, что является нежелательным и практически компенсирует преимущества рекуррентных алгоритмов.
Результаты использования синтезированного алгоритма (сплошная линия, рис.1) показывает, что этот эффект можно частично исключить, применяя более устойчивую процедуру нахождения параметров авторегрессии.
В частности, для приведенных примеров синтезированный алгоритм позволяет сохранять долговременную численную устойчивость вплоть до п= 10000, что лежит в пределах разумной максимальной длины временных рядов различной природы.
Рис. 1 Зависимость отношения р / р ( от числа отсчетов временного ряда
4. ВЫВОДЫ
В работе синтезирован эффективный РНК-алгоритм спектрального оценивания, позволяющий в отличии от известного алгоритма повысить его долговременную численную устойчивость. Это достигается тем, что в синтезированном алгоритме использовалась устойчивая к вычислительным ошибкам процедура триангуляризации корреляционной матрицы (метод Хаусхолдера). Предложен рекуррентный вариант указанной процедуры применительно к РНК-алгоритму спектрального оценивания. Показано, что использование синтезированного алгоритма позволило расширить диапазон числа обусловленности корреляционной матрицы входной последовательности, при котором сохраняется долговременная численная устойчивость алгоритма до 60 db при числе отсчетов входной последовательности равном 1 0000, что соизмеримо с пределами максимальной длины временных рядов различной природы. Таким образом, применение синтезированного алгоритма в задачах спектрального оценивания позволит избежать эффекта искажения СПМ, возникающего в реальных вычислительных устройствах за счет конечной длины машинных слов.
Литература
1. Marple S.L. Jr. "Digital spectral analysis with application"(Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1987)
2. Haykin S. "Adaptive Filter Theory"(Prentice Hall, 1996, 989p)
3. Wilkinson J..H. "The algebraic eigenvalue problem: monographs in Numerical Analysis"(Clarendon Press, Oxford, 1988)
4. Golub G.H., Van Loan C.F. Matrix Computations. - Johns Hopkins University Press, 1993.
5. Householder A.S. "A class of method for inverting matrices". J. Sos. Indust. Appl. Math, 1958, Vol.6, # 2, pp. 189-195.
