# Собственные вектора и значения трехдиагональной матрицы

## Вычисление характеристического полинома

После того, как исходная действительная симметричная матрица была сведена к трехдиагональной форме, одним из возможных способов отыскать ее собственные значения становится прямое вычисление корней характеристического полинома pn(l). Характеристический полином трехдиагональной матрицы для любого значения l может быть вычислен с помощью эффективной процедуры рекурсии. Полиномы меньших порядков, получающиеся в процессе этой рекурсии, образуют последовательность Штурма, которая может послужить для локализации собственных значений на интервалах действительной оси. После локализации могут применяться методы поиска корней, например, метод Ньютона или дихотомия. Соответствующие собственные вектора могут быть найдены обратной итерацией.

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

## Алгоритмы QR и QL

Основная идея алгоритма QR заключается в том, что любая действительная матрица может быть представлена в форме
A = QR,
где Q ортогональная, а R -- верхняя треугольная матрицы. Для матрицы общего вида это разложение строится применением алгоритма преобразования Хаусхолдера с целюь уничтожения элементов столбцов матрицы A ниже диагонали.

Теперь рассмотрим матрицу, сформированную произведением этих факторов в обратном порядке:

A' = RQ.
Поскольку Q ортогональна, получаем R = QTA. Таким образом,
A' = QTAQ.
Видно, что A' является ортогональной трансформацией A.

Можно проверить, что QR-преобразование сохраняет следующие свойства матрицы: симметрию, трехдиагональность и форму Гессенберга (определенную ниже).

Не обязательно выбирать один из факторов A в виде верхней треугольной матрицы: можно точно также выбрать нижнюю треугольную. Это приведет к QL-алгоритму, поскольку разложение A будет записано в форме

A = QL,
где L -- нижняя треугольная матрица.

Вспомним, что редукцию Хаусхолдера мы проводили, начиная с n-ого (последнего) столбца исходной матрицы. Для уменьшения ошибки округления мы советовали расположить ее наибольшие элементы ближе к правому нижнему углу. Если мы после этого хотим диагонализировать получившуюся трехдиагональную матрицу, то алгоритм QL приведет к меньшей ошибке округления, чем QR, таким образом далее мы будем говорить об алгоритме QL.

Алгоритм QL состоит из цепочки ортогональных преобразований:

As = QsLs, As+1 = LsQs (= QsTAsQs).
Следующая (неочевидная!) теорема является основой для алгоритма в применении к матрице A общего вида: (i) Если A имеет собственные значения, различные по модулю |li|, то As стремится к нижней треугольной форме при неограниченном возрастании s. Собственные значения на диагонали при этом появляются в порядке возрастания модуля. (ii) Если A имеет собственное значение |lj| кратности p, то Ak стремится к нижней треугольной форме, за исключением блока размерности p, примыкающего к диагонали, собственные значения этого блока стремятся к lj. Доказательство этой теоремы достаточно длинное, см. [8].

Объем арифметических действий алгоритма QL имеет порядок O(n3) на итерацию для матрицы общего вида, что является недопустимым. Однако для трехдиагональной матрицы этот объем составляет O(n) на итерацию, а для матрицы в форме Гессенберга O(n2), что делает подобные формы матриц высокоэффективными.

В этой секции мы рассматриваем случай действительной, симметричной, трехдиагональной матрицы. Таким образом, все собственные значения li действительны. Согласно теореме, если имеется lj кратности p, то должно быть по крайней мере (p-1) нулей сверху и снизу от диагонали. Таким образом, матрицу можно будет разбить на подматрицы, которые можно будет диагонализировать раздельно, а усложнение диагональных блоков, имеющее место для матриц общего вида, здесь не важно.

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

