Участница:Larisayaryzheva/Решение СЛАУ методом сопряжённых градиентов

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

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

Метод сопряжённых градиентов — это итерационный численный метод, предназначенный для решения систем линейных алгебраических уравнений (СЛАУ) вида [math]\displaystyle{ Ax = b }[/math], где матрица [math]\displaystyle{ A }[/math] является положительно определённой и симметричной . Основная идея алгоритма заключается в минимизации квадратичного функционала, соответствующего данной системе.

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

Необходимо найти решение СЛАУ [math]\displaystyle{ Ax = b }[/math], где матрица [math]\displaystyle{ A }[/math] является положительно определённой и симметричной ([math]\displaystyle{ A^T = A \gt 0 }[/math]). Данная задача эквивалентна минимизации следующего функционала [math]\displaystyle{ f(x) = \frac{1}{2} \langle A x, x \rangle - \langle b, x \rangle }[/math].

Классически данная задача решается с помощью итерационного метода спуска: [math]\displaystyle{ x_{k+1} = x_{k} + α_{k}d_{k} }[/math].

На каждой итерации выбирается направление спуска [math]\displaystyle{ d_{k} }[/math] и шаг [math]\displaystyle{ α_{k} }[/math], где [math]\displaystyle{ d_{k} }[/math] – антиградиент [math]\displaystyle{ f(x) }[/math], а [math]\displaystyle{ \alpha_k = \arg\min_{\alpha \ge 0} f(x_k + \alpha_k d_k) = -\frac{\langle g_k, d_k \rangle}{\langle A d_k, d_k \rangle}, \quad g_k := A x_k - b. }[/math]

Сопряжённые направления

Главная особенность метода сопряжённых градиентов заключается в выборе направлений [math]\displaystyle{ d_{k} }[/math], которые являются A-сопряжёнными, то есть удовлетворяют условию: [math]\displaystyle{ \langle A d_{i}, d_{j} \rangle = 0 }[/math] для всех [math]\displaystyle{ i ≠ j }[/math].

Общая схема метода сопряжённых градиентов

1. Пусть даны [math]\displaystyle{ d_{0},...,d_{k}: \langle A d_{i}, d_{j} \rangle = 0 }[/math] для [math]\displaystyle{ i ≠ j }[/math]

2. [math]\displaystyle{ d_{k+1} = -g_{k+1} + β_{k}d_{k} }[/math], где [math]\displaystyle{ g_{k+1} = A x_k - b }[/math]

3. [math]\displaystyle{ β_{k} }[/math] ищем исходя из условий на сопряжённость направлений: [math]\displaystyle{ 0 = \langle A d_k, d_{k+1} \rangle= -\,\langle A d_k, g_{k+1} \rangle + \beta_k \langle A d_k, d_k \rangle }[/math]

Отсюда [math]\displaystyle{ β_k = \frac{\langle A d_k, g_{k+1} \rangle}{\langle A d_k, d_k \rangle} }[/math]

Инициализация: [math]\displaystyle{ g_0 = A x_0 - b }[/math] [math]\displaystyle{ d_0 = -g_0 }[/math] [math]\displaystyle{ k = 0 }[/math] Цикл: Пока [math]\displaystyle{ \|g_k\| \gt \varepsilon \|b\| }[/math] выполнять: [math]\displaystyle{ \quad 1.\; }[/math] Вычисляем шаг [math]\displaystyle{ \alpha_k = -\dfrac{\langle g_k, d_k \rangle}{\langle A d_k, d_k \rangle} }[/math] [math]\displaystyle{ \quad 2.\; }[/math] Обновляем приближение решения [math]\displaystyle{ x_{k+1} = x_k + \alpha_k d_k }[/math] [math]\displaystyle{ \quad 3.\; }[/math] Обновляем невязку [math]\displaystyle{ g_{k+1} = A x_{k+1} - b }[/math] [math]\displaystyle{ \quad 4.\; }[/math] Вычисляем коэффициент [math]\displaystyle{ \beta_k = \dfrac{\langle g_{k+1}, A d_k \rangle}{\langle A d_k, d_k \rangle} }[/math] [math]\displaystyle{ \quad 5.\; }[/math] Обновляем направление [math]\displaystyle{ d_{k+1} = -g_{k+1} + \beta_k d_k }[/math] [math]\displaystyle{ \quad 6.\; }[/math] Увеличиваем счётчик [math]\displaystyle{ k = k + 1 }[/math] Результат: [math]\displaystyle{ x_k }[/math].

Оптимизация метода сопряжённых градиентов

Выше приведённая схема метода сопряжённых градиентов включает в себя вычисление двух матрично-векторных произведений за один прогон цикла. Данную неэффективность можно устранить, используя следующие факты:

[math]\displaystyle{ 1. }[/math] Метод сопряжённых градиентов построен таким образом, что [math]\displaystyle{ \langle g_{i}, g_{k+1} \rangle = 0 }[/math] и [math]\displaystyle{ \langle g_{k+1}, d_{i} \rangle }[/math] для [math]\displaystyle{ i \lt = k }[/math] (это можно проверить по индукции и использовать свойство сопряжённости направлений)

