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

Основные авторы статьи (разделы 1, 2.4, 2.7): Башлыков Л.А. и Волков Н.И.
Авторы отдельных разделов (2.4, 2.7): Григорьев М.А., Заспа А.Ю., Левин А.Д., Слюсарь Д.Р. и Фролов А.А.
В работе также принимали участие: Антонов А.С. и Коньшин И.Н.

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

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

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

Алгоритм Ланцоша был опубликован венгерским математиком и физиком Корнелием Ланцощем в 1950 году. Этот метод является частным случаем алгоритма Арнольда в случае, если исходная матрица [math]A[/math] - симметрична, и был представлен как итерационный метод вычисления собственных значений симметричной матрицы. Этот метод позволяет за [math]k[/math] итераций вычислять [math]k[/math] приближений собственных значений и собственных векторов исходной матрицы. Хотя алгоритм и был эффективным в вычислительном смысле, но он на некоторое время был предан забвению из-за численной неустойчивости. Только в 1970 Ojalvo и Newman модифицировали алгоритм для использования в арифметике с плавающей точкой. Новый метод получил название алгоритма Ланцоша с полной переортогонализацией.

В данной статье рассматривается исходная версия алгоритма.

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

Математическое описание алгоритма Ланцоща для вычисления собственных значений симметричной матрицы А требует 1) математического описания метода Ланцоща для построения крыловского подпространства и 2) математического описания т.н. процедуры Рэлея-Ритца.

1.2.1 Метод Ланцоща для построения крыловского подпространства

Пусть дано уравнение [math]Ax = b[/math], причём вектор [math]b[/math] известен и имеется способ вычисления [math]Ax[/math].

Вычислим [math]y_{1} = b[/math] , [math]y_{2} = Ab[/math], [math]y_{3} = A^{2}b[/math],...,[math]y_{n} = A^{(n - 1)}b[/math], где [math]n[/math] - порядок матрицы [math]A[/math]. Определим матрицу [math]K = [y_{1},...,y_{n}][/math].

Тогда [math]AK = [Ay_{1},...,Ay_{n - 1}, Ay_{n}] = [y_{2},...,y_{n},A^n{y_{1}}][/math]. В предположении, что матрица [math]K[/math] - невырождена, можно вычислить вектор [math]c = -K^{-1}A^{ny_{1}}[/math]. Поскольку первые [math]n - 1[/math] столбцов в [math]AK[/math] совпадают с последними [math]n - 1[/math] столбцами в [math]K[/math], то справедливо: [math]AK = K[e_{2},e_{3},...,e_{n},-c] = KC[/math], то есть [math]K^{-1}AK = C[/math], где [math]C[/math] - верхняя хессенбергова матрица.

Заменим [math]K[/math] ортогональной матрицей [math]Q[/math] так, что при любом [math]k[/math] линейные оболочки первых [math]k[/math] столбцов в [math]K[/math] и [math]Q[/math] - одно и то же подпространство. Это подпространство называется крыловским и точно определяется как [math]\kappa_{k}(A, b)[/math] - линейная оболочка векторов [math]b[/math],[math]Ab[/math],...,[math]A^{(k - 1)}b[/math]. Суть алгоритма Ланцоша заключается в вычислении стольких первых столбцов в [math]Q[/math], сколько необходимо для получения требуемого приближения к решению [math]Ax = b [/math] ([math]Ax = {\lambda}x[/math]).

Используем QR-разложение матрицы [math]K = QQ[/math]
Тогда [math]K^{-1}AK = (R^{-1}Q^{T})A(QR) = C[/math], откуда [math]Q^{T}AQ = RCR^{-1} = H[/math]

Матрица [math]H[/math] является верхней хессенберговой в силу того, что верхней хессенберговой является матрица [math]C[/math], а матрицы [math]R[/math] и [math]R^{-1}[/math] - верхние треугольные.

Однако алгоритм Ланцоша подразумевает симметричность матрицы [math]A[/math]. Положим [math]Q = [Q_{k}, Q_{u}][/math], где [math]Q_{k} = [q_{1},...,q_{k}][/math] и [math]Q_{u} = [q_{k+1},...,q_{n}][/math]. В итоге получается, что [math]H = Q^{T}AQ = [Q_{k},Q_{u}]^{T}A[Q_{k},Q_{u}][/math] Из симметричности [math]A[/math], таким образом, следует симметричность и трехдиагональность матрицы [math]H = T[/math]. Теперь можно вывести две основные формулы, используемые в методе Ланцоша:

