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

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

Материал из Алговики
Перейти к навигации Перейти к поиску

Авторы статьи:
Волков Н. И. ( группа 623 ) - разделы 1.1 , 1.3 , 1.7 , 1.9 , 2.7
Башлыков. Л. А. ( группа 616 ) - разделы 1.2 , 1.4
Совместно: 1.5 , 1.6 , 1.8 , 1.10

Symbol wait.svgЭта работа ждет рассмотрения преподавателем
Дата последней правки страницы:
21.10.2016
Авторы этой статьи считают, что задание выполнено.



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


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

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

Метод Ланцоша — процедура, применяемая для нахождения [math]k[/math] собственных значений симметричной матрицы произвольного размера путём вычисления их приближений. Суть метода заключается в итеративном построении специальной матрицы [math]T_{j}[/math], а затем вычислении её собственных значений и собственных векторов.

Здесь мы рассмотрим упрощенный вариант алгоритма Ланцоша, который подразумевает отсутствие влияния ошибок округления на вычислительный процесс. На практике используются модификации алгоритма Ланцоша, обладающие устойчивостью и без этого допущения ( например, алгоритм Ланцоша с полной переортогонализацией ).

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

Этот раздел явно требует некоторой доработки.

На вход алгоритму Ланцоша подаётся симметричная матрица порядка [math]n[/math], для которой необходимо найти собственные значения. Для этого из ортонормированных векторов Ланцоша строится матрица [math]Q_{k} = [q_{1}, .. , q{k}][/math] , и в качестве приближенных собственных значений исходной матрицы принимаются числа Ритца, т.е. собственные значения симметричной трёхдиагональной матрицы, вычисляемой по основной формуле метода: [math]T_{j} = Q^T_{j} * A * Q_{j}[/math], где [math]Q_{j}[/math] – ортогональная матрица Ланцоша, [math]A[/math] - исходная матрица, [math]T_{j}[/math] – симметричная трехдиагональная матрица, собственные значения которой — приближения собственных значений матрицы [math]A[/math].

При этом алгоритм может быть доработан для оценки погрешности полученных собственных значений [math]A[/math]. Пусть [math]T_{k} = V * \lambda * V^T[/math] - спектральное разложение матрицы [math]T_{k}[/math]. Тогда столбцы матрицы [math]Q_{k} * V[/math] - приближения к соответствующим собственным векторам, называемым векторами Ритца.

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

Вычислительное ядро алгоритма Ланцоща состоит в получении на каждой итерации очередного промежуточного вектора [math]z[/math], получаемого путём умножения исходной матрицы [math]A[/math] на вектор Ланцоша [math]q_{j}[/math], полученный на предыдущей итерации.

Также при значениях [math]k[/math], сопоставимых с [math]n[/math], процедура вычисления собственных значений и собственных векторов симметричной трёхдиагональной матрицы по некоторому выбранному алгоритму подразумевает значительный объём расчётов и может считаться вычислительным ядром. Однако на практике метод используется прежде всего при малых [math]k[/math].

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

Метод Ланцоша соединяет метод Ланцоша для построения крыловского подпространства в процедурой Рэлея-Ритца, подразумевающей вычисление собственных значений симметричной трёхдиагональной матрицы.

Первая часть алгоритма представляет собой строго последовательные итерации, на каждой из которых вначале строится очередной столбец матрицы [math]Q_{j}[/math] из первых [math]j[/math] ортонормированных векторов Ланцоша и матрица [math]T_{j} = Q^T_{j} * A * Q_{j}[/math] порядка j. Итерационный процесс завершается, когда [math]j = k[/math]

Вторая часть алгоритма представляет собой поиск собственных значений и собственных векторов симметричной трёхдиагональной матрицы [math]T_{j}[/math], что реализуется отдельным алгоритмом ( QR-алгоритм, метод «Разделяй-и-властвуй», итерации Арнольди и.т.д. ), который мы будем рассматривать как отдельную макрооперацию.

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

Коэффициенты [math]alpha_{j}[/math] соответствуют диагональным элементам матрицы [math]T_{j}[/math], коэффициенты [math]beta_{j}[/math] – наддиагональным и поддиагональным элементам той же матрицы.

