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

Question

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

当 A = B 并取向量 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，与1e38.

相比，1e20的误差可以认为是0

总结一下，大多数时候看绝对误差是没有意义的。相反，您必须查看相对错误：

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

这给你大约 1e-17。这确实比 double.

的数值精度小

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

Sparse matrix multiplication in Eigen giving wrong result?

sparse-matrix

eigen

eigen3