• [math]AQ_{j} = \beta_{j - 1}q_{j - 1} + \alpha_{j}q_{j} + \beta_{j}q_{j + 1}[/math] - получается приравниванием столбцов [math]j[/math] в обеих частях равенства [math]AQ = QT[/math].
• [math]q_{j}^{T}Aq_{j} = \alpha_{j}[/math] - получается домножением обеих частей предыдущего соотношения на [math]q_{j}^{T}[/math] и учетом ортогональности [math]Q[/math].

1.2.2 Процедура Рэлея-Ритца

Пусть [math]Q = [Q_{k}, Q_{u}][/math] - произвольная ортогональная матрица порядка [math]n[/math], причём [math]Q_{k}[/math] и [math]Q_{u}[/math] имеют размеры [math]nk[/math] и [math]n(n - k)[/math].

Пусть [math]T = Q^{T}AQ = [Q_{k},Q_{u}]^{T}A[Q_{k},Q_{u}][/math]. Обратим внимание, что в случае [math]k = 1[/math] выражение [math]T_{k} = Q_{k}^{T}AQ_{k}[/math] - отношение Рэлея [math]\rho(Q_{1},A)[/math], то есть число [math](Q_{1}^{T}AQ_{1}) / (Q_{1}^{T}Q_{1})[/math]. Определим теперь суть процедуры Рэлея-Ритца:

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

Существует несколько важных фактов о процедуре Рэлея-Ритца, характеризующих получаемые приближения собственных значений матрицы [math]A[/math].

• Минимум величины [math]||AQ_{k} - Q_{k}R||_{2}[/math] по всем симметричным [math]k^2[/math] матрицам [math]R[/math] достигается при [math]R = T_{k}[/math]. При этом [math]||AQ_{k} - Q_{k}R||_{2} = ||T_{ku}||_{2}[/math]. Пусть [math]T_{k} = V{\lambda}V^{T}[/math] - спектральное разложение матрицы [math]T_{k}[/math]. Минимум величины [math]||AP_{k} - P_{k}D||_{2}[/math], когда [math]P_{k}[/math] пробегает множество [math]nk [/math] матрицы с ортонормированными столбцами, таких, что [math]span(P_{k}) = span(Q_{k})[/math], а [math]D[/math] пробегает множество диагональных [math]k^2[/math] матриц также равен [math]||T_{ku}||_{2}[/math] и достигается для [math]P_{k} = Q_{k}V[/math] и [math]D = \lambda[/math].

Пусть [math]T_{k}[/math], [math]T_{ku}[/math], [math]Q_{k}[/math] - те же самые матрицы. Пусть [math]T_{k} = V{\lambda}V^{T}[/math] - всё то же спектральное разложение, причём [math]V = [v_{1},...,v_{k}][/math] - ортогональная матрица, а [math]\lambda = diag(\theta_{1},...,\theta_{k})[/math]. Тогда:

• Найдутся [math]k[/math] необязательно наибольших собственных значений [math]\alpha_{1},...,\alpha_{k}[/math] матрицы [math]A[/math], такие, что [math]|\theta_{i} - \alpha_{i}| \lt = ||T_{ku}||_{2}[/math] для [math]i = 1,...,k[/math]. Если [math]Q_{k}[/math] вычислена методом Ланцоша,то правая часть этого выражения равняется [math]\beta_{k}[/math], то есть единственному элементу блока [math]T_{ku}[/math], который может быть отличен от нуля. Указанный эдемент находится в правом верхнем углу блока.
• [math]||A(Q_{k}v_{i}) - (Q_{k}v_{i})\theta_{i}||_{2} = ||T_{ku}v_{i}||_{2}[/math]. Таким образом, разность между числом Ритца [math]\theta_{j}[/math] и некоторым собственным значением [math]\alpha[/math] матрицы [math]A[/math] не превышает величины [math]||T_{ku}v_{i}||_{2}[/math], которая может быть много меньше, чем [math]||T_{ku}||_{2}[/math]. Если матрица [math]Q_{k}[/math] вычислена алгоритмом Ланцоша, то эта величина равна [math]\beta_{k}|v_{i}(k)|[/math], где [math]v_{i}(k)[/math] - k-ая компонента вектора [math]v_{i}[/math]. Таким образом, можно дещево вычислить норму невязки - левой части выражения, не выполняя ни одного умножения вектора на матрицы [math]Q_{k}[/math] или [math]A[/math].
• Не имея информации о спектре матрицы [math]T_{u}[/math], нельзя дать какую-либо содержательную оценку для погрешности в векторе Ритца [math]Q_{k}v_{i}[/math]. Если известно, что [math]\theta_{j}[/math] отделено расстоянием, не меньшим [math]g[/math], от прочих собственных значений матриц [math]T_{k}[/math], [math]T_{u}[/math] и матрица [math]Q_{k}[/math] вычислена алгоритмом Ланцоша, то угол [math]\theta[/math] между [math]Q_{k}v_{i}[/math] и точным собственным вектором [math]A[/math] можно оценить формулой [math]\frac{1}{2}\sin{2\theta} \lt = \frac{b_{k}}{g}[/math]

