[непроверенная версия]

Версия 11:51, 30 января 2017

Авторы: Абдулпотиев А.А. Шаблон:Временная статья

Авторы: Абдулпотиев А.А.

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

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

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

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

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

Сам алгоритм объединяет метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца интерпретации собственных значений некоторой вычисляемой матрицы как приближений к собственным значениям исходной матрицы.

Пусть [math]Q = [Q_k,Q_u][/math] - ортогональная матрица порядка [math]n[/math], причем [math]Q_k[/math] и [math]Q_u[/math] имеют соответственно размеры [math]n \times k[/math] и [math]n \times (n-k)[/math]. Столбцы матрицы [math]Q_k[/math] вычисляются методом Ланцоша..

Запишем следующие соотношения:

[math]T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] = \left[ \begin{array}{cc} Q_k^T A Q_k & Q_k^T A Q_u\\ Q_u^T A Q_k & Q_u^T A Q_u \end{array} \right] = \left[ \begin{array}{cc} T_{k} & T_{ku}^T\\ T_{ku} & T_{u} \end{array} \right][/math]

Метод Рэлея-Ритца заключается в интерпретации собственных значений матрицы [math]T_k = Q_k^T A Q_k[/math] как приближенных к собственным значениям матрицы [math]A[/math]. Эти приближения называются числами Ритца. Если [math]A = V \Lambda V^T[/math] есть спектральное разложение матрицы [math]T_k[/math], то столбцы матрицы [math]Q_k V[/math] являются приближениями соответствующих собственных векторов и называются векторами Ритца.

Числа и векторы Ритца являются оптимальными приближениями к собственным значениям и собственным векторам матрицы [math]A[/math].

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

Алгоритм Ланцоша с полной переортогонализацией. На каждой итерации запускается повторный процесс ортогонализации Грамма-Шмидта [math]z = z - \textstyle \sum_{i=1}^{j-1} (z^T q_i)q_i \displaystyle[/math]
Алгоритм Ланцоша с выборочной ортогонализацией. Согласно теореме Пейджа векторы [math]q_k[/math] теряют ортогональность весьма систематическим образом, а именно, приобретая большие компоненты в направлениях уже сошедшихся векторов Ритца (Именно это приводит к появлению повторных копий сошедшихся чисел Ритца). Поэтому можно выборочно ортогонализировать векторы [math]q_k[/math] вместо того, чтобы на каждом шаге ортогонализировать вектор ко всем ранее сошедшимся векторам, как это делается при полной ортогонализации. Таким образом достигается (почти) ортогональности векторов Ланцоша, затрачивая очень небольшая дополнительную работу.

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

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

[math] \begin{align} q_1 = & b / \Vert b \Vert_2, \; \beta_0 = 0,\; q_0 = 0\\ for \; & j = 1 \; to \; k \\ & z = A\,q_j\\ & \alpha_j = q_j^T z\\ & z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\ & t = j-1\\ & for \; i=1 \; to \; t \\ & \; \; \; if \; \beta_t \|v_i (t)\| \le \sqrt{\epsilon} \Vert T_t \Vert_2 \; then\\ & \; \; \; \; \; \; z = z - (y_{i,t} ^ T z) y_{i,t} \;/*\; y_{i,t} = Q_t v_i \;*/ \\ & end \; for\\ & \beta_j = \Vert z \Vert_2\\ & if \; \beta_j = 0 \; then\\ & \; \; \; Stop\; the\; algorithm\\ & q_{j+1} = z / \beta_j\\ & /*Calculate\; eigenvalues\; and\; eigenvektors\; of\; T_j\; */\\ end \; & for \end{align} [/math]

Где [math]v_i[/math] - это столбцы ортогональной матрицы [math]V[/math] из спектрального разложения [math]T_t = V \Lambda V^T[/math], а [math]y_{t,i} = Q_t v_i[/math] - столбцы матрицы [math]Q_t V[/math] - векторы Ритца.

