Метод Якоби

Практика №8

Метод Якоби состоит из цепочки ортогональных преобразований подобия матрицы. Каждое преобразование (ротация Якоби) -- это плоский поворот с целью обнуления одного из внедиагональных элементов матрицы. Последовательные преобразования не сохраняют уже установленные нулевые элементы, но вместе с тем внедиагональные элементы становятся меньше и меньше до тех пор, пока матрица не станет диагональной с точностью до машинного нуля. Накопление в процессе преобразований произведения трансформационных матриц дает матрицу собственных векторов, в то время как диагональные элементы являются собственными значениями.

Метод Якоби всегда является робустным для действительных симметричных матриц. Для матриц размера, скажем больше 10, этот алгоритм более медленный, чем метод QR (см. ниже). Однако алгоритм Якоби значительно проще других более эффективных методов, таким образом он рекомендуется в случаях, когда время выполнения не является критическим.

Базовая матрица ротации Якоби P имеет вид:

Все диагональные элементы матрицы ротации Якоби равны 1, за исключением двух элементов в строках (и столбцах) p и q. Все внедиагональные элементы равны нулю, за исключением двух элементов, равных s и (-s). Числа c и s являются косинусом и синусом угла поворота f, так что c²+s²=1.

Ротация Якоби преобразует матрицу A согласно формуле:

A' = P _pq^T AP _pq.

Идея метода Якоби -- в том, чтобы постараться обнулить внедиагональные элементы серией ротаций. Для того, чтобы обнулилось a'_pq, угол ротации должен быть следующим:

q = cot(2f) = (c² - s²)/(2cs) = (a_qq- a_pp)/(2a_pq).

Полагая t=s/c, для определения t прийдем к следующему уравнению:

t² + 2tq - 1 = 0.

Меньший по модулю корень этого уравнения соответствует вращению на угол, меньший p/4; выбор этого угла на каждой из стадий дает наиболее устойчивый алгоритм. Используя формулу для решения квадратного уравнения с дискриминантом в знаменателе, меньший корень определяется как:

t = sign (q)/(|q|+ sqrt (q²+1)).

Если величина q такова, что q² даст переполнение, полагаем t = 1/(2q). По значению t определяются c и s:

c = 1/ sqrt (t²+1), s = tc.

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

a'_pp = a_pp - ta_pq
a'_rp = a_rp - s(a_rq + ga_rp)
a'_rq = a_rq + s(a_rp - ga_rq),

где t (тангенс половины угла поворота) определяется, как t = s/(1+c).

Сходимость метода Якоби можно оценить, рассматривая сумму квадратов внедиагональных элементов S; уравнения преобразований дают для этой суммы после трансформации следующее соотношение:

S = S_{r! =s}(|a_rs|²); S' = S - 2|a_pq|².

Поскольку преобразования ортогональные, сумма диагональных элементов при этом возрастает соответственно на |a_pq|². Сумма же квадратов внедиагональных элементов монотонно убывает. Поскольку она ограничена снизу нулем и поскольку мы можем каждый раз выбирать какой хотим элемент a_pq для обнуления, сумма будет стремиться к нулю.

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

D = V ^T AV, где V = P ₁ P ₂ P ₃...,

а P _i -- матрицы ротации Якоби.

Столбцы матрицы V являются собственными векторами. Они вычисляются рекуррентно, на каждой стадии по формулам:

V' = VP _i,

где первоначально V -- единичная матрица. Для преобразования элементов V используются формулы:

v'_rs = v_rs (s! =p, s! =q)
v'_rp = cv_rp - sv_rq
v'_rq = sv_rp + cv_rq

И так далее

на 14 шаге получим матрицу