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

Эта работа прошла предварительную проверку Дата последней правки страницы: 23.11.2016 Данная работа соответствует формальным критериям. Проверено ASA.

Метод Якоби вычисления собственных значений симметричной матрицы ЗП
Последовательный алгоритм
Последовательная сложность	[math]O(n^3)[/math]
Объём входных данных	[math]n^2[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(n^2)[/math]
Ширина ярусно-параллельной формы	[math]O(n)[/math]

Авторы: Зиновьев В.С. (написание кода, информационный граф, общее и мат описание алгоритма), Пронин А.М. (написание кода, информационный граф, вычислительное ядро алгоритма, макроструктура алгоритма). В остальные разделы статьи вклад равноценный.

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

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

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

Метод Якоби приводит симметричную действительную матрицу к диагональной форме с помощью некоторой (в принципе, бесконечной) последовательности плоских вращений в разных плоскостях исходного n-мерного векторного пространства. Метод Якоби состоит из некоторой последовательности шагов, на каждом из которых обнуляется один из внедиагональных элементов (обычно наибольший по модулю), в результате чего старая матрица [math]A_{i}[/math] преобразуется в некоторую новую матрицу [math]A_{i+1}[/math]. Делается это с помощью преобразования подобия, использующего матрицу плоского вращения [math]J_{i}[/math] (подробное описание матрицы плоского вращения приведено ниже в разделе "Математическое описание алгоритма"), которая изменяет лишь две координаты [math]n[/math]-мерных векторов, а именно, координаты с номерами [math]j[/math] и [math]k[/math]. При таком шаге матрица [math]A_{i+1}[/math] отличается от матрицы [math]A_{i}[/math] лишь элементами двух строк с номерами [math]j[/math], [math]k[/math] и двух столбцов с теми же номерами ([math]j[/math] и [math]k[/math]). Все же другие элементы матрицы [math]A[/math] остаются неизменными. Число шагов, за которое матрица [math]A[/math] приводится таким способом к диагональной форме, в общем случае является бесконечным. (это свойство метода Якоби более подробно рассмотрено в разделе "Свойства алгоритма") В результате метод Якоби представляет собой некоторый бесконечный итерационный процесс, на каждом шаге которого преобразующаяся матрица становится всё более и более близкой к диагональной форме (в смысле уменьшения суммы квадратов всех своих внедиагональных элементов). На практике обычно выбирается заранее заданная погрешность [math]\varepsilon[/math] и процесс приведения исходной матрицы к диагональному виду прекращается, когда наибольший по модулю внедиагональный элемент матрицы становится меньше [math]\varepsilon[/math]. Таким образом, применяя на практике метод Якоби, мы приводим исходную матрицу к диагональной форме с точностью до некоторой вычислительной погрешности [math]\varepsilon[/math] в результате конечного числа шагов. Повышение этой точности требует увеличения числа шагов метода Якоби. В конечном итоге у полученной в результате этого процесса матрицы на диагонали будут находиться все собственные значения (с определённой погрешностью) исходной симметрической матрицы.

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[/math]                         [math]k[/math]

[math] J_i = R(j,k,\theta) \equiv \begin{matrix} \\ \\ \\ j \\ \\ k \\ \\ \\ \\ \end{matrix} [/math] [math] \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 Вычислительное ядро алгоритма

Вычислительным ядром отдельной итерации является вычисление элементов матрицы [math]A[/math] при применении поворота Якоби. Конкретно, изменяются [math]2[/math] строки и [math]2[/math] столбца, что означает изменение [math]4n[/math] элементов. Отметим, что вычисление нового значения одного элемента требует [math]2[/math] операции умножения и одну операцию сложения, что позволяет утверждать, что основной вклад в вычислительное ядро вносят умножения. Однако в целом, следует отметить, что число операций в алгоритме во многом определяется общим числом итераций, а оно квадратично, в то время как отдельная итерация имеет линейную сложность. Более того, так как алгоритм является итерационным и останавливается при достижении некоторой точности, число итераций может изменяться в зависимости от желаемой точности.

1.4 Макроструктура алгоритма

На макроуровне имеет смысл выделять следующие блоки алгоритма:

1) Блок вычисления параметров вращения на текущей итерации (выбор обнуляемого элемента и расчет угла вращения).