Согласно теореме Пейджа, [math]y_{t,i}^T q_{t+1} = O(\epsilon \Vert A \Vert_2) / \beta_t \Vert v_i(t) \Vert[/math], то есть компонента [math]y_{t,i}^T q_{t+1}[/math] вычисленного вектора Ланцоша [math]q_{t+1}[/math] в направлении вектора Ритца [math]y_{t,i} = Q_t v_i[/math] обратно пропорциональна величине [math]\beta_t \Vert v_i(t) \Vert[/math], являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца [math]q_{t,i}[/math], и следует произвести процедуру переортогонализации.

Величина погрешности [math]\beta_t \Vert v_i(t) \Vert[/math] считается малой, если она меньше, чем [math]\sqrt{\epsilon} \Vert A \Vert_2[/math], так как в этом случае, согласно теореме Пейджа, компонента [math]\Vert y_{i,t}^T q_{t+1} \Vert = \Vert y_{i,t}^T z \Vert / \Vert z \Vert_2[/math] скорее всего превосходит уровень [math]\sqrt{\epsilon}[/math] (на практике используется [math]\Vert T_t \Vert_2[/math] вместо [math]\Vert A \Vert_2[/math], так как первое известно, а второе не всегда).

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

У описанного алгоритма два ядра, на которые приходится основное время работы алгоритма:

Вычисление вектора [math]z = A q_j: O (n^2)[/math] операций умножения.
Выборочная переортогонализация [math]z = z - (y_{i,t} ^ T z) y_{i,t}[/math], где [math]i = 1\dots t[/math].

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

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

Каждая итерация первого этапа может быть разложена на следующие макроединицы:

Умножение матрицы на вектор
Вычисление скалярного произведения векторов
Умножение векторов на числа и вычитание векторов
При выполнении условия необходимости переортогонализации – перемножение и вычитание векторов
Вычисление нормы вектора
Умножение вектора на число

Второй этап, процедура Рэлея-Ритца, заключается в нахождении собственных значений и собственных векторов матрицы и интерпретации их как приближенных значений собственных чисел и собственных векторов исходной матрицы, он может быть выполнен любым способом, например, QR-итерацией.

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

На вход алгоритму подаются исходная матрица [math]A[/math], вектор начального приближения [math]b[/math] и число [math]k[/math]. Далее последовательно выполняются итерационно для всех [math]j: 1 \le j \le k[/math] следующие шаги:

Умножается матрица на вектор для нахождения вектора [math]z[/math]
Вычисляется коэффициент [math]\Alpha_j[/math] посредством скалярного перемножения векторов
Вычисляется новое значение вектора [math]z[/math]
Производятся оценка погрешностей [math]\beta_t \Vert v_i(t) \Vert[/math], где [math]t = j - 1[/math] и выборочная переортогонализация
Вычисление очередного значения [math]\beta_j[/math] - норма вектора [math]z[/math]
Вычисление и добавление к матрице [math]Q[/math] нового столбца
Вычисление собственных значений и собственных векторов матрицы [math]T_j[/math]

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

Вычисление вектора [math]z[/math] - [math]n^2[/math] операций умножения и [math]n(n-1)[/math] операций сложения
Вычисление скалярного произведения векторов - [math]n[/math] операций умножения и [math](n-1)[/math] сложений
Вычисление вектора [math]z[/math] - [math]2n[/math] умножений и [math]2n[/math] сложений
Оценка погрешностей - [math]kn[/math] операций умножения
Выборочная переортогонализация - [math]O(2nj)[/math] умножений и [math]O(2j(n-1) + n)[/math] сложений
Вычисление нормы вектора - [math]n[/math] операций умножения и [math](n-1)[/math] сложений
Вычисление нового столбца матрицы [math]Q[/math] - [math]n[/math] операций умножения

Итого на каждой итерации:

Умножений: [math]n^2 + n + 2n + kn + O(2nj) + n + n = n^2 + 5n + kn + O(2nj)[/math]
Сложений: [math]n(n-1) + n-1 + 2n + O(2j(n-1) + n) + n-1 = n^2 + 3n + O(2j(n-1) + n) -2[/math]

Итого за все время выполнения алгоритма:

Умножений: [math]kn^2 + 5kn + k^2n + O(k^2n + kn)[/math]
Сложений: [math]kn^2 + 3kn + O((k^2+k)(n-1)+kn) = kn^2 + 3kn + O(k^2(n-1)+2kn-k)[/math]

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

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

