本征:稀疏矩阵 pruned() 不会删除低于阈值的条目

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

我想将 Eigen 用于稀疏矩阵乘法,在每次迭代中我想删除低于特定阈值的所有条目。在我看来,Eigen 只删除恰好等于零的元素。

我是运行 Eigen 3.3.7,用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;

在上面的 S.nonzeros() returns 16 中,而不是像 S2.nonzeros() 中预期的 12。

输出为:

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

sparseView 还有第二个参数 reference。最后,两者的乘积将决定阈值,因此您应该使用以下行:

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

得到你似乎想要的。

执行评估的实际代码在 MathFunctions.h

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

其中双精度类型的默认 prec 是(当前)1e-12

如果您阅读 sparseViewdoc,您会发现第一个参数不是绝对阈值,而是 参考 非零(或者如果你愿意的话,也可以是预期的震级)。然后第二个可选参数是相对阈值。这与 pruned() 的逻辑相同。如果你想要一个绝对阈值,那么你可以这样做:

S = D.sparseView(1,0.5);
S = D.sparseView(0.5,1);