Участник:Grinch96/QR-факторизация методом Грама-Шмидта с последующей реортогонализацией
| QR-факторизация | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]\displaystyle{ O(n^3) }[/math] |
| Объём входных данных | [math]\displaystyle{ n^2 }[/math] |
| Объём выходных данных | [math]\displaystyle{ 2\cdot n^2 }[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]\displaystyle{ O(n) }[/math] |
| Ширина ярусно-параллельной формы | [math]\displaystyle{ O(n^2) }[/math] |
Автор: Г. А. Балыбердин
Свойства и структура алгоритма
Общее описание алгоритма
Пусть [math]\displaystyle{ A = (a_1, ..., a_n) }[/math] - вещественная матрица [math]\displaystyle{ n × n }[/math], определитель которой не равен [math]\displaystyle{ 0 }[/math]. Во многих приложениях требуется решать линейную систему [math]\displaystyle{ Ax=b }[/math] c плохо-обусловленной матрицей [math]\displaystyle{ A }[/math]. При решении данной задачи, чтобы не увеличивать число обусловленности матрицы [math]\displaystyle{ A }[/math], ее можно представить в виде [math]\displaystyle{ A=QR }[/math], где матрица [math]\displaystyle{ Q \in \mathbb{R}^{n×n} }[/math] состоит из ортонормированных строк, а матрица [math]\displaystyle{ R \in \mathbb{R}^{n×n} }[/math] является верхнетреугольной. В итоге мы получаем так называемую [math]\displaystyle{ QR }[/math]-факторизацию.
Чтобы построить [math]\displaystyle{ QR }[/math]-факторизацию можно воспользоваться процессом ортогонализации Грама-Шмидта[1], однако в условиях машинной арифметики, матрица [math]\displaystyle{ Q }[/math] может получиться далекой от ортогональной.Чтобы этого избежать, на определенных итерациях нужно проводить процесс реортогонализации [2].
Математическое описание алгоритма
Исходные данные: квадратная матрица [math]\displaystyle{ A }[/math] порядка [math]\displaystyle{ n }[/math] (элементы [math]\displaystyle{ a_{ij} }[/math]), определитель которой не равен [math]\displaystyle{ 0 }[/math].
Вычисляемые данные: верхнетреугольная матрица [math]\displaystyle{ R }[/math] порядка [math]\displaystyle{ n }[/math] (элементы [math]\displaystyle{ r_{ij} }[/math]), унитарная матрица [math]\displaystyle{ Q }[/math] порядка [math]\displaystyle{ n }[/math] (элементы [math]\displaystyle{ q_{ij} }[/math]); выполняется соотношение [math]\displaystyle{ QR=A }[/math].
Формулы процесса ортогонализации:
[math]\displaystyle{ \begin{align} & p_{j} = a_{j} - \sum_{i=1}^{j-1} q_{i} (q_{i}^{T} a_{j}), \\ & q_{j} = p_{j} / ||p_{j}||_{2}, \\ & r_{ij} = q_{i}^{T} a_{j},\quad i = 1, \ldots , j - 1.\\ \end{align} }[/math]
Здесь [math]\displaystyle{ q_{i}, \, a_{i} }[/math] столбцы матриц [math]\displaystyle{ Q, \, A }[/math] соответственно.
На [math]\displaystyle{ m }[/math]-ом шаге получаем [math]\displaystyle{ p_{m} \in \mathbb{R}^{n} }[/math]. Если [math]\displaystyle{ ||a_{m}||_{2}/||p_{m}||_{2} \gt k }[/math], запускается процесс реортогонализации. [math]\displaystyle{ k }[/math] произвольная и задает точность ортогональности матрицы [math]\displaystyle{ Q }[/math]. В посвященной методу реортогонализации статье [3] [math]\displaystyle{ k }[/math] рекомендуется брать равным [math]\displaystyle{ 10 }[/math], однако приводятся примеры алгоритмов, где [math]\displaystyle{ k=\sqrt{2} }[/math] или [math]\displaystyle{ k=\sqrt{5} }[/math]. Мы исследуем [math]\displaystyle{ k=10 }[/math].
Суть процесса реортогонализации - это процесс ортогонализации Грама-Шмидта, примененный повторно к вектору [math]\displaystyle{ p_{m} }[/math] на [math]\displaystyle{ m }[/math]-оm шаге алгоритма. Более формально:
[math]\displaystyle{ \begin{align} & \tilde{p}_{m} = p_{m} - \sum_{i=1}^{m-1} q_{i} (q_{i}^{T} p_{m}), \\ & q_{m} = \tilde{p}_{m} / ||\tilde{p}_{m}||_{2}, \\ & \tilde{r}_{mj} = q_{j}^{T} p_{m}, \quad j = 1, \ldots , m - 1, \\ & r_{mj} = \tilde{r}_{mj} + r_{mj}^{0}, \quad j = 1, \ldots , m - 1, \\ & r_{mm} = ||\tilde{p}_{m}||_{2}. \end{align} }[/math]
Здесь [math]\displaystyle{ \tilde{q}_{i} }[/math] столбцы матриц [math]\displaystyle{ \tilde{Q}_{m} }[/math], [math]\displaystyle{ r_{mj}^{0} }[/math] - элементы матрицы [math]\displaystyle{ R }[/math], вычисленные до реортогонализации.
Вычислительное ядро алгоритма
Вычислительное ядро алгоритма, если не требуется реортогонализация, состоит из вычисления:
- [math]\displaystyle{ \dfrac{n(n+1)}{2} }[/math] скалярных произведений векторов, включая вычисление длины вектора;
- [math]\displaystyle{ \dfrac{(n-1)n}{2} }[/math] умножений векторов на число.
На каждом шаге реортогонализации добавляется [math]\displaystyle{ j-1 }[/math] вычислений скалярных произведений векторов на число и столько же умножений вектора на число. Здесь реортогонализация производится на [math]\displaystyle{ j }[/math]-ом шаге. Однако сколько раз алгоритм будет реортогонализовывать вектора зависит от свойств матрицы.
Макроструктура алгоритма
Как записано и в описании ядра алгоритма, основную часть процесса составляют скалярные произведения векторов и умножения вектора на число:
[math]\displaystyle{ \begin{align} & p_{j} = a_{j} - \sum_{i=1}^{j-1} q_{i} (q_{i}^{T} a_{j}), \\ & q_{j} = p_{j} / ||p_{j}||_{2}. \end{align} }[/math]
как если алгоритм реортогонализует повторно вектора, так и без реортогонализаций.
Схема реализации последовательного алгоритма
Последовательность исполнения процесса следующая:
1. [math]\displaystyle{ q_1 = \frac{a_1}{||a_1||_2} }[/math]
Далее для всех [math]\displaystyle{ i = 2, n }[/math]:
2. [math]\displaystyle{ r_{ij} = q_{j}^{T} a_{i},\quad j = 1, \ldots , i - 1 }[/math]
3. [math]\displaystyle{ p_{i} = a_{i} - \sum_{j=1}^{i-1} q_{j} r_{ij} }[/math]
Если [math]\displaystyle{ ||a_{i}||_{2}/||p_{i}||_{2} \gt k }[/math], выполнить 4, 5, 6. Иначе перейти сразу к 7:
4. [math]\displaystyle{ r_{ij}^{0} = q_{j}^{T} p_{i},\quad j = 1, \ldots , i - 1 }[/math]
5. [math]\displaystyle{ p_{i} = p_{i} - \sum_{j=1}^{i-1} q_{j} r_{ij}^{0}, \quad j = 1, \ldots , i - 1 }[/math]
6. [math]\displaystyle{ r_{ij} = r_{ij}^{0} + r_{ij},\quad j = 1, \ldots , i - 1 }[/math]
7. [math]\displaystyle{ r_{ii} = ||p_{i}||_{2} }[/math]
8. [math]\displaystyle{ q_{i} = p_{i} / r_{ii} }[/math]
Здесь (как и везде) под второй нормой подразумевается:
[math]\displaystyle{ ||x||_2 = \sqrt{x_1^2 + x_2^{2} + \ldots + x_n^2} }[/math]
Последовательная сложность алгоритма
Для факторизации матрицы порядка n процессом Грама-Шмидта в последовательном варианте без реортогонализации требуется:
- [math]\displaystyle{ n }[/math] извлечений квадратного корня,
- [math]\displaystyle{ n }[/math] делений,
- [math]\displaystyle{ (n-1)(n^2+2n-1) }[/math] сложений (вычитаний),
- [math]\displaystyle{ n(n^2+2n-2) }[/math] умножений.
Если реортогонализация работает на каждом шаге (чего никогда не случается, однако предположим худший случай), то количество всех операций удвоится. Таким образом, умножения и сложения (вычитания) составляют основную часть алгоритма.
При классификации по последовательной сложности, процесс ортогонализации Грама-Шмидта с реортогонализацией относится к алгоритмам с кубической сложностью.
Информационный граф
Опишем вершины информационного графа.
1. [math]\displaystyle{ a_i, \; i = 1 \ldots n }[/math] — столбцы исходной матрицы [math]\displaystyle{ A }[/math]; подаются на вход программе.
2. [math]\displaystyle{ N }[/math] — операция взятия нормы вектора.
3. [math]\displaystyle{ div }[/math] — операция деления вектора на число.
4. [math]\displaystyle{ Sc }[/math] — операция взятия скалярного произведения двух векторов.
5. [math]\displaystyle{ F }[/math] — операция вычитания из вектора произведение числа на вектор.
6. [math]\displaystyle{ + }[/math] — операция суммирования двух чисел.
7. [math]\displaystyle{ r_{ij} }[/math], [math]\displaystyle{ q_i, \; i = 1 \ldots n, \; j = 1 \ldots n }[/math] — элементы матрицы [math]\displaystyle{ R }[/math] и столбцы матрицы [math]\displaystyle{ Q }[/math] соответственно; выходные данные программы.
Информационный граф построен для матрицы размерности [math]\displaystyle{ 3 }[/math] при выполнении реортогонализации для третьего вектора.
Ресурс параллелизма алгоритма
Входные и выходные данные алгоритма
Свойства алгоритма
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
Локальность данных и вычислений
Возможные способы и особенности параллельной реализации алгоритма
Масштабируемость алгоритма и его реализации
Динамические характеристики и эффективность реализации алгоритма
Выводы для классов архитектур
Существующие реализации алгоритма
Литература
<references \>
[[Категория:]] [[Категория:]]
- ↑ Тыртышников Е.Е. Методы численного анализа // М.: 2006. 83 с.
- ↑ Luc Giraud, Julien Langou, Miroslav Rozložník, Jasper van den Eshof. Rounding error analysis of the classical Gram-Schmidt orthogonalization process // Springer-Verlag, 2005.
- ↑ Luc Giraud and Jilien Langou. A robust criterion for the modified Gram–Schmidt algorithm with selective reorthogonalization // Society for Industrial and Applied Mathematic, 2003.