Ниже представлен граф одной итерации алгоритма с разделением на [math]n[/math] потоков.

Рисунок 1. Детальный граф алгоритма Ланцоша с выборочной переортогонализацией.

Oper - простая операция сложения/вычитания/умножения/деления с числами с плавающей точкой.
'?' - проверка условия на необходимость переортогонализации.

На вход итерации подается матрица [math]A[/math] и вектор [math]q_j[/math], на выходе получаем новый вектор [math]q_{j+1}[/math].

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

На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов:

Вычисление вектора [math]z[/math] - умножение матрицы размера [math]n \times n[/math] на вектор длины [math]n[/math] - [math]n[/math] ярусов сложений, [math]n[/math] операций умножений в каждом
Вычисление коэффициента [math]\alpha_j[/math] - скалярное произведение векторов - [math]n[/math] ярусов сложений, 1 операция умножения в каждом
Вычисление нового значения вектора [math]z[/math]: 2 яруса умножений длины [math]n[/math] (умножение вектора на число), 2 яруса вычитаний длины [math]n[/math]
Выборочная переортогонализация: последовательно для всех [math]i \le k[/math] выполняется [math]j[/math] ярусов сложений, [math]n[/math] операций умножений в каждом, [math]n[/math] ярусов сложений, [math]j[/math] операций умножений в каждом, один ярус вычитаний размера [math]n[/math]
Вычисление [math]\beta_j[/math] - скалярное произведение векторов - [math]n[/math] ярусов сложений, 1 операция умножения в каждом
Вычисление [math]q_j[/math] - деление вектора на число - один ярус делений размера [math]n[/math]

Сложность алгоритма по ширине ЯПФ - [math]O(n)[/math].

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

Входные данные: вещественная матрица [math]A[/math] размера [math]n \times n[/math], вектор [math]b[/math] длины [math]n[/math], число [math]k[/math]

Объем входных данных: [math]n^2 + n + 1[/math]

Выходные данные: матрица [math]Q[/math] размера [math]n \times k[/math]

Объем выходных данных: [math]nk[/math]

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

Соотношение последовательной и параллельной сложности алгоритма равно [math]n[/math], что говорит о том, что параллельная реализация теоретически должна давать ощутимый выигрыш
Вычислительная мощность алгоритма - отношение числа операций к суммарному объему входных и выходных данных - примерно равно [math]2k[/math]
Алгоритм устойчив, так как операция переортогонализации как раз направлена на исключение проблемы округления вычисляемых значений
Алгоритм не является детерминированным, так как, во-первых, используется для разряженных матриц, и, во-вторых, является итерационным с выходом по точности - число операций заранее не определено
Алгоритм содержит длинные дуги, так как на протяжении всего процесса его выполнения нужно хранить исходную матрицу [math]A[/math] и матрицу [math]Q[/math] - они нужны на каждой итерации

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

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

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

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

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

Исследование масштабируемости алгоритма проводилось на суперкомпьютере IBM Blue Gene/P ВМК МГУ.

Для тестирования была выбрана реализация BLZPACK. BLZPACK написан на FORTRAN 77 в 2001 году и использует MPI для параллелизма.

Для запуска BLZPACK необходима программа-драйвер. Использовался один из встроенных драйверов DRVMPI (файл drv/double/drvmpi.f) c измененными параметрами -- они меняются в зависимости от размера матрицы и количества требуемых пар собственных значений/векторов. Параметры устанавливаются в начале файла drvmpi.f.

Поскольку в мануале (doc/manual.ps) нету полного описания этих параметров, прокомментируем некоторые из них. Например, для матрицы размера 7500х7500 использовались следующие значения:

      INTEGER          LEIG
      PARAMETER        (LEIG   =     30)

      INTEGER          LISTOR
      PARAMETER        (LISTOR =  15000)

      INTEGER          LN
      PARAMETER        (LN     =    7500)

      INTEGER          LRSTOR
      PARAMETER        (LRSTOR =  10000000)

      INTEGER          NCUV
      PARAMETER        (NCUV   =      1)

      INTEGER          NNZMAX
      PARAMETER        (NNZMAX =  28128750)

      DOUBLE PRECISION ONE
      PARAMETER        (ONE    =  1.0D0)

      DOUBLE PRECISION ZERO
      PARAMETER        (ZERO   =  0.0D0)