В алгоритме Ланцоша из ортонормированных векторов строится матрица [math]Q_{k} = [q_{1},...,q_{k}][/math] и в качестве приближенных собственных значений [math]A[/math] принимаются числа Ритца - собственные значения симметричной трёхдиагональной матрицы [math]T_{k} = Q_{k}^{T}AQ_{k}[/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}AQ_{j}[/math] порядка [math]j[/math]. Итерационный процесс завершается, когда [math]j = k[/math]. Вторая часть алгоритма представляет собой поиск собственных значений и собственных векторов симметричной трёхдиагональной матрицы [math]T_{j}[/math], что реализуется отдельным алгоритмом, на практике используется QR-алгоритм.

Особенность методов крыловского подпространства состоит в предположении, что матрица [math]A[/math] доступна в виде "черного ящика", а именно подпрограммы, вычисляющей по заданному вектору [math]z[/math] произведение [math]y = Az[/math] (а также, возможно, произведение [math]y = A^{T}z[/math], если матрица [math]A[/math] несимметрична). Другими словами, прямой доступ к элементам матрицы и их изменение не используются, т.к.

• Умножение матрицы на вектор - наиболее дешевая нетривиальная операция, которую можно проделать с разреженной матрицей: в случае разреженной [math]A[/math], содержащей [math]m[/math] ненулевых элементов, для матрично-векторного умножения нужны [math]m[/math] скалярных умножений и не более [math]m[/math] сложений.
• [math]A[/math] может быть не представлена в виде матрицы явно, а доступна именно через подпрограмму для вычисления произведений [math]Ax[/math].

Таким образом, как QR-алгоритм поиска собственных значений, собственных векторов матрицы [math]T_{j}[/math] и оценки погрешности в них результирующей трехдиагональной матрицы, так и умножение матрицы на вектор внутри каждой итерации могут быть рассмотрены в качестве отдельных макроопераций. Проблемой является то, что именно эти манипуляции порождают основную вычислительную сложность алгоритма. Поскольку матрица, к которой применяется QR-алгоритм, имеет порядок [math]k[/math] и на практике бывает мала, сделаем допущение о плотности матрицы [math]A[/math]. Тогда умножение этой матрицы на вектор внутри каждой итерации нужно будет реализовывать в рамках самого алгоритма. Еще одним важным замечанием является момент вычисления собственных значений матрицы [math]T_{j}[/math]. Её собственные значения, полученные на итерации [math]p[/math], никак не используются на итерации [math]p + 1[/math], поэтому целесообразно вынести эту процедуру за основной цикл. Однако на практике в силу малых размеров матрицы [math]T_{j}[/math] вычисление соответствующих величин производится непосредственно в теле основного цикла, что позволяет сразу оценить достигнутую точность вычислений.

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

