Участник:Egor Luckjanow 403/Алгоритм Грамма-Шмидта: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
(Автора добавили, ибо мерзко как-то, когда у статьи автора нету...)
Строка 12: Строка 12:
 
<math>...</math><br>
 
<math>...</math><br>
 
Мы применяем этот метод к столбцам исходной матрицы A и получаем разложение A = QR, где Q составлена из столбцов q_1,...,q_n, полученных в результате применения метода Грама-Шмидта, а R -- матрица коэффициентов, которая будет верхнетреугольной в силу вида метода.
 
Мы применяем этот метод к столбцам исходной матрицы A и получаем разложение A = QR, где Q составлена из столбцов q_1,...,q_n, полученных в результате применения метода Грама-Шмидта, а R -- матрица коэффициентов, которая будет верхнетреугольной в силу вида метода.
 +
=== Вычислительное ядро алгоритма ===
 +
Ядро алгоритма -- последовательное нахождение Q и R, как в математическом описании.
 +
=== Макроструктура алгоритма ===
 +
Макрооперации: скалярное произведение ddot; умножение матрицы на вектор с прибавлением другого вектора daxpy.
 +
=== Схема реализации последовательного алгоритма ===
 +
<source lang="cpp">
 +
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;
 +
}
 +
}
 +
}
 +
</source>
 +
=== Последовательная сложность алгоритма ===
 +
Сложность последовательного алгоритма -- <math>O(n^3)</math>.
 +
=== Информационный граф ===
 +
=== Ресурс параллелизма алгоритма ===
 +
Параллельная сложность (в терминах операций сложения и умножения) -- <math>O(n^2)</math>.
 +
=== Входные и выходные данные алгоритма ===
 +
Входные данные -- размер матрицы и сама матрица.
 +
Выходные данные -- матрицы <math>Q</math> и <math>R</math> из QR-разложения поданной на вход матрицы.
 +
=== Свойства алгоритма ===
 +
Алгоритм Грама-Шмидта обладает свойством потери ортогональности. Поэтому для корректной его реализации необходимо время от времени делать реортогонализацию, а также обеспечить максимально возможную точность вычислений с плавающей точкой.
 +
Можно реализовать алгоритм для работы с рациональными матрицами. Тогда проблема с потерей ортогонализации пропадает (в виду идеальной точности вычислений в поле частных колец вида <math>\mathbb{Z}/<m></math> на компьютере).
 +
== Программная реализация алгоритма ==
 +
=== Особенности реализации последовательного алгоритма ===
 +
Проблема с потерей ортогональности уже описана в свойствах алгоритма.
 +
Следует также отметить большое количество пересылок.
 +
=== Локальность данных и вычислений ===
 +
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
Метод Грама-Шмидта плохо приспособлен для параллельной реализации. Единственный способ, известный автору -- отдать каждому процессу по группе столбцов <math>q_i</math> для вычисления.
 +
Число пересылок тогда составляет <math>n(n+1)/2</math>.
 
== Литература ==
 
== Литература ==
 
#Е. Е. Тыртышников, Методы численного анализа
 
#Е. Е. Тыртышников, Методы численного анализа

Версия 13:11, 6 декабря 2021

Автор: Егор Лукьянов, студент 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.7 Информационный граф

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.2 Локальность данных и вычислений

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

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

3 Литература

  1. Е. Е. Тыртышников, Методы численного анализа
  2. Ф. Р. Гантмахер, Теория матриц
  3. С. Ленг, Алгебра
  4. Е. Е. Тыртышников, Основы алгебры
Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Egor_Luckjanow_403/Алгоритм_Грамма-Шмидта&oldid=28722