Участник:Galkina/Метод Якоби вычисления собственных значений симметричной матрицы

Метод Якоби вычисления собственных значений симметричной матрицы
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(n^3) }[/math]
Объём входных данных	[math]\displaystyle{ \frac{n (n + 1)}{2} }[/math]
Объём выходных данных	[math]\displaystyle{ n }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ O(n^2) }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ O(n) }[/math]

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

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

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

Метод Якоби — итерационный алгоритм для вычисления собственных значений и собственных векторов вещественной симметричной матрицы. Карл Густав Якоб Якоби, в честь которого назван этот метод, предложил его в 1846 году, хотя использоваться он начал только в 1950-х годах с появлением компьютеров. Суть алгоритма заключается в том, чтобы по заданной симметрической матрице [math]\displaystyle{ A = A_0 }[/math] построить последовательность ортогонально подобных матриц [math]\displaystyle{ A_1,A_2,\dotsc,A_m }[/math], сходящуюся к диагональной матрице, на диагонали которой стоят собственные значения. Для построения этой последовательности применяется специально подобранная матрица поворота [math]\displaystyle{ J_i }[/math], такая что норма наддиагональной части [math]\displaystyle{ of\!f(A)=\sqrt{\sum\limits_{1 \le j \lt k \le n} a_{jk}^2} }[/math] уменьшается при каждом повороте матрицы [math]\displaystyle{ A }[/math]. Вычисление останавливается, когда угол поворота становится близок к нулю, либо когда максимальный внедиагональный элемент становится меньше заранее выбранного порогового значения.

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

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

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

Матрица [math]\displaystyle{ A_{i+1} }[/math] получается из [math]\displaystyle{ A_i }[/math] по формуле [math]\displaystyle{ A_{i+1}={J_i}^TA_iJ_i }[/math], где [math]\displaystyle{ J_i }[/math] — ортогональная матрица, называемая вращением Якоби. При подходящем выборе [math]\displaystyle{ J_i }[/math] матрица [math]\displaystyle{ A_m }[/math] для больших m будет близка к диагональной матрице [math]\displaystyle{ \Lambda }[/math].

Матрица [math]\displaystyle{ J_i }[/math] выбирается так, чтобы сделать нулями пару внедиагональных элементов матрицы [math]\displaystyle{ A_{i+1} }[/math].

                                                 [math]\displaystyle{ j }[/math]                           [math]\displaystyle{ k }[/math]

[math]\displaystyle{ J_i = \begin{matrix} \\ \\ \\ j \\ \\ k \\ \\ \\ \\ \end{matrix} }[/math] [math]\displaystyle{ \begin{bmatrix} & 1 & & & & & & & & \\ & & 1 & & & & & & & \\ & & & \ddots & & & & & & \\ & & & & \cos(\theta) & & -\sin(\theta) & & & \\ & & & & & \ddots & & & & \\ & & & & \sin(\theta) & & \cos(\theta) & & & \\ & & & & & & & \ddots & & \\ & & & & & & & & 1 & \\ & & & & & & & & & 1 \\ \end{bmatrix} }[/math]

Обозначим [math]\displaystyle{ s = \sin \theta }[/math] и [math]\displaystyle{ c = \cos \theta }[/math]. Тогда матрица [math]\displaystyle{ A_{i+1} }[/math] состоит из следующих элементов:

[math]\displaystyle{ \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]\displaystyle{ \theta }[/math] так, чтобы [math]\displaystyle{ a_{jk}^{(i+1)} = 0 }[/math] и [math]\displaystyle{ a_{kj}^{(i+1)} = 0 }[/math]. Отсюда получим равенство

[math]\displaystyle{ \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{tg}(2\theta) \equiv \tau }[/math].

Положим [math]\displaystyle{ t = \frac{s}{c} = \operatorname{ctg}(\theta) }[/math] и заметим, что [math]\displaystyle{ t^2 - 2t\tau + 1 = 0 }[/math]. Решая квадратное уравнение, находим [math]\displaystyle{ t = \frac{\operatorname{sign}(\tau)}{|\tau| + \sqrt{1+\tau^2}}, c = \frac{1}{\sqrt{1+t^2}}, s = tc }[/math].

Выбор параметров [math]\displaystyle{ j }[/math] и [math]\displaystyle{ k }[/math] производится путем построчного циклического обхода внедиагональных элементов матрицы [math]\displaystyle{ A }[/math].

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

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