Вычислить столбец q1 матрицы Qk q1 = b / ||b||
q0 - нулевой вектор
beta_0 = 0
foreach (j = 1 .. k ) 
{
     z = A * qj
     alpha_j = qj * z 
     z = z – alpha_j * qj – beta_(j-1) * q(j – 1)
     beta_j = ||z||
     if ( beta_j == 0 )
     {
         break // Все ненулевые собственные значения найдены. Выход может быть досрочным.
     }
     else 
     {
         q(j + 1) = z/beta_j
     }
}
Вычислить собственные значения и собственные векторы симметричной трехдиагональной матрицы Tj наиболее подходящим методом.
Полученный вектор Lj есть искомые собственные значения.

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

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

  • Умножение квадратной матрицы порядка [math]n[/math] на вектор длины [math]n[/math] и наоборот требует по [math]n^2[/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[/math] QR-алгоритмом требует [math]O(k^2)[/math] операций, вычисление также и собственных векторов требует примерно [math]6 * k ^ 3[/math], то есть [math]O(k^3)[/math] операций; при использовании метода «Разделяй-и-властвуй» для вычисления собственных значений и векторов аналогичной матрицы в среднем затрачивается [math]O(k^{2.3})[/math] операций.

Таким образом, для выполнения одной итерации метода Ланцоща требуется 1 операция вычисления квадратного корня, [math]n[/math] умножений, а также по [math]n ^ 2 + n + 2 * n + n[/math] сложений и умножений, а также суммарно [math]O(k ^ 3)[/math] операций, требуемых для поиска собственных значений и собственных векторов матрицы [math]T_{j}[/math]. То есть последовательная сложность алгоритма Ланцоша составляет [math]O(n ^ 2) + O(k ^ 3)[/math]. Это, очевидно, меньше, чем [math]O(n ^ 3)[/math] операций, требуемых в алгоритмах вычисления всех собственных значений произвольных симметричных матриц, в чём и заключается выгода метода Ланцоша.

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

Информационный граф алгоритма требует следующего замечания: подразумевается, что, например, умножение квадратной матрицы на вектор распараллеливается на [math]n[/math] потоков, где [math]n[/math] - порядок матрицы, а не на [math]n * \log{n}[/math] потоков, чего можно было бы достичь, используя суммирование сдваиванием. Это допущение аналогично допущениям из статей от матрично-векторных операциях. Это же допущение будет использовано и при оценке ресурса параллелизма алгоритма. Приводится граф без отображения входных и выходных данных.

Граф алгоритма для матрицы А порядка 5 и k = 3. Op - одна итерация алгоритма, Eigen - процедура вычисления собственных значений
Внутренний граф итерации алгоритма для матрицы A порядка 5. Vector sum - нахождение суммы элементов вектора-поэлементного произведения векторов z и qj. Op - вычитания из z векторов qj и q(j - 1), домноженных на коэффициенты. Norm - вычисление нормы вектора z.

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

Каждая итерация алгоритма Ланцоща выполняется строго последовательно. Однако внутри итерации возможно распараллеливание алгоритма.

  • Умножение квадратной матрицы порядка [math]n[/math] на вектор длины [math]n[/math] требует последовательного выполнения [math]n[/math] ярусов умножений и сложений.
  • Все остальные операции в рамках одной итерации, ведущие к вычислению матрицы [math]Tj[/math] выполняются строго последовательно. При этом вычисление значений векторов может быть выполнено за 1 ярус умножений / сложений, а вычисление чисел - за [math]\log{n}[/math] таких операций.
  • Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма и рассматривается в соответствующей статье.

Таким образом, при классисифкации по ширине ЯПФ алгоритм обладает линейной сложностью. В случае классификации по высоте ЯПФ алгоритм также имеет линейную сложность.

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

Входные данные: квадратная матрица [math]A[/math] ( элементы [math]a_{ij}[/math] ), вектор [math]b[/math], число [math]k[/math]
Объём входных данных: [math]n^2 + n + 1[/math]
Выходные данные: вектор [math]L[/math], матрица [math]E[/math] ( элементы [math]e_{pq}[/math] )
Объём выходных данных: [math]k + k * n[/math]

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

В случае неограниченности ресурсов соотношение последовательной и параллельной сложности алгоритма Ланцоща по сумме всех итераций представляет собой [math]k ^ 3 + k * n ^ 2[/math] к [math]X + k * n[/math], где [math]X[/math] - параллельная сложность алгоритма, используемого для вычисления собственных значений. В случае использования метода вычисления собственных значений с линейной сложностью при классисифкации по высоте ЯПФ это соотношение будет равно [math](k ^ 2 + n ^ 2) / (k + n)[/math].

Алгоритм Ланцоша не является полностью детерминированным в том смысле, что возможно выполнение числа итераций алгоритма, меньшего, чем заданное ( из-за того, что вычислены все ненулевые собственные значения матрицы [math]A[/math] ).

Важное свойство метода Ланцоша состоит в том, что первыми в матрице [math]T_{j}[/math] появляются собственные значения с максимальной величиной по модулю. Таким образом, метод особенно хорошо подходит для вычисления собственных значений матрицы [math]A[/math], находящихся на краях её спектра.

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

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

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

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

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

Для проверки масштабируемости алгоритма из реализации была исключена непосредственно задача поиска собственных значений матрицы [math]T_{j}[/math] на каждой итерации, так как это является предметом отдельных статей ( QR-алгоритм, метод «Разделяй-и-властвуй», итерации Арнольди и.т.д. ). Таким образом, проверяемая на масштабируемость задача сводится к итерационному вычислению коэффициентов матрицы [math]T_{j}[/math]. Матрица [math]A[/math] в рассматриваемой реализации хранится построчно и построчно же разделяется между процессорами.

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

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

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

Существует несколько как последовательных, так и параллельных реализаций алгоритма Ланцоша, включенных в программные библиотеки для поиска собственных значений матриц. Свойства этих реализаций можно увидеть в таблицах ниже. Первая таблица — относительно недавние реализации, вторая - «музейные экспонаты». Следует уделить особенное внимание столбцам «тип метода» и «параллелизация». В первом из них значение только N соответствует описанному варианту без реортогонализации, F – полной реортогонализации, P – частичной реортогонализации, S – выборочной реортогонализации. Во втором значние «none» соответствует отсутствию параллельной реализации, M – параллельной реализации посредством MPI, O – параллельной реализации посредством OpenMP. Их приведение целесообразно, так как на практике алгоритм Ланцоша без переортогонализации неустойчив. Также указываются библиотеки, в которых реализованы более глубокие модификации метода Ланцоша, с указанием изменений в графе «тип метода».

Таблица 1 - «современные» реализации.
Название библиотеки Язык программирования Дата появления, версия библиотки Тип метода Параллелизация
BLKLAN C/Matlab 2003 P none
BLZPACK F77 2000, 04/00 P + S M
IETL C++ 2006, 2.2 N none
SLEPc C/F77 2009, 3.0.0 All M
TRLAN F90 2006 Dynamic thick-restart M
PROPACK F77/Matlab 2005, 2.1 / 1.1 P, finds SVD O
IRBLEIGS Matlab 2000, 1.0 Indefinitie symmetric none
Таблица 2 - «исторические» реализации.
Название библиотеки Язык программирования Дата появления, версия библиотки Тип метода Параллелизация
ARPACK F77 1995, 2.1 Arnoldi iterations, impliicit restart M
ARPACK++ C++ 1998, 1.1 Arnoldi iterations, impliicit restart none
LANCZOS F77 1992 N none
LANZ F77 1991, 1.0 P none
LASO F77 1983, 2 S none
NAPACK F77 1987 N none
QMRPACK F77 1996 Designed for nonsymmetric matrices ( lookahead ) none
SVDPACK C/F77 1992 P, finds SVD none
Underwood F77 1975 F, block version none

3 Литература

Деммель Дж. Вычислительная линейная алгебра. Теория и приложения - М.: Мир, 2001. Параллельные вычисления (Воеводин В.В., Воеводин Вл.В.) - Спб, изд-во "БХВ-Петербург", 2002 A Survey of Software for Sparse Eigenvalue Problems - V. Hern´andez, J. E. Rom´an, A. Tom´as, V. Vidal