Об оценке спектра линейного оператора Текст научной статьи по специальности «Математика»

УДК 519.71

Вестник СПбГУ. Сер. 1. 2012. Вып. 1

ОБ ОЦЕНКЕ СПЕКТРА ЛИНЕЙНОГО ОПЕРАТОРА*

С. М. Ермаков1, К. О. Видяева2

1. С.-Петербургский государственный университет, д-р физ.-мат. наук, профессор, ermakov@math.spbu.ru

2. С.-Петербургский государственный университет, аспирант, karina.vidyaeva@gmail.com

В работе описаны некоторые новые подходы к анализу спектра линейных операторов. Эти подходы основаны на обобщении известного метода Крылова [2] вычисления коэффициентов минимального многочлена матрицы. Обобщение, полученное в работе, позволяет иметь дело с матрицами бесконечного порядка, что важно, например, в теории массового обслуживания и ряда других приложений. Полученные результаты имеют также связь с БЭЛ-разложениями и методами рекуррентного сдвига («Гусеница» [3, 4]), используемыми при анализе временных рядов.

Дальнейшее изложение будет опираться на лемму 1, имеющую достаточно общий характер.

1°. Гассмотрим последовательность вещественных чисел

т

& = * = 1,2,..., (1)

к= 1

где ад; и Ьгк —числа (м.б. комплексные), Т может быть бесконечным, и в этом случае предполагается, что ряды (1) абсолютно сходятся и векторы а = и Ь =

принадлежат I2. Все Ък различны, но могут быть равными по модулю:

М > ы > ••• > \Ът\- (2)

Имеет место следующая

Лемма 1. Пусть I = 1,2, ■ ■ ■, имеют вид (1) и для некоторого Ь, 1 ^ Ь ^ Т, константы I = 0,1,..., Ь, и возрастающей последовательности номеров Ь, Ь = 1,2,..., существует решение системы уравнений

ь

= 3=0,1,...,Ь-1. (3)

1=0

Тогда если

(1) к = 1,... ,Т, удовлетворяет неравенствам (2), но \Ьь\ > |Ьь+1|, \Ьь\ > О,

(2) ак^0, к=1,...,Ь,

(3) для некоторого I = 1о, 0 ^ 1о ^ Ь, с^ = 1, то существует предел

Итсг(4)=сг, / = 0,1,...,Ь, (4)

* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00769-а). © С. М. Ермаков, К. О. Видяева, 2012

и многочлен

ь

Q(x) = J2cixl

1=0

имеет своими корнями числа Ъ]~, k = 1,..., L,

Q(x) = с (х — bk), с = const. (5)

к= 1

L

Доказательство. Обозначим <34(ж) = с\ х'' и будем считать 1о = Ь (с^ = 1),

1=0

что не умаляет общность. Подставляя в (3) выражение (4) для ^+1+], получаем

Е ^ Е = Е Е = о. «».....

1=0 к= 1 к=1 г=о

Или

L TL

fc=i fc=L+i г=о

Для каждого j разделим j-e равенство на Имеем

ht+j l т ht+j

„.Л— л —1,2,...,L. (6)

й=1 1 г=о й=ь+1

Заметим, что С^^Ък) ограничены |<3((Ьй)| < С1 < сю, /г = 1,2,..., в силу предположения о существовании решения системы (3) и условия < \Ъь\, к = Ь + 1, Ь + 2,..., т

и Ы < с2 < со (абсолютная сходимость ряда при Т —>• оо). Поэтому

й=Ь+1

fc=L+l

< с\ ■ со • Qt+3, а = тах

4 ' 4 Ь+1<к<Т

Таким образом, совокупность равенств (6), рассматриваемая как система линейных уравнений относительно a,i~Qt(bk), к = 1,..., L, имеет в качестве правой части вектор, стремящийся к нулевому вектору с ростом t, и определитель, отличный от нуля при каждом t (определитель Вандермонда). В силу сказанного a]~Qt{bk) —>• 0,

оо

к = 1,..., L, и, в соответствии с (2), lim Qt(bk) = 0, что и доказывает лемму. □

оо

Замечание 1. При L = Т второе слагаемое в равенствах (6), j = 0,..., L — 1, отсутствует. Полученная система является однородной, и Q(x) может быть вычислено при t = 1 (метод Крылова). Предположим теперь, что в гильбертовом пространстве Н определён оператор 1С, и для X, Y из Н мы умеем вычислять скалярные произведения £t = (1СгХ,У). Если X и Y имеют ненулевые проекции на все инвариантные подпространства Hfc оператора 1С (предполагается, что Н = |J Щ, Щ р| Н; = 0, к, I = 1, 2,...,

к

спектр 1С дискретен), то £t могут быть представлены в виде

