Участник:Galkina/Метод Якоби вычисления собственных значений симметричной матрицы
![]() | Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Frolov и ASA. |
Метод Якоби вычисления собственных значений симметричной матрицы | |
Последовательный алгоритм | |
Последовательная сложность | O(n^3) |
Объём входных данных | \frac{n (n + 1)}{2} |
Объём выходных данных | n |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | O(n^2) |
Ширина ярусно-параллельной формы | O(n) |
Авторы описания: А.С.Галкина (входные и выходные данные, математическое описание алгоритма и др.), И.А.Плахов (ресурс параллелизма алгоритма,последовательная сложность алгоритма и др.) Вклад авторов считать равноценным.
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Метод Якоби — итерационный алгоритм для вычисления собственных значений и собственных векторов вещественной симметричной матрицы. Карл Густав Якоб Якоби, в честь которого назван этот метод, предложил его в 1846 году, хотя использоваться он начал только в 1950-х годах с появлением компьютеров. Суть алгоритма заключается в том, чтобы по заданной симметрической матрице A = A_0 построить последовательность ортогонально подобных матриц A_1,A_2,\dotsc,A_m, сходящуюся к диагональной матрице, на диагонали которой стоят собственные значения. Для построения этой последовательности применяется специально подобранная матрица поворота J_i, такая что норма наддиагональной части of\!f(A)=\sqrt{\sum\limits_{1 \le j \lt k \le n} a_{jk}^2} уменьшается при каждом повороте матрицы A. Вычисление останавливается, когда угол поворота становится близок к нулю, либо когда максимальный внедиагональный элемент становится меньше заранее выбранного порогового значения.
1.2 Математическое описание алгоритма
Исходные данные: симметрическая матрица A (элементы a_{ij}).
Вычисляемые данные: диагональная матрица \Lambda (элементы \lambda_{ij}).
Матрица A_{i+1} получается из A_i по формуле A_{i+1}={J_i}^TA_iJ_i, где J_i — ортогональная матрица, называемая вращением Якоби. При подходящем выборе J_i матрица A_m для больших m будет близка к диагональной матрице \Lambda.
Матрица J_i выбирается так, чтобы сделать нулями пару внедиагональных элементов матрицы A_{i+1} [1].
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} состоит из следующих элементов:
- \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.
Выбор параметров j и k производится путем построчного циклического обхода внедиагональных элементов матрицы A.
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).
Эту процедуру, в свою очередь, можно разделить на две логические части:
- Определение угла поворота \theta по элементам матрицы a_{jj}, a_{kk} и a_{jk};
- Поворот матрицы A (изменяются лишь строки и столбцы, соответствующие индексам j и k).
1.5 Схема реализации последовательного алгоритма
Алгоритм можно описать следующим образом [1]:
1. Выбрать пару индексов j,k 2. Обратиться к процедуре Jakobi-Rotation[math](A,j,k)[/math] Если [math]A[/math] не достаточно близка к диагональной матрице, перейти к шагу 1.
Существует несколько способов выбора пары j,k. Наиболее простой и быстрый способ — построчный циклический обход внедиагональных элементов матрицы A:
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(A,j,k) — это следующий алгоритм:
Если [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]
1.6 Последовательная сложность алгоритма
Для осуществления одной итерации метода Якоби для матрицы размера (n\times n) требуется выполнить:
Для вычислительного ядра —
- 4n+4 умножений,
- 2n+2 сложений.
Для остальной части алгоритма —
- 3 умножения,
- 3 деления,
- 5 сложений (вычитаний),
- 2 извлечения квадратного корня.
Умножения и сложения (вычитания) составляют основную часть алгоритма.
Так как выбор индексов j и k осуществляется путем перебора внедиагональных элементов матрицы, для полного прохода потребуется \frac{n(n-1)}{2} итераций.
Таким образом, при классификации по последовательной сложности метод Якоби вычисления собственных значений симметричной матрицы относится к алгоритмам с кубической сложностью.
1.7 Информационный граф
Граф алгоритма состоит из трёх групп вершин.
Первой группе вершин (J) соответствует вычисление значений c и s.
Второй группе вершин (А) соответствует вычисление значений элементов a_{jm}^{(i+1)} = a_{mj}^{(i+1)} и a_{km}^{(i+1)} = a_{mk}^{(i+1)}, m \ne j,k.
Третьей группе вершин (B) соответствует вычисление значений элементов a_{jj}^{(i+1)} и a_{kk}^{(i+1)} .
1.8 Ресурс параллелизма алгоритма
Для выполнения одной итерации требуется последовательно выполнить следующие действия:
- вычислить \tau (2 операции сложения, 1 операция деления)
- вычислить t (1 операция сравнения, 1 операция взятия модуля, 2 операции сложения, 1 операция умножения, 1 операция деления, 1 операция извлечения квадратного корня)
- вычислить 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 сложения).
Суммарное количество итераций зависит от входных данных и заданного значения погрешности. Один "цикл" (полный проход по внедиагональным элементам) осуществляется за N = \frac{n(n-1)}{2} итераций. Асимптотически метод сходится квадратично [2]. Как правило, для сходимости метода Якоби требуется от 5 до 10 циклов, что хуже, чем у конкурирующих алгоритмов [1].
Таким образом, при классификации по высоте ЯПФ метод Якоби относится к алгоритмам со сложностью O(n^2) , а при классификации по ширине ЯПФ — к алгоритмам со сложностью O(n).
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 Свойства алгоритма
Метод Якоби является самым медленным из имеющихся алгоритмов вычисления собственных значений симметричной матрицы. Тем не менее, он способен вычислять малые собственные числа и отвечающие им собственные векторы с гораздо большей точностью, чем другие методы [1]. Кроме того, он не предполагает первоначального приведения матрицы к трехдиагональной форме.
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов для метода Якоби является линейным.
Вычислительная мощность алгоритма линейна.
Метод Якоби не детерминирован, так как является итерационным алгоритмом с выходом по точности: число итераций зависит от входных данных и порогового значения.
Дуги информационного графа, исходящие из вершин, соответствующих вычислениям значений c и s, образуют пучки рассылок линейной мощности.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
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.
На следующих рисунках приведены графики времени выполнения и производительности выбранной реализации алгоритма в зависимости от изменяемых параметров запуска
На Рис.4 видно, что при увеличении количества процессов производительность сначала увеличивается, а затем уменьшается. Время работы алгоритма резко увеличивается при увеличении размеров матрицы, что обусловлено большим количеством сообщений, обрабатываемых корневым процессом. Полученные результаты говорят о плохой масштабируемости алгоритма.
Реализация алгоритма доступна по ссылке.
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
Существует библиотека JACOBI_EIGENVALUE, реализующая последовательный вариант метода Якоби на языках программирования С, С++, FORTRAN77, FORTRAN90, MATLAB, Python. Данная библиотека распространяется по лицензии GNU LESSER GENERAL PUBLIC LICENSE.
Существует также параллельная реализация метода на платформе CUDA.
3 Литература
<references \>
- ↑ Перейти обратно: 1,0 1,1 1,2 1,3 Деммель Дж. Вычислительная линейная алгебра. - М.: МИР, 2001 - С. 244-248
- ↑ Wilkinson J. H. The Algebraic Eigenvalue Problem // Oxford University Press, UK, 1965, P. 270