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

2006_ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА_Сер. 1_Вып. 1

МАТЕМАТИКА

УДК 517.977

А. Е. Барабанов, В. И. Ситников

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

Введение

Классическая задача расчета оптимальных в среднеквадратичном смысле фильтров для стационарных процессов имеет известное решение методом пространства состояний в форме Шеннона или спектральным методом Винера—Колмогорова [1]. При расчете параметров фильтров в случае больших размерностей возникают значительные вычислительные погрешности, а также проблемы с вычислительными ресурсами — памятью и быстродействием. Так, в ряде задач обработки звукового сигнала, в том числе для подавления акустического эха, количество существенных параметров спектральных плотностей становится равным нескольким тысячам.

В спектральном методе наибольшую сложность имеет процедура факторизации многочлена, положительно определенного на единичной окружности. В данной работе предлагается быстрый метод факторизации с одновременным улучшением точности вычислений. Известно, что операция факторизации тесно связана с решением соответствующих уравнений Риккати [2] и [3]. Наиболее полной работой, охватывающей все, что было сделано в области теории Риккати до 1991 года, считается [4]. В [5] можно найти описание теории дискретного уравнения Риккати.

Вычислительные проблемы факторизации

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

с(г) = спх-п + ап-1г-п+1 + ... + вхх-1 + со + с\х + ... + спхп.

Далее будем предполагать, что значения с(г) на единичной окружности вещественны и положительны. Требуется найти многочлен р(г) = ро + Р1^ + ...+ рпгп с корнями вне единичного круга, удовлетворяющий уравнению с(г) = р(г)р(г-1).

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

© А. Е. Барабанов, В. И. Ситников, 2006

Пример. В качестве примера возьмем квазимногочлен порядка п = 70, вычислим его корни, а затем обратно вычислим коэффициенты квазимногочлена по корням. На рис. 1 сплошной линией обозначены значения коэффициентов исходного многочлена, а пунктирной линией — результат вычислений, проведенных при помощи пакета программ Ма!;ЬаЪ 6.0.

0 10 20 30 40 50 60 70

Рис. 1. Коэффициенты исходного и восстановленного по корням многочленов.

Нетрудно видеть, что расхождение результатов уже для п = 70 существенно. В задачах обработки акустического эха, где порядки достигают нескольких тысяч, результаты вычислений по этому способу становятся непригодными.

Сведение к уравнению Риккати

Введем обозначения для коэффициентов факторизумого квазимногочлена с(г) со свободным членом с0, для коэффициентов искомого квазимногочлена р(г) со свободным членом ро и для вспомогательных матриц е и 7:

Р =

(Р\

\Рп

/сЛ

\сп)

е\ =

0

7 =

/0 1 о ...

V0 •••

0

. . 1

0 0/

Будем считать, что сп = 0, так как в противном случае можно понизить степень. Покажем связь между факторизацией и уравнением Риккати.

Теорема 1. Между множеством .многочленов р(г), удовлетворяющих уравнению факторизации р(г)р(г—1) = с(г) и условию р(0) > 0, и множеством решений Н уравнения Риккати

со

— eT He

существует взаимно-однозначное соответствие. По матрице Н коэффициенты многочлена р(г) определяются формулами

ро = \/ со — еТ Не, р = —(с — 1Не).

Ро

По коэффициентам многочлена р(г) матрица Н определяется равенством

п — 1

Н = £ 7кррТ(7к)Т.

к=0

При этом р(г)/ро есть характеристический многочлен матрицы

А = (с-ТЯе)^

со — еТ Не '

то есть р(г) = ро д.вЬ(1 — Лг).

Доказательство. Пусть Ь(г) = (1,г,..., гп—1), так что многочлены выражаются через Ь(г) и векторы коэффициентов:

р(г) = ро + гЬ(г)р, с(г)=со + гЬ(г)с + сТ г—1ЬТ (г—1).

Из определений матрицы 7 и вектора е следует, что Ь(г)(1 — 7г) = еТ, и поэтому

Ь(г) = еТ (I — 7г) — 1 = еТ + Ь(г)7г.

Для любой матрицы Н

Ь(г)(Н — 7Н7Т )ЬТ (г—1) = (еТ + гЬ(г)7)Н (е + 7Т ЬТ (г—1)г—1) —

— Ь(г)7Н7Т ЬТ (г—1) = = еТ Не + гЬ(г)7Не + еТ 7Т ЬТ (г—1)г—1.

Пусть задан факторизующий многочлен р(г). Определим матрицу Н по формуле в формулировке теоремы и докажем, что она удовлетворяет уравнению Риккати. Из определения следует, что

ТТ = рр

H - JHJT - —T