т

6 = í = 1,2,...,

к=1

где Ад — собственные числа /С, и результаты леммы очевидным образом применимы для приближенного вычисления X¡..

Алгоритм 1. Вычисление Ai,..., A¿. Пусть L < Т и вычислены с достаточной степенью точности числа £tj£t+b • • • • • • Для достаточно больших значений t.

Считаем t настолько большим, что ¿р близки к своим предельным значениям Cj. Составляем систему из L уравнений относительно cf.

co€t + ci£t+i Н-----Ь cL£t+L = О,

co6+i + H-----cL^t+L+1 = 0,

:

SoCí+l-l + Clít+L+Í H-----h Cl£í+2L-1 = 0.

Поскольку Cj должны быть известны с точностью до общего множителя, фиксируем Ъь = 1.

Решив полученную систему, получим приближенное значение чисел Cj. Затем следует решить уравнение XL + iA¿_1 + ... + со = 0. При L = Т получаем точное значение Cj. Нужно отметить, что при решении конкретной задачи мы можем заранее не знать, каким следует выбрать L и правильно ли выбраны начальные векторы X и Y (a,j ф^ 0). Окончательные выводы приходится делать на основании численного эксперимента.

Замечание 2. Если К, = А является матричным оператором, то вопрос о выборе векторов X и У, порождающих представление векторов (1),

& = (А*Х,¥), A=\\aij\\lj=v

достаточно полно изучен в линейной алгебре [2]. Линейная оболочка векторов АгХ образует подпространства Крылова, a Q(x) при L = Т является минимальным многочленом матрицы А. Практически компоненты X и Y следует получать с помощью датчика случайных чисел.

Алгоритм 2. Пусть ji,... ,js, s < N, —различные числа из множества 1,... ,Т.

S

Предположим, что нам известны А..., Ajg и R(x) = П (ж — Aj'k)- Тогда выбором

к=1

начального вектора R(X) мы можем обратить в нуль коэффициенты разложения (1), стоящие перед найденными собственными числами (ортогонализация). Такой подход позволяет в некоторых случаях вычислять оставшиеся собственные числа с большей точностью.

Действительно, пусть £t = (Аг(К(А)Х),У), где R(X)—многочлен, корнями которого являются собственные числа А^,..., Ajs, R(x) = xs + а\х+ • • • + as. Тогда применение матричного оператора R(A) к X уничтожает всю информацию о собственных числах с номерами ji,... ,js. Алгоритм состоит в применении алгоритма 1 к модифицированному начальному вектору X.

2°. Метод Монте-Карло. Многие прикладные задачи теории массового обслуживания, расчета ядерных реакторов и др. включают вычисление итераций линейных операторов с помощью метода Монте-Карло. Особый интерес представляет информация о первых собственных числах, а метод Монте-Карло позволяет оценивать

(КгХ, У) при больших Ь со сравнительно малой трудоемкостью. Мы далее опишем алгоритм метода Монте-Карло для матричного оператора А. Обычно алгоритм вычисления связан с моделированием цепи Маркова с п состояниями, удовлетворяющей условиям согласования. Условия согласования для цепи, определяемой стохастической матрицей перехода V = Нл,¿11^=1, таковы:

Рг^ > 0, если а^ ф 0, 1,3 = 1 ,...,п. (8)

Алгоритм основан на методах моделирования дискретных распределений.

Считаем заданными векторы X = {х\,..., хп)', У = (у1,...,уп) и матрицу ^ = и предполагаем, что имеет место представление (1) (штрих обозначает

здесь и далее транспонирование). Опишем один из алгоритмов метода Монте-Карло для оценивания = (АгХ,У) = (Х'(А'У,У). Зададим стохастическую матрицу Р =

п

Ц^Ц, согласованную с матрицей А'. Обычно полагают^ = \а'^\/в® = \а'ц\-

¿=1

Заметим теперь, что для любого вектора Z = ..., гп) и вектора р° = (р®, • • • ,Рп)> с ним согласованного, случайный вектор

„/

-1

PÏPij

j=i

является несмещенной оценкой вектора ZA', где i выбирается случайным образом в соответствии с распределением р°, a j—в соответствии с распределением

(рп,... ,pin).

В общем случае вектор с компонентами

zioai0i1 ■ ■ ■ а'н-гн .

Vh = —-> 4 = 1

Pi0Pioil ■ ■ -Pit-lit

