Участник:Зиновьев Владимир/Метод Якоби вычисления собственных значений симметричной матрицы ЗП: различия между версиями
Строка 89: | Строка 89: | ||
\frac{a_{jj} - a{kk}}{2a_{jk}} = \frac{c^2-s^2}{2sc} = ctg(2\theta) = r. | \frac{a_{jj} - a{kk}}{2a_{jk}} = \frac{c^2-s^2}{2sc} = ctg(2\theta) = r. | ||
</math> | </math> | ||
− | Пусть <math>tg(\theta) = t = \frac{s}{c}</math> . Заметим, что выполняется квадратное уравнение <math>t^2 + 2rt - 1 = 0</math> . Решая уравнение, получаем <math>t = \frac{sign | + | Пусть <math>tg(\theta) = t = \frac{s}{c}</math> . Заметим, что выполняется квадратное уравнение <math>t^2 + 2rt - 1 = 0</math> . Решая уравнение, получаем <math>t = \frac{sign(r)}{|r| + \sqrt{1+r^2}}</math> . Следовательно <math>c = \frac{1}{\sqrt{1+t^2}}</math>, <math>s = ct</math> |
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === |
Версия 00:24, 12 октября 2016
Разложение Холецкого | |
Последовательный алгоритм | |
Последовательная сложность | [math][/math] |
Объём входных данных | [math][/math] |
Объём выходных данных | [math][/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math][/math] |
Ширина ярусно-параллельной формы | [math][/math] |
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 3 Литература
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]
Пусть [math]c = \cos(\theta)[/math] и [math]s = sin(\theta)[/math], [math]a^{i}_{lk} = a_{lk}[/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] \frac{a_{jj} - a{kk}}{2a_{jk}} = \frac{c^2-s^2}{2sc} = ctg(2\theta) = r. [/math]
Пусть [math]tg(\theta) = t = \frac{s}{c}[/math] . Заметим, что выполняется квадратное уравнение [math]t^2 + 2rt - 1 = 0[/math] . Решая уравнение, получаем [math]t = \frac{sign(r)}{|r| + \sqrt{1+r^2}}[/math] . Следовательно [math]c = \frac{1}{\sqrt{1+t^2}}[/math], [math]s = ct[/math]
1.3 Вычислительное ядро алгоритма
1.4 Макроструктура алгоритма
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.7 Информационный граф
1.8 Ресурс параллелизма алгоритма
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].