Search code examples
c++eigenmatrix-multiplicationeigen3

Eigen: sparse matrix pruned() does not remove entries below threshold


I want to use Eigen for sparse matrix multiplications, where in each iteration I want to remove all entries below a certain threshold. It seems to me that Eigen only removes elements exactly equal zero.

I am running Eigen 3.3.7, compiling with g++.

#include <Eigen/Sparse>
#include <Eigen/Dense>
#include <iostream>

using namespace Eigen;
typedef SparseMatrix<double> CscMat;            
typedef SparseMatrix<double,RowMajor> CsrMat;    

int N = 4;
CsrMat S, S2;

MatrixXd D(N, N), D2(N,N);
D << 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16;
D *= 0.1;

S = D.sparseView(0.5);
std::cout << D  << std::endl;
std::cout << S.nonZeros()  << std::endl;

D2 = D;
D2 = (D2.array() < 0.5).select(0, D2);
S2 = D2.sparseView();
std::cout << D  << std::endl;
std::cout << S2.nonZeros() << std::endl;

In the above S.nonzeros() returns 16, instead of expected 12 like in S2.nonzeros().

The output is:

0.1 0.2 0.3 0.4
0.5 0.6 0.7 0.8
0.9   1 1.1 1.2
1.3 1.4 1.5 1.6
16

  0   0   0   0
0.5 0.6 0.7 0.8
0.9   1 1.1 1.2
1.3 1.4 1.5 1.6
12

Solution

  • There is a second parameter to sparseView which is reference. In the end, the product of the two will determine the threshold, so you should use the line:

    S = D.sparseView(0.5, 1.0 - std::numeric_limits<double>::epsilon());
    

    to obtain what you seem to want.

    The actual code that does the evaluation is in MathFunctions.h

    static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y, 
                                         const RealScalar& prec)
    {
       return numext::abs(x) <= numext::abs(y) * prec;
    }
    

    where the default prec of type double is (currently) 1e-12.