так как Jn = 0. Подставим это равенство в уравнение факторизации:

c(z) = p(z)p(z-1) = (ро + zL(z)p)(po + pT LT (z-1)z-1) =

= p0 + zL(z)ppo + popT LT (z-1)z-1 + L(z)ppT LT (z-1) = = pi + zL(z)ppo + popTLT(z-1)z-1 + L(z)(H - JHJT)LT(z-1) = = p2o + eT He + zL(z)(ppo + JHe) + (ppo + JHe)T LT (z-1)z-1.

Приравнивая коэффициенты при одинаковых степенях, находим, что

со = p0 + eT He, с = ppo + JHe.

Если po = 0, то порядок квазимногочлена c(z) = p(z)p(z-1) не больше п—1, что противоречит предположению cn = 0. Поэтому po = 0, вектор p явно выражается через H и

Н = JHJT + ррт = JHJT + (c-JHe)(c-JHe)^

со — eT He

что совпадает с уравнением Риккати.

Пусть H — некоторое решение уравнения Риккати. Определим вектор r и многочлен r(z) равенствами

ттт , N zL(z)r

f = с — JHe, r(z) = l +-Чгтг-

со — eT He

Вектор С удовлетворяет уравнению

T

гп rrT

Н - JHJT =

со — eT He

Докажем, что функция R(z)=r(z)r(z-1 )(co—eTHe) равна c(z):

T

К{г)=с0-егНе + гЬ{г)г + гтЬ1\г-1)г-1 + Ь{г)-'—— Ьт {г-1) =

со — ет Не

= с0 — ет Не + гЬ(г)г + тт Ьт (г-1)г-1 + Ь(г)(Н — JHJT )Ьт (г-1) = = с0 — ет Не + гЬ(г)(т + JHe) + (г + JHe)T Ьт (г-1)г-1 + ет Не = = с0 + гЬ(г)с + ст Ьт (г-1)г-1 = с(г).

По условию, с(г) > 0 при \г\ = 1. Однако, если \г\ = 1, то с(г) = \т(г)\2(со—етНе). Следовательно, со — етНе > 0 и искомый факторизующий многочлен определяется равенством р(г) = (со — е т Не)12 т(г).

Наконец, матрица А из формулировки теоремы имеет стандартную форму Фробе-ниуса А = J — р-1рет, и поэтому коэффициенты ее характеристического многочлена стоят в первом столбце. □

Следствие 1. Пусть Н —произвольное решение уравнения Риккати

со

eTHe

Тогда eTHe < со и H > 0.

В [7] уравнение Риккати (1) решается с использованием сигнум-функции, что требует значительных вычислительных затрат. Ниже предлагается более экономный метод простой итерации.

Монотонность оператора Риккати

Пусть A, C, D —заданные матрицы. Введем обозначение для оператора Риккати Ric(H) = AHA* + (D - AHC)(R - C*HC)-1(D - AHC)*,

заданного на множестве <1от(Ше) всех самосопряженных матриц Н, для которых матрица Я — С * НС невырождена. Этому оператору соответствует коэффициент усиления

К (Н) = (Б — АНС )(Я — С* НС)-1.

Теорема 2. Пусть Н1,Н2 € ¿от(Ше). Определим АК = К(Нг) — К(Н2), Аг = А — КС*, А2 = А — К2С*. Тогда

Ше(Н2) — Вле(Н\) = А1(Н2 — Н1)А*1 + АК (Я — С *Н2С )АК * = = А2(Н2 — Н\)А2 — АК (Я — С *Н1С )АК *.

Доказательство. Введем обозначения Кг = К(Нг), К2 = К(Н2). Из определения К(■) следует, что

К (Н )(Я — С *НС)К (Н )* = К (Н )(Б — АНС )* = (Б — АНС )К (Н )*

для любой матрицы Н € <1от(Я1е). Поэтому

—Кг(Я — С*НС)К* = Кг(Я — С*НС)К* — 2Яе [(Б — АН С)К*] = = К (Я — С*Н2С)К* — КС* (Н2 — Нг )СК* — — 2Яе [(Б — АН2С)К*] — 2ЯеА(Н2 — Нг)СК*.

Подставим правую часть в выражение для АШе=Яле(Н2)—Яле(Н1):

АШе = А(Н2 — Нг )А* + К2(Я — С * Н2С )К* — Кг (Я — С *НгС )К * = = А(Н2 — Нг )А* — 2ЯеА(Н2 — Нг )СК* + Кг (Я — С*Н2С)К* + + К2(Я — С* Н2С)К** — 2ЯеК2(Я — С* Н2С)К* + Кг (Я — С* Н2С)К* = = (А — КгС*)(Н2 — Нг )(А 1 — Кг С* )* + АК (Я — С*Н2 С )АК *,

