Масштабирование применяется к несимметричным (неэрмитовым) матрицам A общего вида и состоит в таких преобразованиях матрицы, которые приводят к уменьшению ее нормы. Уменьшение нормы целесообразно осуществлять потому, что погрешность результатов при вычислении собственных значений пропорциональна евклидовой норме исходной матрицы. Непосредственно перед процессом масштабирования производится выделение "изолированных" собственных значений исходной матрицы с помощью предварительной перестановки строк и столбцов, т.е. вычисляется матрица
| T X Y |
PAP = | 0 B Z |
| 0 0 R |
Здесь P - матрица перестановок, матрицы T и R представляют собой верхние треугольные матрицы, диагональные элементы которых и образуют "изолированные" собственные значения матрицы A. Эти матрицы образуются в тех случаях, когда матрица A0, являющаяся внедиагональной частью матрицы A, имеет нулевые строки или столбцы (тогда соответствующие диагональные элементы матрицы A представляют собой ее собственные значения). Матрица B является квадратной, причем соответствующая ей матрица B0 не имеет нулевых строк или столбцов. Матрицы X, Y и Z являются прямоугольными. Заметим, что если необходимо определить только собственные значения матрицы A, то дальнейшие вычисления можно проводить только для матрицы B, а это значит, что простыми перестановками строк и столбцов можно снизить размерность решаемой задачи.
Затем строится такая вещественная невырожденная диагональная матрица D, что матрица D - 1 B D является масштабированной (сбалансированной), т.е. суммы абсолютных величин элементов соответственных ее строк и столбцов почти равны. Диагональные элементы матриц D является степенями основания системы счисления, используемой в арифметике с плавающей запятой на данной машине. После описанной процедуры масштабирования исходная матрица примет вид:
| T XD Y |
| 0 D-1 BD D-1 Z |
| 0 0 R |
Следует иметь в виду, что масштабирование представляет собой процесс, позволяющий уменьшить ошибки, возникающие только при вычислениях, выполняемых на последующих этапах решения проблемы собственных значений, но не может уменьшить погрешности (в результатах), связанные с неточностью задания элементов исходной матрицы.
Пусть A - произвольная матрица порядка n. QR - алгоритм определяется следующими соотношениями:
Ak = Qk Rk , Ak+1 = Rk Qk = Q*k Ak Qk , k = 1, ... ,
где A1 = A, Qk - унитарная, Rk - правая треугольная. На каждом шаге строится разложение текущей матрицы в произведение унитарной и правой треугольной, которые затем перемножаются в обратном порядке. Все матрицы A2, ..., Ak унитарно подобны исходной матрице. Если все собственные значения матрицы различны по модулю, то можно показать, что последовательность матриц {Ak} сходится к верхней треугольной матрице. В общем случае последовательность {Ak} сходится по форме к правой квазитреугольной матрице.
Одно из достоинств QR - алгоритма состоит в том, что он не меняет правую почти треугольную форму исходной матрицы. Поэтому обычно матрицу, для которой определяются собственные значения, предварительно приводят к верхней почти треугольной форме. Тогда и все матрицы Ak будут иметь тот же вид. Это предварительное преобразование очень эффективно, так как один шаг QR - алгоритма для почти треугольной матрицы выполняется примерно в n/2 раз быстрее, чем для полной матрицы.
Сходимость QR - алгоритма зависит от отношений модулей собственных значений и в случае плохого разделения их может быть медленной. Сходимость процесса может быть ускорена использованием сдвигов матрицы на каждом шаге.
QR - алгоритм со сдвигами определяется следующими соотношениями
A1 = A ,
Ak - vk I = Qk Rk , Ak+1 = Rk Qk + vk I , k = 1, ... .
Числа v1, ..., vk называются сдвигами.
Пусть матрица A является правой почти треугольной. Тогда
по - прежнему все матрицы Ak имеют такой же вид,
и все матрицы Ak унитарно подобны исходной
матрице A. Пусть выполнены k шагов QR - алгоритма
со сдвигом.
Предположим, что собственные значения
λ1 ( k ), ..., λn ( k )
занумерованы таким образом, что выполняются неравенства
k k
П | λ1(k) - vm | ≥ П | λ2(k) - vm | ≥ ...
m =1 m =1
k
≥ П | λn(k) - vm |
m =1
(верхний индекс k указывает на то, что нумерация, вообще говоря, зависит от k). Тогда можно показать, что поддиагональные элементы a (k)i+1, i матрицы Ak есть величины такого порядка:
a(k)i+1, i =
k
= O ( k2 (s-1) П [ ( λ(k)i+1 - vm ) / ( λ i(k) - vm ) ] ) ,
m =1
где s - максимальный порядок канонического ящика Жордана для матрицы A. Выбирая сдвиги подходящим образом, можно получить ускорение сходимости QR - алгоритма. Наиболее часто используются две стратегии сдвигов: сдвиг 1 - го порядка, определяемый по последнему диагональному элементу матрицы выполняемого шага, и сдвиг 2 - го порядка, определяемый по собственному значению подматрицы 2 - го порядка, расположенной в правом нижнем углу матрицы выполняемого шага. При этом скорость сходимости QR - алгоритма становится асимптотически квадратичной, а в отдельных случаях и кубической. Но ни одна из этих двух стратегий не гарантирует глобальную сходимость метода для случая матрицы общего вида.
Для вещественных матриц общего вида используется QR - алгоритм с двойным сдвигом, один шаг которого объединяет два шага стандартного QR - алгоритма соответственно с двумя вещественными сдвигами или с парой комплексно сопряженных сдвигов, что позволяет избежать использования комплексной арифметики.
QR - алгоритм эффективно модифицируется для симметричной трехдиагональной матрицы.
Если для вещественной симметричной трехдиагональной матрицы с ненулевыми внедиагональными элементами в качестве сдвига vk брать собственное значение матрицы 2 - го порядка, стоящей в нижнем углу матрицы Ak, которое ближе к последнему диагональному элементу (если оба собственных значения находятся на одинаковом расстоянии от этого элемента, то в качестве vk берется меньшее из них), то точный QR - алгоритм со сдвигами всегда сходится.
Скорость сходимости при этом всегда не хуже квадратичной и почти всегда лучше кубической. На практике для вычисления каждого собственного значения обычно бывает достаточно 1 - 5 итераций. Однако возможна и медленная сходимость, для предотвращения которой в соответствующие подпрограммы вставлено ограничение на число итераций. При этом, если для вычисления какого - либо собственного значения требуется более 30 итераций, то происходит выход из подпрограммы и выдается соответствующая диагностика.
Перечислим основные достоинства QR - алгоритма: алгоритм сохраняет почти треугольную форму исходной матрицы, что ведет к значительному сокращению обьема вычислений, алгоритм эффективно модифицируется в случае эрмитовых матриц, алгоритм устойчив к влиянию ошибок округления.
Пусть A - произвольная матрица. QL - алгоритм строит последовательность унитарных матриц Qk и левых треугольных матриц Lk по следующим формулам:
Ak = Qk Lk , Ak+1 = Lk Qk , k = 1, ... .
Свойства QL - агоритма аналогичны свойствам QR - алгорима. Различие состоит в том, что матрицы Ak в QL - алгоритме сходятся по форме к левой квазитреугольной матрице. Для QL - алгоритма могут быть применены аналогичные приемы ускорения. QR - алгоритм удобнее применять в тех случаях, когда наибольшие элементы матрицы расположены в верхнем левом углу, а QL - алгоритм - когда в нижнем правом.
Метод бисекций - эффективный метод решения частичной проблемы собственных значений вещественной симметричной трехдиагональной матрицы.
Пусть A - вещественная симметричная трехдиагональная матрица порядка n; α1, ..., αn - ее диагональные элементы, и β2, ..., βn - ее поддиагональные элементы; λ1 ≤ ... ≤ λn - ее собственные значения. Пусть λ - произвольное вещественное число. Для главных миноров матрицы A - λI имеют место следующие простые рекуррентные формулы:
p0(λ) = 1 , p1(λ) = α1 - λ ,
pi(λ) = (αi - λ) pi -1(λ) - β2i pi -2(λ) , i = 2, ..., n .
Известно, что число совпадений знаков в последовательности p0 (λ), ..., pn (λ) (если какое - либо pi (λ) равно нулю, то ему приписывается произвольный знак) равно числу собственных значений матрицы A больших λ. Это свойство позволяет, беря различные значения λ, теоретически локализовать любое собственное значение матрицы A с любой точностью. В самом деле, пусть, например, требуется определить k - oе собственное значение λk и пусть известно, что λk ∈ (a0, b0]. Взяв λ = (a0 + b0)/2, определяем, какому из полуинтервалов (a0, λ] или (λ, b0] принадлежит λk, границы полученного полуинтервала обозначаем (a1, b1], делим его пополам и т.д. Таким образом строится последовательность вложенных полуинтервалов, содержащих искомое собственное значение, причем длина каждого следующего в 2 раза меньше предыдущего.
При реальных вычислениях обычно вычисляют не сами pi, а их отношения. Каждое собственное значение метода бисекции может быть вычислено с точностью kε || A ||, где ε - отностельная машинная точность, || A || - спектральная норма, причем k не зависит от порядка матрицы. Наличие близких или кратных собственных значений не оказывает никакого влияния на реализацию метода.
Метод бисекции может быть непосредственно применен и к симметричной ленточной матрице, если ширина ленты невелика и требуется вычислить лишь несколько собственных значений. Как и для трехдиагональной матрицы вычисляются не сами главные миноры, а их отношения. Для этого можно использовать метод Гаусса с перестановками.
Различные модификации метода Ланцоша [ 5, 6 ] являются эффективными методами решения различных спектральных задач для больших симметричных разреженных или регулярной структуры матриц.
Приведем описание простого метода Ланцоша. Для заданной симметричной матрицы A порядка n и начального вектора v, v ≠ 0, строится последовательность векторов v1, ..., vj и последовательность трехдиагональных матриц T1, ..., Tj по следующим формулам:
v1 = v / (vT v)1/2 , β1 = 0 ,
ui = Avi - vi -1 βi , αi = uiT vi ,
T1 = α1 , i = 1 ,
| | 0 |
| | . |
Ti = | Ti - 1 | 0 | , i > 1 , i = 1, ..., j ,
| _______| βi |
| 0 ... 0 βi αi |
wi +1 = ui - vi αi , βi +1 = (wTi +1 wi +1)1/2 , vi +1 = wi +1/ βi +1
Метод Ланцоша может быть рассмотрен как реализация для спектральной задачи Az = λ z метода проекций Ритца - Галеркина на последовательность подпространств Крылова
K( j ) = span {v1, Av1, ..., Aj -1v1} , j = 1, ..., N ,
при этом числа Ритца совпадают с собственными значениями θi ( j ), i = 1, ..., j, матрицы Tj, получаемой на j - ом шаге метода Ланцоша, а векторы Ритца равны Vj si ( j ) , где si ( j ) - соответствующий θi ( j ) собственный вектор матрицы Tj, а Vj - матрица, столбцами которой являются векторы Ланцоша v1, ..., vj.
Метод Ланцоша удобен, во - первых, тем, что не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому метод Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены. Во - вторых, часто уже при небольшом ( j ~ 2 √N) числе шагов метода крайние пары Ритца ( θi ( j ), Vj si ( j )) оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны. В - третьих, для оценки точности по невязке получаемых на j - ом шаге приближений сами векторы Ритца yi ( j ) = Vj si ( j ) вычислять не требуется, т.к. выполнено соотношение
|| Ayi( j ) - θi( j ) yi( j ) || = βj +1 | {si( j )}j | ,
где {s}j - j - ая компонента вектора s.
К сожалению, простой метод Ланцоша не является устойчивым к влиянию ошибок округления. Но его модификация - метод Ланцоша с выборочной ортогонализацией - устойчива к влиянию ошибок округления и отличается от простого метода только тем, что на некоторых шагах метода перед вычислением βi + 1 и vi + 1 проводится ортогонализация вектора wi + 1 относительно некоторых векторов Ритца. Эта модификация используется для вычисления нескольких крайних собственных пар матриц рассматриваемого типа.
Иногда метод Ланцоша приходится использовать итерационным образом. Дело в том, что если требуется вычислить собственные векторы, то необходимо хранить вычисляемые векторы Ланцоша. Если работа при этом производится только с оперативной памятью, то при большом порядке исходной матрицы (а именно для таких матриц и предназначен метод Ланцоша) число шагов j, которые могут быть выполнены до того момента, как оперативная память будет исчерпана, будет очень незначительно, а следовательно, искомые собственные пары могут не успеть сойтись с требуемой точностью. В этом случае метод Ланцоша можно использовать итерационным образом: после выполнения заданного числа шагов процесс прерывается и, если еще не все требуемые пары Ритца сошлись с требуемой точностью, то формируется новый начальный вектор, являющийся линейной комбинацией векторов Ритца, еще не сошедшихся с заданной точностью, и процесс начинается заново.
Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого метода Ланцоша имеет смысл использовать блочный (ленточный) метод Ланцоша. В блочном методе Ланцоша в отличие от простого метода вместо одного начального вектора берется ортонормированная система из m, m > 1, векторов.
Аналогом векторов ui, vi, wi
простого метода в блочном являются прямоугольные матрицы размера
N * m, аналогом скаляров
βi, αi -
квадратные матрицы порядка m, аналогом трехдиагональной матрицы
Ti - ленточная матрица с полушириной ленты m.
Такие ленточные матрицы могут иметь собственные значения кратности
вплоть до m. Поэтому блочный метод с m начальными
вектрами используется, когда известно, что среди искомых
собственных значений многие имеют кратность m. При
использовании блочного метода все копии кратного собственного
значения будут найдены на одном шаге, в простом же методе Ланцоша
они определяются последовательно.
Блочный метод, как и простой, может быть дополнен выборочной
ортогонализацией.
Простой метод Ланцоша, дополненный специальным анализом получаемых результатов, используется для вычисления всех собственных значений без учета их кратности. Обычно при этом требуется выполнить 2N ÷ 5N шагов метода.
Пусть A - эрмитова матрица,
λ1, ..., λn -
собственные значения матрицы A,
z1, ..., zn -
соответствующие собственные векторы.
Пусть λj -
известное приближение к λj.
Рассмотрим решение системы уравнений
( A - λj I ) u1 = b ,
где b - нормированный начальный вектор. Если
n
b = ∑ γi zi ,
i =1
то точное решение рассматриваемой системы запишется в виде
n
ui = ∑ { γi / (λi - λj) } zi .
i =1
Обозначим
a = min | λi - λj | , ε = | λj - λj | . i ≠ j
Если ε << a, то сравнивая векторы b и u1, видим, что разложение u1 содержит значительно большую часть вектора zj, чем разложение вектора b. Для того чтобы u1 было хорошим приближением к zj, необходимо, чтобы γj не было малым. Решая теперь систему
( A - λj I ) u2 = u1 , получим n u2 = ∑ { γi / (λi - λj)2 } zi, i =1
т.е. разложение u2 содержит еще большую часть вектора zj, чем разложение u1.
Продолжая этот процесс, будем получать векторы все более близкие к zj, если только γj не было очень малым. На k - м шаге процесса получим
n
uk = ∑ { γi / (λi - λj)k } .
i =1
Перечислим некоторые трудности, возникающие при использовании метода обратных итераций:
| 1. | выбор начального вектора; |
| 2. | так как λj - приближенное собственное значение, то решаемая система линейных уравнений будет плохо обусловлена; |
| 3. | в случае кратных или близких собственных значений метод обратных итераций не обеспечивает ортогональность вычисляемых собственных векторов, поэтому требуется проводить их переортогонализацию. |