aij(s) ~ (li / lj)s.
Хотя li < lj, сходимость может быть медленной при близости этих значений. Сходимость может быть ускорена методом сдвига: если k -- любая константа, то (A - k1) имеет собственные значения (li - k). Если разложить
As - ks1 = QsLs, так что As+1 = LsQs + ks1 = QsTAsQs,
то скорость сходимости будет определяться отношением
(li - ks)/(lj - ks).
Самое главное здесь выбрать на каждой стадии такое ks, чтобы максимально ускорить сходимость. На начальном этапе хорошим выбором является значение ks, близкое к l1, наименьшему собственному значению. В этом случае первая строка внедиагональных элементов будет быстро сведена к нулю. Однако, величина l1 априорно неизвестна. На практике, очень эффективной стратегией (хотя и не доказано, что она оптимальна) является вычисление собственных значений начального минора размером (2 x 2) матрицы A. После этого ks полагается равной тому собственному значению, которое ближе к a11.

Обобщая, представим, что уже найдены (r - 1) собственных значений A. Тогда матрицу можно уменьшить, обнулив первые (r - 1) столбцов и строк:

A = | ( ( ( ( (
(
(
(
(
(
(
0
...
...
0
)
)
)
)
)
)
)
)
)
)
)
| ... |  |  |  |  |  |  |  |  |
|  | 0 |  |  |  |  |  |  |  |
... |  |  | dr | er |  |  | ... |  |  |
... |  |  | er | dr+1 |  |  |  |  |  |
|  |  |  |  | ... |  | 0 |  |  |
|  |  |  |  |  | dn-1 | en-1 |  |  |
0 |  |  | ... |  | 0 | en-1 | dn |  |  |

Значение ks выбирается равным тому собственному значению ведущего минора размером (2 x 2), которое ближе к dr. Можно показать, что скорость сходимости алгоритма при выборе этой стратегии обычно кубическая (и в худшем случае квадратичная при вырожденных собственных значениях). Эта исключительно быстрая сходимость делает алгоритм весьма привлекательным.

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

Как уже говорилось выше, QL-декомпозиция матрицы общего вида производится последовательностью преобразований Хаусхолдера. Однако для трехдиагональной матрицы значительно более эффективно применять плоские ротации Якоби Ppq. Элементы a12, a23, a34, ..., an-1,n обнуляются серией ротаций P12, P23, P34, ..., Pn-1,n. По симметрии, будут обнулены также поддиагональные элементы a21, a32, ..., an,n-1. Таким образом, каждая матрица трансформации Qs будет представлять из себя произведение плоских ротаций:

QsT = P1(s)P2(s)...Pn-1(s),
где Pi обнуляет ai,i+1. Заметим, что в написанном выше выражении участвует именно QsT, а не Qs, поскольку было определено L = QTA.

## QL-алгоритм с неявными сдвигами

Алгоритм, описанный выше, может быть очень эффективным. Однако когда элементы A сильно различаются по порядку величины, вычитание большого ks из диагональных элементов может привести к потере точности для собственных значений с небольшой величиной. Эта трудность обходится применением QL-алгоритма с неявными сдвигами. Неявный QL-алгоритм математически эквивалентен исходному, однако при вычислениях не используется явное вычитание ks1 из A.

Алгоритм основан на следующей лемме. Пусть A -- симметричная невырожденная матрица и B = QTAQ, где Q -- ортогональна, а B -- трехдиагональна с положительными внедиагональными элементами. Тогда Q и B полностью определяются на основании только последней строки QT. Доказательство. Пусть qiT обозначает i-ую строку матрицы QT, а qi -- i-ый столбец матрицы Q. Соотношение BQT = QTA может быть записано как

( ( ( ( (
(
(
(
b1
g1
)
)
)
)
)
)
)
)
*
(
(
(
(
(
(
(
(
q1T
)
)
)
)
)
)
)
)
=
(
(
(
(
(
(
(
(
q1T
)
)
)
)
)
)
)
)
A
a2 | b2 | g2 |  |  |  |  | q2T | q2T |  |  |  |  |  |  |  |  |
|  |  | ... |  |  |  | ... | ... |  |  |  |  |  |  |  |  |
|  |  |  | an-1 | bn-1 | gn-1 | qn-1T | qn-1T |  |  |  |  |  |  |  |  |
|  |  |  |  | an | bn | qnT | qnT |  |  |  |  |  |  |  |  |

n-ая строка этого матричного уравнения выглядит как