[math]\displaystyle{ 2. }[/math] [math]\displaystyle{ g_{k+1} = Ax_{k+1} - b = A(x_{k} + α_{k}d_{k}) - b = g_{k} + α_{k}Ad_{k} }[/math]

Отсюда [math]\displaystyle{ Ad_{k} = α_{k}^{-1}(g_{k+1} - g_{k}) }[/math]

[math]\displaystyle{ 3. }[/math] Значит
[math]\displaystyle{ \beta_k = \dfrac{\langle g_{k+1} - g_k,\; g_{k+1} \rangle}{\langle g_{k+1} - g_k,\; d_k \rangle}= \dfrac{\| g_{k+1} \|^2}{-\langle g_k,\; d_k \rangle} }[/math]
[math]\displaystyle{ = \dfrac{\| g_{k+1} \|^2}{-\langle -g_k + \beta_{k-1} d_{k-1},\; g_k \rangle}= \dfrac{\| g_{k+1} \|^2}{\| g_k \|^2} }[/math]
[math]\displaystyle{ 4. }[/math] Аналогично [math]\displaystyle{ \alpha_k = \dfrac{\| g_k \|^2}{\langle A d_k,\; d_k \rangle} }[/math]

Общая схема эффективного метода сопряжённых градиентов

Теперь можем записать эффективную версию метода сопряжённых градиентов, использующий всего одно матрично-векторное произведение.

Инициализация: [math]\displaystyle{ g_0 = A x_0 - b }[/math] [math]\displaystyle{ d_0 = -g_0 }[/math] [math]\displaystyle{ k = 0 }[/math] Цикл: Пока [math]\displaystyle{ \|g_k\| \gt \varepsilon \|b\| }[/math] выполнять: [math]\displaystyle{ \quad 1.\; }[/math] Вычисляем шаг [math]\displaystyle{ \alpha_k = \dfrac{\| g_k \|^2}{\langle A d_k,\; d_k \rangle} }[/math] [math]\displaystyle{ \quad 2.\; }[/math] Обновляем приближение решения [math]\displaystyle{ x_{k+1} = x_k + \alpha_k d_k }[/math] [math]\displaystyle{ \quad 3.\; }[/math] Обновляем невязку [math]\displaystyle{ g_{k+1} = g_{k} + α_{k}Ad_{k} }[/math] [math]\displaystyle{ \quad 4.\; }[/math] Вычисляем коэффициент [math]\displaystyle{ \beta_k = \dfrac{\| g_{k+1} \|^2}{\| g_k \|^2} }[/math] [math]\displaystyle{ \quad 5.\; }[/math] Обновляем направление [math]\displaystyle{ d_{k+1} = -g_{k+1} + \beta_k d_k }[/math] [math]\displaystyle{ \quad 6.\; }[/math] Увеличиваем счётчик [math]\displaystyle{ k = k + 1 }[/math] Результат: [math]\displaystyle{ x_k }[/math].

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

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

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

Алгоритм метода сопряжённых градиентов представляет собой итерационный процесс, включающий следующие макрооперации:

1. Вычисление невязки текущего приближения

2. Вычисление матрично-векторного произведения

3. Вычисление скалярных произведений и норм векторов

4. Обновление текущего приближения и направлений поиска

5. Проверка условия останова

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

Подробная схема реализации метода сопряжённых градиента и эффективного метода сопряженных градиентов представлена в соответствующих разделах Общая схема метода сопряжённых градиентов и Общая схема эффективного метода сопряжённых градиентов.

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

На каждой итерации метода сопряжённых градиентов выполняются линейно-алгебраические операции. Их стоимостная оценка даёт итоговую асимптотическую сложность итерации.

• Сложение двух векторов размерности [math]\displaystyle{ n }[/math] требует [math]\displaystyle{ n }[/math] операций сложения: сложность [math]\displaystyle{ O(n) }[/math].
• Скалярное произведение двух векторов требует [math]\displaystyle{ n }[/math] умножений и [math]\displaystyle{ n }[/math] сложений: сложность [math]\displaystyle{ O(n) }[/math].
• Вычисление нормы вектора [math]\displaystyle{ \|a\| }[/math] эквивалентно вычислению скалярного произведения [math]\displaystyle{ \langle a, a\rangle }[/math] с последующим извлечением корня: сложность [math]\displaystyle{ O(n) }[/math].
• Матрично-векторное произведение [math]\displaystyle{ A x }[/math] для плотной матрицы размера [math]\displaystyle{ n \times n }[/math] требует [math]\displaystyle{ n^2 }[/math] умножений и [math]\displaystyle{ n^2 }[/math] сложений: сложность [math]\displaystyle{ O(n^2) }[/math].

