← Все примеры

Метод GMRES

Обобщённый метод минимальных невязок (Generalized Minimal Residual) — итерационный метод решения системы линейных уравнений $Ax = b$ для произвольных (в том числе несимметричных) матриц. В отличие от метода сопряжённых градиентов, GMRES не требует, чтобы матрица была симметричной и положительно определённой.

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

GMRES ищет приближённое решение в подпространстве Крылова:

$K_m = \text{span}{r_0,; Ar_0,; A^2 r_0,; \ldots,; A^{m-1} r_0}$

где $r_0 = b - Ax_0$ — начальная невязка. На каждом шаге строится ортонормированный базис этого подпространства (процесс Арнольди), и в нём минимизируется норма невязки.

Если за $m$ шагов решение не найдено, выполняется рестарт: текущее приближение становится новым начальным, и процесс повторяется.

Разбор

Параметры и объявление переменных

алг main
нач
  цел n, m, maxRestart
  вещ tol
  n := 64
  m := 20
  maxRestart := 20
  tol := 0.0000000001

• n — размер системы (матрица $64 \times 64$)
• m — размер подпространства Крылова (каждые $m$ шагов — рестарт)
• maxRestart — максимальное число рестартов
• tol — допуск: если норма невязки $|r| \leq \text{tol}$, решение найдено

Массивы

  вещ таб x[0:n-1]
  вещ таб b[0:n-1]
  вещ таб r[0:n-1]
  вещ таб w[0:n-1]
  вещ таб V[0:m,0:n-1]
  вещ таб H[0:m,0:m-1]
  вещ таб cs[0:m-1]
  вещ таб sn[0:m-1]
  вещ таб g[0:m]
  вещ таб y[0:m-1]
  вещ таб tmpv[0:n-1]

• x — вектор решения, b — правая часть, r — невязка
• V — матрица базисных векторов подпространства Крылова ($m+1$ вектор длины $n$)
• H — верхняя матрица Хессенберга размера $(m+1) \times m$
• cs, sn — косинусы и синусы вращений Гивенса
• g — правая часть малой задачи наименьших квадратов
• y — решение малой треугольной системы

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

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

Правая часть $b = 1$, начальное приближение $x = 0$. Матрица $A$ — матрица Гильберта: $A_{ij} = \frac{1}{i + j + 1}$.

Внешний цикл (рестарты)

  it := 0
  нц пока it < maxRestart
    compute_residual(n, x, b, r, w)
    beta := norm2(n, r)
    rnorm := beta
    если rnorm <= tol то выход все

На каждом рестарте вычисляем невязку $r = b - Ax$ и её норму $\beta = |r|_2$. Если невязка уже достаточно мала — выходим.

Начало процесса Арнольди

    нц для i от 0 до n-1
      V[0,i] := r[i] / beta
    кц
    g[0] := beta
    нц для j от 1 до m
      g[j] := 0.0
    кц

Первый базисный вектор: $v_0 = r / \beta$. Вектор $g$ инициализируется как $(\beta, 0, 0, \ldots, 0)$ — правая часть для задачи наименьших квадратов.

Внутренний цикл Арнольди

    нц для j от 0 до m-1
      jArn := j
      | w = A * v_j
      нц для i от 0 до n-1
        w[i] := 0.0
      кц
      нц для k от 0 до n-1
        tmpv[k] := V[j,k]
      кц
      matvec(n, tmpv, w)

На каждом шаге $j$ вычисляем $w = A \cdot v_j$ — умножаем матрицу на текущий базисный вектор.

Ортогонализация (Грама-Шмидта)

      нц для i от 0 до j
        нц для k от 0 до n-1
          tmpv[k] := V[i,k]
        кц
        hi := dot(n, w, tmpv)
        H[i,j] := hi
        нц для k от 0 до n-1
          w[k] := w[k] - hi * V[i,k]
        кц
      кц
      hip1 := norm2(n, w)
      H[j+1,j] := hip1

Для каждого предыдущего базисного вектора $v_i$ вычисляем проекцию $h_{ij} = \langle w, v_i \rangle$ и вычитаем её из $w$. В результате $w$ становится ортогонален всем предыдущим векторам. Элемент $h_{j+1,j} = |w|$ — длина остатка после ортогонализации.

