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

Метод Якоби вычисления собственных значений симметричной матрицы
Последовательный алгоритм
Последовательная сложность	[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 Общее описание алгоритма

Метод Якоби вычисления собственных значений симметричной матрицы — итерационный алгоритм, названный в честь предложившего его для конкретной матрицы 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]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] не слишком мал
   [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 \\
      A &= J_i^T\cdot A\cdot J_i, \qquad
      J_i = R(j,k,\theta),\ c = \cos \theta,\ s = \sin \theta 
     \end{align}
[/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. Off>eps — сравнение матрицы с диагональной: вычисление нормы наддиагональных элементов [math]of\!f(A)=\sqrt{\sum\limits_{1 \le j \lt k \le n} a_{jk}^2}[/math] и сравнение с [math]\epsilon[/math].

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

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

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

2. 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]
  • Последующие — результаты вычислений вершины четвертой группы.

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

3. [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];
  • Последующие — результаты вычислений вершины четвертой группы.

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

4. [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. Граф алгоритма с отображением входных и выходных данных.

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

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

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

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

1. Поворот матрицы [math]A[/math] с параллельным выполнением [math]2(n-2)[/math] операций вычисления элементов матрицы [math]a_{jl}^{(i+1)} = a_{lj}^{(i+1)}[/math]   и   [math]a_{kl}^{(i+1)} = a_{lk}^{(i+1)}[/math],   [math]m \ne j,k[/math], а также [math]2[/math] операций вычисления элементов [math]a_{jj}^{(i+1)} [/math]   и   [math]a_{kk}^{(i+1)} [/math].

Суммарное количество итераций зависит от входных данных и заданного значения погрешности [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 циклов, что хуже, чем у конкурирующих алгоритмов). Его преимущество заключается в том, что он не предполагает первоначального приведения матрицы к трехдиагональной форме. Также метод до сих пор используется на практике для вычисления малых собственных чисел и отвечающих им собственных векторов, так как он делает это с гораздо большей точностью, чем другие методы.

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

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

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

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

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

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

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