Участник:Dlitvinov/Алгоритм Ланцоша для арифметики с плавающей точкой с выборочной ортогонализацией: различия между версиями
Dlitvinov (обсуждение | вклад) (Новая страница: «Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 611 группа - все описан…») |
ASA (обсуждение | вклад) |
||
(не показано 18 промежуточных версий 3 участников) | |||
Строка 1: | Строка 1: | ||
− | + | Авторы статьи: Литвинов Д.А., 611, Шер А.В., 628. Авторы внесли равнозначный вклад в написание статьи. | |
{{algorithm | {{algorithm | ||
| name = Алгоритм Ланцоша для арифметики с плавающей точкой с выборочной ортогонализацией | | name = Алгоритм Ланцоша для арифметики с плавающей точкой с выборочной ортогонализацией | ||
− | | serial_complexity = <math> | + | | serial_complexity = <math>O(kn^2 + k^2n)</math> |
− | | pf_height = <math> | + | | pf_height = <math>O(kn + k^2)</math> |
− | | pf_width = <math> | + | | pf_width = <math>O(n)</math> |
− | | input_data = <math> | + | | input_data = <math>n^2 + n + 1</math> |
− | | output_data = <math> | + | | output_data = <math>nk</math> |
}} | }} | ||
Строка 28: | Строка 28: | ||
\end{array} \right]</math> | \end{array} \right]</math> | ||
− | Метод Рэлея-Ритца заключается в интерпретации собственных значений матрицы <math>T_k = Q_k^T A Q_k</math> как приближенных к собственным значениям матрицы <math>A</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>\Vert A P_k - P_k D \Vert_2</math> достигается при <math>P_k = Q_k V</math> и <math>D = \Lambda</math> и равен <math>\Vert T_ku \Vert_2</math>. | Числа и векторы Ритца являются ''оптимальными'' приближениями к собственным значениям и собственным векторам матрицы <math>A</math>: минимум величины <math>\Vert A P_k - P_k D \Vert_2</math> достигается при <math>P_k = Q_k V</math> и <math>D = \Lambda</math> и равен <math>\Vert T_ku \Vert_2</math>. | ||
Строка 36: | Строка 36: | ||
== Математическое описание алгоритма == | == Математическое описание алгоритма == | ||
− | Пусть <math>A</math> - заданная | + | Пусть <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}\\ | ||
+ | & for \; i \le k \\ | ||
+ | & \; \; \; if \; \beta_k \|v_i (k)\| \le \sqrt{\epsilon} \Vert T_k \Vert_2\\ | ||
+ | & \; \; \; \; \; \; z = z - (y_{i,k} ^ T z) y_{i,k}\\ | ||
+ | & end \; for\\ | ||
+ | & \beta_j = \Vert z \Vert_2\\ | ||
+ | & if \; \beta_j = 0 \; then\\ | ||
+ | & \; \; \; Quit\; the\; algorithm\\ | ||
+ | & q_{j+1} = z / \beta_j\\ | ||
+ | end \; & for | ||
+ | \end{align} | ||
+ | </math> | ||
+ | |||
+ | Где <math>v_i</math> - это столбцы ортогональной матрицы <math>V</math> из спектрального разложения <math>T_k = V \Lambda V^T</math>, а <math>y_{k,i} = Q_k v_i</math> - столбцы матрицы <math>Q_k V</math> - векторы Ритца. | ||
+ | |||
+ | Согласно теореме Пейджа, <math>y_{k,i}^T q_{k+1} = O(\epsilon \Vert A \Vert_2) / \beta_k \Vert v_i(k) \Vert</math>, то есть компонента <math>y_{k,i}^T q_{k+1}</math> вычисленного вектора Ланцоша <math>q_{k+1}</math> в направлении вектора Ритца <math>y_{k,i} = Q_k v_i</math> обратно пропорциональна величине <math>\beta_k \Vert v_i(k) \Vert</math>, являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца <math>q_{k,i}</math>, и следует произвести процедуру переортогонализации. | ||
+ | |||
+ | Величина погрешности <math>\beta_k \Vert v_i(k) \Vert</math> считается малой, если она меньше, чем <math>\sqrt{\epsilon} \Vert A \Vert_2</math>, так как в этом случае, согласно теореме Пейджа, компонента <math>\Vert y_{i,k}^T q_{k+1} \Vert = \Vert y_{i,k}^T z \Vert / \Vert z \Vert_2</math> скорее всего превосходит уровень <math>\sqrt{\epsilon}</math> (на практике используется <math>\Vert T_k \Vert_2</math> вместо <math>\Vert A \Vert_2</math>, так как первое известно, а второе не всегда). | ||
+ | |||
+ | == Вычислительное ядро алгоритма == | ||
+ | |||
+ | У описанного алгоритма два ядра, на которые приходится основное время работы алгоритма: | ||
+ | |||
+ | * Вычисление вектора <math>z = A q_j : O(n^2)</math> операций умножения | ||
+ | * Выборочная переортогонализация <math>z = z - (y_{i,k} ^ T z) y_{i,k} : O(n k^2)</math> операций умножения (на практике коэффициент маленький, что обеспечивает бОльшую скорость выполнения алгоритма, чем в случае с полной переортогонализацией) | ||
+ | |||
+ | == Макроструктура алгоритма == | ||
+ | |||
+ | Макроструктура алгоритма Ланцоша с выборочной переортогонализацией представляется в виде совокупности метода Ланцоша с выборочной переортогонализацией построения крыловского подпространства и процедуры Рэлея-Ритца. | ||
+ | |||
+ | Каждая итерация первого этапа может быть разложена на следующие макроединицы: | ||
+ | * Умножение матрицы на вектор | ||
+ | * Вычисление скалярного произведения векторов | ||
+ | * Умножение векторов на числа и вычитание векторов | ||
+ | * При выполнении условия необходимости переортогонализации – перемножение и вычитание векторов | ||
+ | * Вычисление нормы вектора | ||
+ | * Умножение вектора на число | ||
+ | |||
+ | Второй этап, процедура Рэлея-Ритца, заключается в нахождении собственных значений и собственных векторов матрицы и интерпретации их как приближенных значений собственных чисел и собственных векторов исходной матрицы, он может быть выполнен любым способом, например, QR-итерацией. | ||
+ | |||
+ | == Схема реализации последовательного алгоритма == | ||
+ | |||
+ | На вход алгоритму подаются исходная матрица <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_k \Vert v_i(k) \Vert</math> и выборочная переортогонализация | ||
+ | # Вычисление очередного значения <math>\beta_j</math> - норма вектора <math>z</math> | ||
+ | # Вычисление и добавление к матрице <math>Q</math> нового столбца | ||
+ | # Вычисление собственных значений и собственных векторов матрицы <math>T_j</math> | ||
+ | |||
+ | == Последовательная сложность алгоритма == | ||
+ | |||
+ | * Вычисление вектора <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>. | ||
+ | |||
+ | == Информационный граф == | ||
+ | |||
+ | Ниже представлен граф одной итерации алгоритма с разделением на <math>n</math> потоков. | ||
+ | |||
+ | [[Файл:Algorithm_scheme_(v._3).svg|1350px|thumb|center| Рисунок 1. Детальный граф алгоритма Ланцоша с выборочной переортогонализацией.<br><br> | ||
+ | Oper - простая операция сложения/вычитания/умножения/деления с числами с плавающей точкой.<br> | ||
+ | '?' - проверка условия на необходимость переортогонализации.<br><br> | ||
+ | На вход итерации подается матрица <math>A</math> и вектор <math>q_j</math>, на выходе получаем новый вектор <math>q_{j+1}</math>.]] | ||
+ | |||
+ | |||
+ | == Ресурс параллелизма алгоритма == | ||
+ | |||
+ | На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов: | ||
+ | |||
+ | * Вычисление вектора <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>. | ||
+ | |||
+ | == Входные и выходные данные алгоритма == | ||
+ | |||
+ | Входные данные: вещественная матрица <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> | ||
+ | |||
+ | == Свойства алгоритма == | ||
+ | |||
+ | * Соотношение последовательной и параллельной сложности алгоритма равно <math>n</math>, что говорит о том, что параллельная реализация теоретически должна давать ощутимый выигрыш | ||
+ | * Вычислительная мощность алгоритма - отношение числа операций к суммарному объему входных и выходных данных - примерно равно <math>2k</math> | ||
+ | * Алгоритм устойчив, так как операция переортогонализации как раз направлена на исключение проблемы округления вычисляемых значений | ||
+ | * Алгоритм не является детерминированным, так как, во-первых, используется для разряженных матриц, и, во-вторых, является итерационным с выходом по точности - число операций заранее не определено | ||
+ | * Алгоритм содержит длинные дуги, так как на протяжении всего процесса его выполнения нужно хранить исходную матрицу <math>A</math> и матрицу <math>Q</math> - они нужны на каждой итерации | ||
+ | |||
+ | = Программная реализация алгоритма = | ||
+ | |||
+ | == Особенности реализации последовательного алгоритма == | ||
+ | |||
+ | == Локальность данных и вычислений == | ||
+ | |||
+ | == Возможные способы и особенности параллельной реализации алгоритма == | ||
+ | |||
+ | == Масштабируемость алгоритма и его реализации == | ||
+ | |||
+ | Исследование масштабируемости алгоритма проводилось на суперкомпьютере | ||
+ | [//hpc.cmc.msu.ru/ IBM Blue Gene/P ВМК МГУ.] | ||
+ | |||
+ | Для тестирования была выбрана реализация | ||
+ | [//crd-legacy.lbl.gov/~osni/marques.html#BLZPACK BLZPACK]. BLZPACK написан | ||
+ | на FORTRAN 77 в 2001 году и использует MPI для параллелизма. | ||
+ | |||
+ | Для запуска BLZPACK необходима программа-драйвер. Использовался один из | ||
+ | встроенных драйверов DRVMPI (файл <tt>drv/double/drvmpi.f</tt>) c измененными | ||
+ | параметрами -- они меняются в зависимости от размера матрицы и количества | ||
+ | требуемых пар собственных значений/векторов. Параметры устанавливаются в начале | ||
+ | файла drvmpi.f. | ||
+ | |||
+ | Поскольку в мануале (<tt>doc/manual.ps</tt>) нету полного описания этих параметров, | ||
+ | прокомментируем некоторые из них. Например, для матрицы размера 7500х7500 | ||
+ | использовались следующие значения: | ||
+ | |||
+ | <pre> | ||
+ | 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) | ||
+ | </pre> | ||
+ | |||
+ | *<tt>LEIG</tt> -- размер массива для хранения полученных собственных | ||
+ | векторов/значений. | ||
+ | *<tt>LISTOR</tt> -- размер массива <tt>ISTOR</tt>. Формула определения его минимального размера приведена в мануале. | ||
+ | *<tt>LN</tt> -- порядок матрицы <tt>A</tt>, т.е. <tt>n</tt>. | ||
+ | *<tt>LRSTOR</tt> -- размер массива <tt>RSTOR</tt>. Формула определения его минимального размера приведена в мануале. | ||
+ | *<tt>NNZMAX</tt> -- количество строк во входном файле с матрицей А. Драйвер <tt>DRVMPI</tt> прочитывает симметричную верхнетреугольную матрицу из файла в координатном формате (строка, столбец, значение); таким образом, для плотной матрицы <tt>NNZMAX</tt> будет равен <math>(n^2 + n)/2</math> | ||
+ | |||
+ | Сам запуск <tt>DRVMPI</tt> конфигурируется в файле <tt>drvsp1.dat</tt> -- там задаётся количество требуемых | ||
+ | собственных значений/векторов, имя файла со входной матрицей и некоторые другие параметры. | ||
+ | |||
+ | Входные матрицы генерировались следующим простым скриптом на Python 2.7: | ||
+ | |||
+ | <pre> | ||
+ | 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 | ||
+ | |||
+ | </pre> | ||
+ | |||
+ | Был использован компилятор GNU Fortran 77 (g77), точнее, его обертка | ||
+ | <tt>mpif77</tt> с поддержкой MPI. На Blue Gene/P используется собственная | ||
+ | реализация MPI от IBM, основанная на MPICH. Для сборки BLZPACK на Blue Gene/P | ||
+ | нужно заменить имя компилятора (ключ FC) в файле <tt>sys/MACROS/Linux</tt> на | ||
+ | <tt>mpif77</tt> и запустить <tt>./creator -f -mpi</tt> в корне проекта. После | ||
+ | этого можно собрать драйвер, зайдя в директорию <tt>drv/data</tt> и выполнив | ||
+ | <tt>make drvmpi.x</tt>. | ||
+ | |||
+ | Использовался SMP-режим работы Blue Gene/P, когда на каждом вычислительном | ||
+ | узле выполняется всего один MPI-процесс. | ||
+ | |||
+ | Пример запуска программы: | ||
+ | <pre> | ||
+ | mpisubmit.bg -n 64 -w 00:15:00 -m smp --stdout eigen.out --stderr eigen.err drvmpi.x | ||
+ | </pre> | ||
+ | |||
+ | Тестировались следующие конфигурации: | ||
+ | * Порядок матриц 1000, 2500, 5000, 7500 | ||
+ | * Количество процессоров 4, 8, 16, 32, 64, 128 | ||
+ | |||
+ | У матрицы размера 1000 производился поиск 100 пар собственных векторов/значений, | ||
+ | у остальных -- 30 пар. Результаты приведены в таблице 1: | ||
+ | |||
+ | {| class="wikitable" style="border: none; background: none;" | ||
+ | |+ align="bottom" style="caption-side: bottom" | Табл. 1. Результаты тестирования. В ячейках таблицы -- время работы алгоритма в секундах. | ||
+ | ! colspan="2" rowspan="2" style="border: none; background: none;"|[[File:Pfeil_SO.svg|none|20px]] | ||
+ | ! colspan="6"| Число процессоров | ||
+ | |- | ||
+ | ! 4 !! 8 !! 16 !! 32 !! 64 !! 128 | ||
+ | |- | ||
+ | ! rowspan="4"| Размер матрицы | ||
+ | ! 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> | ||
+ | |||
+ | == Динамические характеристики и эффективность реализации алгоритма == | ||
+ | |||
+ | == Выводы для классов архитектур == | ||
+ | |||
+ | == Существующие реализации алгоритма == | ||
+ | |||
+ | Неплохой обзор существующих реализаций находится здесь: | ||
+ | [http://slepc.upv.es/documentation/reports/str6.pdf A Survey of Software for Sparse Eigenvalue Problems] | ||
+ | |||
+ | Ниже будут перечислены только те реализации, которые поддерживают именно выборочную ортогонализацию. | ||
+ | |||
+ | * [http://crd-legacy.lbl.gov/~osni/ BLZPACK] FORTRAN 77, MPI. 2000 год. | ||
+ | * [http://www.cas.mcmaster.ca/~qiao/BLKLAN BLKLAN] C, MATLAB. Нет поддержки параллелизма. 2003 год. | ||
+ | * [http://sun.stanford.edu/~rmunk/PROPACK/ PROPACK] FORTRAN 77, MATLAB. OpenMP. Предназначен для нахождения SVD разложения, но алгоритм Ланцоша может быть вызван и напрямую. 2005 год. | ||
+ | * [http://www.netlib.org/lanz/ LANZ] FORTRAN 77. Нет поддержки параллелизма. 1991 год. | ||
+ | * [http://www.netlib.org/svdpack/ SVDPACK] C, FORTRAN 77. Нет поддержки параллелизма. Предназначен для нахождения SVD разложения. 1992 год. | ||
+ | |||
+ | = Литература = | ||
+ | |||
+ | * Дж. Деммель, Вычислительная линейная алгебра. Теория и приложения. Пер. с англ. - М.: Мир, 2001. |
Текущая версия на 16:01, 3 февраля 2017
Авторы статьи: Литвинов Д.А., 611, Шер А.В., 628. Авторы внесли равнозначный вклад в написание статьи.
Алгоритм Ланцоша для арифметики с плавающей точкой с выборочной ортогонализацией | |
Последовательный алгоритм | |
Последовательная сложность | [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 Общее описание алгоритма
- 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 Общее описание алгоритма
Алгоритм Ланцоша применяется для поиска небольшого числа собственных значений большой разреженной симметричной матрицы и объединяет метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца интерпретации собственных значений некоторой вычисляемой матрицы как приближений к собственным значениям заданной матрицы.
Пусть [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]\Vert A P_k - P_k D \Vert_2[/math] достигается при [math]P_k = Q_k V[/math] и [math]D = \Lambda[/math] и равен [math]\Vert T_ku \Vert_2[/math].
При применении метода Ланцоша в арифметике с плавающей точкой ортогональность в вычисленных векторах Ланцоша [math]q_k[/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}\\ & for \; i \le k \\ & \; \; \; if \; \beta_k \|v_i (k)\| \le \sqrt{\epsilon} \Vert T_k \Vert_2\\ & \; \; \; \; \; \; z = z - (y_{i,k} ^ T z) y_{i,k}\\ & end \; for\\ & \beta_j = \Vert z \Vert_2\\ & if \; \beta_j = 0 \; then\\ & \; \; \; Quit\; the\; algorithm\\ & q_{j+1} = z / \beta_j\\ end \; & for \end{align} [/math]
Где [math]v_i[/math] - это столбцы ортогональной матрицы [math]V[/math] из спектрального разложения [math]T_k = V \Lambda V^T[/math], а [math]y_{k,i} = Q_k v_i[/math] - столбцы матрицы [math]Q_k V[/math] - векторы Ритца.
Согласно теореме Пейджа, [math]y_{k,i}^T q_{k+1} = O(\epsilon \Vert A \Vert_2) / \beta_k \Vert v_i(k) \Vert[/math], то есть компонента [math]y_{k,i}^T q_{k+1}[/math] вычисленного вектора Ланцоша [math]q_{k+1}[/math] в направлении вектора Ритца [math]y_{k,i} = Q_k v_i[/math] обратно пропорциональна величине [math]\beta_k \Vert v_i(k) \Vert[/math], являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца [math]q_{k,i}[/math], и следует произвести процедуру переортогонализации.
Величина погрешности [math]\beta_k \Vert v_i(k) \Vert[/math] считается малой, если она меньше, чем [math]\sqrt{\epsilon} \Vert A \Vert_2[/math], так как в этом случае, согласно теореме Пейджа, компонента [math]\Vert y_{i,k}^T q_{k+1} \Vert = \Vert y_{i,k}^T z \Vert / \Vert z \Vert_2[/math] скорее всего превосходит уровень [math]\sqrt{\epsilon}[/math] (на практике используется [math]\Vert T_k \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,k} ^ T z) y_{i,k} : O(n k^2)[/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_k \Vert v_i(k) \Vert[/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.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:
Число процессоров | |||||||
---|---|---|---|---|---|---|---|
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.