Метод Якоби (вращений) для симметричных матриц с циклическим исключением и барьерами

Авторы описания: А.С.Галкина (входные и выходные данные, математическое описание алгоритма и др.), И.А.Плахов (ресурс параллелизма алгоритма, последовательная сложность алгоритма и др.), А.В.Фролов (общая редакция и правка страницы).

Метод Якоби — итерационный алгоритм для вычисления собственных значений и собственных векторов вещественной симметричной матрицы. Карл Густав Якоб Якоби, в честь которого назван этот метод, предложил его в 1846 году. Однако использоваться метод начал только в 1950-х годах с появлением компьютеров.

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Суть алгоритма заключается в том, чтобы для заданной симметрической матрице $A = A^{(0)}$ построить последовательность ортогонально подобных матриц $A^{(1)},A^{(2)},\dotsc,A^{(m)}$, сходящуюся к диагональной матрице, на диагонали которой стоят собственные значения $A$. Для построения этой последовательности применяется специально подобранная матрица вращения $J_i$, такая что норма наддиагональной части $\| A^{(i)} \|_{off} =\sqrt{\sum\limits_{1 \le j \lt k \le n} (a_{jk}^{(i)})^2}$ уменьшается при каждом двустороннем вращении матрицы $A^{(i+1)}={J_i}^TA^{(i)}J_i$.

1.2 Математическое описание алгоритма

Исходные данные: симметрическая матрица $A$ (элементы $a_{ij}$).

Вычисляемые данные: диагональная матрица $\Lambda$ (элементы $\lambda_{ij}$).

Норма наддиагональной части $\| A^{(i)} \|_{off} =\sqrt{\sum\limits_{1 \le j \lt k \le n} (a_{jk}^{(i)})^2}$ уменьшается при каждом двустороннем вращении матрицы $A^{(i+1)}={J_i}^TA^{(i)}J_i$.

Это достигается обнулением элемента матрицы $A^{(i)}$ в матрице $A^{(i+1)}$. Если он расположен в j-й строке и k-м столбце, то

                                                       $j$                             $k$

$J_i = \begin{matrix} \\ \\ \\ j \\ \\ k \\ \\ \\ \\ \end{matrix}$ $\begin{bmatrix} & 1 & & & & & & & & \\ & & 1 & & & & & & & \\ & & & \ddots & & & & & & \\ & & & & \cos(\theta) & & -\sin(\theta) & & & \\ & & & & & \ddots & & & & \\ & & & & \sin(\theta) & & \cos(\theta) & & & \\ & & & & & & & \ddots & & \\ & & & & & & & & 1 & \\ & & & & & & & & & 1 \\ \end{bmatrix}$

Если обозначить $s = \sin \theta$ и $c = \cos \theta$, то матрица $A^{(i+1)}$ состоит из следующих элементов, отличающихся от элементов $A^{(i)}$:

$\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}$

Можно выбрать $\theta$ так, чтобы $a_{jk}^{(i+1)} = 0$ и $a_{kj}^{(i+1)} = 0$. Отсюда следует равенство

$\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$.

Если $a_{jj}^{(i)} = a_{kk}^{(i)}$, то выбирается $\theta = \frac{\pi}{4}$, в противном случае

вводится $t = \frac{s}{c} = \operatorname{tg}(\theta)$ и тогда $t^2 - 2t\tau + 1 = 0$. Решение квадратного уравнения даёт $t = \frac{\operatorname{sign}(\tau)}{|\tau| + \sqrt{1+\tau^2}}, c = \frac{1}{\sqrt{1+t^2}}, s = tc$.

Вычисление останавливается, когда выполняются критерии близости к диагональной матрице. Это малость максимального по абсолютной величине внедиагонального элемента матрицы.

1.3 Вычислительное ядро алгоритма

Рассматривая отдельную итерацию, можно считать, что вычислительное ядро составляют множественные вычисления элементов матрицы $a_{jm}^{(i+1)} = a_{mj}^{(i+1)}$   и   $a_{km}^{(i+1)} = a_{mk}^{(i+1)}$,   $m \ne j,k$   в процессе применения матрицы поворота $J_i$ к матрице $A$:

$\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}$

каждое из которых повторяется по $(n-2)$ раза, а также вычисление элементов $a_{jj}^{(i+1)}$   и   $a_{kk}^{(i+1)}$:

$\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}$

1.4 Макроструктура алгоритма

Основную часть метода составляет процедура применения вращения к матрице $A$ (здесь и далее верхние индексы, содержащие номер матрицы, опускаются), которая в дальнейшем будет обозначена как Jakobi-Rotation$(A,j,k)$.

