我的块GMRES重启实现异常：收敛失败问题排查

架构老林

2026-6-12

块GMRES算法实现故障排查

需求背景

我正在实现块GMRES算法（求解线性方程组$Ax=B$，其中$B$是$n×r$矩阵，$r\ll n$），目标是先完成Python版本（可与Scipy实现对比），之后移植到C++。

现有代码

正常运行的重启GMRES（标量版本）

import numpy as np

def gmres_restart(A, b, x0=None, tol=1e-5, max_iters=100, restart_dim=30):
    n = len(b)
    if x0 is None:
        x = np.zeros(n)
    else:
        x = x0.copy()

    r = b - A @ x
    beta = np.linalg.norm(r)
    if beta < tol:
        return x

    for outer_it in range(0, max_iters, restart_dim):
        # Arnoldi process
        V = np.zeros((n, restart_dim + 1))
        H = np.zeros((restart_dim + 1, restart_dim))
        V[:, 0] = r / beta

        for j in range(restart_dim):
            w = A @ V[:, j]
            for i in range(j + 1):
                H[i, j] = np.dot(V[:, i], w)
                w -= H[i, j] * V[:, i]
            H[j + 1, j] = np.linalg.norm(w)

            if H[j + 1, j] < 1e-14:
                break

            V[:, j + 1] = w / H[j + 1, j]

        # Solve least squares:
        e1 = np.zeros(j + 2)
        e1[0] = beta
        H_small = H[:j + 2, :j + 1]
        y, *_ = np.linalg.lstsq(H_small, e1, rcond=None)

        # Update solution
        x += V[:, :j + 1] @ y

        # New residual
        r = b - A @ x
        beta = np.linalg.norm(r)
        print(f"GMRES iter {outer_it + restart_dim}, residual = {beta:.2e}")
        if beta < tol:
            print("Converged.")
            return x

    print("Did not converge.")
    return x

块Arnoldi过程实现（基于Y. Saad《NUMERICAL METHODS FOR LARGE EIGENVALUE PROBLEMS（第二版）》算法6.8）

def block_arnoldi(A, V0, m):
    n, s = V0.shape
    V = np.zeros((n, s * (m + 1)))
    H = np.zeros((s * (m + 1), s * m))

    # Normalize the initial block
    Q, _ = np.linalg.qr(V0, mode='reduced')
    V[:, :s] = Q

    for j in range(m):
        # Compute W = A * V_j
        Vj = V[:, j * s : (j + 1) * s]
        Wj = A @ Vj

        # Orthogonalize against previous blocks
        for i in range(j + 1):
            Vi = V[:, i * s : (i + 1) * s]
            Hij = Vi.T @ Wj
            H[i * s : (i + 1) * s, j * s : (j + 1) * s] = Hij
            Wj -= Vi @ Hij

        # QR decomposition to get next block of orthonormal vectors
        Q, R = np.linalg.qr(Wj, mode='reduced')
        V[:, (j + 1) * s : (j + 2) * s] = Q
        H[(j + 1) * s : (j + 2) * s, j * s : (j + 1) * s] = R

    return V, H

待修复的块GMRES重启算法

def block_gmres_restart(A, B, X0=None, tol=1e-5, max_iters=100, restart_dim=30):
    n, s = B.shape
    if X0 is None:
        X = np.zeros((n, s))
    else:
        X = X0.copy()

    R = B - A @ X
    beta = np.linalg.norm(R, ord='fro')

    if beta < tol:
        return X

    for outer in range(0, max_iters, restart_dim):
        V, H = block_arnoldi(A, R, restart_dim)

        # Solve least-squares
        e1 = np.zeros(((restart_dim + 1) * s, s))
        e1[:s, :] = np.linalg.norm(R, axis=0).reshape(-1, 1) * np.eye(s)
        H_small = H[:(restart_dim + 1) * s, :restart_dim * s]
        Y, *_ = np.linalg.lstsq(H_small, e1, rcond=None)

        # Update solution
        X += V[:, :restart_dim * s] @ Y

        # New residual
        R = B - A @ X
        res_norm = np.linalg.norm(R, ord='fro')
        print(f"Block GMRES iter {outer + restart_dim}, residual = {res_norm:.2e}")

        if res_norm < tol:
            print("Converged.")
            return X

    print("Did not converge.")
    return X

问题现象

当$r>1$时，块GMRES完全无法工作；
当$r=1$时，部分场景下无法收敛（例如$n=30$时残差反而增大），但标量GMRES可正常收敛；

测试案例：

np.random.seed(42)  # fix randomness

