Блочный GMRES: как избежать стагнации — QR для начальной невязки

Разбор правильной постановки малой МНК-задачи в блочном GMRES для нескольких правых частей: зачем нужен QR-разбор начальной невязки и чем опасна диагональ норм.

Блочный GMRES обманчиво близок к варианту с одной правой частью: та же идея крыловского подпространства, тот же шаг наименьших квадратов, похожая стратегия перезапуска. Подвох в том, как начальная невязка встраивается в малую задачу МНК. Если этот элемент задан неверно, метод для нескольких правых частей может стагнировать или даже расходиться. Ниже — минимальный рабочий разбор, который показывает ловушку и точную правку.

Постановка задачи и минимальная опорная версия GMRES

Для ориентира начнём с компактного GMRES с перезапусками для одной правой части. Он использует процесс Арнольди и малую задачу наименьших квадратов, как это описано в классических источниках.

import numpy as np
def gmres_cycle(M, rhs, x_init=None, eps=1e-5, max_steps=100, krylov_dim=30):
    n = len(rhs)
    sol = np.zeros(n) if x_init is None else x_init.copy()
    res = rhs - M @ sol
    beta0 = np.linalg.norm(res)
    if beta0 < eps:
        return sol
    for outer_k in range(0, max_steps, krylov_dim):
        basis = np.zeros((n, krylov_dim + 1))
        hess = np.zeros((krylov_dim + 1, krylov_dim))
        basis[:, 0] = res / beta0
        for j in range(krylov_dim):
            wv = M @ basis[:, j]
            for i in range(j + 1):
                hess[i, j] = np.dot(basis[:, i], wv)
                wv -= hess[i, j] * basis[:, i]
            hess[j + 1, j] = np.linalg.norm(wv)
            if hess[j + 1, j] < 1e-14:
                break
            basis[:, j + 1] = wv / hess[j + 1, j]
        g = np.zeros(j + 2)
        g[0] = beta0
        Hsub = hess[:j + 2, :j + 1]
        coef, *_ = np.linalg.lstsq(Hsub, g, rcond=None)
        sol += basis[:, :j + 1] @ coef
        res = rhs - M @ sol
        beta0 = np.linalg.norm(res)
        if beta0 < eps:
            return sol
    return sol

Наивный блочный GMRES и в чём он ошибается

Прямолинейный «порт» на несколько правых частей заменяет вектор невязки на блочную невязку, а Арнольди — на блочный Арнольди. Задача МНК решается в норме Фробениуса. Два фрагмента ниже задают процедуру блочного Арнольди и блочный GMRES с перезапуском, повторяющий структуру варианта с одной правой частью.

def arnoldi_block(M, W0, mstep):
    n, blk = W0.shape
    W = np.zeros((n, blk * (mstep + 1)))
    T = np.zeros((blk * (mstep + 1), blk * mstep))
    Q, _ = np.linalg.qr(W0, mode='reduced')
    W[:, :blk] = Q
    for j in range(mstep):
        Bj = W[:, j * blk : (j + 1) * blk]
        Zj = M @ Bj
        for i in range(j + 1):
            Bi = W[:, i * blk : (i + 1) * blk]
            Tij = Bi.T @ Zj
            T[i * blk : (i + 1) * blk, j * blk : (j + 1) * blk] = Tij
            Zj -= Bi @ Tij
        Qj, Rj = np.linalg.qr(Zj, mode='reduced')
        W[:, (j + 1) * blk : (j + 2) * blk] = Qj
        T[(j + 1) * blk : (j + 2) * blk, j * blk : (j + 1) * blk] = Rj
    return W, T