Нормировка нового базисного вектора

      если hip1 <> 0.0 то
        нц для k от 0 до n-1
          V[j+1,k] := w[k] / hip1
        кц
      иначе
        нц для k от 0 до n-1
          V[j+1,k] := 0.0
        кц
      все

Если $|w| \neq 0$, нормируем и получаем новый базисный вектор $v_{j+1} = w / |w|$.

Вращения Гивенса

      нц для i от 0 до j-1
        tmp :=  cs[i] * H[i,j] + sn[i] * H[i+1,j]
        H[i+1,j] := -sn[i] * H[i,j] + cs[i] * H[i+1,j]
        H[i,j] := tmp
      кц

      denom := sqrt(H[j,j]*H[j,j] + H[j+1,j]*H[j+1,j])
      если denom <> 0.0 то
        cs[j] := H[j,j] / denom
        sn[j] := H[j+1,j] / denom
      иначе
        cs[j] := 1.0
        sn[j] := 0.0
      все
      H[j,j] :=  cs[j] * H[j,j] + sn[j] * H[j+1,j]
      H[j+1,j] := 0.0

Вращения Гивенса — это ортогональные преобразования, которые обнуляют поддиагональные элементы матрицы $H$, приводя её к верхнетреугольному виду. Каждое вращение определяется парой (косинус $c$, синус $s$):

$c = \frac{h_{jj}}{\sqrt{h_{jj}^2 + h_{j+1,j}^2}}, \quad s = \frac{h_{j+1,j}}{\sqrt{h_{jj}^2 + h_{j+1,j}^2}}$

Сначала применяем все предыдущие вращения к новому столбцу $H$, затем вычисляем новое вращение.

Обновление правой части и проверка невязки

      tmp :=  cs[j] * g[j]
      g[j+1] := -sn[j] * g[j]
      g[j] := tmp

      rnorm := g[j+1]; если rnorm < 0.0 то rnorm := -rnorm все
      если rnorm <= tol то
        jArn := j
        выход
      все
    кц

Применяем то же вращение к вектору $g$. Элемент $|g_{j+1}|$ — это оценка нормы невязки. Если она меньше допуска, досрочно выходим из внутреннего цикла.

Обратная подстановка

    нц для i от jArn до 0 шаг -1
      sum := g[i]
      нц для k от i+1 до jArn
        sum := sum - H[i,k] * y[k]
      кц
      если H[i,i] <> 0.0 то
        y[i] := sum / H[i,i]
      иначе
        y[i] := 0.0
      все
    кц

Решаем верхнетреугольную систему $Hy = g$ обратной подстановкой — снизу вверх.

Обновление решения

    нц для k от 0 до n-1
      sum := 0.0
      нц для i от 0 до jArn
        sum := sum + V[i,k] * y[i]
      кц
      x[k] := x[k] + sum
    кц

Обновляем приближение: $x := x_0 + V \cdot y$, где $V$ — матрица базисных векторов, $y$ — решение малой системы.

Вывод результата

  compute_residual(n, x, b, r, w)
  вещ res_l2_sq, res_l2, res_inf, vabs
  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)

  вывод "gmres: n=", n, ", m=", m, нс
  вывод "restarts=", it, нс
  вывод "||r||_2= ", res_l2, нс
  вывод "||r||_inf= ", res_inf, нс
кон

После завершения выводим параметры задачи, число рестартов и нормы финальной невязки ($L_2$ и $L_\infty$).

Вспомогательные алгоритмы

