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

Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Frolov и 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.2, 1.3), А. Кухтинов (разделы 1.4, 1.5, 1.9). Вклад в остальные разделы считаем равным.

1 ЧАСТЬ. Свойства и структура алгоритмов

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

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

Суть алгоритма — приведение симметрической матрицы к диагональной с собственными значениями на диагонали путем вращения.

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

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

Вычисляемые данные: диагональная матрица [math]\Lambda=\{\lambda_{ii}\}[/math]. Диагональные элементы [math]\Lambda[/math] — собственные значения матрицы [math]A[/math], столбцы — приближенные собственные вектора матрицы [math]A[/math].

По заданной матрице [math]A=A_0[/math] строится последовательность ортогонально подобных матриц [math]A_1, A_2,\ldots, A_m[/math], сходящихся к [math]\Lambda[/math].

[math]A_{i+1}={J_i}^TA_iJ_i[/math], где [math] J_i[/math] — ортогональная матрица, называемая вращением Якоби.

Матрица [math]J_i[/math] выбирается так, чтобы обнулить элементы [math](j,k)[/math] и [math](k,j)[/math] матрицы [math]A_{i+1}[/math]:

                                                                           [math]j[/math]                         [math]k[/math]

[math] \begin{align} J_i = R(j,k,\theta) = \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} \end{align} [/math]

Для угла вращения [math]\theta[/math]:

[math]\begin{bmatrix}a^{(i+1)}_{jj} & a^{(i+1)}_{jk} \\ a^{(i+1)}_{kj} & a^{(i+1)}_{kk}\end{bmatrix}= \begin{bmatrix}\cos\theta & -\sin\theta \\ \sin\theta &\cos\theta\end{bmatrix}^{T} \begin{bmatrix}a^{(i)}_{jj} & a^{(i)}_{jk} \\ a^{(i)}_{kj} & a^{(i)}_{kk}\end{bmatrix} \begin{bmatrix}\cos\theta & -\sin\theta \\ \sin\theta &\cos\theta\end{bmatrix}= \begin{bmatrix}\lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix},[/math] где [math]\lambda_1[/math] и [math]\lambda_2[/math] — собственные значения подматрицы [math]\begin{bmatrix}a^{(i)}_{jj} & a^{(i)}_{jk} \\ a^{(i)}_{kj} & a^{(i)}_{kk}\end{bmatrix}[/math].

Вычислим [math]s = \sin \theta[/math] и [math]c = \cos \theta[/math], исходя из предыдущего матричного равенства:

[math]\begin{bmatrix}\lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}= \begin{bmatrix}a_{jj}c^2+a_{kk}s^2+2sca_{jk} & sc(a_{kk} - a_{jj}) + a_{jk}(c^2-s^2) \\ sc(a_{kk} - a_{jj}) + a_{jk}(c^2-s^2) & a_{jj}s^2+a_{kk}c^2 - 2sca_{jk} \end{bmatrix}[/math].

Приравняем внедиагональный элемент нулю: [math]\dfrac{a_{jj}-a_{kk}}{2a_{jk}}=\dfrac{c^2-s^2}{2sc}=\dfrac{\cos2\theta}{\sin2\theta}=\operatorname{ctg}2\theta\equiv\tau[/math].

При [math]a_{jj}=a_{kk}[/math], [math]\theta=\dfrac{\pi}{4}[/math]. Тогда [math]s=\dfrac{\sqrt{2}}{2}[/math], [math]c=\dfrac{\sqrt{2}}{2}[/math]. В противном случае переходим к вычислениям [math]s, c[/math].

Положим [math]tg\theta=t=\dfrac{s}{c}[/math]. Тогда [math]t^2+2t\tau-1=0[/math].

Решая квадратное уравнение, получим: [math]t=\dfrac{\operatorname{sign}\tau}{\left|\tau\right|+\sqrt{1+\tau^2}}, c=\dfrac{1}{\sqrt{1+t^2}}, s=t\cdot c.[/math]

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