что доказывает первое равенство. Второе доказывается аналогично. □

Следствие 2. Если Н2>Нг и Я- С* Н2С > 0; то Шс(Н2) > Д»с(#1). Таким образом, оператор Яле возрастает на множестве Н, для которых С * НС < Я. Для нахождения решения уравнения Риккати Н = Ше(Н) применим метод простой итерации:

Нк+1 = Ше(Нк), к > 0,

с некоторой начальной матрицей Но.

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

Теорема 3. Пусть уравнение Риккати Н = Я1е(Н) имеет решение Н и Я — С * Н С > 0. Пусть матрица А устойчива. Тогда метод простой итерации с начальным данным Но = 0 монотонно сходится к минимальному неотрицательному решению уравнения Риккати.

Доказательство. Докажем по индукции, что Нк < Н для любого к. При к = 0 это следует из неравенства

Н = ^2 Ак(Б — АНС)(Я — С*НС)-1 (Б — АННС)(Ак)* > 0 = Н0.

к=0

По следствию 2 из неравенства Нк < Н следует Нк+1 = Шс(Нк) < Шс(Н) = Н, что доказывает индукционный шаг.

Таким образом, Нк < И и поэтому Я — С*НкС > Я — С*НС > 0 для любого к. Докажем монотонность (Нк) по индукции. При к = 0 это очевидно, так как 0 = Но < ПЯ-1П* = Н\. Далее, по следствию 2 из условия Нк < Нк+\ следует, что Нк+\ = Шс(Нк) < Ягс(Нк+1) = Нк+2, что доказывает индукционный шаг.

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

Пусть Нт1п — минимальное неотрицательное решение уравнения Риккати. Тогда Нт1п < Н и, поэтому, Я — С*НттС > 0. Заменим в проведенном доказательстве теоремы Н на Нт1п. Тогда получим, что Нж < Нт1п. Но Нж > 0, поэтому Нж = Нт1п. □

Стабилизирующее решение

Решение Н уравнения Риккати называется стабилизирующим, если матрица А — К(Н)С устойчива. Из теоремы 3 не следует, что предельная матрица Нж является стабилизирующим решением уравнения Риккати. В задаче факторизации многочленов матрица А = . устойчива, так все ее собственные числа равны нулю. Кроме того, Я есть число. В следующем утверждении будем предполагать, что размерность Я равна 1.

Теорема 4. Пусть у уравнения Риккати Н = Ягс(Н) существуют такие решения Нт\п, Нтах, что .матрица Ат[п = А — К(Нт[п)С* устойчива, а .матрица Атах = А — К(Нтах)С* — антиустойчива. Пусть, кроме того, Я — С*Нт\пС > 0.

Тогда для любого решения Н этого уравнения Риккати Нт[п < Н < Нтах и Я — С*НС > 0.

Доказательство. Пусть Н — некоторое решение уравнения Риккати Н = Ше(Н). Введем обозначения ДНт;п = Н — Нт;п, ДКт;п = К(Н) — К(Нт;п). Тогда по теореме 2

ДНтш = Ат^НтАп + ДКт1п(Я — С* НС)ДК**п,п.

Из устойчивости матрицы Ат^ следует сходимость ряда в следующем решении последнего уравнения:

ДНтШ = £ Атш[ДКтШ(Я—с* НС)ДКтш](Атш )к.

к=0

Допустим, что Я — С*НС < 0. Тогда ДНт;п = Н — Нт;п < 0. Следовательно, Я — С*НС > Я — С*Нт\пС > 0, что противоречит предположению. Следовательно, Я — С*НС > 0 и, поэтому, ДНт;п = Н — Нт;п > 0.

Таким образом, доказано, что для любого решения Н уравнения Риккати Я — С * НС > 0 и Н > Нт1п. В частности, это верно и для решения Н = Нтах.

Докажем, что Н < Нтах для произвольного решения Н уравнения Н = Я1с(Н). Введем обозначения ДНтах = Нтах — Н, ДКтах = К(Нтах) — К(Н). В силу второго равенства из теоремы 2

ДНтах = АтахДНтахАтах — ДКтах(Я — С* НС)ДК*т&х-

Все собственные числа матрицы Amax лежат вне единичного круга, поэтому A^^ ^ 0 экспоненциально при к ^ Отсюда, последнее уравнение имеет решение

AHmax = £ Am1x[AKmax(R - C* HC)AK*max](A*max)-k. k=l

Поскольку R — С*HC > 0, получаем AHmax = Hmax — Н > 0. □

Следствие 3. Пусть выполнены условия теоремы 4- Тогда в итерационной процедуре Hk+i = Ric(Hk) при Ho = 0 предельная матрица Нж = limHk является минимальным и стабилизирующим решением уравнения Риккати-

Экономный алгоритм факторизации

В соответствии с теоремой 1 уравнение Риккати для задачи факторизации многочлена удовлетворяет условиям теорем 3 и 4. По этим теоремам метод простой итерации

Ж°) = О, .>0,

c0 — e*H (k> e

монотонно сходится к решению H уравнения Риккати, а многочлен с коэффициентами

Ро = \/ со — е*Не, р = — (с — JHe)

Ро

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

p(z )p(z-1) = c(z).

Структура матриц этого уравнения позволяет упростить вычислительный алгоритм. Введем обозначения (H(kдля элементов матрицы Hи (p)"=1 для эле-менетов вектора p(k):

Рок) = Veo ~ е*Нке, р{к) = -¡-(с - JF(fc)e), к > 0.

p0

Матрица J является оператором сдвига вектора вверх, поэтому

{п(к ■■■ нП 0\

H (к+1) =

н (к) Hn2 0

2n

Нг!П 0

+

■■ 0 0/

Поскольку вектор e есть первый орт,

(к) Po ) =

co - Hk,

' (Pk)2

(к) (к) Pl Pn-1 (k) (k) P 1 P n

(k) (k) (k) (k)

Pl Pn-l Pl Pn X

•• (к) (k)

. P 2 P n

P(k) P(k) (Лк))2

Pn-1Pn (Pn ) /

p(k) =

(к) Po '

Í ci - Hk \

c H(k) cn 1 Hn1

1

n

Подставляя в это выражение первые столбцы матриц Нпри г < п, получим следующие рекуррентные формулы для расчета коэффициентов факторизующего многочлена:

p(k+1) _

ро _

n-1 „

со - £ (Р+^

i=0

Л ( n-1-j \

(k+1) Л V^ (k-i) (k-i) I 1 о

о = ^ТТ) \ ci- Ъ Pi+1 Pi+1+j ) > з = 1, 2,..., П.

P j = "TlxТТ I с:

Начальные данные нулевые: р^ = 0 при 1 < ] < п. Поэтому при к < п в суммах присутствуют только к слагаемых.

Последовательность коэффициентов рО^ убывающая, так как матрицы Н «, а вместе с ними и элементы = со — (рО^)2, возрастают. Алгоритм останавливает работу, когда изменение коэффициента рО^ оказывается в пределах машинной точности вычислений.

Количество итераций алгоритма зависит от близости нулей функции с(-) к единичной окружности. Проведенные эксперименты показали, что если расстояние между нулями с(-) и единичной окружностью не меньше 0.04, то количество итераций превосходит п не более, чем на п/2. Поэтому вычислительная сложность всего алгоритма оказывается пропорциональной 3п2/4. Относительно небольшое количество операций повышает одновременно точность результата расчета. В приведенном ранее примере погрешность оценок коэффициентов, определенных по итеративному алгоритму, не превосходит 10-3.

Заключение

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

Summary

A. E. Barabanov, V. I. Sitnikov. High degree quasi-polynomial factorization.

Computational problems of polynomial spectral factorization are discussed. Factorization is reduced to a Riccati equation of a special type. It is proved that the Riccati operator is monotonous and the iterative procedure converges to the stabilizing solution. A memory reduced computational algorithm is presented.

Литература

1. Фомин В. Н. Рекуррентное оценивание и адаптивная фильтрация. М.: Наука, 1984. С. 167.

2. Brian D. O., Anderson, Konrad L. Hitz, and N.D. Diem. Recursive algorithm for spectral factorization // IEEE transactions on circuits and systems. Vol. Cas-21. N6. 1974.

3. Якубович В. А. Частотная теорема в теории управления / Сиб. мат. журнал. №2. 1973.

4. Bittanti S., Laub A. J., and Willems J. C. The Riccati Equation. Springer Verlag, Berlin. 1991.

5. Malinen J. Discrete Time ffi-Algebraic Riccati Equations. Doctoral dissertation. Technical Report A428, Helsinki University of Technology. Espoo, Finland. 2000

6. Jan C. Willems. Last Squares Stationary Optimal Control and the Algebraic Riccati Equation // IEEE Transactions on Automatic Control. Vol. AC-16. N6. 1971.

7. Алиев Ф.А., Бордюг Б. А., Ларин В. Б. Дискретное обобщенное уравнение Риккати и факторизация матричных полиномов // Автоматика. 1990. №4. С39-46.

Статья поступила в редакцию 6 сентября 2005 г.