Уровень алгоритма

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

Материал из Алговики
(разн.) ← Предыдущая | Текущая версия (разн.) | Следующая → (разн.)
Перейти к навигации Перейти к поиску
Symbol wait.svgЭта работа прошла предварительную проверку
Дата последней правки страницы: 08.02.2017
Данная работа соответствует формальным критериям.
Проверено ASA.



Алгоритм Ланцоша для точной арифметики (без ортогонализации)
Последовательный алгоритм
Последовательная сложность [math]O(kn^2)[/math]
Объём входных данных [math]\frac{n(n+1)}{2} + 1[/math]
Объём выходных данных [math]kn+k[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы [math]O(k*log(n))[/math]
Ширина ярусно-параллельной формы [math]O(n^2)[/math]


Авторы: А.В. Шостик.

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

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

Алгоритм Ланцоша - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача [math]Ax = \lambda x[/math]). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

Алгоритм был предложен в 1950 году венгерским физиком и математиком Корнелием Ланцошем. Основан на понятии построения крыловского подпространства и процедуре Рэлея-Ритца: строит последовательность трехдиагональных матриц с тем свойством, что экстремальные собственные значения для каждой последующей трехдиагональной матрицы дают все более точные оценки экстремальных собственных значений для А.

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)[1].

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

Исходные данные:

  • симметрическая матрица [math]A=A^T[/math] размерности [math]n[/math], элементы матрицы [math]a_{ij}=a_{ji}[/math]
  • начальный вектор [math]v = {v_1, v_2, ..., v_n}, v \neq \theta[/math]

Так как матрица является симметрической, достаточно хранить [math]\frac{n(n+1)}{2}[/math] ее элементов.

Для исходной матрицы [math]A[/math] и вектора [math]v[/math] строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов [math][v, Av, A^2v, ..., A^{k-1}v][/math], называемых векторами Ланцоша. Каждый из векторов Ланцоша ортонормируется: [math]q_k=\frac{A^kv}{\|A^kv\|}[/math], и на шаге [math]k[/math] алгоритма имеется Крыловская матрица [math]Q = {q_0, ..., q_{k-1}}[/math] размерности [math]n \times k[/math].

Соответствующий текущей итерации [math]k[/math] вектор Ланцоша считается [math]k[/math]-ым приближением собственного вектора исходной матрицы [math]A[/math]. В качестве приближенных собственных значений матрицы [math]A[/math] берутся собственные значения симметричной трехдиагональной матрицы [math]T_k = Q^T_k A Q[/math] (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует [math]O(k^3)[/math] операций.[2].

Выходные данные:

  • совокупность собственных векторов матрицы [math]T_k[/math] (т.е. [math]k[/math]-ых приближений собственных векторов матрицы [math]A[/math])
  • совокупность собственных значений матрицы [math]T_k[/math] (т.е. [math]k[/math]-ых приближений собственных значений матрицы [math]A[/math])

1.3 Вычислительное ядро алгоритма

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

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

Итерация алгоритма:

  • вычисление нормы вектора,
  • деление вектора на число,
  • умножение матрицы на вектор,
  • линейная комбинация векторов (умножение на число и сложение).

Финальный расчет:

  • вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

Псевдокод алгоритма[3]: 

Вход  : размерность матрицы [math]n[/math], элементы симметричной матрицы [math]A[/math], количество итераций [math]k[/math]
Выход : собственные вектора матрицы [math]T_k[/math], матрица собственных значений

[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]procedure(T_k)[/math];

[math]procedure()[/math] - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

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

Глобально алгоритм состоит из инициализации, [math]k[/math] итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:

  • вычисление нормы вектора - [math](n-1)[/math] сложений, [math]n[/math] умножений, операция извлечения корня
  • деление вектора на число - [math]n[/math] операций
  • умножение матрицы на вектор - [math]n^2[/math] операций
  • скалярное произведение векторов - [math]n[/math] операций
  • линейная комбинация векторов - [math]2n[/math] умножений и вычитаний
  • вычисление собственных векторов и собственных значений матрицы [math]T_k[/math] методом "Разделяй и властвуй" - [math]O(k^3)[/math] операций.

Инициализация:

  • вычисление нормы вектора,
  • деление вектора на число.

Одна итерация:

  • умножение матрицы на вектор,
  • скалярное произведение векторов,
  • линейная комбинация векторов,
  • вычисление нормы вектора,
  • деление вектора на число.

Финальный расчет:

  • вычисление собственных векторов и собственных значений матрицы

Итого последовательная сложность алгоритма составляет [math]O(3n+k(n^2+n+2n+2n+3n))+O(k^3)[/math].

Так как на практике порядок матрицы [math]n[/math] намного превышает количество итераций [math]k[/math], сложность алгоритма Ланцоша можно рассматривать [math]O(kn^2)[/math].

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

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

Рисунок 1. Информационный граф алгоритма Ланцоша.
[math]\|.\|[/math] — операция вычисления нормы,
[math]/[/math] — операция деления,
[math]A[/math] — операция умножения матрицы на вектор,
[math]*[/math] — операция скалярного произведения,
[math]L[/math] — операция вычисления линейной комбинации векторов,
▽ - проверка условия [math]\mathsf{(\beta_i = 0)}[/math] и выход из цикла в случае выполнения условия

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

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

1.9 Входные и выходные данные алгоритма

Входные данные: симметрическая матрица [math]A=A^T[/math] размерности [math]n[/math], количество итераций [math]k[/math].

Объем входных данных: [math]\frac{n(n+1)}{2} + 1[/math].

Выходные данные: [math]k[/math] собственных векторов, собственных значений матрицы на [math]k[/math]-ом приближении.

Объем выходных данных: [math]kn+k[/math].

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

  • Сложность последовательной реализации алгоритма [math]O(kn^2)[/math].
  • Сложность параллельного алгоритма Ланцоша по высоте ЯПФ [math]O(k\log n)[/math], по ширине ЯПФ [math]O(n^2)[/math].
    • В умножении матрицы на вектор сложение n элементов выполняется в [math]\log_2 n[/math] ярусов шириной [math]\frac{n}{2}, \frac{n}{4}, ... , 1[/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]k[/math] если все собственные значения матрицы вычисляются раньше.
  • Для найденных собственных значений справедливо: [math]\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})[/math], т.е. в первую очередь находятся максимальные по модулю собственные значения.
  • Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

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

