我的 GMRES 实施有什么问题？答案

【问题标题】：What's wrong with my GMRES implementation?我的 GMRES 实施有什么问题？
【发布时间】：2016-10-24 01:19:50
【问题描述】：

我正在尝试在 Jupyter Notebook 中实现 GMRES，即（以防你不知道）：

这是我的代码：

import numpy as np

def GMRes(A, b, x0, e, nmax_iter, restart=None):
    r = b - np.asarray(np.dot(A, x0)).reshape(-1)

    x = []
    q = [0] * (nmax_iter)

    x.append(r)

    q[0] = r / np.linalg.norm(r)

    h = np.zeros((nmax_iter + 1, nmax_iter))

    for k in range(nmax_iter):
        y = np.asarray(np.dot(A, q[k])).reshape(-1)

        for j in range(k):
            h[j, k] = np.dot(q[j], y)
            y = y - h[j, k] * q[j]
        h[k + 1, k] = np.linalg.norm(y)
        if (h[k + 1, k] != 0 and k != nmax_iter - 1):
            q[k + 1] = y / h[k + 1, k]

        b = np.zeros(nmax_iter + 1)
        b[0] = np.linalg.norm(r)

        result = np.linalg.lstsq(h, b)[0]

        x.append(np.dot(np.asarray(q).transpose(), result) + x0)

    return x

在我看来应该是正确的，但是当我执行时：

A = np.matrix('1 1; 3 -4')
b = np.array([3, 2])
x0 = np.array([1, 2])

e = 0
nmax_iter = 5

x = GMRes(A, b, x0, e, nmax_iter)

print(x)

注意：目前e 什么都不做。

我明白了：

[array([0, 7]), array([ 1.,  2.]), array([ 1.35945946,  0.56216216]), array([ 1.73194463,  0.80759216]), array([ 2.01712479,  0.96133459]), array([ 2.01621042,  0.95180204])]

x[k] 应该接近(32/7, -11/7)，因为这是结果，但它却接近(2, 1)，我做错了什么？

【问题讨论】：

标签： python numpy scipy numerical-methods jupyter-notebook

【解决方案1】：

我认为算法给出了正确的结果。

您正在尝试解决 Ax=b 其中：

如果您尝试手动求解，则您尝试求解的矩阵运算等效于可以使用替换求解的系统。

如果你尝试解决它，你会发现解决方案是：

您的算法给出的解决方案相同。

您可以使用 scipy 中已有的 GMRES 实现来仔细检查：

import scipy.sparse.linalg as spla
import numpy as np

A = np.matrix('1 1; 3 -4')
b = np.array([3, 2])
x0 = np.array([1, 2])
spla.gmres(A,b,x0)

哪些输出

array([ 2.,  1.])

【讨论】：

我是世界上最愚蠢的人，浪费了很多时间，因为我解错了方程 xD...谢谢！

【解决方案2】：

请注意，此算法正在收敛到正确的结果，但收敛速度太慢。 GMRES 迭代收敛的最大次数不应超过矩阵 A 的维数。如果矩阵 A 的维数为 n，则第 (n+1) 个 Arnoldi 向量应为零，例如我们应该能够用 n 个 Arnoldi 向量完全跨越 Krylov 空间。我将只应用以下补丁，并且事情应该可以正常工作：

-    for k in range(nmax_iter):
+    for k in range(min(nmax_iter, A.shape[0])):
         y = np.asarray(np.dot(A, q[k])).reshape(-1)

-        for j in range(k):
+        for j in range(k + 1):

那么解向量序列应该是：

 [array([ 1.        ,  0.35294118]), array([ 2.,  1.])]

例如从 n = 2 开始，我们在两次迭代中收敛。

【讨论】：