Участник:Egor Luckjanow 403/Алгоритм Грамма-Шмидта

Автор: Егор Лукьянов, студент 403 группы ВМК МГУ

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

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

Алгоритм Грама-Шмидта предназначен для построения QR-разложения матрицы.

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

В основе работы алгоритма лежит метод Грама-Шмидта для построения ортонормированного базиса на основе уже существующего:
[math]q_1 = \frac{a_1}{||a_1||}[/math]
[math]p_2 = a_2 - (a_2, q_1)q_1[/math]
[math]q_2 = \frac{p_2}{||p_2||}[/math]
[math]p_3 = a_3 - (a_3, q_1)q_1 - (a_3, q_2)q_2[/math]
[math]q_3 = \frac{p_3}{||p_3||}[/math]
[math]...[/math]
Мы применяем этот метод к столбцам исходной матрицы A и получаем разложение A = QR, где Q составлена из столбцов q_1,...,q_n, полученных в результате применения метода Грама-Шмидта, а R -- матрица коэффициентов, которая будет верхнетреугольной в силу вида метода.

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

Ядро алгоритма -- последовательное нахождение Q и R, как в математическом описании.

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

Макрооперации: скалярное произведение ddot; умножение матрицы на вектор с прибавлением другого вектора daxpy.

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

void QR(int N, double* A, double* Q, double* R)
{
	int one = 1;
	for (int i = 0; i < N; ++i) {
		R[i * N + i] = 1;
		//int one = 1;
		for (int k = 1; k < i + 1; ++k) {
			R[i - k + i * N] = ddot(&N, A + i * N, &one, Q + (i - k) * N, &one) / ddot(&N, Q + (i - k) * N, &one, Q + (i - k) * N, &one);
		}
		for (int j = 0; j < N; ++j) {
			Q[j + i * N] = A[j + i * N];
			/*for (int k = 0; k < i; ++k) {
				Q[j + i * N] -= R[i - k - 1 + i * N] * Q[j + (i - k - 1) * N];
			}*/
		}
		for (int k = 0; k < i; ++k) {
			double tmp = -R[i - k - 1 + i * N];
			daxpy(&N, &tmp, Q + (i - k - 1) * N, &one, Q + i * N, &one);
		}
		/*
		double norm = dnrm2(&N, Q + i * N, &one);
		for (int j = 0; j < N; ++j) {
			Q[j + i * N] /= norm;
		}
		for (int k = 0; k < i + 1; ++k) {
			R[k + i * N] *= norm;
		}*/
		for (int k = i + 1; k < N; ++k) {
			R[k + i * N] = 0;
		}
	}
	for (int i = 0; i < N; ++i) {
		double norm = dnrm2(&N, Q + i * N, &one);
		for (int j = 0; j < N; ++j) {
			Q[j + i * N] /= norm;
		}
		for (int j = 0; j < N; ++j) {
			R[i + j * N] *= norm;
		}
	}
}

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

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

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

Параллельная сложность (в терминах операций сложения и умножения) -- [math]O(n^2)[/math].

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

Входные данные -- размер матрицы и сама матрица. Выходные данные -- матрицы [math]Q[/math] и [math]R[/math] из QR-разложения поданной на вход матрицы.

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

Алгоритм Грамма-Шмидта обладает свойством потери ортогональности. Поэтому для корректной его реализации необходимо время от времени делать реортогонализацию, а также обеспечить максимально возможную точность вычислений с плавающей точкой. Можно реализовать алгоритм для работы с рациональными матрицами. Тогда проблема с потерей ортогонализации пропадает (в виду идеальной точности вычислений в поле частных колец вида [math]\mathbb{Z}/\lt m\gt [/math] на компьютере).

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

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

Проблема с потерей ортогональности уже описана в свойствах алгоритма. Следует также отметить большое количество пересылок.

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

Метод Грамма-Шмидта плохо приспособлен для параллельной реализации. Единственный способ, известный автору -- отдать каждому процессу по группе столбцов [math]q_i[/math] для вычисления. Число пересылок тогда составляет [math]n(n+1)/2[/math].

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

Алгоритм увеличивает свою работу при добавлении независимых вычислительных ядер. Алгоритм обладает сильной масштабируемостью.

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

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

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

Реализация архитектурно независима. Зависимы только макрооперации, которые предполагаются уже реализованными.

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

В силу потери ортогональности алгоритм Грамма-Шмидта пользуется небольшой популярностью в сравнении с альтернативными методами QR-разложения, так что популярных параллельных реализаций нет. Альтернативный (не Грамм-Шмидт) метод представлен, например, в пакете LAPACK и является одним из лучших по производительности.

3 Литература

1. Е. Е. Тыртышников, Методы численного анализа
2. Ф. Р. Гантмахер, Теория матриц
3. С. Ленг, Алгебра
4. Е. Е. Тыртышников, Основы алгебры