где ¿о —^ «1 —^ • • • —'Н отрезок траектории цепи Маркова с начальным распределением р° и переходной матрицей Р, является несмещенной оценкой для Z(A')t, а случайная величина £t = ЩгУн оценивает (Z(A')t,Y) = (АгХ,У), если положить Z = X'. Трудоемкость вычисления £t оценивается как 0(_/Vlog2n) [1], где N — число независимо вычисленных оценок для £t. Это делает описанный выше алгоритм эффективным при больших (ив.

При использовании метода Монте-Карло вместо решения системы (8) лучше пользоваться нормальными уравнениями метода наименьших квадратов, который для некоторого, выбранного заранее, натурального то может выглядеть следующим образом:

L— 1 fc+m 2

Е Е ® + • • • + 4%+l+ï) =min, (9)

j=0 t=k

qt > 0, то > L. Для простоты далее полагаем qt = 1.

Если N —> оо и случайной ошибкой можно пренебречь, то имеет место результат, вполне аналогичный лемме 1. А именно

Лемма 2. В условиях леммы 1 система нормальных уравнений

I к^т

= р = 0,1,... ,L — 1, (10)

j=о t=k

имеет при к —> оо своим решением определенные с точностью до постоянного множителя коэффициенты многочлена

ь

<2 (ж) = с (х - АД

з=1

где с — константа.

Доказательство. Действительно,

т

¿=1

(а^ + • • • + атХ'т+р)

г=к

(I +

¿=0

+ «2Л2 ( Е СГЛ2 ]+••• + О-тХ'т ( Е С? и=0 / Ь'=°

0, р = 0,...,Ь-1. (И)

Введя обозначение = ^ с- А^, преобразуем нашу систему в систему из Ь

уравнении вида

к-\-т к-\-т

а^(Х1) ]Г А'[ (а1Х\+р + .. . + аТХ'т+Р) + .. , + атд(Ат) ]Г А^ (а1Х\+р +.. . + атХ'т+Р) =0; г=к г=к

при Ь = Т имеем =0, г = 1,..., Ь; если Ь < Т, то

к^т

агЯ(Х1) Е Х{ (а1Х\+р + • • • + аТХ*+р) + • • • +

г=к

к-\-т

+аьЯ(Хь) Е АI (а1Х[+р + • • • + аТХ*+р) + • • • +

г=к

к-\-т

+атд(Ат) Е (а1Х[+р + • • • + аТХ*+р) =0, р = 0,...,Ь-1.

г=к

Достаточно до множить обе части р- го равенства на 1 / \Хк^гт+р \ и повторить рассуждения, проведенные при доказательстве леммы 1. □

Вопрос о поведении ошибки коэффициентов с^ при наличии случайной погрешности у весьма сложен. Здесь мы имеем проблему, аналогичную проблемам «ошибка в переменных» в регрессионных задачах. Как известно [6], эти проблемы далеки от полного исследования. В нашем случае ситуация облегчается, так как мы можем распоряжаться выбором числа испытаний N, и вопрос о поведении погрешности может быть решен экспериментально.

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

Первая группа результатов относится к случаю, когда случайная ошибка отсутствует и применима лемма 1.

Рассматривалась матрица восьмого порядка с вещественными собственными числами. Вычисления производились при различных значениях параметров Ь и то, где Ь — длина окна, а то — степень матрицы, то есть = (АтХ, У). Приведена таблица результатов (см. табл. 1) с наилучшими значениями параметров. В левой колонке записаны точные собственные числа исследуемой матрицы, в последующих колонках — результаты применения алгоритма 1 с параметрами Ь и то.

Таблица 1

Ь; т числа 8; 1 7; 1 6; 2 5; 5 4; 8 3; 14 2; 25 1; 14

30.1947 30.1947 30.1947 30.1947 30.1947 30.1947 30.1947 30.1947 30.1947

13.2645 13.2645 13.2645 13.2645 13.2645 13.2645 13.2645 13.2641 -

-8.2408 -8.2408 -8.2408 -8.241 -8.2405 -8.2414 -8.2455 - -

6.0134 6.0134 6.0498 5.9697 5.9532 5.8677 - - -

-3.7765 -3.7765 -3.6329 -3.5384 -4.0049 - - - -

3.7439 3.7439 4.3031 3.1941 - - - - -

-2.5051 -2.5051 -1.4452 - - - - - -

1.3058 1.3058 - - - - - - -

Как мы видим, если длина окна Ь равна размерности матрицы, то все собственные числа матрицы находятся без погрешности уже при то = 1 (см. замечание 1).

В таблице 2 приведены результаты вычисления собственных чисел, найденных по алгоритму 2, то есть с учетом полученного значения первого собственного числа.

Таблица 2

Ь; т числа 8; 1 7; 1 6; 2 5; 5 4; 8 3; 14 2; 25

30.1947 30.1947 30.1947 30.1947 30.1947 30.1947 30.1947 30.1947

13.2645 13.2645 13.2645 13.2645 13.2645 13.2645 13.2645 13.2642

-8.2408 -8.2408 -8.2408 -8.2408 -8.2407 -8.2410 -8.2438 -

6.0134 6.0134 6.0612 5.9981 5.9872 6.1038 - -

-3.7765 -3.7765 -3.6432 -3.66 -4.0111 - - -

3.7439 3.7439 4.3031 3.2293 - - - -

-2.5051 -2.5051 -1.5317 - - - - -

1.3058 1.3058 - - - - - -

Как показывают исследование, в большинстве случаев учет значения Ах улучшает точность вычисления последующих собственных чисел.

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

Собственные числа матриц с близкими собственными числами также находились достаточно хорошо, точность вычислений не изменялась.

Также приводится пример, когда значения вычислялись методом Монте-Карло. Определение второго собственного числа требовало достаточно большого числа испытаний порядка N = 10000 и более. Для удобства вычислений матрица нормировалась, все собственный числа были уменьшены в 13 раз. Точными значениями были Ах = 2.322; А2 = 1.0203. Полученные в процессе моделирования случайные оценки не являются несмещенными. Величина смещения имеет порядок 0{1/Ы). При

L; т числа 9; 2 8; 4 7; 6

1.0000 1.0000 1.0000 1.0000

—0.1466+0.0862Í —0.1465+0.0862Í —0.1465+0.0862Í —0.1465+0.0862Í

—0.1466—0.0862Í —0.1465—0.0862Í —0.1465—0.0862Í —0.1465—0.0862Í

-0.0420+0.1387Í -0.0418+0.1387Í —0.0418+0.1386Í —0.0418+0.1387Í

—0.0420—0.1387Í -0.0418-0.1387Í —0.0418—0.1386Í —0.0418—0.1387Í

0.1397 0.1398 0.1398 0.1397

—0.0561+0.0595Í —0.0495+0.0343Í —0.0633+0.0399Í -0.0426

—0.0561—0.0595Í —0.0495—0.0343Í —0.0633—0.0399Í -

0.0486+0.0236Í 0.0208 - -

0.0486—0.0236Í - - -

6,8 5,10 4,13 3,16

1.0000 1.0000 1.0000 1.0000

—0.1465+0.0862Í —0.1418+0.0882Í —0.1459+0.0730Í —0.1556+0.0898Í

—0.1465—0.0862Í —0.1418—0.0882Í —0.1459—0.0730Í -0.1556—0.0898i

—0.0418+0.1387Í 0.0059+0.1533Í -0.0047 -

—0.0418—0.1387Í 0.0059—0.1533Í - -

0.1398 - - -

- - - -

больших N эмпирическая дисперсия может служить показателем погрешности. Результаты вычислений первых двух собственных чисел и величины дисперсии второго собственного числа приведены в табл. 4.

Таблица 4

L = 9 т = 1 2.3225 1.0274 2.3329 0.9627 2.3315 1.1249 2.3221 0.9188 2.2893 0.9554 ?2 0.0053

L = 8 т = 1 2.3229 0.9740 1.8076 1.1065 2.3240 1.1574 2.3230 1.0392 2.3232 1.1368 ?2 0.0046

L = 7 т= 2 2.3234 1.5309 2.3182 1.6366 2.3219 1.1866 2.3196 1.5032 2.3218 0.9735 ?2 0.0611

L = 6 т = 5 2.3217 1.1715 2.3221 1.4639 2.3259 0.8086 2.3198 1.8559 2.3264 1.7143 ?2 0.1425

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

Литература

1. Ермаков С. М. Метод Монте-Карло в вычислительной математике. Вводный курс. СПб.: Невский Диалект; М.: Бином, 2009. 192 с.

2. Гантмахер Ф. Р. Теория матриц. Изд. 4-е, дополненное. Москва: Наука, 1988. 548 с.

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

4. Ермаков С. M., Котова Л. Ю. О выборе базисных функций в регрессионном анализе // Сб. работ кафедры статистического моделирования СПбГУ. 1999. С. 3-43.

5. Capon J. High resolution frequency-wavenumber spectrum analysis. Proc IEEE. Vol. 57, N8. 1969. P. 1408-1418.

6. Fuller W. A. Measurement error models. New York: Willy, 1987. 440 p.

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