Умножение матрицы на вектор (матрица Гильберта $A_{ij} = \frac{1}{i+j+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 = \sqrt{\sum 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)
кон

Скалярное произведение $\langle a, b \rangle = \sum 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
кон

Вычисление невязки $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

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

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

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

  | векторы/матрицы
  вещ таб x[0:n-1]      | решение
  вещ таб b[0:n-1]      | правая часть
  вещ таб r[0:n-1]      | текущая невязка
  вещ таб w[0:n-1]      | рабочий вектор (A*v_j)
  вещ таб V[0:m,0:n-1]  | базис Крылова (m+1 векторов длины n)
  вещ таб H[0:m,0:m-1]  | верх. Гессенберг (размер (m+1) x m)
  вещ таб cs[0:m-1]     | косинусы Гивенса
  вещ таб sn[0:m-1]     | синусы Гивенса
  вещ таб g[0:m]        | правая часть для НУ
  вещ таб y[0:m-1]      | решение малого треугольника
  вещ таб tmpv[0:n-1]   | временный вектор для копирования строк V

  | служебные
  цел i, j, k, it, jArn
  вещ sum, absv, beta, rnorm, tmp, denom, hi, hip1, rho

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

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

  | главный цикл рестартов
  it := 0
  нц пока it < maxRestart
    compute_residual(n, x, b, r, w)
    beta := norm2(n, r)
    rnorm := beta
    если rnorm <= tol то выход все

    | v0 = r / beta
    нц для i от 0 до n-1
      V[0,i] := r[i] / beta
    кц
    | g := [beta, 0, ..., 0]
    g[0] := beta
    нц для j от 1 до m
      g[j] := 0.0
    кц

    | внутренний цикл Арнольди с поворотами Гивенса
    jArn := m - 1     | по умолчанию берём последний шаг (если не выйдем раньше)
    нц для j от 0 до m-1
      jArn := j
      | w = A * v_j
      нц для i от 0 до n-1
        w[i] := 0.0
      кц
      | копируем строку V[j,*] во временный вектор и умножаем
      нц для k от 0 до n-1
        tmpv[k] := V[j,k]
      кц
      matvec(n, tmpv, w)

      | ортогонализация по предыдущим v_i
      нц для i от 0 до j
        | скопировать V[i,*] в tmpv для вызова dot
        нц для k от 0 до n-1
          tmpv[k] := V[i,k]
        кц
        hi := dot(n, w, tmpv)
        H[i,j] := hi
        нц для k от 0 до n-1
          w[k] := w[k] - hi * V[i,k]
        кц
      кц
      hip1 := norm2(n, w)
      H[j+1,j] := hip1
      если hip1 <> 0.0 то
        нц для k от 0 до n-1
          V[j+1,k] := w[k] / hip1
        кц
      иначе
        нц для k от 0 до n-1
          V[j+1,k] := 0.0
        кц
      все

      | применяем прежние повороты к столбцу j матрицы H
      нц для i от 0 до j-1
        tmp :=  cs[i] * H[i,j] + sn[i] * H[i+1,j]
        H[i+1,j] := -sn[i] * H[i,j] + cs[i] * H[i+1,j]
        H[i,j] := tmp
      кц

      | новый поворот для (H[j,j], H[j+1,j])
      denom := sqrt(H[j,j]*H[j,j] + H[j+1,j]*H[j+1,j])
      если denom <> 0.0 то
        cs[j] := H[j,j] / denom
        sn[j] := H[j+1,j] / denom
      иначе
        cs[j] := 1.0
        sn[j] := 0.0
      все
      H[j,j] :=  cs[j] * H[j,j] + sn[j] * H[j+1,j]
      H[j+1,j] := 0.0

      | обновляем g = Q^T * e1*beta на лету
      tmp :=  cs[j] * g[j]
      g[j+1] := -sn[j] * g[j]
      g[j] := tmp

      | оценка невязки
      rnorm := g[j+1]; если rnorm < 0.0 то rnorm := -rnorm все
      если rnorm <= tol то
        jArn := j
        выход
      все
    кц
    | если не вышли по допуску внутри — jArn уже равен m-1 по умолчанию
    | решаем верхнетреугольную систему R*y = g (R = H[0..jArn, 0..jArn])
    нц для i от 0 до jArn
      y[i] := 0.0
    кц
    нц для i от jArn до 0 шаг -1
      sum := g[i]
      нц для k от i+1 до jArn
        sum := sum - H[i,k] * y[k]
      кц
      если H[i,i] <> 0.0 то
        y[i] := sum / H[i,i]
      иначе
        y[i] := 0.0
      все
    кц

    | x := x + V_{0..jArn} * y
    нц для k от 0 до n-1
      sum := 0.0
      нц для i от 0 до jArn
        sum := sum + V[i,k] * y[i]
      кц
      x[k] := x[k] + sum
    кц

    | проверим норму после обновления
    compute_residual(n, x, b, r, w)
    rnorm := norm2(n, r)
    если rnorm <= tol то выход все

    it := it + 1
  кц
  | финальная невязка и её нормы
  compute_residual(n, x, b, r, w)
  вещ res_l2_sq, res_l2, res_inf, vabs
  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)

  вывод "gmres: n=", n, ", m=", m, нс
  вывод "restarts=", 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]
  кц
кон

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