📖 Документация Qumir

← Вернуться в Playground

← Все примеры

Метод сопряжённых градиентов (CG)

Метод сопряжённых градиентов -- итерационный алгоритм решения систем линейных алгебраических уравнений (СЛАУ) -- набора уравнений вида $Ax = b$, где матрица $A$ симметрична и положительно определена. В отличие от прямого метода Гаусса, CG не хранит и не модифицирует матрицу, а лишь умножает её на вектор, что делает его эффективным для больших разреженных систем.

Задача

Найти вектор $x$, удовлетворяющий системе:

$Ax = b$

где $A$ -- симметричная положительно определённая матрица. В примере используется матрица Гильберта:

$A_{i,j} = \frac{1}{i + j + 1}, \quad i,j = 0, \ldots, n-1$

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

Идея алгоритма

  1. Задать начальное приближение $x_0 = 0$. Вычислить начальную невязку $r_0 = b - Ax_0 = b$ и начальное направление спуска $p_0 = r_0$.
  2. На каждой итерации:
    • Вычислить произведение $Ap = A \cdot p$.
    • Найти длину шага: $\alpha = \frac{\langle r, r \rangle}{\langle p, Ap \rangle}$.
    • Обновить решение: $x := x + \alpha \cdot p$.
    • Обновить невязку: $r := r - \alpha \cdot Ap$.
    • Вычислить коэффициент: $\beta = \frac{\langle r_{\text{new}}, r_{\text{new}} \rangle}{\langle r_{\text{old}}, r_{\text{old}} \rangle}$.
    • Обновить направление: $p := r + \beta \cdot p$.
  3. Остановиться, когда $|r|_2 \le \varepsilon$ или достигнуто максимальное число итераций.

Разбор

Параметры и данные

Задаём размер $n = 64$, максимум итераций и допуск. Объявляем рабочие векторы: $x$ (решение), $b$ (правая часть), $r$ (невязка), $p$ (направление спуска), $Ap$ (результат умножения $A \cdot p$):

  | параметры
  цел n, maxIter
  вещ tol
  n := 64              | размер (берите больше, если нужно)
  maxIter := 500       | максимум итераций CG
  tol := 0.0000000001  | допуск по норме невязки

  | векторы
  вещ таб x[0:n-1]      | решение
  вещ таб b[0:n-1]      | правая часть
  вещ таб r[0:n-1]      | текущая невязка r = b - A x
  вещ таб p[0:n-1]      | направление спуска
  вещ таб Ap[0:n-1]     | A * p

Инициализация

Устанавливаем $b_i = 1$, $x_i = 0$, вычисляем начальную невязку $r = b - Ax$ (при $x = 0$ это просто $r = b$), копируем $r$ в $p$ и считаем начальное скалярное произведение $\langle r, r \rangle$:

  | b := 1, x := 0
  нц для i от 0 до n-1
    b[i] := 1.0
    x[i] := 0.0
  кц

  | начальная невязка и начальное направление
  compute_residual(n, x, b, r, Ap)
  нц для i от 0 до n-1
    p[i] := r[i]
  кц
  rsold := dot(n, r, r)

Основной цикл CG

Каждая итерация состоит из одного умножения матрицы на вектор (matvec), двух скалярных произведений (dot) и двух обновлений векторов. Шаг $\alpha$ определяет, как далеко двигаться вдоль направления $p$, а коэффициент $\beta$ обеспечивает сопряжённость нового направления:

  it := 0
  нц пока it < maxIter
    | Ap := A * p
    matvec(n, p, Ap)
    pAp := dot(n, p, Ap)
    если pAp = 0.0 то выход все
    alpha := rsold / pAp

    | x := x + alpha * p; r := r - alpha * Ap
    нц для i от 0 до n-1
      x[i] := x[i] + alpha * p[i]
      r[i] := r[i] - alpha * Ap[i]
    кц

    rnorm := norm2(n, r)
    если rnorm <= tol то выход все

    rsnew := dot(n, r, r)
    если rsold <> 0.0 то
      | p := r + (rsnew/rsold) * p
      вещ beta
      beta := rsnew / rsold
      нц для i от 0 до n-1
        p[i] := r[i] + beta * p[i]
      кц
    иначе
      | rsold == 0: возьмём p := r
      нц для i от 0 до n-1
        p[i] := r[i]
      кц
    все

    rsold := rsnew
    it := it + 1
  кц