2)Блок, реализующий непосредственно вращение Якоби.

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[/math] применяются вращения с вычисленными параметрами:

  [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)Если матрица [math]A[/math] не достаточно близка к диагональной, происходит переход к шагу 1.

Вообще говоря, несмотря на то, что поворот, обеспечивающий обнуление наибольшего внедиагонального элемента, обеспечивает схождение метода за наименьшее число итераций, поиск данного наибольшего элемента добавляет порядка [math]n^2[/math] операций сравнения, а возможно и присваивания на каждой итерации. Поэтому с вычислительной точки зрения выгоднее не выбирать наибольший элемент, а по очереди перебирать все наддиагональные элементы, большие заданного порога.

В последующих разделах рассматривается именно такой вариант алгоритма.

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

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

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

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

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

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

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

Если выбирать не наибольший элемент матрицы, а все по очереди, порядок числа итераций сохраняется, но исчезают операции сравнения. Таким образом, итоговая сложность алгоритма будет составлять [math]O(N^3)[/math].

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

Аналитически информационный граф данного алгоритма можно описать следующим образом:

1)Последовательно идущие макроярусы вычисления параметров текущей итерации и вращения Якоби с этими параметрами. Число данных ярусов варьируется в зависимости от входных данных, но порядок составляет [math]O(n^2)[/math].

2)Вычисление параметров текущей итерации содержит конечное число операций умножения, сложения и вычисления квадратного корня.

3)Макроярус вращения Якоби содержит вычисление [math]2n+2[/math] элементов матрицы и [math]2n-2[/math] операций присваивания. На данном ярусе ширина графа достигает максимального значения. Вначале вычисляются значения в строках [math]j[/math] и [math]k[/math] (это соответствует умножению матрицы [math]R^T[/math] на матрицу [math]A[/math]). Далее, в силу того, что при умножении результата на матрицу [math]A[/math] изменяются только столбцы [math]j[/math] и [math]k[/math] и итоговая матрица симметрична, можно не пересчитывать значения в столбцах, а получить их простым присваиванием по формуле: [math]a_{pk}=a_{kp}[/math], [math]a_{pj}=a_{jp}[/math], [math]p \neq k,j[/math]. Для элементов [math]a_{jj}[/math] и [math]a_{kk}[/math] необходимо произвести вычисления, аналогичные вычислениям для элементов строк.

Визуальное представление графа алгоритма представлено на схемах ниже:

Рисунок 1. Граф алгоритма на макроуровне. Rarams - вычисление параметров, Rotate - получение [math]A_{i+1}[/math] из [math]A_{i}[/math].

Рисунок 2. Структура блока Params.

Рисунок 3. Структура блока Rotate.

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

Как видно из представления алгоритма в ярусно-параллельной форме ( п.1.7), максимальная ширина ярусов составляет [math]2n[/math] вершин. Для выполнения данного алгоритма в параллельном варианте следует выполнить [math]O(\frac{n(n+1)}{2})[/math] итераций. Рассмотрим отдельную итерацию с точки зрения ярусно-параллельной формы алгоритма.

1)Выполнение 8 ярусов с конечным числом операций сложения, умножения, взятия квадратного корня и нахождения знака.

2)Выполнение 2 ярусов вычисления элементов в изменяемых строках и столбцах.

На данных двух ярусах осуществляется [math]4n[/math] умножений и [math]2n[/math] сложений. Данный вид параллелизма относится к координатному. Таким образом, при наличии бесконечного числа процессоров, итоговая сложность алгоритма составляет [math]O(\frac{n(n+1)}{2})[/math].

