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

Авторы: Хакимов А. С.

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

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

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 1950 году. Его эффективность обусловлена экономией памяти для хранения матриц и экономией вычислительных ресурсов. Алгоритм итерационный и использует степенной метод для поиска наибольших собственных значений и векторов матриц. Основной недостаток алгоритма заключается в накоплении ошибок округления, для решения которых появились методы поддержания ортогонализации т.н. векторов Ланцоша. Здесь мы рассмотрим выборочный метод поддержания ортогонализации, который существенно экономит процессорное время.

На вход алгоритма подаётся [math]A = A^T[/math],

[math] A = \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1\ n-1} & a_{1\ n} \\ a_{12} & a_{22} & a_{23} & \cdots & a_{2\ n-1} & a_{2\ n} \\ a_{13} & a_{23} & a_{33} & \cdots & a_{3\ n-1} & a_{3\ n} \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ a_{1\ n-1} & \cdots & \cdots & a_{n-2\ n-1} & a_{n-1\ n-1} & a_{n-1\ n} \\ a_{1\ n} & \cdots & \cdots & a_{n-2\ n} & a_{n-1\ n} & a_{n\ n} \\ \end{pmatrix} [/math] [math] ,\, \;[/math]

случайный вектор [math]b[/math], как первое приближение собственного вектора матрицы и [math]k [/math] - количество собственных значений и собственных векторов, которые требуется найти.

Матрица [math]Q_j = [q_1, q_2, \dots, q_j][/math] размерности [math]n \times j[/math] строится на каждой итерации и состоит из ортонормированных векторов Ланцоша. А в качестве приближенных собственных значений берутся числа Ритца [math]\theta_i [/math], - собственные значения симметричной трехдиагональной матрицы [math]T_j = Q^T_j A Q_j[/math] размерности [math]j \times j[/math].

[math] T_j = \begin{pmatrix} \alpha_1 & \beta_1 \\ \beta_1 & \alpha_2 & \beta_2 \\ & \beta_2 & \ddots & \ddots \\ & & \ddots & \ddots & \beta_{j-1} \\ & & & \beta_{j-1} & \alpha_j \end{pmatrix}\; (2). [/math]

Однако, векторы [math]q_j [/math] теряют ортогональность вследствие приобретения больших компонент в направлениях векторов Ритца [math]y_{i,j} = Q_j v_i [/math], отвечающих сошедшимся числам Ритца [math] \theta_i [/math]. Поэтому чтобы построить [math]q_j [/math], предлагается на каждом шаге следить за оценками погрешностей [math]\beta_{t}|v_i(t)|, i = 1 \dots t, t = j - 1 [/math], где [math]v_i(t) [/math] - [math]t[/math]-я компонента собственного вектора [math]v_i [/math]. И когда какая-то оценка становится слишком малой, проводить ортогонализацию вектора Ланцоша [math]z [/math]. Величина [math]\beta_{t}|v_i(t)| [/math] считается малой, если она меньше, чем [math]\sqrt{\varepsilon}||T_{t}|| [/math], где [math]\varepsilon[/math] - доступная машинная точность чисел.

После следует вычисление собственных значений [math] \theta_j [/math] и собственных векторов [math]v_j [/math] полученной трехдиагональной матрицы [math]T_j[/math], например, с помощью метода "разделяй и властвуй"

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

[math] q_{1} = b_{j}/\|b\|_2, \beta_0=0, q_0=0[/math]
[math] for\, j=1\,\, to\, \, k[/math]
    [math]z=Aq_j,[/math]
    [math]\alpha_j=q_j^Tz,[/math]
    [math]z=z-\alpha_jq_j-\beta_{j-1}q_{j-1},[/math]
    Провести выборочную ортогонализацию по отношению
    к сошедшимся векторам Ритца
    [math]for\, i \leqslant k,[/math] таких, что [math]\, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{k}\| \, [/math]
    [math]z = z-(y^T_{i,k},z)y_{i,k}[/math]
    [math]end \, for[/math]
    [math]\beta_{j}=\|z\|_2[/math]
    [math]q_{j+1}=z/\beta_{j}, [/math]
    Вычислить собственные значения и собственные векторы
    [math]\, \, T_{j} \, \, [/math]и оценки погрешности в них 
[math]end \, for[/math]

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

• Умножение матрицы на вектор, [math]z=Aq_j, [/math].

• Ортогонализация по отношению к сошедшимся векторам [math]z = z-(y^T_{i,k},z)y_{i,k} [/math] для [math]i = 1 \dots t[/math]

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

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

• Умножение матрицы на вектор. Состоит из операций умножения вектора на число и сложения векторов.

[math]z=Aq_j,[/math].

• Вычисление вектора [math]q_{j+1}[/math] как скалярное произведение двух других векторов.

[math]\alpha_j=q_j^Tz,[/math]

[math]z=z-\alpha_jq_j-\beta_{j-1}q_{j-1},[/math]

[math]q_{j+1}=z/\|z\|_2[/math]

• Ортогонализация вектора Ланцоша с помощью скалярного произведения:

[math]z = z-(y^T_{i,t},z)y_{i,t}[/math]

• Алгоритм "разделяй и властвуй" для вычисления собственных значений матрицы.

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

Алгоритм выполняется в следующей последовательности:

[math]1.\, \, \beta_0 = 0,\; q_0 = 0[/math] #Инициализируем векторы

[math]2.\, \, \|b\|_2 = \sqrt{\sum\limits_{j=1}^{n} b_j^2}[/math] #Вычисляем норму вектора начального приближения.

[math]3.\, \, q_{1_{j}} = \frac{b_{j}}{\|b\|_2}, \; j = 1,\, \dots\,, n[/math] #Нормализуем вектор начального приближения.

[math]4.\, \,[/math] Начинается цикл, в котором не более чем [math]k[/math]. На [math]i[/math]-й итерации производятся следующие вычисления:

[math]4.1.\, \, z_j = \sum\limits_{m=1}^{n} a_{jm} q_{i_m}, \; j = 1,\,\dots\,, n[/math] #Считаем результат применения линейного оператора [math]A[/math] к вектору [math]q_i[/math].

[math]4.2.\, \, \alpha_i = \sum\limits_{j=1}^{n}q_{i_j} z_j[/math] #Получаем результат скалярного произведения векторов [math]q_i[/math] и [math]z[/math].

[math]4.3.\, \, z_j = z_j - \alpha_i q_{i_j} - \beta_{i-1}q_{i-1_j}, \, j = 1,\,\dots\,, n[/math] #Вычисляем линейную комбинацию векторов.

[math]4.4.\, \,[/math] Начинается цикл, в котором производится выборочная ортогонализация векторов Ланцоша

[math]4.4.1.\, \,[/math] Вычисляется вектор Ритца [math]y_{i,k_d} = \sum\limits_{x=1}^{k} q_{x_d} v_{i_x}, \; d = 1,\,\dots\,, n,[/math], где  [math]y_{i,k_d}[/math] -  [math]d[/math]-я компонента вектора Ритца [math]y,[/math]
[math]4.4.2.\, \,[/math] Производится скалярное умножение [math]y_{i,k}^T[/math]  и [math] z: \; \; \gamma = \sum\limits_{d=1}^{n}y_{i,k_d} z_d,[/math]
[math]4.4.3.\, \,[/math] Ищется новое собственное значение [math] \theta_j [/math] и собственный вектор [math] v_j [/math] для полученной матрицы [math] T_j,[/math]

[math]4.5.\, \, \beta_i = \|z\|_2 = \sqrt{\sum\limits_{j=1}^{n} z_j^2}[/math] #Считаем норму вектора [math]z[/math].

[math]4.6.[/math] Проверка равенства [math]\beta_i == 0[/math] # Если норма оказалась равной нулю, то завершаем итерации и переходим к вычислению собственных векторов и собственных значений полученной матрицы. В обратном случае, продолжаем выполнения итераций.

[math]4.7.\, \, q_{i+1_j} = \frac{z_j}{\beta_i}, \; j = 1,\, \dots \,, n[/math] #Нормируем вектор [math]z[/math].

[math]4.8.\,[/math] Если выполнили [math]k[/math] итераций, то завершаем выполнение итераций, переходим к следующему шагу. Иначе начинаем последующую итерацию цикла.

[math]5.[/math] Вычисляем собственные значения и собственные вектора полученной матрицы [math]T_k[/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]

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

Здесь:
[math]Ax[/math]  — перемножение матрицы и вектора,
[math]\|x\|[/math]  — макрооперация взятия нормы,
[math]\alpha x[/math]  — умножение скаляра на вектор,
[math]\cdot [/math]  — операция скалярного умножения векторов,
[math]\mathsf{F}[/math]  — операция [math]z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i[/math],
[math]\mathsf{f}[/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]

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

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

Исследование алгоритма на масштабируемость производилось на суперкомпьютере "Ломоносов" с использованием технологии MPI. Для тестирования была выбрана вот эта реализация алгоритма.

Набор параметров, которые изменялись при тестировании:

• Число процессоров от 16 до 128.
• Размерность матрицы от 5000 до 50000 с шагом 5000.

Табл. 1. Результаты тестирования. В ячейках таблицы -- время работы алгоритма в секундах.

		Число процессоров
		16	32	64	128
Размер матрицы	5000	1,07	0,71	0,62	0,54
	10000	2,39	1,94	1,17	0,79
	15000	4,35	3,38	1,89	1,59
	20000	7,63	5,19	3,78	2,35
	15000	12,41	7,33	4,32	3,44
	30000	16,68	9,36	5,02	4,68
	35000	22,58	11,92	8,04	6,20
	40000	29,31	15,90	8,47	7,43
	45000	37,04	19,42	10,55	9,04
	50000	45,64	23,39	16,09	10,25

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

Реализация в проекте Vienna VL

Реализация в проекте IETL

Реализация в проекте ARPACK

Lanczos Plus Plus

ED Lanczos

Paralleles Rechnen 2

Cuda-Arnoldi