n = 30
A = np.diag(np.linspace(1, 100, n)) + 0.001 * np.random.randn(n, n)
A = (A + A.T) / 2  # Ensure symmetry

b = np.random.randn(n, 1)

# 标量GMRES正常收敛
x = gmres_restart(A, b[:, 0], restart_dim=10)
# 块GMRES残差增大，无法收敛
x = block_gmres_restart(A, b, restart_dim=10)

故障原因分析与修复

核心错误：初始残差的最小二乘右端项构造错误

在块GMRES中，初始残差$R_0 = B - AX_0$需要先做QR分解得到$R_0 = V_0 R_0^{{qr}$（其中$V_0$是正交块），这一步你在`block_arnoldi`中已经完成，但在构造最小二乘问题的右端项时，错误地用了**列范数乘单位阵**，而非QR分解得到的$R_0}{qr}$。

具体修复步骤

修正初始残差的处理：在block_gmres_restart中，先对初始残差做QR分解，保留分解后的$R_0^{qr}$，并将正交块$V_0$传入block_arnoldi：

def block_gmres_restart(A, B, X0=None, tol=1e-5, max_iters=100, restart_dim=30):
    n, s = B.shape
    if X0 is None:
        X = np.zeros((n, s))
    else:
        X = X0.copy()

    R = B - A @ X
    beta = np.linalg.norm(R, ord='fro')
    if beta < tol:
        return X
    
    # 对初始残差做QR分解，得到正交块V0和上三角阵R0_qr
    V0, R0_qr = np.linalg.qr(R, mode='reduced')

    for outer in range(0, max_iters, restart_dim):
        # 传入正交块V0而非原始残差R
        V, H = block_arnoldi(A, V0, restart_dim)

        # 构造最小二乘右端项：用QR分解得到的R0_qr
        e1 = np.zeros(((restart_dim + 1)*s, s))
        e1[:s, :] = R0_qr
        H_small = H[:(restart_dim + 1)*s, :restart_dim*s]
        Y, *_ = np.linalg.lstsq(H_small, e1, rcond=None)

        # 更新解
        X += V[:, :restart_dim*s] @ Y

        # 计算新残差并重新做QR分解，为下一轮重启准备
        R = B - A @ X
        res_norm = np.linalg.norm(R, ord='fro')
        print(f"Block GMRES iter {outer + restart_dim}, residual = {res_norm:.2e}")
        if res_norm < tol:
            print("Converged.")
            return X
        # 更新V0为新残差的正交块
        V0, R0_qr = np.linalg.qr(R, mode='reduced')

    print("Did not converge.")
    return X

修正块Arnoldi的冗余QR分解：由于我们已经在block_gmres_restart中对初始残差做了QR分解，block_arnoldi中可以去掉对V0的重复QR分解，避免不必要的计算：

def block_arnoldi(A, V0, m):
    n, s = V0.shape
    V = np.zeros((n, s * (m + 1)))
    H = np.zeros((s * (m + 1), s * m))

    # 直接使用传入的正交块V0，无需重复QR
    V[:, :s] = V0

    for j in range(m):
        Vj = V[:, j * s : (j + 1) * s]
        Wj = A @ Vj

        for i in range(j + 1):
            Vi = V[:, i * s : (i + 1) * s]
            Hij = Vi.T @ Wj
            H[i * s : (i + 1) * s, j * s : (j + 1) * s] = Hij
            Wj -= Vi @ Hij

        Q, R = np.linalg.qr(Wj, mode='reduced')
        V[:, (j + 1) * s : (j + 2) * s] = Q
        H[(j + 1) * s : (j + 2) * s, j * s : (j + 1) * s] = R

    return V, H

额外说明

块GMRES的最小二乘问题本质是$\min_Y | R_0^{qr} - H Y |_F$，其中$H$是块Arnoldi生成的块Hessenberg矩阵，只有右端项正确，才能得到正确的更新量$Y$；
重启时需要对新残差重新做QR分解，为下一轮块Arnoldi提供正交初始块，保证迭代的稳定性；
符号问题：解更新的X += ...是正确的，与标量GMRES逻辑一致，无需修改为减号。

内容的提问来源于stack exchange，提问作者Pierre Beaujean

我的块GMRES重启实现异常：收敛失败问题排查

需求背景 #

现有代码 #

正常运行的重启GMRES（标量版本） #

块Arnoldi过程实现（基于Y. Saad《NUMERICAL METHODS FOR LARGE EIGENVALUE PROBLEMS（第二版）》算法6.8） #

待修复的块GMRES重启算法 #

问题现象 #

故障原因分析与修复 #

核心错误：初始残差的最小二乘右端项构造错误 #

具体修复步骤 #

额外说明 #