Sparse matrix multiplication in Eigen giving wrong result?

Question

I am using Eigen in a project of mine, and I am running into a strange issue. I have complex sparse matrices A and B (1500x1500 or larger), and am multiplying them together with coefficients.

When A = B, and taking vector x of ones, I expect that

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

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

etc. and I do get this result for all these cases. (A.isApprox(B) evaluates to 1 and (A-B).norm() = 0).

However, when I multiply the matrices by doubles, as in

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

I get a nonzero result, which doesn't make sense to me, as scalars should commute with the matrices. In fact, if I do,

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

I get zero. Any time the coefficients are interspersed in the matrix manipulation, I get a nonzero result.

I am not using any compiler optimizations, etc.

What am I doing wrong here?

Edit: I have worked up a simple example. Maybe I'm missing something dumb, but here it is. This gives me an error of 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;
}

'''

These might just be some small numerical errors due to the usual accuracy limits of floating point arithmetics, What does "nonzero result" mean here quantitatively? If the norm of the result is much smaller than one, i.e., nominally zero, then there is nothing to worry about. — RHertel
By nonzero, I mean that some entries come back as greater than 0.1. I will check for numerical errors (and maybe it’s just that there are so many small errors that compound). Thanks for the suggestion. — ar24
Hi @RHertel, I have included an example in my edit. If you change the integer Mp, the result changes quite a bit. The only differences between the matrices N and D are the coefficients, i.e. N = c_1*c_2*...*c_MpX ^ Mp, and D = d_1*d_2...*d_Mp*X^Mp. — ar24

ggael ggael · Accepted Answer · 2019-11-18T15:24:29

The problem is that without a meaningful reference value defining what is the expected order of magnitude of a non-zero value, it is impossible to conclude whether 1e20 is a huge or a tiny value.

In your case, the norm of the matrices N and D are about 1e20 and 1e18 respectively, and the norm of N*D is about 1e38. Given that the relative precision of double is about 1e-16, an error of 1e20 can be considered as 0 compared to 1e38.

To summarize, it is most of the time meaningless to look at the absolute error. Instead, you have to look at the relative error:

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

which gives you about 1e-17. This is indeed smaller that the numerical precision of double.

Sparse matrix multiplication in Eigen giving wrong result?

1 Answers