Алгоритм итерационный, на каждой итерации выполняется вычисление очередного столбца матрицы [math]Q[/math], а также вычисление элементов очередной пары ( диагонального и околодиагонального ) элементов матрицы [math]T[/math] с использованием значений столбца матрицы [math]Q[/math], полученного на предыдущей итерации и значения околодиагонального элемента матрицы [math]T[/math] с предыдущей итерации. Вектор [math]z[/math] носит вспомогательный характер. Коэффициенты [math]\alpha_{j}[/math] соответствуют диагональным элементам матрицы [math]T_{j}[/math], коэффициенты [math]\beta_{j}[/math] – наддиагональным и поддиагональным элементам той же матрицы. Далее приводится пример псевдокода. По окончании итераций вычислеяются собственные значения и собственные векторы матрицы [math]T_{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]6k^3[/math], то есть [math]O(k^3)[/math] операций; при использовании метода «Разделяй-и-властвуй» для вычисления собственных значений и векторов аналогичной матрицы в среднем затрачивается [math]O(k^{2})[/math] операций.

Таким образом, для выполнения одной итерации метода Ланцоща требуется 1 операция вычисления квадратного корня, [math]n[/math] умножений, а также по [math]n ^ 2 + n + 2n + 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] потоков, чего можно было бы достичь, используя суммирование сдваиванием. Это допущение аналогично допущению из статьи Умножение матрицы на вектор и Это же допущение будет использовано и при оценке ресурса параллелизма алгоритма. Приводится граф без отображения входных и выходных данных.

Рис. 1: Граф алгоритма для матрицы А порядка 5 и k = 3. Op - одна итерация алгоритма, Eigen - процедура вычисления собственных значений

