CC BY
38
УДК 519.71
Вестник СПбГУ. Сер. 1. 2012. Вып. 1
ОБ ОЦЕНКЕ СПЕКТРА ЛИНЕЙНОГО ОПЕРАТОРА*
С. М. Ермаков1, К. О. Видяева2
1. С.-Петербургский государственный университет, д-р физ.-мат. наук, профессор, [email protected]
2. С.-Петербургский государственный университет, аспирант, [email protected]
В работе описаны некоторые новые подходы к анализу спектра линейных операторов. Эти подходы основаны на обобщении известного метода Крылова [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 г.