Вычисление и вывод невязки

После завершения итераций заново вычисляем точную невязку $r = b - Ax$ и выводим её нормы:

  | финальная невязка и её нормы
  compute_residual(n, x, b, r, Ap)
  вещ res_l2_sq, res_l2, res_inf
  res_l2_sq := 0.0
  res_inf := 0.0
  нц для i от 0 до n-1
    vabs := r[i]; если vabs < 0.0 то vabs := -vabs все
    res_l2_sq := res_l2_sq + r[i] * r[i]
    если vabs > res_inf то res_inf := vabs все
  кц
  res_l2 := sqrt(res_l2_sq)

  вывод "cg: n=", n, ", maxIter=", maxIter, нс
  вывод "iters=", it, нс
  вывод "||r||_2= ", res_l2, нс
  вывод "||r||_inf= ", res_inf, нс

Вспомогательные функции

Программа использует четыре вспомогательных подпрограммы:

matvec -- умножение матрицы Гильберта на вектор ($w_i = \sum_j \frac{v_j}{i+j+1}$), сложность $O(n^2)$:

алг matvec(цел n, вещ таб v[0:n-1], рез вещ таб w[0:n-1])
нач
  цел i, j
  вещ sum
  нц для i от 0 до n-1
    sum := 0.0
    нц для j от 0 до n-1
      sum := sum + (1.0 / (i + j + 1.0)) * v[j]
    кц
    w[i] := sum
  кц
кон

norm2 -- евклидова норма $|v|_2 = \sqrt{\sum_i v_i^2}$:

алг вещ norm2(цел n, вещ таб v[0:n-1])
нач
  цел i
  вещ sum
  sum := 0.0
  нц для i от 0 до n-1
    sum := sum + v[i] * v[i]
  кц
  знач := sqrt(sum)
кон

dot -- скалярное произведение $\langle a, b \rangle = \sum_i a_i \cdot b_i$:

алг вещ dot(цел n, вещ таб a[0:n-1], вещ таб b[0:n-1])
нач
  цел i
  вещ sum
  sum := 0.0
  нц для i от 0 до n-1
    sum := sum + a[i] * b[i]
  кц
  знач := sum
кон

compute_residual -- вычисление невязки $r = b - Ax$:

алг compute_residual(цел n, вещ таб x[0:n-1], вещ таб b[0:n-1], рез вещ таб r[0:n-1], вещ таб w[0:n-1])
нач
  цел i
  matvec(n, x, w)
  нц для i от 0 до n-1
    r[i] := b[i] - w[i]
  кц
кон

Каждая итерация CG стоит $O(n^2)$ (одно умножение матрицы на вектор). Общее число итераций зависит от числа обусловленности матрицы $A$: в точной арифметике метод сходится за $n$ шагов, но на практике для плохо обусловленных матриц может потребоваться значительно больше.

Полная программа