Эту процедуру, в свою очередь, можно разделить на две логические части:

1. Определение угла поворота $\theta$ по элементам матрицы $a_{jj}$, $a_{kk}$ и $a_{jk}$;
2. Поворот матрицы $A$ (изменяются лишь строки и столбцы, соответствующие индексам $j$ и  $k$).

1.5 Схема реализации последовательного алгоритма

В данной версии метода Якоби используется наиболее простой способ выбора параметров j, k — построчный циклический обход внедиагональных элементов матрицы $A$:

repeat
    for j=1 to n-1
        for k=j+1 to n
            выполнить процедуру Jakobi-Rotation (A,j,k) 
        end for
    end for
пока A не достаточно близка к диагональной матрице

Процедура Jakobi-Rotation (A,j,k) — это следующий алгоритм:

Если |a(j,k)| достаточно мал, вычисление заканчивается. В противном случае выполняются следующие действия:
     1. Если  a(j,j)==a(k,k), то угол theta = pi/4 .
        В остальных случаях находим параметры tau, t, c, s:
          
          tau = (a(j,j)-a(k,k)/2/a(j,k) 
          t = sign(tau)/(|tau|+sqrt(1+tau^2)) 
          c = 1/sqrt(1+t^2) 
          s = ct
          
     2. Выполняется поворот матрицы (изменяются лишь строки и столбцы, соответствующие индексам k и j):
          A = R*(j,k,theta)AR(j,k,theta), c = cos (theta), s = sin (theta)

В отличие от метода без барьеров в данной версии метода критерии малости меняются по мере приближения к его окончанию. С этим и связано название версии - "с барьерами", поскольку их выбор - довольно сложная процедура. Этим обусловлено то, что данное описание - не одного алгоритма, а группы алгоритмов, различающихся выбором последовательности барьеров.

1.6 Последовательная сложность алгоритма

Для осуществления одной итерации метода Якоби для матрицы размера $(n\times n)$ требуется выполнить:

Для вычислительного ядра —

• $4n+4$ умножений,
• $2n+2$ сложений.

Для остальной части алгоритма —

• $3$ умножения,
• $3$ деления,
• $5$ сложений (вычитаний),
• $2$ извлечения квадратного корня.

Умножения и сложения (вычитания) составляют основную часть алгоритма.

1.8 Ресурс параллелизма алгоритма

Для выполнения одной итерации требуется последовательно выполнить следующие действия:

1. вычислить $\tau$ (2 операции сложения, 1 операция деления)
2. вычислить $t$ (1 операция сравнения, 1 операция взятия модуля, 2 операции сложения, 1 операция умножения, 1 операция деления, 1 операция извлечения квадратного корня)
3. вычислить $c$ (1 операция сложения, 1 операция умножения, 1 операция деления, 1 операция извлечения квадратного корня) и $s$ (1 операция умножения)

После этого выполняется ярус, отвечающий за применение поворота к матрице $A$ с параллельным выполнением $2(n-2)$ операций вычисления элементов матрицы $a_{jm}^{(i+1)} = a_{mj}^{(i+1)}$   и   $a_{km}^{(i+1)} = a_{mk}^{(i+1)}$,   $m \ne j,k$ , а также 2 операций вычисления элементов $a_{jj}^{(i+1)}$   и   $a_{kk}^{(i+1)}$. Наибольшее количество операций содержится в вычислении последних двух элементов (по 6 умножений и 3 сложения).

1.9 Входные и выходные данные алгоритма

Входные данные: матрица $A$ (элементы $a_{ij}$). Дополнительные ограничения:

• $A$ – симметрическая матрица, т. е. $a_{ij}= a_{ji}, i, j = 1, \ldots, n$.

Объём входных данных: $\frac{n (n + 1)}{2}$ (в силу симметричности достаточно хранить только диагональ и над/поддиагональные элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом.

Выходные данные: вектор собственных значений $\Lambda$ (элементы $\lambda_{ii}$).

Объём выходных данных: $n$.

1.10 Свойства алгоритма

Метод Якоби является самым медленным из имеющихся алгоритмов вычисления собственных значений симметричной матрицы. Тем не менее, он способен вычислять малые собственные числа и отвечающие им собственные векторы с гораздо большей точностью, чем другие методы . Кроме того, он не предполагает первоначального приведения матрицы к трехдиагональной форме.

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным.

Вычислительная мощность алгоритма линейна.

Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности: число итераций зависит от входных данных и порогового значения.