Исследование проводилось на суперкомпьютере "Ломоносов"[1].

Исследовалась реализация[4] на языке C++ с использованием MPI. Компилятор GCC 5.2.1; версия MPI 1.8.4;

Для компиляции и запуска использовались следующие команды: 

$ mpic++ -std=c++0x -O2 <CPP_FILE> -lm -static-libstdc++ -o <EXE_FILE>
$ mpirun -np <N_PROCESSORS> <EXE_FILE>
  • Число процессоров: 2, 4, 8, 16, 32, 64, 128.
  • Размерность матрицы [math]N[/math] от 5 000 до 50 000, шаг 5 000.
  • Количество итераций алгоритма Ланцоша для всех экспериментов 200.

Результаты:

Время выполнения алгоритма
2 4 8 16 32 48 64 96 128
5000 1,5243 1,0951 0,9334 0,6631 0,7456 0,7038 0,7989 0,9985 0,985
10 000 5,667 4,2625 2,7502 3,6762 3,4435 3,6192 3,5194 3,8486 3,2956
15 000 12,5207 6,5995 7,2897 8,5086 8,539 9,0018 8,9427 15,0132 7,7088
20 000 23,8143 15,851 10,6313 13,7241 12,8051 14,3183 13,2946 14,0864 13,7499
25 000 37,5419 17,9436 15,5451 18,3944 19,6076 18,5444 19,1347 21,8366 21,4454
30 000 55,014 26,2702 22,3943 26,4935 25,5618 29,6784 28,0853 27,8872 31,2415
35 000 71,3969 37,4811 30,2274 34,0041 38,1028 38,4912 38,0174 37,1407 40,4463
40 000 100,326 47,4248 40,2295 44,0332 50,0015 52,9428 46,2961 51,4372 49,7504
45 000 142,93 66,4829 52,1382 55,4683 62,3042 65,0885 60,0512 57,8891 64,3764

Ниже приведен график зависимости времени выполнения программы от числа процессоров и размера исходной матрицы.

Рисунок 2. Время работы программы в зависимости от количества процессоров и размера матрицы.
[math]Procs[/math] — количество процессоров,
[math]Size[/math] — размерность матрицы [math]N[/math],
[math]Time[/math] — время работы алгоритма в сек.

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

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

Алгоритм Ланцоша для точной арифметики без переортогонализации на практике используется гораздо реже, чем алгоритм Ланцоша с частичной или полной ортогонализацией, но, тем не менее, существует множество реализаций данного алгоритма, как последовательных, так и параллельных. Ниже перечислены наиболее известные и обозреваемые в литературе [5][6][7] проекты на базе данного алгоритма.

  • UNDERWOOD (1977) [8]. Без переортогонализации и с частичной реортогонализацией. Язык Fortran 77.
  • PROPACK (1998) [9]. Частичная реортогонализация. Язык Fortran 77, Matlab. Присутствует реализация с распараллеливанием на базе OpenMP.
  • BLKLAN (2004) [10]. Без реортогонализации, с частичной реортогонализацией. Язык C, Matlab.
  • IETL (2006) [11]. Только без реортогонализации. Язык C++.
  • SLEPc (2009) [12]. Без реортогонализации, с выборочной реортогонализацией, с частичной реортогонализацией, с полной реортогонализацией. Язык C. Присутствует реализация с распараллеливанием на базе MPI.
  • Lanczos Plus Plus (2009) [13]. Пример реализации алгоритма Ланцоша в применении к задаче расчета взаимодействия сильно кореллированных электронов. С частичной переортогонализацией, с полной переортогонализацией. Язык C++.

3 Литература

  1. ↑ Деммель Д. Вычислительная линейная алгебра
  2. ↑ https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm
  3. ↑ Деммель Д. Вычислительная линейная алгебра
  4. ↑ http://pastebin.com/3W4uUVGj
  5. ↑ Hernandez V. et al. A survey of software for sparse eigenvalue problems //Universidad Politecnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-6, – 2009.
  6. ↑ Hernandez V. et al. Lanczos methods in SLEPc //Universidad Politécnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-5. – 2006. – С. 136.
  7. ↑ http://www.grycap.upv.es/slepc
  8. ↑ Golub, G. H. and R. Underwood. The Block Lanczos Method for Computing Eigenvalues. In Mathematical Software III (edited by J. Rice), pp. 364{377. Academic Press, New York, NY, USA
  9. ↑ http://www.daimi.au.dk/PB/537
  10. ↑ Liu, G., W. Xu, and S. Qiao. Block Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS-04-07-SQ, Department of Computing and Software, McMaster University, Hamilton, Ontario, Canada
  11. ↑ http://www.comp-phys.org/software/ietl/
  12. ↑ Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems //ACM Transactions on Mathematical Software (TOMS).– Т. 31. – №. 3. – С. 351-362.
  13. ↑ https://github.com/g1257/LanczosPlusPlus
Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Shostix/Алгоритм_Ланцоша_для_точной_арифметики_(без_переортогонализации)&oldid=22468