Классический метод Якоби (вращений) для симметричных матриц с выбором по всей матрице: различия между версиями
[досмотренная версия] | [досмотренная версия] |
Frolov (обсуждение | вклад) м (→Литература) |
Frolov (обсуждение | вклад) |
||
Строка 97: | Строка 97: | ||
# Определение угла поворота <math>\theta</math> по элементам матрицы <math>a_{jj}</math>, <math>a_{kk}</math> и <math>a_{jk}</math>; | # Определение угла поворота <math>\theta</math> по элементам матрицы <math>a_{jj}</math>, <math>a_{kk}</math> и <math>a_{jk}</math>; | ||
# Поворот матрицы <math>A</math> (изменяются лишь строки и столбцы, соответствующие индексам <math>j</math> и <math>k</math>). | # Поворот матрицы <math>A</math> (изменяются лишь строки и столбцы, соответствующие индексам <math>j</math> и <math>k</math>). | ||
+ | |||
+ | === Схема реализации последовательного алгоритма === | ||
+ | Алгоритм можно описать следующим образом: | ||
+ | |||
+ | 1. Выбрать пару индексов j,k так, чтобы элемент в этой позиции был максимальным по абсолютной величине из всех внедиагональных | ||
+ | 2. Обратиться к процедуре '''''Jakobi-Rotation<math>(A,j,k)</math>''''' | ||
+ | Если <math>A</math> не достаточно близка к диагональной матрице, перейти к шагу 1. | ||
+ | |||
+ | Процедура '''''Jakobi-Rotation<math>(A,j,k)</math>''''' — это следующий алгоритм: | ||
+ | Если <math>|a_{jk}|</math> достаточно мал, вычисление заканчивается. В противном случае выполняются следующие действия: | ||
+ | 1. Если <math>a_{jj}==a_{kk}</math>, то угол <math>\theta=\frac{\pi}{4}</math>. | ||
+ | В остальных случаях находим параметры <math>\tau,\ t,\ c,\ s</math>: | ||
+ | <math>\begin{align} | ||
+ | \tau &= \frac{a_{jj}-a_{kk}}{2\,a_{jk}} \\ | ||
+ | t &= \frac{sign(\tau)}{|\tau|+\sqrt{1+\tau^2}} \\ | ||
+ | c &= \frac{1}{\sqrt{1+t^2}} \\ | ||
+ | s &= c\,t | ||
+ | \end{align}</math> | ||
+ | 2. Выполняется поворот матрицы (изменяются лишь строки и столбцы, соответствующие индексам <math>i</math> и <math>j</math>): | ||
+ | <math>A = R^T(j,k,\theta)\cdot A\cdot R(j,k,\theta), \qquad J_i = R(j,k,\theta),\ c = \cos \theta,\ s = \sin \theta </math> | ||
== Литература == | == Литература == |
Версия 19:38, 22 декабря 2017
Авторы описания: А.С.Галкина, А.В.Фролов (общее редактирование и правка).
Метод Якоби — итерационный алгоритм для вычисления собственных значений и собственных векторов вещественной симметричной матрицы. Карл Густав Якоб Якоби, в честь которого назван этот метод, предложил его в 1846 году. Однако использоваться метод начал только в 1950-х годах с появлением компьютеров.
Содержание
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Суть алгоритма заключается в том, чтобы для заданной симметрической матрице [math]A = A^{(0)}[/math] построить последовательность ортогонально подобных матриц [math]A^{(1)},A^{(2)},\dotsc,A^{(m)}[/math], сходящуюся к диагональной матрице, на диагонали которой стоят собственные значения [math]A[/math]. Для построения этой последовательности применяется специально подобранная матрица вращения [math]J_i[/math], такая что норма наддиагональной части [math]\| A^{(i)} \|_{off} =\sqrt{\sum\limits_{1 \le j \lt k \le n} (a_{jk}^{(i)})^2}[/math] уменьшается при каждом двустороннем вращении матрицы [math]A^{(i+1)}={J_i}^TA^{(i)}J_i[/math].
1.2 Математическое описание алгоритма
Исходные данные: симметрическая матрица [math]A[/math] (элементы [math]a_{ij}[/math]).
Вычисляемые данные: диагональная матрица [math]\Lambda[/math] (элементы [math]\lambda_{ij}[/math]).
Норма наддиагональной части [math]\| A^{(i)} \|_{off} =\sqrt{\sum\limits_{1 \le j \lt k \le n} (a_{jk}^{(i)})^2}[/math] уменьшается при каждом двустороннем вращении матрицы [math]A^{(i+1)}={J_i}^TA^{(i)}J_i[/math].
Это достигается выбором максимального по абсолютной величине элемента матрицы [math]A^{(i)}[/math] и его обнулением[1] в матрице [math]A^{(i+1)}[/math]. Если он расположен в j-й строке и k-м столбце, то
[math]j[/math] [math]k[/math]
[math] J_i = \begin{matrix} \\ \\ \\ j \\ \\ k \\ \\ \\ \\ \end{matrix} [/math] [math] \begin{bmatrix} & 1 & & & & & & & & \\ & & 1 & & & & & & & \\ & & & \ddots & & & & & & \\ & & & & \cos(\theta) & & -\sin(\theta) & & & \\ & & & & & \ddots & & & & \\ & & & & \sin(\theta) & & \cos(\theta) & & & \\ & & & & & & & \ddots & & \\ & & & & & & & & 1 & \\ & & & & & & & & & 1 \\ \end{bmatrix} [/math]
Если обозначить [math]s = \sin \theta[/math] и [math]c = \cos \theta[/math], то матрица [math]A^{(i+1)}[/math] состоит из следующих элементов, отличающихся от элементов [math]A^{(i)}[/math]:
- [math]\begin{align} a_{jj}^{(i+1)} &= c^2\, a_{jj}^{(i)} - 2\, s c \,a_{jk}^{(i)} + s^2\, a_{kk}^{(i)} \\ a_{kk}^{(i+1)} &= s^2 \,a_{jj}^{(i)} + 2 s c\, a_{jk}^{(i)} + c^2 \, a_{kk}^{(i)} \\ a_{jk}^{(i+1)} &= a_{kj}^{(i+1)} = (c^2 - s^2 ) \, a_{jk}^{(i)} + s c \, (a_{kk}^{(i)} - a_{jj}^{(i)} ) \\ a_{jm}^{(i+1)} &= a_{mj}^{(i+1)} = c \, a_{jm}^{(i)} - s \, a_{km}^{(i)} & m \ne j,k \\ a_{km}^{(i+1)} &= a_{mk}^{(i+1)} = s \, a_{jm}^{(i)} + c \, a_{km}^{(i)} & m \ne j,k \\ a_{ml}^{(i+1)} &= a_{ml}^{(i)} &m,l \ne j,k \end{align}[/math]
Можно выбрать [math]\theta[/math] так, чтобы [math]a_{jk}^{(i+1)} = 0[/math] и [math]a_{kj}^{(i+1)} = 0[/math]. Отсюда следует равенство
- [math] \frac{a_{jj}^{(i)} - a_{kk}^{(i)}}{2 a_{jk}^{(i)}} = \frac{c^2 - s^2}{2sc} = \frac{\cos (2\theta)}{\sin (2\theta)} = \operatorname{ctg}(2\theta) \equiv \tau [/math].
Если [math] a_{jj}^{(i)} = a_{kk}^{(i)}[/math], то выбирается [math]\theta = \frac{\pi}{4}[/math], в противном случае
вводится [math]t = \frac{s}{c} = \operatorname{tg}(\theta)[/math] и тогда [math]t^2 - 2t\tau + 1 = 0[/math]. Решение квадратного уравнения даёт [math]t = \frac{\operatorname{sign}(\tau)}{|\tau| + \sqrt{1+\tau^2}}, c = \frac{1}{\sqrt{1+t^2}}, s = tc[/math].
Вычисление останавливается, когда выполняются критерии близости к диагональной матрице. Это малость максимального по абсолютной величине внедиагонального элемента матрицы.
1.3 Вычислительное ядро алгоритма
Рассматривая отдельную итерацию, можно считать, что вычислительное ядро составляют множественные вычисления элементов матрицы [math]a_{jm}^{(i+1)} = a_{mj}^{(i+1)}[/math] и [math]a_{km}^{(i+1)} = a_{mk}^{(i+1)}[/math], [math]m \ne j,k[/math] в процессе применения матрицы поворота [math]J_i[/math] к матрице [math]A[/math]:
- [math]\begin{align} a_{jm}^{(i+1)} &= a_{mj}^{(i+1)} = c \, a_{jm}^{(i)} - s \, a_{km}^{(i)} & m \ne j,k \\ a_{km}^{(i+1)} &= a_{mk}^{(i+1)} = s \, a_{jm}^{(i)} + c \, a_{km}^{(i)} & m \ne j,k, \end{align}[/math]
каждое из которых повторяется по [math] (n-2) [/math] раза, а также вычисление элементов [math]a_{jj}^{(i+1)} [/math] и [math]a_{kk}^{(i+1)} [/math]:
- [math]\begin{align} a_{jj}^{(i+1)} &= c^2\, a_{jj}^{(i)} - 2\, s c \,a_{jk}^{(i)} + s^2\, a_{kk}^{(i)} \\ a_{kk}^{(i+1)} &= s^2 \,a_{jj}^{(i)} + 2 s c\, a_{jk}^{(i)} + c^2 \, a_{kk}^{(i)} \end{align}[/math]
1.4 Макроструктура алгоритма
Основную часть метода составляет процедура применения вращения к матрице [math]A[/math] (здесь и далее верхние индексы, содержащие номер матрицы, опускаются), которая в дальнейшем будет обозначена как Jakobi-Rotation[math](A,j,k)[/math].
Эту процедуру, в свою очередь, можно разделить на две логические части:
- Определение угла поворота [math]\theta[/math] по элементам матрицы [math]a_{jj}[/math], [math]a_{kk}[/math] и [math]a_{jk}[/math];
- Поворот матрицы [math]A[/math] (изменяются лишь строки и столбцы, соответствующие индексам [math]j[/math] и [math]k[/math]).
1.5 Схема реализации последовательного алгоритма
Алгоритм можно описать следующим образом:
1. Выбрать пару индексов j,k так, чтобы элемент в этой позиции был максимальным по абсолютной величине из всех внедиагональных 2. Обратиться к процедуре Jakobi-Rotation[math](A,j,k)[/math] Если [math]A[/math] не достаточно близка к диагональной матрице, перейти к шагу 1.
Процедура Jakobi-Rotation[math](A,j,k)[/math] — это следующий алгоритм:
Если [math]|a_{jk}|[/math] достаточно мал, вычисление заканчивается. В противном случае выполняются следующие действия: 1. Если [math]a_{jj}==a_{kk}[/math], то угол [math]\theta=\frac{\pi}{4}[/math]. В остальных случаях находим параметры [math]\tau,\ t,\ c,\ s[/math]: [math]\begin{align} \tau &= \frac{a_{jj}-a_{kk}}{2\,a_{jk}} \\ t &= \frac{sign(\tau)}{|\tau|+\sqrt{1+\tau^2}} \\ c &= \frac{1}{\sqrt{1+t^2}} \\ s &= c\,t \end{align}[/math] 2. Выполняется поворот матрицы (изменяются лишь строки и столбцы, соответствующие индексам [math]i[/math] и [math]j[/math]): [math]A = R^T(j,k,\theta)\cdot A\cdot R(j,k,\theta), \qquad J_i = R(j,k,\theta),\ c = \cos \theta,\ s = \sin \theta [/math]
2 Литература
- ↑ Деммель Дж. Вычислительная линейная алгебра. - М.: МИР, 2001.