如何解决大条件数的 Ax=b

问题描述

我正在处理高度病态的矩阵 (cond> 10^25)。我已经尝试了大部分 scipy 方法：splu、gmers solve 和来自 pyamg 的 pyamg.krylov.steepest_descent。

也是迭代细化方法：

def iterRef(A,b,tol=1e-5,maxiter=100,verbose=True):
    """
    Solve the equation a x = b for x using Iterative Refinement method.

    :param A: [(M,M) array_like] A square matrix.
    :param b: [(M,) array like]Right-hand side matrix in a x = b.
    :param maxIter: [int] max number of iteration for convergence.
    :param t:  #!
    :return x: (M,) or (M,N) ndarray
        Solution to the system a x = b. Shape of the return matches the shape of b.

    **Reference:**
        Burden,R.L. and Faires,J.D.,2011. Numerical analysis.
    """

    
    # declarations
    n = len(b)
    xx = np.zeros_like(b)
    r = np.zeros_like(b)

    lu = splu(A)
    x = lu.solve(b)
    res = np.sum(A.dot(x)-b)
    print("res :: A * x - b = {:e}".format(res))

    # check if converged
    if np.abs(res) < tol:
        print("IR ::: A * x - b = {:e}".format(res))
        return x

    k = 1                                       # step 1
    while (k <= maxiter):                       # step 2
        r = b - A.dot(x)                        # step 3
        y = lu.solve(r)                         # step 4
        xx = copy(x + y)                        # step 5

        if k == 1:                              # step 6
            # t = 16
            # COND = np.linalg.norm(y)/np.linalg.norm(xx) * 10**t
            # print("cond is ::::",COND)
            COND = np.linalg.cond(A.toarray())

        norm_ = np.linalg.norm(x-xx) # * np.linalg.norm(x) * 1e10
        print("iteration {:3d},norm = {:e}".format(k,norm_))

        if norm_ < tol:                         # step 7

            if verbose:
                print("Conditional number of matrix A is: {:e}".format(COND))
                print("The procedure was successful.")
                print("IR: A * x - b = {:e}".format(np.sum(A.dot(xx)-b)))
                print(f"number of iteration is : {k:d}")
                print(" ")
            return xx

        k += 1                                 # step 8
        x = copy(xx)                           # step 9

    print("Max iteration exceeded.")
    print("The procedure was not successful.")
    print("Conditional number of matrix A is: {:e}".format(COND))
    print(" ")

    return None

准确率取决于条件数，对于更大的条件，准确率会消失。

我尝试了 eigen3 的一些方法：Catalogue of decompositions offered by Eigen：

VectorXd x = A.partialPivLu().solve(b);
VectorXd x = A.fullPivLu().solve(b);
VectorXd x = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b);

这是 python 和 eigen 中完整示例的 github 链接。是否可以使用具有四倍精度或更高精度的特征模块？有什么例子吗？我不确定我们可以用 scipy 模块做到这一点。

C++ 文件的输出为：

condition number is :3.16172e+28
The relative error for                             partialPivLu  is        0.0000000000000011:
The relative error for                                fullPivLu  is        0.0157277957331164:
The relative error for                            householderQr  is        0.0000000000000020:
The relative error for                      colPivHouseholderQr  is        0.0000000000577384:
The relative error for                     fullPivHouseholderQr  is        0.0157277957331138:
The relative error for          completeOrthogonalDecomposition  is        0.0157277957331138:
The relative error for                                      llt  is                      -nan:
The relative error for                                     ldlt  is        0.0000000001045718:
The relative error for                                     ldlt  is 396613336624311706845184.0000000000000000:
The relative error for                                   bdcSvd  is        0.0157277957329992:
The relative error for                                jacobiSvd  is        0.0157277957528586:

我将 python 和 C++ 代码放在附加的链接中以防万一。

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

eigen3 scipy scipy solver