anqn-1T + bnqnT = qnTA.
Поскольку матрица Q ортогональна,
qnTqm = dnm.
Если теперь умножим обе части уравнения на qn, то найдем:
bn = qnTAqn,
которая однозначно известна при заданном qn. Кроме того, из уравнения следует:
anqn-1T = zn-1T, где обозначено zn-1T = qnTA - bnqnT,
последняя величина также известна. Таким образом,
an2 = zn-1Tzn-1, или an = |zn-1|, а также qn-1T = zn-1T/an.
(Поскольку согласно гипотезе an ненулевое). Аналогично можно по индукции доказать, что если нам известны qn, qn-1, ..., qn-j, а также все a, b и g вплоть до уровня (n - j), то мы можем определить все значения на уровне (n - (j+1)).

Для практического применения этой леммы предположим, что мы как-то нашли трехдиагональную матрицу A"s+1, такую что

A"s+1 = Q"sTA"sQ"s,
где Q"sT ортогональна и ее последняя строка совпадает с последней строкой матрицы QsT из оригинального QL-алгоритма. Тогда будем иметь Q"s = Qs и A"s+1 = As+1.

Далее, из оригинального алгоритма видно, что последняя строка QsT совпадает с последней строкой Pn-1(s). Но вспомним, что Pn-1(s) служит для плоской ротации с целью обнуления элемента (n-1,n) матрицы A - ks1. Простое вычисление показывает, что параметры c и s этой матрицы должны иметь вид:

c = (dn - ks)/sqrt(en2 + (dn - ks)2), s = (- en-1)/sqrt(en2 + (dn - ks)2).
Матрица Pn-1(s)AsPn-1(s)T является трехдиагональной с двумя дополнительными элементами:

( ( ( ( (
(
(
(
...
)
)
)
)
)
)
)
)
| ... |  |  |  |  |  |
| x | x | x |  |  |  |
|  | x | x | x | x |  |
|  |  | x | x | x |  |
|  |  | x | x | x |  |

Она должна быть приведена к трехдиагональной форме с помощью ортогональной матрицы, имеющей последнюю строку вида [0,0, ..., 0,1]; таким образом, последняя строка Q"sT останется равной строке Pn-1(s). Это можно сделать последовательностью преобразований Хаусхолдера или Гивенса. Для матрицы указанного вида преобразование Гивенса эффективнее. Производим поворот в плоскости (n-2,n-1) для обнуления элемента (n-2,n) (по симметрии, будет также обнулен элемент (n, n-2)). После этого трехдиагональная форма сохраняется за исключением дополнительных элементов (n-3, n-1) и (n-1, n-3). Эти в свою очередь обнуляются поворотом в плоскости (n-3, n-2) и т.д. Таким образом, требуется последовательность (n - 2) трансформаций Гивенса. В результате будем иметь

QsT = Q"sT = P"1(s)P"2(s)...P"n-2(s)Pn-1(s),
где P" -- матрицы трансформаций Гивенса, а Pn-1 -- плоская ротация из оригинального алгоритма. По значению QsT находим значение A на новой итерации. Заметим, что сдвиг ks производится неявно при вычислении параметров c и s.

Помещенная ниже программа tqli ("Tridiagonal QL Implicit"), алгоритмически основанная на опубликованных в [2, 3], на практике работает очень хорошо. Число итераций для первых собственных значений должно быть около 4 - 5, но помимо этого на этих итерациях также уменьшаются внедиагональные элементы правого нижнего угла. Таким образом, дальнейшие собственные значения выясняются значительно быстрее. Среднее число итераций на каждое собственное значение обычно составляет 1.3 - 1.6. Число действий на одной итерации имеет порядок O(n) с достаточно большим множителем, скажем ~20n. Таким образом, общее число действий по диагонализации оценивается как ~20n(1.5)n ~ 30n2 операций. Если требуются собственные вектора, то число дополнительных операций будет значительно больше; их количество оценивается как 3n3.

## Программа tqli

`#include <math.h>`
`#include "nrutil.h" /* Здесь определяются некоторые утилиты типа выделения памяти */`

