Участник:Dlitvinov/Алгоритм Ланцоша для арифметики с плавающей точкой с выборочной ортогонализацией
Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 611 группа - все описанные разделы.
Алгоритм Ланцоша для арифметики с плавающей точкой с выборочной ортогонализацией | |
Последовательный алгоритм | |
Последовательная сложность | [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 Масштабируемость алгоритма и его реализации
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
3 Литература
- Дж. Деммель, Вычислительная линейная алгебра. Теория и приложения. Пер. с англ. - М.: Мир, 2001.