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

Разложение Холецкого
Последовательный алгоритм
Последовательная сложность	[math][/math]
Объём входных данных	[math][/math]
Объём выходных данных	[math][/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math][/math]
Ширина ярусно-параллельной формы	[math][/math]

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

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

1.1.1 Симметричность и положительная определённость матрицы

Симметричность матрицы позволяет хранить и вычислять только чуть больше половины её элементов, что почти вдвое экономит как необходимые для вычислений объёмы памяти, так и количество операций в сравнении с, например, разложением по методу Гаусса.

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

Пусть [math]A = A_0[/math] — симметричная матрица. В методе Якоби для [math]A_0[/math] строится последовательность ортогонально подобных матриц [math]A_1, A_2,...[/math], которая в конечном счёте сходится к диагональной матрице, у которой диагональными элементами являются все собственные значения матрицы [math]A[/math]. Матрица [math]A_{i}[/math] получается из матрицы [math]A_{i-1}[/math] по формуле

  [math]A_{i}=J^{T}_iA_{i-1}J_i ,[/math]

где [math]J_i[/math] - ортогональная матрица, называемая вращением Якоби. Таким образом,

  [math]A_{m}=J^{T}_{m-1}A_{m-1}J_{m-1}[/math]
       [math]=J^{T}_{m-1}J^{T}_{m-2}A_{m-2}J_{m-2}J_{m-1}[/math]
       [math]=J^{T}_{m-1}...J^{T}_{0}A_{0}J_{0}...J_{m-1}[/math]
       [math]=J^{T}AJ[/math]

При подходящем выборе вращений [math]J_i[/math] матрица [math]A_m[/math] для больших [math]m[/math] будет близкой к диагональной матрице [math]\Lambda[/math]. Таким образом, [math]\Lambda \approx J^{T}AJ[/math] или [math]J\Lambda J^{T} \approx A[/math], следовательно, столбцы матрицы [math]J[/math] являются собственными значениями матрицы [math]A[/math]. Матрица [math]J_i[/math] выбирается таким образом, чтобы сделать нулями пару внедиагональных элементов матрицы [math]A_{i+1}=J^{T}_{i}A_{i}J_{i}[/math]:

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

Угол вращения [math]\theta[/math] выбирается так, чтобы элементы [math](j,k)[/math] и [math](k,j)[/math] в [math]A_{i+1}[/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]

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

1)Находится максимальный наддиагональный элемент [math]a_{jk}[/math] матрицы [math]A[/math].

2)Вычисляются параметры вращения Якоби:

  [math]\tau= \frac{(a_{jj}-a_{kk})}{2a_{jk}}[/math]
  [math]t=sign(\tau)(|\tau|+\sqrt{1+\tau^2})[/math]
  [math]c= \frac{1}{1+\sqrt{1+\tau^2}}[/math]
  [math]s=c*t[/math]

3)К текущей матрице А применяются вращения с вычисленными параметрами:

  [math]A=R^T(j,k,\theta)*A*R(j,k,\theta)[/math]

Стоит отметить, что вычисленные на предыдущем шаге [math]c[/math] и [math]s[/math] есть [math]\cos {\theta}[/math] и [math]\sin {\theta}[/math] соответственно.

4)Если матрица А не достаточно близка к диагональной, происходит переход к шагу 1.

Стоит отметить, что за [math]\frac{n*(n-1)}{2}[/math] итераций алгоритма гарантируется получение результата, т.к. сумма квадратов наддиагональных элементов матрицы А на каждой итерации уменьшается на [math]a_{jk}^2[/math].

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

Сначала оценим сложность одной итерации:

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

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

3)[math]4*n[/math] операций умножения и [math]2*n[/math] операций сложения.

4)операция сравнения

В наиболее неудачном случае количество итераций равняется числу наддиагональных элементов, т.е. [math]\frac{n*(n-1)}{2}[/math]. Таким образом, суммарная сложность алгоритма составляет [math]O(N^4)[/math] сравнений, и [math]O(N^3)[/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]\frac{n (n + 1)}{2}[/math] по строкам своей нижней части, хранение в виде так называемого ступенчатого массива (такой вариант двумерного массива, в котором каждая строка имеет свой размер).

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