Участник:Danyanya/Алгоритм Ланцоша для точной арифметики (без переортогонализации): различия между версиями
Danyanya (обсуждение | вклад) |
ASA (обсуждение | вклад) |
||
(не показаны 53 промежуточные версии 4 участников) | |||
Строка 1: | Строка 1: | ||
+ | {{Assignment|ASA|Konshin}} | ||
+ | |||
{{algorithm | {{algorithm | ||
| name = Алгоритм Ланцоша для точной арифметики (без переортогонализации) | | name = Алгоритм Ланцоша для точной арифметики (без переортогонализации) | ||
| serial_complexity = <math>O(kn^2)</math> | | serial_complexity = <math>O(kn^2)</math> | ||
− | | pf_height = <math>O(k)</math> | + | | pf_height = <math>O(k \log(n))</math> |
| pf_width = <math>O(n^2)</math> | | pf_width = <math>O(n^2)</math> | ||
| input_data = <math>\frac{n(n + 1)}{2}</math> | | input_data = <math>\frac{n(n + 1)}{2}</math> | ||
Строка 8: | Строка 10: | ||
}} | }} | ||
− | Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]] | + | Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]] (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 2.4, 2.7), [[Участник:m.grigoriev|М.А.Григорьев]] (1.2, 1.7, 1.8, 1.9, 2.4, 2.7, 3) |
+ | |||
== Свойства и структура алгоритма == | == Свойства и структура алгоритма == | ||
Строка 14: | Строка 17: | ||
=== Общее описание алгоритма === | === Общее описание алгоритма === | ||
− | ''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году | + | ''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>. |
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>. | Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>. | ||
− | В данной статье рассмотрен упрощенный вариант алгоритма Ланцоша, | + | В данной статье рассмотрен упрощенный вариант алгоритма Ланцоша, подразумевающий отсутствие влияния ошибок округления на вычислительный процесс. |
− | Данный алгоритм является неустойчивым, вследствие чего на практике применяется модифицированный алгоритм Ланцоша с полной переортогонализацией. | + | Данный алгоритм является неустойчивым, вследствие чего на практике применяется модифицированный [[Алгоритм Ланцоша с полной переортогонализацией | алгоритм Ланцоша с полной переортогонализацией]] предложенный в 1970 Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref>. |
=== Математическое описание алгоритма === | === Математическое описание алгоритма === | ||
Строка 26: | Строка 29: | ||
{{Шаблон: ASymmetric}} | {{Шаблон: ASymmetric}} | ||
− | Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца. Иными словами, из оргонормированных векторов Ланцоша | + | Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца[https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Ritz_method]. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>. |
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>: | В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>: | ||
Строка 43: | Строка 46: | ||
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === | ||
− | + | Вычислительное ядро на каждой итерации алгоритма Ланцоша состоит в вычислении очередного промежуточного вектора <math>z</math> произведением исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации: | |
:<math>z = Aq_i</math> | :<math>z = Aq_i</math> | ||
− | + | ||
+ | При больших <math>k</math>, сопоставимых с <math>n</math>, к вычислительному ядру можно отнести вычисление собственных значений и собственных векторов трёхдиагональной матрицы. Однако на практике <math>k</math> много меньше <math>n</math>. | ||
+ | |||
=== Макроструктура алгоритма === | === Макроструктура алгоритма === | ||
− | + | Исходя из предложенной последовательной реализации метода, макрооперациями в алгоритме являются: | |
− | * Процедура | + | * Процедура итеративного построения трехдиагональной симметричной матрицы, включающая: |
− | ** умножение матрицы на вектор; | + | ** умножение матрицы на вектор (состоит из умножения вектора на число и сложения векторов); |
** скалярное произведение векторов; | ** скалярное произведение векторов; | ||
− | ** | + | ** линейная комбинация векторов (сложение/умножение на вещественные числа); |
− | * Вычисление | + | ** вычисление нормы (скалярное произведение векторов и вычисление квадратного корня); |
+ | * Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы. | ||
+ | |||
+ | При этом, первая макрооперация (построение трехдиагональной матрицы) выполняется строго последовательно. Единственное, что можно распараллелить - это умножение матрицы на вектор. Однако это операция достаточно легковесная, если оперировать, например, разреженной матрицей. | ||
+ | |||
+ | Вычисление собственных значений полученной матрицы на практике выполняется с помощью QR-алгоритма и выходит за рамки описанного алгоритма. | ||
=== Схема реализации последовательного алгоритма === | === Схема реализации последовательного алгоритма === | ||
Строка 62: | Строка 72: | ||
'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. | '''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. | ||
− | Алгоритм на псевдокоде: | + | Алгоритм <ref>Деммель Д. Вычислительная линейная алгебра</ref> на псевдокоде: |
<math> | <math> | ||
Строка 92: | Строка 102: | ||
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций. | * Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций. | ||
− | Таким образом, в худшем случае, алгоритм Ланцоша имеет сложность <math> O( | + | Суммарное число операций в алгоритме без учета вычисления собственных значений в трехдиагональной памяти: |
+ | * <math> kn^2 + n(4k + 1) </math> операций умножения; | ||
+ | * <math> k(n^2+3n-2) + n-1</math> операций сложения/вычитания; | ||
+ | * <math> n(k + 1) </math> операций деления; | ||
+ | * <math> k + 1 </math> операций вычислений квадратного корня. | ||
+ | |||
+ | Таким образом, в худшем случае, алгоритм Ланцоша имеет сложность <math> O(kn^2) </math> и односится к квадратичным алгоритмам по <math>n</math> при малых <math>k</math>. | ||
=== Информационный граф === | === Информационный граф === | ||
Информационный граф алгоритма можно разбить на две части: | Информационный граф алгоритма можно разбить на две части: | ||
− | * | + | * Граф алгоритма с отображением входных и выходных данных. |
− | [[Файл: | + | * Граф линейного оператора. |
+ | [[Файл:mgrigoriev2.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма.<br/> | ||
+ | <math>\|x\|</math> — вычисление нормы, <br/> | ||
+ | <math>/</math> — операция деления, <br/> | ||
+ | <math>\mathsf{F}</math> — вычисление линейной комбинации векторов, <br/> | ||
+ | <math>*</math> — операция умножения,<br/> | ||
+ | <math>\mathsf{Ax}</math> — линейный оператор (Рисунок 2),<br/> | ||
+ | ▽ - условие выхода <math>\mathsf{(\beta_i = 0)}</math>]] | ||
− | * | + | [[Файл:mgrigoriev3.png|thumb|center|600px|Рисунок 2. Информационный граф линейного оператора <math>\mathsf{Ax}</math>. <br/> <math>+</math>, <math>*</math> - операции сложения и умножения.]] |
− | [[Файл: | + | [[Файл:mgrigoriev1.png|600px|thumb|center|Рисунок 3. Граф алгоритма с указанием входных и выходных данных.<br/> |
+ | <math>\|x\|</math> — вычисление нормы, <br/> | ||
+ | <math>/</math> — операция деления, <br/> | ||
+ | <math>\mathsf{F}</math> — вычисление линейной комбинации векторов, <br/> | ||
+ | <math>*</math> — операция умножения,<br/> | ||
+ | <math>\mathsf{Ax}</math> — линейный оператор (Рисунок 2),<br/> | ||
+ | ▽ - условие выхода <math>\mathsf{(\beta_i = 0)}</math>,<br/> | ||
+ | <math>\mathsf{A, b, q_0 = 0, \beta_0 = 0}</math> — входные данные, <br/> | ||
+ | <math>\mathsf{\alpha_i, \beta_i}</math> — выходной результат.]] | ||
+ | |||
+ | Из информационного графа видно, что паралеллизация алгоритма возможна при выполнении линенйного оператора <math>Ax</math> на каждой итерации алгоритма. Затем необходима синхронизация между процессами. | ||
=== Ресурс параллелизма алгоритма === | === Ресурс параллелизма алгоритма === | ||
Строка 108: | Строка 141: | ||
Важно заметить, что итерации алгоритма выполняются строго последовательно, а распараллеливание возможно только внутри итераций. | Важно заметить, что итерации алгоритма выполняются строго последовательно, а распараллеливание возможно только внутри итераций. | ||
− | * Умножение <math>A^{ | + | * Умножение <math>A^{n * n}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложенийж |
+ | ** При этом сложение элементов вектора длины <math>n</math> можно выполнить за <math>\log(n)</math> [https://algowiki-project.org/ru/%D0%9D%D0%B0%D1%85%D0%BE%D0%B6%D0%B4%D0%B5%D0%BD%D0%B8%D0%B5_%D1%81%D1%83%D0%BC%D0%BC%D1%8B_%D1%8D%D0%BB%D0%B5%D0%BC%D0%B5%D0%BD%D1%82%D0%BE%D0%B2_%D0%BC%D0%B0%D1%81%D1%81%D0%B8%D0%B2%D0%B0_%D1%81%D0%B4%D0%B2%D0%B0%D0%B8%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5%D0%BC]; | ||
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус): | * Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус): | ||
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма. | * Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма. | ||
− | Исходя из вышеизложенного, алгоритм Ланцоша обладает | + | Исходя из вышеизложенного, алгоритм Ланцоша обладает <math>O(k\log n)</math>-ой сложностью по высоте ЯПФ и <math>O(n^2)</math>-ой по ширине ЯПФ. |
=== Входные и выходные данные алгоритма === | === Входные и выходные данные алгоритма === | ||
Строка 118: | Строка 152: | ||
'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>. | '''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>. | ||
− | '''Объём входных данных:''' <math>n | + | '''Объём входных данных:''' <math>n(n + 1) + 1 </math>. |
'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>. | '''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>. | ||
− | '''Объём выходных данных:''' <math>k | + | '''Объём выходных данных:''' <math>k(n + 1)</math>. |
=== Свойства алгоритма === | === Свойства алгоритма === | ||
+ | * Сложность последовательного алгоритма <math>O(kn^2)</math>; | ||
+ | * По классификации по высоте ЯПФ Алгоритм Ланцоша является алгоритмом со сложностью <math>O(k\log n)</math>, при классификации по ширине ЯПФ его сложность <math>O(n^2)</math>; | ||
+ | ** Таким образом, отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>; | ||
+ | * Вычислительная мощность алгоритма Ланцоша без переортогонализации из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. При <math>k</math> много меньше <math>n</math> вычислительная мощность ≈ <math> 2k</math>; | ||
* Алгоритм Ланцоша без переортогонализации не является детерминированным из-за того, что возможно выполнение меньшего числа итераций алгоритма, из-за того, что все собственные значения уже вычислены; | * Алгоритм Ланцоша без переортогонализации не является детерминированным из-за того, что возможно выполнение меньшего числа итераций алгоритма, из-за того, что все собственные значения уже вычислены; | ||
− | |||
* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения); | * Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения); | ||
− | |||
* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются; | * Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются; | ||
− | |||
* Нестабильность алгоритма (эффект ''ложной сходимости'') присуще плавающей арифметике, из-за ошибок округления в которой на очередном этапе может быть построен линейно зависимый от исходных новый вектор, что повлечет невозможность дальнейшего приближения собственных значений. | * Нестабильность алгоритма (эффект ''ложной сходимости'') присуще плавающей арифметике, из-за ошибок округления в которой на очередном этапе может быть построен линейно зависимый от исходных новый вектор, что повлечет невозможность дальнейшего приближения собственных значений. | ||
Строка 144: | Строка 179: | ||
=== Масштабируемость алгоритма и его реализации === | === Масштабируемость алгоритма и его реализации === | ||
− | + | Для исследования масшабируемости алгоритма была написана реализация[https://github.com/danyanya/ss16-task1] на языке C++ с использованием MPI. Реализация была протестирована на суперкомпьютере Ломоносов[http://parallel.ru/cluster/]. | |
+ | |||
+ | Были исследованы: | ||
+ | * время выполнения программы в зависимости от размера входных данных (матрицы); | ||
+ | * время выполнения в зависимости от размера параллельных нод. | ||
+ | |||
+ | Дополнительно было измерено время работы исходя из оптимизационных флагов компилятора GCC [https://gcc.gnu.org]. | ||
+ | |||
+ | Параметры запуска алгоритма: | ||
+ | * размер матрицы от 20000 до 175000 с шагом 2500; | ||
+ | * количество процессоров от 8 до 128 с шагом 8. | ||
+ | |||
+ | Программа запускалась на суперкомьютере "Ломоносов" со следующими характеристиками: | ||
+ | * Компилятор GCC 5.2.1; | ||
+ | * Версия MPI 1.8.4; | ||
+ | * Сборка проводилась командой: <pre>mpic++ -std=c++0x -O2 <CPP_FILE> -lm -static-libstdc++</pre> | ||
+ | |||
+ | Полученная эффективность реализации варьируется в пределах от 2% (на маленьких входных данных) до 45% (на больших входных данных и максимальном числе нодов). | ||
+ | |||
+ | На рисунке 4 представлена эффективность программы при отсутствии флага оптимизации O2. | ||
+ | |||
+ | [[Файл:Grigorev_lanczos_eff_no_02.png|thumb|center|600px|Рисунок 4. Эффективность параллельной реализации алгоритма Ланцоща на MPI.]] | ||
+ | |||
+ | На рисунке 5 представлена эффективность алгоритма с оптимизацией компилятора GCC (-O2). | ||
+ | |||
+ | [[Файл:Grigorev_lanczos_eff_o2.png|thumb|center|600px|Рисунок 5. Эффективность параллельной реализации алгоритма Ланцоща на MPI с использованием оптимизации компилятора -O2.]] | ||
+ | |||
+ | |||
+ | Как видно из графиков выше, средняя эффективность минимальна при малых входных данных (на размере матрицы в 10000). | ||
+ | В среднем данный показатель держится между 9% и 14%. | ||
+ | |||
+ | При этом, оптимизация компилятора MPIC++ позволяет получить выигрыш в эффективности более чем в 3 раза (45,5% против 14,7%). Однако, как можно заметить из графиков, при оптимизации алгоритм ведет себя более нестабильно, скорее всего из-за значительной траты времени на работу с памятью (более планые участки свидетельствует о том, что необходимые данные нашлись в кеше). | ||
=== Динамические характеристики и эффективность реализации алгоритма === | === Динамические характеристики и эффективность реализации алгоритма === | ||
Строка 154: | Строка 220: | ||
В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них: | В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них: | ||
− | 1. The IETL Project | + | 1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref> |
+ | |||
+ | 2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref> | ||
− | + | 3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref> | |
− | + | 4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref> | |
− | |||
− | ... | + | Также существуют официальные реализации на других языках (например R). Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL. |
+ | |||
+ | Для проверки масшабируемости, алгоритм был реализован на языке C++ с применением MPI [https://github.com/danyanya/ss16-task1]. | ||
== Литература == | == Литература == | ||
− | |||
− | |||
− |
Текущая версия на 16:59, 19 декабря 2016
Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Konshin и ASA. |
Алгоритм Ланцоша для точной арифметики (без переортогонализации) | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(kn^2)[/math] |
Объём входных данных | [math]\frac{n(n + 1)}{2}[/math] |
Объём выходных данных | [math]k(n + 1)[/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math]O(k \log(n))[/math] |
Ширина ярусно-параллельной формы | [math]O(n^2)[/math] |
Основные авторы описания: Д.Р.Слюсарь (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 2.4, 2.7), М.А.Григорьев (1.2, 1.7, 1.8, 1.9, 2.4, 2.7, 3)
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Алгоритм Ланцоша поиска собственных значений был опубликован Корнелием Ланцошем в 1950 году [1]. Этот итерационный алгоритм применим только к эрмитовым матрицам [math]A[/math]. Метод позволяет за [math]k[/math] итераций вычислять [math]k[/math]-ое приближение собственных значений и собственных векторов исходной матрицы [math]A[/math].
В данной статье рассмотрен упрощенный вариант алгоритма Ланцоша, подразумевающий отсутствие влияния ошибок округления на вычислительный процесс.
Данный алгоритм является неустойчивым, вследствие чего на практике применяется модифицированный алгоритм Ланцоша с полной переортогонализацией предложенный в 1970 Ojalvo и Newman[2].
1.2 Математическое описание алгоритма
На вход алгоритма подается эрмитова матрица [math]A = A^\dagger[/math] (в вещественном случае матрица симметрична) ,
- [math] A = \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1\ n-1} & a_{1\ n} \\ a_{12} & a_{22} & a_{23} & \cdots & a_{2\ n-1} & a_{2\ n} \\ a_{13} & a_{23} & a_{33} & \cdots & a_{3\ n-1} & a_{3\ n} \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ a_{1\ n-1} & \cdots & \cdots & a_{n-2\ n-1} & a_{n-1\ n-1} & a_{n-1\ n} \\ a_{1\ n} & \cdots & \cdots & a_{n-2\ n} & a_{n-1\ n} & a_{n\ n} \\ \end{pmatrix} [/math]
Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца[1]. Иными словами, из оргонормированных векторов Ланцоша [3] на каждой итерации строится матрица [math]Q_k = [q_1, q_2, \dots, q_k][/math] размерности [math]n \times k[/math]. В качестве приближенных собственных значений матрицы [math]A[/math] берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы [math]T_k = Q^T_k A Q[/math]:
[math] T_k = \begin{pmatrix} \alpha_1 & \beta_1 & 0 & \dots & 0 \\ \beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\ 0 & \beta_2 & \ddots & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\ 0 & \dots & \dots & \beta_{k-1} & \alpha_k \end{pmatrix} [/math]
На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы [math]T_k[/math], с помощью которых и будет найдены искомые собственные векторы исходной матрицы [math]A[/math].
1.3 Вычислительное ядро алгоритма
Вычислительное ядро на каждой итерации алгоритма Ланцоша состоит в вычислении очередного промежуточного вектора [math]z[/math] произведением исходной матрицы [math]A[/math] на вектор [math]q_i[/math] с предыдущей итерации:
- [math]z = Aq_i[/math]
При больших [math]k[/math], сопоставимых с [math]n[/math], к вычислительному ядру можно отнести вычисление собственных значений и собственных векторов трёхдиагональной матрицы. Однако на практике [math]k[/math] много меньше [math]n[/math].
1.4 Макроструктура алгоритма
Исходя из предложенной последовательной реализации метода, макрооперациями в алгоритме являются:
- Процедура итеративного построения трехдиагональной симметричной матрицы, включающая:
- умножение матрицы на вектор (состоит из умножения вектора на число и сложения векторов);
- скалярное произведение векторов;
- линейная комбинация векторов (сложение/умножение на вещественные числа);
- вычисление нормы (скалярное произведение векторов и вычисление квадратного корня);
- Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.
При этом, первая макрооперация (построение трехдиагональной матрицы) выполняется строго последовательно. Единственное, что можно распараллелить - это умножение матрицы на вектор. Однако это операция достаточно легковесная, если оперировать, например, разреженной матрицей.
Вычисление собственных значений полученной матрицы на практике выполняется с помощью QR-алгоритма и выходит за рамки описанного алгоритма.
1.5 Схема реализации последовательного алгоритма
Исходные данные: симметричная матрица [math]A[/math], случайный вектор [math]b[/math].
Вычисляемые данные: собственные вектора матрицы [math]T_k[/math] являющиеся столбцами матрицы [math]Q_k V[/math], и матрица собственных значений [math]\Lambda[/math], где [math]V, \Lambda[/math] из спектрального разложения [math]T_k = V\Lambda V^T[/math].
Алгоритм [4] на псевдокоде:
[math] \begin{align} q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\ for \; & i = 1 \; to \; k \\ & z = Aq_i\\ & \alpha_i = q^T_i z\\ & z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\ & \beta_i = \|z\|_2\\ & If \; \beta_i == 0 \; then \\ & \; \; \; \; exit\\ & else \\ & \; \; \; \; q_{i+1} = z / \beta_i \\ end \; & for \end{align} [/math]
После этого вычисляются собственные значения и собственные вектора симметричной трехдиагональной матрицы [math]T_k[/math] наиболее удобным образом.
1.6 Последовательная сложность алгоритма
Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
- Умножение квадратной матрицы [math]n * n[/math] на вектор длины [math]n[/math]. Требует [math]n * n[/math] умножений и сложений;
- Скалярное произведение векторов длины [math]n[/math]. Требует [math]n[/math] умножений и сложений;
- Сложение векторов длины [math]n[/math]. Требует [math]n[/math] сложений;
- Умножение вектора длины [math]n[/math] на скаляр. Требует [math]n[/math] умножений;
- Нахождение квадратичной нормы вектора длины [math]n[/math]. Требует [math]n[/math] умножений и сложений + извлечения квадратного корня;
- Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера [math] k \times k [/math]. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует [math]~O(k^{2.3})[/math] операций.
Суммарное число операций в алгоритме без учета вычисления собственных значений в трехдиагональной памяти:
- [math] kn^2 + n(4k + 1) [/math] операций умножения;
- [math] k(n^2+3n-2) + n-1[/math] операций сложения/вычитания;
- [math] n(k + 1) [/math] операций деления;
- [math] k + 1 [/math] операций вычислений квадратного корня.
Таким образом, в худшем случае, алгоритм Ланцоша имеет сложность [math] O(kn^2) [/math] и односится к квадратичным алгоритмам по [math]n[/math] при малых [math]k[/math].
1.7 Информационный граф
Информационный граф алгоритма можно разбить на две части:
- Граф алгоритма с отображением входных и выходных данных.
- Граф линейного оператора.
Из информационного графа видно, что паралеллизация алгоритма возможна при выполнении линенйного оператора [math]Ax[/math] на каждой итерации алгоритма. Затем необходима синхронизация между процессами.
1.8 Ресурс параллелизма алгоритма
Важно заметить, что итерации алгоритма выполняются строго последовательно, а распараллеливание возможно только внутри итераций.
- Умножение [math]A^{n * n}[/math] на вектор длины [math]n[/math] требует [math]n[/math] ярусов умножений и сложенийж
- При этом сложение элементов вектора длины [math]n[/math] можно выполнить за [math]\log(n)[/math] [2];
- Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
- Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.
Исходя из вышеизложенного, алгоритм Ланцоша обладает [math]O(k\log n)[/math]-ой сложностью по высоте ЯПФ и [math]O(n^2)[/math]-ой по ширине ЯПФ.
1.9 Входные и выходные данные алгоритма
Входные данные: симметричная вещественная матрица [math]A[/math], случайный вектор [math]b[/math], число итераций [math]k[/math].
Объём входных данных: [math]n(n + 1) + 1 [/math].
Выходные данные: вектор собственных значений [math]\Lambda[/math], матрица собственных векторов [math]E[/math].
Объём выходных данных: [math]k(n + 1)[/math].
1.10 Свойства алгоритма
- Сложность последовательного алгоритма [math]O(kn^2)[/math];
- По классификации по высоте ЯПФ Алгоритм Ланцоша является алгоритмом со сложностью [math]O(k\log n)[/math], при классификации по ширине ЯПФ его сложность [math]O(n^2)[/math];
- Таким образом, отношение последовательной сложности к параллельной [math]\frac{kn^2}{k \log{n}}[/math];
- Вычислительная мощность алгоритма Ланцоша без переортогонализации из последовательной сложности алгоритма [math]\frac{k(2n^2+8n-1)+3n}{n^2+2k}[/math]. При [math]k[/math] много меньше [math]n[/math] вычислительная мощность ≈ [math] 2k[/math];
- Алгоритм Ланцоша без переортогонализации не является детерминированным из-за того, что возможно выполнение меньшего числа итераций алгоритма, из-за того, что все собственные значения уже вычислены;
- Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы [math]A[/math], находящихся на границе ее спектра (в [math]T_{j}[/math] в первую очередь появляются максимальные по модулю собственные значения);
- Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;
- Нестабильность алгоритма (эффект ложной сходимости) присуще плавающей арифметике, из-за ошибок округления в которой на очередном этапе может быть построен линейно зависимый от исходных новый вектор, что повлечет невозможность дальнейшего приближения собственных значений.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Для исследования масшабируемости алгоритма была написана реализация[3] на языке C++ с использованием MPI. Реализация была протестирована на суперкомпьютере Ломоносов[4].
Были исследованы:
- время выполнения программы в зависимости от размера входных данных (матрицы);
- время выполнения в зависимости от размера параллельных нод.
Дополнительно было измерено время работы исходя из оптимизационных флагов компилятора GCC [5].
Параметры запуска алгоритма:
- размер матрицы от 20000 до 175000 с шагом 2500;
- количество процессоров от 8 до 128 с шагом 8.
Программа запускалась на суперкомьютере "Ломоносов" со следующими характеристиками:
- Компилятор GCC 5.2.1;
- Версия MPI 1.8.4;
- Сборка проводилась командой:
mpic++ -std=c++0x -O2 <CPP_FILE> -lm -static-libstdc++
Полученная эффективность реализации варьируется в пределах от 2% (на маленьких входных данных) до 45% (на больших входных данных и максимальном числе нодов).
На рисунке 4 представлена эффективность программы при отсутствии флага оптимизации O2.
На рисунке 5 представлена эффективность алгоритма с оптимизацией компилятора GCC (-O2).
Как видно из графиков выше, средняя эффективность минимальна при малых входных данных (на размере матрицы в 10000).
В среднем данный показатель держится между 9% и 14%.
При этом, оптимизация компилятора MPIC++ позволяет получить выигрыш в эффективности более чем в 3 раза (45,5% против 14,7%). Однако, как можно заметить из графиков, при оптимизации алгоритм ведет себя более нестабильно, скорее всего из-за значительной траты времени на работу с памятью (более планые участки свидетельствует о том, что необходимые данные нашлись в кеше).
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (ARPACK), так и неофициальных реализаций, выложенных на Github. Среди них:
1. The IETL Project [5]
2. MatLab [6]
3. ARPACK [7]
4. Julia Math [8]
Также существуют официальные реализации на других языках (например R). Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.
Для проверки масшабируемости, алгоритм был реализован на языке C++ с применением MPI [6].
3 Литература
- ↑ Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).
- ↑ Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).
- ↑ https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша
- ↑ Деммель Д. Вычислительная линейная алгебра
- ↑ http://www.comp-phys.org/software/ietl/
- ↑ https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554
- ↑ http://www.caam.rice.edu/software/ARPACK/
- ↑ https://github.com/JuliaMath/IterativeSolvers.jl