[math]\displaystyle{ \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]\displaystyle{ (n-2) }[/math] раза, а также вычисление элементов [math]\displaystyle{ a_{jj}^{(i+1)} }[/math]   и   [math]\displaystyle{ a_{kk}^{(i+1)} }[/math]:

[math]\displaystyle{ \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]

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

Основную часть метода составляет процедура применения вращения к матрице [math]\displaystyle{ A }[/math], которая в дальнейшем будет обозначена как Jakobi-Rotation[math]\displaystyle{ (A,j,k) }[/math].

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

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

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

Алгоритм можно описать следующим образом:

repeat
    выбрать пару индексов j,k
    обратиться к процедуре Jakobi-Rotation[math]\displaystyle{ (A,j,k) }[/math]
пока [math]\displaystyle{ A }[/math] не достаточно близка к диагональной матрице

Существует несколько способов выбора пары [math]\displaystyle{ j,k }[/math]. Наиболее простой и быстрый способ — построчный циклический обход внедиагональных элементов матрицы [math]\displaystyle{ A }[/math]:

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

Процедура Jakobi-Rotation[math]\displaystyle{ (A,j,k) }[/math] — это следующий алгоритм:

if [math]\displaystyle{ |a_{jk}| }[/math] не слишком мал
    [math]\displaystyle{ \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 \\
      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 
     \end{align} }[/math]
end if

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

Для осуществления одной итерации метода Якоби для матрицы размера [math]\displaystyle{ (n\times n) }[/math] требуется выполнить:

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

• [math]\displaystyle{ 4n+4 }[/math] умножений,
• [math]\displaystyle{ 2n+2 }[/math] сложений.

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

• [math]\displaystyle{ 3 }[/math] умножения,
• [math]\displaystyle{ 3 }[/math] деления,
• [math]\displaystyle{ 5 }[/math] сложений (вычитаний),
• [math]\displaystyle{ 2 }[/math] извлечения квадратного корня.

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

Так как выбор индексов [math]\displaystyle{ j }[/math] и [math]\displaystyle{ k }[/math] осуществляется путем перебора внедиагональных элементов матрицы, для полного прохода потребуется [math]\displaystyle{ \frac{n(n-1)}{2} }[/math] итераций.

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

Информационный граф

Граф алгоритма состоит из трёх групп вершин.

Первой группе вершин (J) соответствует вычисление значений c и s.

Второй группе вершин (А) соответствует вычисление значений элементов [math]\displaystyle{ a_{jm}^{(i+1)} = a_{mj}^{(i+1)} }[/math]   и   [math]\displaystyle{ a_{km}^{(i+1)} = a_{mk}^{(i+1)} }[/math],   [math]\displaystyle{ m \ne j,k }[/math].

Третьей группе вершин (B) соответствует вычисление значений элементов [math]\displaystyle{ a_{jj}^{(i+1)} }[/math]   и   [math]\displaystyle{ a_{kk}^{(i+1)} }[/math].

Рисунок 1. Граф одной итерации алгоритма без отображения входных и выходных данных.

Рисунок 2. Внутренний граф вершин J с входными параметрами [math]\displaystyle{ x = a_{jj}^{(i)}, y = a_{kk}^{(i)}, z = a_{jk}^{(i)} }[/math].

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

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

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

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

Суммарное количество итераций зависит от входных данных и заданного значения погрешности. Один "цикл" (полный проход по внедиагональным элементам) осуществляется за [math]\displaystyle{ N = \frac{n(n-1)}{2} }[/math] итераций. Асимптотически метод сходится квадратично. Как правило, для сходимости метода Якоби требуется от 5 до 10 циклов, что хуже, чем у конкурирующих алгоритмов.

Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью [math]\displaystyle{ O(n^2) }[/math] , а при классификации по ширине ЯПФ — к алгоритмам со сложностью [math]\displaystyle{ O(n) }[/math].

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

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

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

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

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

Объём выходных данных: [math]\displaystyle{ n }[/math].

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

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

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

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

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

Дуги информационного графа, исходящие из вершин, соответствующих вычислениям значений [math]\displaystyle{ c }[/math] и [math]\displaystyle{ s }[/math], образуют пучки рассылок линейной мощности.

Программная реализация алгоритма

Существующие реализации алгоритма

Существует библиотека JACOBI_EIGENVALUE, реализующая последовательный вариант метода Якоби на языках программирования С, С++, FORTRAN77, FORTRAN90, MATLAB, Python. Данная библиотека распространяется по лицензии GNU LESSER GENERAL PUBLIC LICENSE.

Существует также параллельная реализация метода на платформе CUDA.

Литература

<references \>