Кроме вышеизложенного, отметим следующую особенность алгоритма. Пусть вращение [math]J_1[/math] изменяет строки и столбцы с номерами [math]j_1[/math] и [math]k_1[/math], а вращение [math]J_2[/math] изменяет строки и столбцы с номерами [math]j_2[/math] и [math]k_2[/math], и при этом выполнено:

[math]j_1 \neq j_2, k_2[/math]

[math]k_1 \neq j_2, k_2[/math]

Тогда вращения [math]J_1[/math] и [math]J_2[/math] перестановочны, и их можно выполнять параллельно. Учитывая перебор обнуляемых элементов поочередно по строкам, данная ситуация достаточно редка. Она не изменяет итоговую сложность алгоритма, но уменьшает высоту ЯПФ на [math]O(n)[/math], а ширину ЯПФ напротив, увеличивает до [math]4n[/math]. Таким образом, итоговый алгоритм имеет ширину ЯПФ [math]O(n)[/math] и высоту [math]O(n^2)[/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].

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

Параллельная сложность выполнения алгоритма составляет [math]O(\frac{n(n+1)}{2}).[/math] Последовательная сложность выполнения алгоритма составляет [math]O(\frac{n^2(n+1)}{2}).[/math] Таким образом, соотношение сложностей является линейным. Более того, и вычислительная мощность данного алгоритма линейна.

В алгоритме используется примерно в два раза больше операций умножения, чем сложения. Отдельно стоит рассмотреть вопрос хранения матрицы [math]A[/math] в процессе работы алгоритма. Хранение лишь нижней или верхней части матрицы в виде одномерного массива возможно; более того, хранение в таком виде обеспечивает экономию памяти в [math]2[/math] раза и защищает от ассиметричности матрицы в результате накопления математических ошибок. Однако, вычисление индексов [math]i[/math] и [math]j[/math] элементов [math]a_{ij}[/math], соответствующих данному номеру [math]k[/math] вообще говоря требует отдельного цикла, что создает дополнительную нагрузку на процессор.

Кроме того, следует отметить следующие особенности алгоритма:

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

2)Если при выборе в качестве обнуляемого элемента [math]a_{ij}[/math] элементы [math]a_{ii} = a_{jj}[/math], то, как видно из формул в пунктах выше, матрицы поворота вырождаются в единичные и матрица на данной итерации не изменяется. Это еще один довод в пользу отказа от выбора наибольшего элемента, т.к. в данном случае произойдет зацикливание алгоритма.

2 Программная реализация алгоритма

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

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

Исследование проводилось на 1U сервере со следующими характеристиками:

• 2 x Intel® Xeon® Processor E5-2697 v2;
• 128GB RAM;

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

• число процессоров [1 : 48] с шагом 1
• размерность матрицы [400:1600] с шагом 200

Использовался компилятор g++ версии 5.4.0 со следующими флагами: -O2, -fopenmp, -std=c++11. Флаг -O2 является опциональным, и означает уровень оптимизации кода при компиляции. Все прочие флаги строго обязательны при компиляции, т.к. в коде используются функции из стандарта С++11.

Рисунок 2. Изменение производительности в зависимости от числа процессоров и размерности матрицы.

Чтобы продемонстрировать важное свойство алгоритма, приведем так же график производительности в зависимости от числа процессов при размере матрицы 1600*1600

Рисунок 2. Изменение производительности в зависимости от числа процессоров при размере матрицы 1600х1600.

За [math]1[/math] принята производительность при одном процессе. Как видно, производительность резко увеличивается при добавлении числа потоков в начале и снижается при дальнейшем добавлении. Это связано с тем, что в одном потоке совершается малое количество арифметических операций. После выполнения операций совершается запись в массив полученных значений, в результате чего возникает конфликт потоков. Поэтому наиболее оптимальное число потоков [math]8-16[/math], что говорит о слабой масштабируемости данного алгоритма.

Параллельная реализация на языке C++

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

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

2)В цифровой библиотеке ieee.org приводится аннотация статьи 2012 года, в которой авторы описывают реализацию данного алгоритма на CUDA.

3 Литература

<references \>