|
Страница 9 из 10
Собственные вектора и значения трехдиагональной матрицы Вычисление характеристического полиномаПосле того, как исходная действительная симметричная матрица была сведена к трехдиагональной форме, одним из возможных способов отыскать ее собственные значения становится прямое вычисление корней характеристического полинома 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. |