`/* QL-алгоритм с неявными сдвигами для определения собственных значений (и собственных`
`векторов) действительной, симметричной, трехдиагональной матрицы. Эта матрица может`
`быть предварительно получена с помощью программы tred2. На входе d[1...n] содержит`
`диагональ исходной матрицы, на выходе - собственные значения. На входе e[1...n]`
`содержит поддиагональные элементы, начиная с e[2]. На выходе массив e разрушается.`
`При необходимости поиска только собственных значений в программе следует`
`закомментировать или удалить инструкции, необходимые только для поиска собственных`
`векторов. Если требуются собственные вектора трехдиагональной матрицы, массив`
`z[1...n][1...n] необходимо инициализировать на входе единичной матрицей. Если`
`требуются собственные вектора матрицы, сведенной к трехдиагональному виду с помощью`
`программы tred2, в массив z требуется загрузить соответствующий выход tred2. В`
`обоих случаях на выходе массив z возвращает матрицу собственных векторов, расположенных`
`по столбцам.`
`*/`

`/* максимальное число итераций */`
`#define MAXITER 30`

`void tqli(float *d, float *e, int n, float **z) {`
`int m,l,iter,i,k;`
`float s,r,p,g,f,dd,c,b;`
`/* удобнее будет перенумеровать элементы e */`
`for(i=2;i<=n;i++) e[i-1]=e[i];`
`e[n]=0.;`
`/* главный цикл идет по строкам матрицы */`
`for(l=1;l<=n;l++) {`
`/* обнуляем счетчик итераций для этой строки */`
`iter=0;`
`/* цикл проводится, пока минор 2х2 в левом верхнем углу начиная со строки l`
`не станет диагональным */`
`do {`
`/* найти малый поддиагональный элемент, дабы расщепить матрицу */`
`for(m=l;m<=n-1;m++) {`
`dd=fabs(d[m])+fabs(d[m+1]);`
`if((float)(fabs(e[m]+dd)==dd) break;`
`}`
`/* операции проводятся, если верхний левый угол 2х2 минора еще не диагональный */`
`if(m!=l) {`
`/* увеличить счетчик итераций и посмотреть, не слишком ли много. Функция`
`nerror завершает программу с диагностикой ошибки. */`
`if(++iter>=MAXITER) nerror("Too many iterations in tqli");`
`/* сформировать сдвиг */`
`g=(d[l+1]-d[l])/(2.*e[l]); r=hypot(1.,g);`
`/* здесь d_m - k_s */`
`if(g>=0.) g+=fabs(r);`
`else g-=fabs(r);`
`g=d[m]-d[l]+e[l]/g;`
`/* инициализация s,c,p */`
`s=c=1.; p=0.;`
`/* плоская ротация оригинального QL алгоритма, сопровождаемая ротациями`
`Гивенса для восстановления трехдиагональной формы */`
`for(i=m-1;i>=l;i--) {`
`f=s*e[i]; b=c*e[i];`
`e[i+1]=r=hypot(f,g);`
`/* что делать при малом или нулевом знаменателе */`
`if(r==0.) {d[i+1]-=p; e[m]=0.; break;}`
`/* основные действия на ротации */`
`s=f/r; c=g/r; g=d[i+1]-p; r=(d[i]-g)*s+2.c*b; d[i+1]=g+(p=s*r); g=c*r-b;`
`/* Содержимое следующего ниже цикла необходимо опустить, если`
`не требуются значения собственных векторов */`
`for(k=1;k<=n;k++) {`
`f=z[k][i+1]; z[k][i+1]=s*z[k][i]+c*f; z[k][i]=c*z[k][i]-s*f;`
`}`
`}`
`/* безусловный переход к новой итерации при нулевом знаменателе и недоведенной`
`до конца последовательности ротаций */`
`if(r==0. && i>=l) continue;`
`/* новые значения на диагонали и "под ней" */`
`d[l]-=p; e[l]=g; e[m]=0.;`
`}`
`} while(m!=l);`
`}`
`}`

К оглавлению
Сравни кредиты