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

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

Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 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 Общее описание алгоритма

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

Пусть [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 Макроструктура алгоритма

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

• Умножение матрицы на вектор
• Сложение, перемножение и умножение на число для векторов
• Вычисление скалярного произведения и норм векторов
• Добавление столбца к матрице

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

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

1. Умножается матрица на вектор для нахождения вектора [math]z[/math]
2. Вычисляется коэффициент [math]\Alpha_j[/math] посредством скалярного перемножения векторов
3. Вычисляется новое значение вектора [math]z[/math]
4. Производятся оценка погрешностей [math]\beta_k \Vert v_i(k) \Vert[/math] и выборочная переортогонализация
5. Вычисление очередного значения [math]\beta_j[/math] - норма вектора [math]z[/math]
6. Вычисление и добавление к матрице [math]Q[/math] нового столбца
7. Вычисление собственных значений и собственных векторов матрицы [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] потоков.

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.7 Существующие реализации алгоритма

• ARPACK library
• KTH ROYAL INSTITUTE OF TECHNOLOGY
• THE NAG LIBRARY

3 Литература

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