LEIG -- размер массива для хранения полученных собственных

векторов/значений.

LISTOR -- размер массива ISTOR. Формула определения его минимального размера приведена в мануале.
LN -- порядок матрицы A, т.е. n.
LRSTOR -- размер массива RSTOR. Формула определения его минимального размера приведена в мануале.
NNZMAX -- количество строк во входном файле с матрицей А. Драйвер DRVMPI прочитывает симметричную верхнетреугольную матрицу из файла в координатном формате (строка, столбец, значение); таким образом, для плотной матрицы NNZMAX будет равен [math](n^2 + n)/2[/math]

Сам запуск DRVMPI конфигурируется в файле drvsp1.dat -- там задаётся количество требуемых собственных значений/векторов, имя файла со входной матрицей и некоторые другие параметры.

Входные матрицы генерировались следующим простым скриптом на Python 2.7:

import random

if __name__ == "__main__":
    n = 7500
    num_rows = (n*n + n) / 2

    with open('matr7500.data', 'w') as f:
        for i in xrange(1, n + 1):
            j = i  # column number
            while j <= n:
                value = random.uniform(-1, 1)
                f.write("    " + str(i) + "    " + str(j) + "    " + '{:.14e}'.format(value) + "\n")
                j += 1

Был использован компилятор GNU Fortran 77 (g77), точнее, его обертка mpif77 с поддержкой MPI. На Blue Gene/P используется собственная реализация MPI от IBM, основанная на MPICH. Для сборки BLZPACK на Blue Gene/P нужно заменить имя компилятора (ключ FC) в файле sys/MACROS/Linux на mpif77 и запустить ./creator -f -mpi в корне проекта. После этого можно собрать драйвер, зайдя в директорию drv/data и выполнив make drvmpi.x.

Использовался SMP-режим работы Blue Gene/P, когда на каждом вычислительном узле выполняется всего один MPI-процесс.

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

mpisubmit.bg -n 64 -w 00:15:00 -m smp  --stdout eigen.out --stderr eigen.err drvmpi.x

Тестировались следующие конфигурации:

Порядок матриц 1000, 2500, 5000, 7500
Количество процессоров 4, 8, 16, 32, 64, 128

У матрицы размера 1000 производился поиск 100 пар собственных векторов/значений, у остальных -- 30 пар. Результаты приведены в таблице 1:

Табл. 1. Результаты тестирования. В ячейках таблицы -- время работы алгоритма в секундах.
		Число процессоров
		4	8	16	32	64	128
Размер матрицы	1000	24	17	14	15	16	21
	2500	48	27	18	13	12	13
	5000	207	92	61	37	24	23
	7500	612	298	153	90	50	42

Заметно, при всех размерах матриц алгоритм перестает масштабироваться начиная с некоторого числа процессоров: для матрицы размера 1000 это 32 узла, для матрицы размера 2500 64 узла и т.д. Это связано с тем, что затраты на коммуникацию между узлами становится больше, чем выгода от распараллеливания вычислений. При меньшем же количестве процессоров алгоритм хорошо масштабируется, соответствуя теоретически вычисленной высоте ярусно-параллельной формы [math]O(kn + k^2)[/math]

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

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

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

Неплохой обзор существующих реализаций находится здесь: A Survey of Software for Sparse Eigenvalue Problems

Ниже будут перечислены только те реализации, которые поддерживают именно выборочную ортогонализацию.

BLZPACK FORTRAN 77, MPI. 2000 год.
BLKLAN C, MATLAB. Нет поддержки параллелизма. 2003 год.
PROPACK FORTRAN 77, MATLAB. OpenMP. Предназначен для нахождения SVD разложения, но алгоритм Ланцоша может быть вызван и напрямую. 2005 год.
LANZ FORTRAN 77. Нет поддержки параллелизма. 1991 год.
SVDPACK C, FORTRAN 77. Нет поддержки параллелизма. Предназначен для нахождения SVD разложения. 1992 год.

3 Литература

Дж. Деммель, Вычислительная линейная алгебра. Теория и приложения. Пер. с англ. - М.: Мир, 2001.

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