def gmres_block_cycle(M, RHS, X_init=None, eps=1e-5, max_steps=100, krylov_dim=30):
    n, blk = RHS.shape
    X = np.zeros((n, blk)) if X_init is None else X_init.copy()
    RES = RHS - M @ X
    if np.linalg.norm(RES, ord='fro') < eps:
        return X
    for outer_k in range(0, max_steps, krylov_dim):
        W, T = arnoldi_block(M, RES, krylov_dim)
        G = np.zeros(((krylov_dim + 1) * blk, blk))
        G[:blk, :] = np.linalg.norm(RES, axis=0).reshape(-1, 1) * np.eye(blk)
        Tsub = T[:(krylov_dim + 1) * blk, :krylov_dim * blk]
        Ycoef, *_ = np.linalg.lstsq(Tsub, G, rcond=None)
        X += W[:, :krylov_dim * blk] @ Ycoef
        RES = RHS - M @ X
        if np.linalg.norm(RES, ord='fro') < eps:
            return X
    return X

Проблема тонкая и не связана с «переворотом знака» обновления. В случае одной правой части правая часть в МНК — это скаляр β, умноженный на первый канонический базисный вектор. Неправильная обобщающая замена на диагональную матрицу из поколоночных норм блочной невязки делает малую систему несогласованной. В итоге решатель минимизирует не ту величину — отсюда рост или застой невязки, как только правых частей становится больше одной.

Суть: QR-разложение блока невязки определяет малую задачу

Для блочного GMRES нужно QR-разложение текущего блока невязки. Вычислите Q0, R0 = qr(R), используйте Q0 как начальный блок в процессе Арнольди и поместите R0 непосредственно в левый верхний блок правой части малой задачи МНК. Это и есть ключевое отличие от случая с одной правой частью, где ту же роль играет скаляр β.

Исправленный блочный GMRES с корректной постановкой МНК

Следующая версия делает именно это. Остальная часть алгоритма остаётся неизменной.

def gmres_block_cycle_fixed(M, RHS, X_init=None, eps=1e-5, max_steps=100, krylov_dim=30):
    n, blk = RHS.shape
    X = np.zeros((n, blk)) if X_init is None else X_init.copy()
    RES = RHS - M @ X
    if np.linalg.norm(RES, ord='fro') < eps:
        return X
    for outer_k in range(0, max_steps, krylov_dim):
        Q0, R0 = np.linalg.qr(RES, mode='reduced')
        W, T = arnoldi_block(M, Q0, krylov_dim)
        G = np.zeros(((krylov_dim + 1) * blk, blk))
        G[:blk, :blk] = R0
        Tsub = T[:(krylov_dim + 1) * blk, :krylov_dim * blk]
        Ycoef, *_ = np.linalg.lstsq(Tsub, G, rcond=None)
        X += W[:, :krylov_dim * blk] @ Ycoef
        RES = RHS - M @ X
        if np.linalg.norm(RES, ord='fro') < eps:
            return X
    return X

Воспроизводимая тестовая заготовка

Небольшая заготовка ниже повторяет описанный выше сценарий. Она полезна для быстрой проверки и сравнения варианта с одной правой частью и блочного варианта на одной и той же матрице.

np.random.seed(42)
n = 20
A = np.diag(np.linspace(1, 100, n)) + 0.001 * np.random.randn(n, n)
A = (A + A.T) / 2
b = np.random.randn(n, 1)
x_ref = gmres_cycle(A, b[:, 0], krylov_dim=10)
X_blk = gmres_block_cycle_fixed(A, b, krylov_dim=10)

Почему этот нюанс важен

Блочный вариант минимизирует норму Фробениуса на блочном крыловском подпространстве. Малая задача МНК должна корректно отражать геометрию невязки при перезапуске. QR-разложение блока невязки фиксирует эту геометрию через R0, и подстановка R0 в правую часть малой системы согласует минимизацию с истинной целью. Замена на диагональные нормы кажется безобидной, но нарушает это соответствие.

Выводы

Если блочный GMRES с несколькими правыми частями даёт сбой, тогда как вариант с одной правой частью работает, проверьте формирование малой задачи МНК при перезапуске. Вычислите Q0, R0 = qr(R), инициализируйте процесс Арнольди блоком Q0 и задайте ведущий блок правой части равным R0. При таком устройстве направление обновления уже корректно, и переворачивать знаки не требуется. Сопоставляя реализации на Python и C++, сохраняйте матричные операции как в коде выше: метод по своей природе итерационный, и явные циклы для шагов Арнольди вполне уместны.