Вычислительное ядро алгоритма составляют вычисления элементов матрицы [math]A[/math] в процессе применения матрицы поворота [math]J[/math] :

1. Вычисления элементов [math]a_{jl}^{(i+1)} = a_{lj}^{(i+1)}[/math] и [math]a_{kl}^{(i+1)} = a_{lk}^{(i+1)}[/math] (каждое [math](n-2)[/math] раз):

[math]\begin{align} a_{jl}^{(i+1)} = a_{lj}^{(i+1)} = c \, a_{jl}^{(i)} - s \, a_{kl}^{(i)}, l \ne j,k;\\ a_{kl}^{(i+1)} = a_{lk}^{(i+1)} = s \, a_{jl}^{(i)} + c \, a_{kl}^{(i)}, l \ne j,k. \end{align}[/math]

2. Вычисления элементов [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 Макроструктура алгоритма

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

repeat
 выбрать пару индексов j, k
 обратиться к процедуре JacobiRotation(A,j,k)
пока А не станет достаточно близка к диагональной матрице

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

В классическом алгоритме Якоби для обнуления на текущей итерации выбирался наибольший элемент [math]a_{jk}[/math]. Поиск наибольшего элемента слишком замедляет процесс вычислений (нужно провести поиск среди [math]\dfrac{n^2-n}{2}[/math] элементов, прежде чем выполнить вращение стоимостью в [math]O(n)[/math] флопов), поэтому для практических вычислений выбор параметров [math]j[/math] и [math]k[/math] производится путем построчного циклического обхода внедиагональных элементов матрицы [math]A[/math].

Таким образом, алгоритм можно описать так:

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

Процедура JacobiRotation(A, j, k):

if [math]|a_{jk}|[/math] не слишком мал
   if [math]a_{jj}=a_{jk}[/math]
      [math]\theta=\dfrac{\pi}{4}[/math]
   end if
   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\cdot t
       \end{align}[/math]
   end else
[math]    
 A = J_i^T\cdot A\cdot J_i, \qquad
 J_i = R(j,k,\theta),\ c = \cos \theta,\ s = \sin \theta     
[/math]
end if

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

Для одной итерации вычислений собственных значений матрицы порядка [math]n[/math] методом Якоби требуется:

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

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

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

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

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

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

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

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

Опишем граф алгоритма аналитически и в виде рисунка.

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

1. Calculate c,s — вычисление [math]c, s[/math] согласно формулам в схеме реализации последовательного алгоритма.

Аргументы функции:

• [math]a^{(i+1)}_{jj}, a^{(i+1)}_{kk}[/math]:
  • Первая итерация — элементы входных данных: диагональные элементы матрицы [math]A[/math];
  • Последующие — результаты вычислений вершины второй группы.
• [math]a^{(i+1)}_{jk}[/math]
  • Первая итерация — элементы входных данных [math]a_{jk}[/math]
  • Последующие — результаты вычислений вершины третьей группы.

Результаты выполнения функции являются промежуточными данными алгоритма.

2. [math]A_{jj}, A_{kk}[/math] — вычисление [math]a^{(i+1)}_{jj}, a^{(i+1)}_{kk}[/math] согласно формулам в разделе Вычислительное ядро алгоритма.

Аргументы функции:

• [math]c, s[/math] — результаты вычислений вершины первой группы.
• [math]a^{(i)}_{jj}, a^{(i)}_{kk}[/math]:
  • Первая итерация — элементы входных данных [math]a_{jj}, a_{kk}[/math];
  • Последующие — результаты вычислений вершины второй группы.
• [math]a^{(i)}_{jk}[/math]:
  • Первая итерация — элементы входных данных [math]a_{jk}[/math];
  • Последующие — результаты вычислений вершины третьей группы.

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

3. [math]A_{jl}, A_{kl}[/math] — вычисление [math]a^{(i+1)}_{jl},a^{(i+1)}_{kl},l\ne j,k[/math] согласно формулам в разделе Вычислительное ядро алгоритма.

Аргументы функции:

• [math]c, s[/math] — результаты вычислений вершины первой группы.
• [math]a^{(i)}_{jl}, a^{(i)}_{kl}[/math]:
  • Первая итерация — элементы входных данных [math]a_{jl},l\ne j,k[/math];
  • Последующие — результаты вычислений вершины третьей группы.

Результат выполнения функции является промежуточным данным алгоритма. (Считается, что алгоритм вычисляет только собственные значения. При вычислении собственных векторов эти данные также были бы выходными)

Описанный граф для четырех итераций представлен на рисунке 1. Отображение дано на макроуровне без входных и выходных данных, дублирующие друг друга дуги показаны как одна. Внутренняя структура каждой вершины на i-той итерации представлена рисунками 2–4.

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

• 

  Рисунок 2. Структура вершины Calculate c,s.

• 

  Рисунок 3. Структура вершины [math]A_{jj}, A_{kk}[/math] (вычисление [math]a_{jj}, a_{kk}[/math]).

• 

  Рисунок 4. Структура вершины [math]A_{jl}, A_{kl}[/math] (вычисление [math]a_{jl},a_{kl},l\ne j,k[/math]). Количество входных и выходных значений, а также количество операций на каждом уровне указано слева.

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

Каждая итерация алгоритма для матрицы порядка [math]n[/math] заключается в последовательном выполнении следующих ярусов:

1. Вычисление [math]c, s[/math] (последовательно всего [math]4[/math] умножений, [math]3[/math] делений, [math]4[/math] сложений (вычитаний), [math]2[/math] операции извлечения квадратного корня, [math]1[/math] операция сравнения, [math]1[/math] операция взятия модуля)
2. Поворот матрицы [math]A[/math] с параллельным вычислением элементов матрицы [math]a_{jl, kl}^{(i+1)}[/math], (высота этой части яруса – 2 операции, ширина – [math]4(n-2)[/math]) и элементов [math]a_{jj,kk}^{(i+1)}[/math] (высота этой части яруса – 4 операции, ширина — 5).

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

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

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

Входные данные:

1. Плотная матрица [math]A=\{a_{ij}\}, A\in\R^{n\times n}[/math]. Дополнительные ограничения:

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

2. Число двойной точности [math]\epsilon[/math], определяющее допустимый порядок малости внедиагональных элементов.

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

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

Выходные данные: вектор собственных значений [math]\lambda_{i}[/math] длины [math]n[/math].

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

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

Из существующих алгоритмов вычисления собственных значений симметричной матрицы метод Якоби является самым медленным (в среднем для сходимости метода Якоби требуется от 5 до 10 циклов, что хуже, чем у конкурирующих алгоритмов). Его преимущество заключается в том, что он не предполагает первоначального приведения матрицы к трехдиагональной форме. Также метод до сих пор используется на практике для вычисления малых собственных чисел и отвечающих им собственных векторов, так как он делает это с гораздо большей точностью, чем другие методы.

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

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

2 ЧАСТЬ. Программная реализация алгоритма

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

Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета. Использовался компилятор gcc версии 5.2 c флагами -fopenmp и -std=c++0x.

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

• число процессоров [1 : 48] с шагом 4;
• размер матрицы [500 : 1500].

Рисунок 5. Зависимость времени выполнения алгоритма от размера матриц и количества процессоров.

Из рис. 6 видно, что при 20-24 процессорах производительность максимальная.

Рисунок 6. Зависимость времени выполнения алгоритма от количества процессоров для фиксированного размера матрицы 500*500.

Реализация на C++ с использованием технологии OpenMP.

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

Библиотека JACOBI_EIGENVALUE содержит последовательную реализацию алгоритма.

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