алг main
нач
  | Conjugate Gradient (CG) для матрицы Гильберта A (A[i,j] = 1/(i+j+1)) и правой части b=1.
  | Один модуль: инициализация → решение → проверка → печать норм невязки.

  | параметры
  цел n, maxIter
  вещ tol
  n := 64              | размер (берите больше, если нужно)
  maxIter := 500       | максимум итераций CG
  tol := 0.0000000001  | допуск по норме невязки

  | векторы
  вещ таб x[0:n-1]      | решение
  вещ таб b[0:n-1]      | правая часть
  вещ таб r[0:n-1]      | текущая невязка r = b - A x
  вещ таб p[0:n-1]      | направление спуска
  вещ таб Ap[0:n-1]     | A * p

  | служебные
  цел i, it
  вещ rsold, rsnew, alpha, pAp, rnorm, vabs

  | b := 1, x := 0
  нц для i от 0 до n-1
    b[i] := 1.0
    x[i] := 0.0
  кц

  | начальная невязка и начальное направление
  compute_residual(n, x, b, r, Ap)
  нц для i от 0 до n-1
    p[i] := r[i]
  кц
  rsold := dot(n, r, r)

  it := 0
  нц пока it < maxIter
    | Ap := A * p
    matvec(n, p, Ap)
    pAp := dot(n, p, Ap)
    если pAp = 0.0 то выход все
    alpha := rsold / pAp

    | x := x + alpha * p; r := r - alpha * Ap
    нц для i от 0 до n-1
      x[i] := x[i] + alpha * p[i]
      r[i] := r[i] - alpha * Ap[i]
    кц

    rnorm := norm2(n, r)
    если rnorm <= tol то выход все

    rsnew := dot(n, r, r)
    если rsold <> 0.0 то
      | p := r + (rsnew/rsold) * p
      вещ beta
      beta := rsnew / rsold
      нц для i от 0 до n-1
        p[i] := r[i] + beta * p[i]
      кц
    иначе
      | rsold == 0: возьмём p := r
      нц для i от 0 до n-1
        p[i] := r[i]
      кц
    все

    rsold := rsnew
    it := it + 1
  кц

  | финальная невязка и её нормы
  compute_residual(n, x, b, r, Ap)
  вещ res_l2_sq, res_l2, res_inf
  res_l2_sq := 0.0
  res_inf := 0.0
  нц для i от 0 до n-1
    vabs := r[i]; если vabs < 0.0 то vabs := -vabs все
    res_l2_sq := res_l2_sq + r[i] * r[i]
    если vabs > res_inf то res_inf := vabs все
  кц
  res_l2 := sqrt(res_l2_sq)

  вывод "cg: n=", n, ", maxIter=", maxIter, нс
  вывод "iters=", it, нс
  вывод "||r||_2= ", res_l2, нс
  вывод "||r||_inf= ", res_inf, нс
кон

| ---- Вспомогательные процедуры/функции ----

| матрично-векторное произведение для матрицы Гильберта: w = A * v,
| A[i,j] = 1/(i+j+1), i,j=0..n-1
алг matvec(цел n, вещ таб v[0:n-1], рез вещ таб w[0:n-1])
нач
  цел i, j
  вещ sum
  нц для i от 0 до n-1
    sum := 0.0
    нц для j от 0 до n-1
      sum := sum + (1.0 / (i + j + 1.0)) * v[j]
    кц
    w[i] := sum
  кц
кон

| евклидова норма ||v||_2
алг вещ norm2(цел n, вещ таб v[0:n-1])
нач
  цел i
  вещ sum
  sum := 0.0
  нц для i от 0 до n-1
    sum := sum + v[i] * v[i]
  кц
  знач := sqrt(sum)
кон

| скалярное произведение <a,b>
алг вещ dot(цел n, вещ таб a[0:n-1], вещ таб b[0:n-1])
нач
  цел i
  вещ sum
  sum := 0.0
  нц для i от 0 до n-1
    sum := sum + a[i] * b[i]
  кц
  знач := sum
кон

| r := b - A*x (использует matvec с временным буфером w)
алг compute_residual(цел n, вещ таб x[0:n-1], вещ таб b[0:n-1], рез вещ таб r[0:n-1], вещ таб w[0:n-1])
нач
  цел i
  matvec(n, x, w)
  нц для i от 0 до n-1
    r[i] := b[i] - w[i]
  кц
кон

▶ Запустить пример