Так как на каждой итерации алгоритма выполняется ровно одно матрично-векторное произведение [math]\displaystyle{ A d_k }[/math], а остальные операции имеют линейную сложность [math]\displaystyle{ O(n) }[/math], общая стоимость одной итерации определяется наиболее дорогой операцией: [math]\displaystyle{ O(n^2). }[/math] Следовательно, последовательная сложность метода при [math]\displaystyle{ k }[/math] итерациях составляет: [math]\displaystyle{ O(k n^2) }[/math]

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

Информационный граф метода сопряжённых градиентов отражает последовательность основных макроопераций внутри одной итерации метода и зависимость данных между ними. Вершины графа соответствуют вычислению матрично-векторного произведения [math]\displaystyle{ q = Ad_k }[/math], вычислению коэффициентов [math]\displaystyle{ α_k }[/math] и [math]\displaystyle{ β_k }[/math], обновлению решения [math]\displaystyle{ x_{k+1} }[/math], невязки [math]\displaystyle{ g_{k+1} }[/math] и направления [math]\displaystyle{ d_{k+1} }[/math]. Граф имеет ярусно-параллельную структуру: отдельные этапы должны выполняться последовательно, однако операции обновления векторов полностью параллелизуемы по компонентам. Итерации метода образуют цикл: после вычисления [math]\displaystyle{ d_{k+1} }[/math] выполняется проверка критерия остановки, и при необходимости управление возвращается к вычислению нового матричного произведения. Ошибка создания миниатюры: Файл не найден

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

Метод сопряжённых градиентов обладает параллельной структурой, так как его вычислительное ядро состоит из операций, которые могут выполняться независимо на элементах векторов или столбцах матрицы. Параллельная сложность алгоритма определяется шириной ярусов его информационного графа - множеством вершин, не имеющих зависимостей между собой. За счёт координатного параллелизма каждый компонент вектора вычисляется независимо, поэтому при размере задачи [math]\displaystyle{ n }[/math] в один момент времени может быть выполнено [math]\displaystyle{ O(n) }[/math] операций, связанных с обновлениями векторов. Наиболее трудоемкой задачей в методе сопряжённых векторов является матричное умножение [math]\displaystyle{ Ad_{k} }[/math]. Данная операция обладает хорошим покоординатным параллелизмом и может быть ускорена. Другим источником параллелизма являются редукции, возникающие при вычислении скалярных произведений и норм векторов. Эти операции имеют сложность [math]\displaystyle{ O(n) }[/math], но их параллельное время составляет [math]\displaystyle{ O(\log n) }[/math]. Последовательные зависимости между итерациями формируют скошенный параллелизм и не позволяют выполнять итерации независимо.

В итоге общая параллельная сложность метода равна [math]\displaystyle{ O(k\log n) }[/math], где k – число операций, а максимальный потенциальный ресурс параллелизма на одной итерации составляет [math]\displaystyle{ O(n) }[/math].

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

Входными данными метода сопряжённых градиентов являются матрица системы [math]\displaystyle{ A \in \mathbb{R}^{n \times n} }[/math], вектор правой части [math]\displaystyle{ b \in \mathbb{R}^n }[/math] и начальное приближение решения [math]\displaystyle{ x_0 }[/math]. Матрица [math]\displaystyle{ A }[/math] должна быть симметричной и положительно определённой, что гарантирует сходимость метода и корректность вычисления сопряжённых направлений. В практических задачах матрица может быть как плотной, так и разреженной. Структурные свойства [math]\displaystyle{ A }[/math] (диагональное преобладание, блочная структура) могут дополнительно уменьшать вычислительные затраты и повышать эффективность параллельной реализации.

Вектор правой части [math]\displaystyle{ b }[/math] и начальное приближение [math]\displaystyle{ x_0 }[/math] являются плотными векторами размерности [math]\displaystyle{ n }[/math]. Начальное приближение может быть выбрано произвольно. На вход также подаётся параметр остановки: допустимая погрешность [math]\displaystyle{ \varepsilon }[/math].

Выходными данными алгоритма являются приближённое решение системы [math]\displaystyle{ x_k }[/math]. При достижении заданной точности [math]\displaystyle{ \|r_k\| \le \varepsilon\|b\| }[/math] метод считается сошедшимся.

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

Метод сопряжённых градиентов обладает рядом особенностей, важных при его практической реализации. Алгоритм включает регулярную последовательность операций, из которых основная вычислительная сложность приходится на вычисление матричного произведения [math]\displaystyle{ Ad_{k} }[/math]. Остальные операции имеют координатный параллелизм и могут эффективно вычисляться на параллельных архитектурах. Несмотря на то что последовательная сложность одной итерации составляет [math]\displaystyle{ O(n^2) }[/math], большая часть вычислений хорошо векторизуется, а редукции, возникающие при подсчёте норм и скалярных произведений, имеют логарифмическую параллельную сложность.

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

Литература

Метод сопряжённых градиентов // BaseGroup Labs. https://basegroup.ru/community/articles/conjugate

Метод сопряжённых градиентов (дополнительный материал к курсу «Методы оптимизации»). ВМК, 2018. http://www.machinelearning.ru/wiki/images/4/4d/MOMO18_Extra3.pdf