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

Эта работа прошла предварительную проверку Дата последней правки страницы: 23.11.2016 Данная работа соответствует формальным критериям. Проверено ASA.

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

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

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

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

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

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

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

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

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

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

                                                 [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]\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].

Выбор параметров [math]j[/math] и [math]k[/math] производится путем построчного циклического обхода внедиагональных элементов матрицы [math]A[/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].

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

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

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

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

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

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

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

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

if [math]|a_{jk}|[/math] не слишком мал
     if [math]a_{jj}==a_{kk}[/math]
          [math]\theta=\frac{\Pi}{4}[/math]
     else
          [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]
     end if
     [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]    
end if

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

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

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

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

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

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

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

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

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

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

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

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

Второй группе вершин (А) соответствует вычисление значений элементов [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].

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

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

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

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

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

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

После этого выполняется ярус, отвечающий за применение поворота к матрице [math]A[/math] с параллельным выполнением [math]2(n-2)[/math] операций вычисления элементов матрицы [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] , а также 2 операций вычисления элементов [math]a_{jj}^{(i+1)} [/math]   и   [math]a_{kk}^{(i+1)} [/math]. Наибольшее количество операций содержится в вычислении последних двух элементов (по 6 умножений и 3 сложения).

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.4 Масштабируемость алгоритма и его реализации

Исследование масштабируемости проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского Государственного Университета имени М.В. Ломоносова. Алгоритм реализован на языке C++ с использованием средств MPI. Компиляция производилась с помощью команды mpicxx с использованием компилятора Intel с OpenMPI и библиотеки MKL. Были подключены следующие модули: openmpi/1.5.5-icc, intel/13.1.0, mkl/4.0.2.146. Программа была запущена в очередях test и regular4.

Значения изменяемых параметров запуска реализации алгоритма:

• число процессов: 2, 4, 8, 16, 32, 64;
• число элементов в матрице: 900, 3600, 8100, 14400, 22500, 32400, 44100, 57600, 72900, 90000.

На следующих рисунках приведены графики времени выполнения и производительности выбранной реализации алгоритма в зависимости от изменяемых параметров запуска

Рисунок 3. Изменение времени выполнения в зависимости от числа процессов и размера входной матрицы

Рисунок 4. Изменение производительности в зависимости от числа процессов и размера входной матрицы

На Рис.4 видно, что при увеличении количества процессов производительность сначала увеличивается, а затем уменьшается. Время работы алгоритма резко увеличивается при увеличении размеров матрицы, что обусловлено большим количеством сообщений, обрабатываемых корневым процессом. Полученные результаты говорят о плохой масштабируемости алгоритма.

Реализация алгоритма доступна по ссылке.

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

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

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

3 Литература

<references \>