GMRES方法求解Ax=b时获取残差向量的问题求助

问题解析 #

你现在得到的ylist里的每个元素是残差的2-范数（标量），不是n维残差向量。scipy的gmres函数默认的callback参数，每次迭代传入的是当前残差的2-范数（单个数值），而非完整的残差向量。另外你的代码里gmres_counter类的append方法缺少self参数，属于Python类方法的基础错误，会导致调用该方法时抛出异常。

修正方案 #

1. 正确获取残差范数（修正现有代码） #

如果你的需求是跟踪残差范数的变化，只需修正类方法的self参数问题，代码如下：

import numpy as np
from scipy.sparse.linalg import gmres
from matplotlib import pyplot as plt

class gmres_counter(object):
    def __init__(self, disp=True):
        self._disp = disp
        self.niter = 0
        self.callbacks = []
        self.internal_list = []
    def append(self, elem):  # 添加self参数，修复类方法语法错误
        self.internal_list.append(elem)
    def getList(self):
        return self.internal_list
    def __call__(self, rk=None):
       self.callbacks.append(rk)
       self.internal_list.append(rk)
       self.niter += 1

n=14
P = np.random.normal(size=[n, n])
b = np.ones(n)

counter = gmres_counter()
x, info = gmres(P, b, callback=counter)

ylist=counter.getList()
print("残差范数列表（长度为迭代次数）：", ylist)

此时ylist中的每个元素是对应迭代步的残差范数，列表长度等于GMRES的迭代次数（通常远小于n）。

2. 获取完整残差向量 #

如果需要获取每个迭代步的n维残差向量r = b - P@x，可以利用scipy 1.7.0及以上版本支持的callback_type='x'参数，让callback接收当前的x估计值，再计算残差向量：

import numpy as np
from scipy.sparse.linalg import gmres
from matplotlib import pyplot as plt

class gmres_residual_tracker(object):
    def __init__(self, A, b):
        self.A = A
        self.b = b
        self.residuals = []  # 存储每一步的残差向量
        self.niter = 0
    def __call__(self, xk=None):
        # 计算当前迭代步的残差向量
        rk = self.b - self.A @ xk
        self.residuals.append(rk)
        self.niter += 1

n=14
P = np.random.normal(size=[n, n])
b = np.ones(n)

tracker = gmres_residual_tracker(P, b)
# 指定callback_type='x'，让callback接收当前的x估计值
x, info = gmres(P, b, callback=tracker, callback_type='x')

residual_vectors = tracker.residuals
print("第1个残差向量：", residual_vectors[0])
print("残差向量数量（等于迭代次数）：", len(residual_vectors))

此时residual_vectors中的每个元素是n维numpy数组（残差向量），列表长度等于迭代次数。

基础知识点补充 #

• Python类的实例方法必须以self作为第一个参数，否则无法访问实例的属性（比如你的append方法之前缺少self，属于基础语法错误）。
• scipy的gmres默认callback传入的是残差范数，只有指定callback_type='x'时，才会传入当前的x估计值；如果需要兼容更早版本的scipy，可以手动在迭代后根据x的历史记录计算残差，但效率不如直接使用callback_type='x'。

内容的提问来源于stack exchange，提问作者Noah F