将高斯消元的结果与 numpy.linalg.solve 的输出进行比较

Question

在下面的代码中，我为一般方形线性系统 Ax=b 实现了高斯消元和部分旋转。我已经测试了我的代码，它产生了正确的输出。我用它来求解 Ax=b，其中 A 是一个随机 100x100 矩阵，b 是一个随机 100x1 向量。**

但是现在我正在寻求一些帮助，以便将我的解决方案与使用 numpy.linalg.solve 获得的解决方案进行比较。如何将此比较添加到我的代码中？

import numpy as np
def GEPP(A, b, doPricing = True):
    '''
    Gaussian elimination with partial pivoting.
    input: A is an n x n numpy matrix
           b is an n x 1 numpy array
    output: x is the solution of Ax=b 
        with the entries permuted in 
        accordance with the pivoting 
        done by the algorithm
    post-condition: A and b have been modified.
    '''
    n = len(A)
    if b.size != n:

        raise ValueError("Invalid argument: incompatible sizes between"+

                     "A & b.", b.size, n)

    # k represents the current pivot row. Since GE traverses the matrix in the 

    # upper right triangle, we also use k for indicating the k-th diagonal 

    # column index.

    # Elimination

    for k in range(n-1):

        if doPricing:

            # Pivot

            maxindex = abs(A[k:,k]).argmax() + k

            if A[maxindex, k] == 0:


                raise ValueError("Matrix is singular.")

            # Swap

            if maxindex != k:

                A[[k,maxindex]] = A[[maxindex, k]]

                b[[k,maxindex]] = b[[maxindex, k]]

        else:

            if A[k, k] == 0:

                raise ValueError("Pivot element is zero. Try setting doPricing to True.")

       #Eliminate

       for row in range(k+1, n):

           multiplier = A[row,k]/A[k,k]

           A[row, k:] = A[row, k:] - multiplier*A[k, k:]

           b[row] = b[row] - multiplier*b[k]

    # Back Substitution

    x = np.zeros(n)

    for k in range(n-1, -1, -1):

        x[k] = (b[k] - np.dot(A[k,k+1:],x[k+1:]))/A[k,k]

    return x



if __name__ == "__main__":
    A = np.round(np.random.rand(100, 100)*10)
    b =  np.round(np.random.rand(100)*10)
    print (GEPP(np.copy(A), np.copy(b), doPricing = False))

Answer 1

要比较这两种解决方案，请使用np.allclose，它会验证两个数组是否在一些合理的余量内逐个元素地一致。人们不应该期望在浮点运算中完全相等。

my_solution = GEPP(np.copy(A), np.copy(b), doPricing=False)
numpy_solution = np.linalg.solve(A, b)
if np.allclose(my_solution, numpy_solution):
    print("Yay")
else:
    print("NOOOOO")

这会打印“耶”。