D = VTAV, где V = P1P2P3...,
а Pi -- матрицы ротации Якоби. Столбцы матрицы V являются собственными векторами. Они вычисляются рекуррентно, на каждой стадии по формулам:
V' = VPi ,
где первоначально V -- единичная матрица. Для преобразования элементов V используются формулы:
v'rs = vrs (s!=p, s!=q)
v'rp = cvrp - svrq
v'rq = svrp + cvrq
Для уменьшения ошибки округления вышеприведенные формулы следует переписать в терминах g (как выше для элементов матрицы A).
Единственным остающимся вопросом теперь является стратегия, по которой нужно перебирать уничтожаемые элементы. Оригинальный алгоритм Якоби от 1846 года на каждой стадии искал в верхней треугольной области наибольший по модулю элемент и обнулял его. Это весьма рациональная стратегия для ручного счета, но она неприемлема для компьютера, поскольку делает число операций на каждой стадии ротации Якоби порядка N2 вместо N.
Для наших целей лучшей стратегией является циклический метод Якоби, где элементы удаляются в строгом порядке. Например, можно просто проходить по строкам: P12 , P13 , ... P1n ; P23 ,...,P2n ;... Можно показать, что скорость сходимости будет в общем случае квадратичной, как для оригинального, так и для циклического метода Якоби. Назовем одно такое множество последовательных n(n-1)/2 преобразований Якоби проходом.
Программа, помещенная ниже, содержит две тонкости:
На первых трех проходах мы проводим ротацию относительно индекса (pq) только тогда, когда |apq|>e для некоторого порогового значения e = S0/(5n2), где S0 -- сумма модулей внедиагональных элементов верхней треугольной матрицы.
После четырех проходов в том случае, когда |apq|<<|app| и |apq|<<|aqq|, полагаем |apq|=0 и пропускаем ротацию. Критерием сравнения является превышение диагональных элементов в 10D+2 раз, где D -- число значащих десятичных цифр.
В приведенной ниже программе симметричная матрица a размером (n x n) на входе загружена в массив a[1...n][1...n]. На выходе наддиагональные элементы a разрушаются, но диагональные и поддиагональные остаются на месте с тем, чтобы сохранить информацию об исходной симметричной матрице. Вектор d[1...n] возвращает собственные значения a. Во время вычисления он содержит текущую диагональ. Матрица v[1...n][1...n] возвращает набор нормализованных собственных векторов, относящихся к d[k] для ее k-го столбца. Параметр nrot -- число ротаций Якоби, необходимых для достижения сходимости.
Типичные матрицы требуют от 6 до 10 проходов для достижения сходимости, или от 3n2 до 5n2 ротаций. Каждая ротация требует порядка 4n операций умножения и сложения, так что общие затраты имеют порядок от 12n3 до 20n3 операций. При вычислении собственных векторов наряду с собственными значениями, число операций на ротации возрастает с 4n до 6n, что означает превышение только на 50 процентов.
Do'stlaringiz bilan baham: |