Рис. 2: Внутренний граф итерации алгоритма для матрицы 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]T_{j}[/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 + kn[/math]

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

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

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

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

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

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

Проведём исследование масштабируемости параллельной реализации алгоритма Ланцоша без переортогонализации согласно методике AlgoWiki. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета.

2.4.1 Реализация 1

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

Рис. 3: Производительность алгоритма для указанных значений размера задачи и числа MPI процессов (по 4 OpenMP нити каждый). Потенциальная максимальная производительность, которая может быть достигнута - 13.6 GFLOPS на 1 MPI процесс. Тогда, например, эффективность вычислений на 100 процессах при размере задачи 4000[math]{\cdot}[/math]4000 составляет приблизительно 14.33%

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

• Число MPI процессов [1 : 200]
• Число OpenMP нитей [1 : 4] (Запуск с 1 OpenMP нитью производился только для 1 MPI процесса, ускорение замерялось относительно этого запуска).
• Размерность матрицы [400 : 4000]

Эффективность выполнения реализации алгоритма:

• Минимальная эффективность - <1% (неактуально, так как в исследовании достигнут предел масштабируемости алгоритма).
• Максимальная эффективность - ~30% (получается на 2-8 процессов при ~640000 ячеек матрицы на процесс; в случае 8 процессов это входная матрица 1600[math]{\cdot}[/math]1600).

Оценка масштабируемости:

• По числу процессов - при увеличении числа процессов эффективность уменьшается на всей области рассмотренных значений, но не слишком интенсивно. Это связано с небольшим объёмом пересылок данных, реализованных через достаточно эффективные операции редукции.
• По размеру задачи - при увеличении размера задачи эффективность вычислений вначале кратковременно возрастает, но затем начинает относительно равномерно убывать на всей рассматриваемой области.
• По двум направлениям - при рассмотрении увеличения, как вычислительной сложности, так и числа процессов по всей рассмотренной области значений наибольшая эффективность наблюдается при соответствии числа процессов и размера задачи. По этой "диагонали" эффективность вначале кратковременно возрастает, затем начинает убывать.

Интересно, что для разного числа процессоров положительные скачки производительности (зачастую - довольно заметные) имеют место на разных объёмах входных данных.

Далее приводятся графики ускорения программы для разных объёмов входных данных. Ускорение считалось относительно варианта с 1 MPI процессом без дальнейшнего распараллеливания посредством OpenMP. На оси абсцисс отмечается число MPI процессов, число OpenMP нитей во всех случаях равно четырём.

Рис. 4: Ускорение для размера задачи 2000[math]{\cdot}[/math]2000

Рис. 5: Ускорение для размера задачи 3000[math]{\cdot}[/math]3000

Рис. 6: Ускорение для размера задачи 4000[math]{\cdot}[/math]4000

Рис. 6: Сравнение ускорения для различных размеров задачи

Используемая реализация алгоритма Реализация алгоритма гибридная, использует как MPI версии 2.2 (использование MPI_IN_PLACE), так и OpenMP. Указанный код компилировался с помощью компилятора IBM (mpixlcxx_r) с опцией -qsmp=omp (она включает в себя оптимизацию -O2). Реализация алгоритма полностью собственная и не использует встроенные библиотеки. Каждая "quadcore node" - это 1 узел системы BG/P, содержащий 4 микропроцессорных ядра IBM PowerPC 450 c суммарной пиковой производительностью 13.6 GFLOPS на узел и 2 GB общей памяти с пропускной способностью 13.6 GB/sec.

2.4.2 Реализация 2

Для исследования масштабируемости рассматриваемого алгоритма Ланцоша без переортогонализации автором была реализована программа на языке С++ с использованием технологии MPI. Дальнейшее исследование было проведено на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса МГУ.
Важно отметить, что реализованная программа выполняет только первую часть алгоритма - формирование трёхдиагональной матрицы [math]T_j[/math]. Вторая часть алгоритма (задача поиска собственных значений матрицы [math]T_j[/math]) может быть решена несколькими методами, каждый из которых обладает своей эффективностью распараллеливания и характерными свойствами. Подробно этот вопрос уже обсуждался в статье ранее, поэтому не будем ещё раз заострять на этом внимание и останавливаться.

Сборка осуществлялась со следующими параметрами:

• Компилятор intel/15.0.090
• Openmpi/1.8.4-icc

В серии проведённых расчётов рассматривались следующие числовые значения параметров:

• Размер [math]n[/math] задачи и, соответственно, матрицы [math]A[/math], задействованной в операции умножения на вектор: [2000 : 30000] с шагом 1000
• Число [math]p[/math] процессоров: 1, 2, 4, 16, 32, 48, 64, 96, 128, 192, 256
• Число итераций в алгоритме: [math]k=350[/math]

Число [math]k[/math] взято таким вопреки информации в начале статьи, где говорится, что значение [math]k[/math] редко превышает [math]10^2[/math] (хоть мы и имеем право брать его любым по нашему усмотрению). Это сделано исключительно для лучшей визуализации пиков на графиках в данном разделе и наглядности получаемых результатов, так как при меньших на порядок значениях параметра [math]k[/math] времена выполнения программ были слишком малы даже при минимальных числах процессоров и максимальных размерах [math]n[/math] задачи.
На рис. 2 изображён график зависимости времени выполнения программы от числа процессоров и размера [math]n[/math] задачи.

Рис. 2 График зависимости времени выполнения программы.

Рис. 3 График зависимости эффективности распараллеливания программы.

Многочисленные тесты показали при малых количествах ядер уменьшение времени выполнения программы почти в 2 раза при увеличении числа ядер в 2 раза на всём диапазоне размеров матриц, что и отражает сильный излом на первом графике. Изменение значений времени можно легко оценить по предоставленной цветовой шкале. Был построен график эффективности распараллеливания (parallel efficiency) в зависимости от числа процессоров и размера задачи, изображённый на рис. 3. Получены следующие результаты численных экспериментов:

• Средняя эффективность распараллеливания программы составила: 42.5576% ;
• Минимальная эффективность реализации: 0.674% (число процессоров 256 и минимальный размер матрицы 2000);
• Максимальная эффективность реализации: 77.773% (число процессоров 2, размер матрицы 20000).

Нужно понимать, что это лишь экстремумы, соответствующие двум экспериментам, которые сами по себе не дают полной картины и не позволяют однозначно судить о параллельных качествах алгоритма и реализующей его программы при столь большом числе параметров и проводимых расчётов.
На основании рис. 3 можно отметить следующие значимые и важные факты:

• Имеется ярко выраженный спад эффективности при минимальном рассматриваемом размере матрицы [math]n=2000[/math] вдоль всей оси OY, отвечающей за число процессоров;
• Присутствуют два пика эффективности (жёлтого цвета на рис. 3): при числе [math]p[/math] процессоров 2 и 128. Затем при дальнейшем увеличении числа процессоров с 128 до 256 наблюдается постепенное уменьшение эффективности реализации. Это связано с увеличением накладных расходов и затрат времени на обмены сообщениями между процессами.

2.4.3 Реализация 3

Для исследования масшабируемости алгоритма была написана реализация[1] на языке C++ с использованием MPI. Реализация была протестирована на суперкомпьютере Ломоносов[2].

Были исследованы:

• время выполнения программы в зависимости от размера входных данных (матрицы);
• время выполнения в зависимости от размера параллельных нод.

Дополнительно было измерено время работы исходя из оптимизационных флагов компилятора GCC [3].

Параметры запуска алгоритма:

• размер матрицы от 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.

Рисунок 4. Эффективность параллельной реализации алгоритма Ланцоща на MPI.

На рисунке 5 представлена эффективность алгоритма с оптимизацией компилятора GCC (-O2).

Рисунок 5. Эффективность параллельной реализации алгоритма Ланцоща на MPI с использованием оптимизации компилятора -O2.

Как видно из графиков выше, средняя эффективность минимальна при малых входных данных (на размере матрицы в 10000). В среднем данный показатель держится между 9% и 14%.

При этом, оптимизация компилятора MPIC++ позволяет получить выигрыш в эффективности более чем в 3 раза (45,5% против 14,7%). Однако, как можно заметить из графиков, при оптимизации алгоритм ведет себя более нестабильно, скорее всего из-за значительной траты времени на работу с памятью (более плавные участки свидетельствует о том, что необходимые данные нашлись в кеше).

2.4.4 Реализация 4

Для исследования масшабируемости алгоритма была написана реализация[4] на языке C с использованием MPI. Реализация была протестирована на суперкомпьютере Ломоносов[5].

Сборка осуществлялась со следующими параметрами:

• gcc-5.2.0
• openmpi-1.8.4
• аргументы компилятора: -std=c11 -Ofast

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

• число процессоров [32 : 256] с шагом 16;
• размер матрицы [5000 : 200000] с шагом 5000.

В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:

• минимальная эффективность реализации 1.92%;
• максимальная эффективность реализации 59.47%.

График полученного распределения эффективности:

Рисунок 4. Параллельная реализация алгоритма Ланцоша. Распределение производительности.

На следующих рисунках приведены графики производительности и эффективности данной реализации в зависимости от изменяемых параметров запуска.

Рисунок 5. Параллельная реализация алгоритма Ланцоша. Изменение производительности в зависимости от числа процессоров и размера матрицы.

Рисунок 6. Параллельная реализация алгоритма Ланцоша. Изменение эффективности в зависимости от числа процессоров и размера матрицы.

Средняя эффективность для матриц с размером больше 25000 составила 52.7%.

2.4.5 Реализация 5

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

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

Для компиляции программы использовались компиляторы intel/15.0.090 и Openmpi/1.8.4-icc.

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

• в качестве размера входной матрицы подавались значения в диапазоне [2000:30000] c шагом 2000.
• Число процессоров варьировалось от 1 до 256.
• В качестве числа, отвечающего за количество выполняемых методом итераций бралось значение 100.

Ниже на рисунке изображен трехмерный график, показывающий зависимость времени выполнения программы от входных данных/

График зависимости времени выполнения программы от входных данных

Также была исследована эффективность распараллеливания:

• Средняя эффективность выполнения алгоритма составила порядка 12%
• Наивысший показатель равен 43%
• Минимальное значение составило менее 1%

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 nonsymmetrical matrices (lookahead)	none
SVDPACK	C/F77	1992	P, finds SVD	none
Underwood	F77	1975	F, block version	none

В настоящее время алгоритм Ланцоша поиска собственных значений квадратной симметричной матрицы включён в несколько библиотек и реализован на различных языках, среди которых такие наиболее распространенные как C, C++, FORTRAN77/90, MathLab. Ссылки на наиболее проверенные и часто используемые реализации алгоритма:

• Так, например, в языке MathLab во встроенном пакете ARPACK найти собственные значения методом Ланцоша можно при помощи вызова функции eigs().

• На языке С существует библиотека PRIMME, название которой расшифровывается как PReconditioned Iterative MultiMethod Eigensolver (итерационные методы поиска собственных значений с предусловием).

• Библиотека NAG Numerical Library также содержит обширный набор функций, позволяющих проводить математический анализ, среди которых есть и реализация алгоритма Ланцоша. Библиотека доступна в трех пакетах: NAG C Library, NAG Fortran Library и NAG Library for .NET, совместна со многими языками и с основными операционными системами (Windows, Linux and OS X, а также Solaris, AIX и HP-UX).

Заслуживают внимания также реализации в проектах IETL и ARPACK .

3 Литература

<references \>