本征中的稀疏矩阵乘法给出错误的结果？

Question

我正在我的项目中使用Eigen，但遇到了一个奇怪的问题。我有复杂的稀疏矩阵A和B（1500x1500或更大），并将它们与系数相乘。

[当A = B，并取1的向量x时，我期望

[(A-B)*x = 0，(A*B-B*A)*x = 0，

(A*A*B*B - B*B*A*A)*x = 0，

等我确实在所有这些情况下都得到了这个结果。 （A.isApprox(B)的值为1，(A-B).norm() = 0）。

但是，当我将矩阵乘以两倍时，如

(c1*A*c2*A*d1*B*d2*B - d1*B*d2*B*c1*A*c2*A)*x，

我得到一个非零的结果，这对我来说没有意义，因为标量应该与矩阵相通。实际上，如果我这样做，

(c1*c2*d1*d2*A*A*B*B - d1*d2*c1*c2*B*B*A*A)*x

我得到零。每当系数散布在矩阵运算中时，我都会得到一个非零的结果。

我没有使用任何编译器优化等。

我在这里做错了什么？

编辑：我已经做了一个简单的例子。也许我想念一些愚蠢的东西，但是确实如此。这给了我10 ^ 20的错误。

'''

#include <iostream>
#include <cmath>
#include <vector>
#include <Eigen/Sparse>
#include <complex>

typedef std::complex<double> Scalar;
typedef Eigen::SparseMatrix<Scalar, Eigen::RowMajor> SpMat;
typedef Eigen::Triplet<Scalar> trip;

int main(int argc, const char * argv[]) {

    double k0 = M_PI;
    double dz = 0.01;
    double nz = 1500;

    std::vector<double> rhos(nz), atten(nz), cp(nz);

    for(int i = 0; i < nz; ++i){
        if(i < 750){
            rhos[i] = 1.5;
            cp[i] = 2500;
            atten[i] = 0.5;
        }
        else{
            rhos[i] = 1;
            cp[i] = 1500;
            atten[i] = 0;
        }
    }

    Scalar ci, eta, n, rho,  drhodz;
    Scalar t1, t2, t3, t4;

    ci = Scalar(0,1);
    eta = 1.0/(40.0*M_PI*std::log10(std::exp(1.0)));

    int Mp = 6;

    std::vector<std::vector<trip> > mat_entries_N(Mp), mat_entries_D(Mp);

    for(int i = 0; i < nz; ++i){
        n = 1500./cp[i] * (1.+ ci * eta * atten[i]);
        rho = rhos[i];

        if(i > 0 && i < nz-1){
            drhodz = (rhos[i+1]-rhos[i-1])/(2*dz);
        }
        else if(i == 0){
            drhodz = (rhos[i+1]-rhos[i])/(dz);
        }
        else if(i == nz-1){
            drhodz = (rhos[i]-rhos[i-1])/(dz);
        }

        t1 = (n*n - 1.);

        t2 = 1./(k0*k0)*(-2./(dz * dz));

        t3 = 1./(k0*k0)*(drhodz/rho*2.*dz);

        t4 = 1./(k0*k0)*(1/(dz * dz));

        /* MATRICES N AND D ARE IDENTICAL EXCEPT FOR COEFFICIENT*/

        double c,d;
        for(int mp = 0; mp < Mp; ++mp){
            c = std::pow(std::sin((mp+1)*M_PI/(2*Mp+1)),2);
            d = std::pow(std::cos((mp+1)*M_PI/(2*Mp+1)),2);
            mat_entries_N[mp].push_back(trip(i,i,(c*(t1 + t2))));
            mat_entries_D[mp].push_back(trip(i,i,(d*(t1 + t2))));
            if(i < nz - 1){
                mat_entries_N[mp].push_back(trip(i,i+1,(c*(-t3 + t4))));
                mat_entries_D[mp].push_back(trip(i,i+1,(d*(-t3 + t4))));
            }
            if(i > 0){
                mat_entries_N[mp].push_back(trip(i,i-1,(c*(t3 + t4))));
                mat_entries_D[mp].push_back(trip(i,i-1,(d*(t3 + t4))));
            }
        }
    }

    SpMat N(nz,nz), D(nz,nz);

    SpMat identity(nz, nz);
    std::vector<trip> idcoeffs;
    for(int i = 0; i < nz; ++i){
        idcoeffs.push_back(trip(i,i,1));
    }
    identity.setFromTriplets(idcoeffs.begin(), idcoeffs.end());

    SpMat temp(nz,nz);
    N = identity;
    D = identity;
    for(int mp = 0; mp < Mp; ++mp){
        temp.setFromTriplets(mat_entries_N[mp].begin(), mat_entries_N[mp].end());
        N = (temp*N).eval();
        temp.setFromTriplets(mat_entries_D[mp].begin(), mat_entries_D[mp].end());
        D = (temp*D).eval();
    }


    std::cout << (N*D - D*N).norm() << std::endl;

    return 0;
}

'''

Answer 1

问题是，如果没有有意义的参考值来定义非零值的期望大小顺序，就无法得出1e20是大还是小值的结论。

在您的情况下，矩阵N和D的范数分别约为1e20和1e18，并且N*D的范数约为1e38。假设double的相对精度约为1e-16，则与1e20相比，1e38的误差可以被认为是0。

总而言之，大多数情况下查看绝对误差是没有意义的。相反，您必须查看相对错误：

std::cout << (N*D - D*N).norm()/(N*D).norm() << std::endl;
为您提供1e-17。确实比double的数值精度小。