在 Eigen 中将低于某个阈值的数字四舍五入为零

Question

我在我已经开发了一段时间的科学应用程序中广泛使用 Eigen。由于我正在实施数值方法，因此低于某个阈值（例如 1e-15）的数字不是兴趣点，并且会减慢计算速度并增加错误率。

因此，我想将低于该阈值的数字四舍五入为 0。我可以用 for 循环来完成它，但是用 for-if 循环锤打多个相对较大的矩阵（每个矩阵 2M 单元或以上）是昂贵的并且因为我需要而减慢了我的速度多做几次。

使用 Eigen 库是否有更有效的方法？

换句话说，我试图在我的计算管道中消除低于特定阈值的数字。

Answer 1

Eigen 有一个名为 UnaryExpr 的方法，它将给定的函数指针应用于矩阵中的每个系数（它也有稀疏和数组变体）。

将测试其性能并相应地更新此答案。

Answer 2

写你想要的最短的方法是

void foo(Eigen::VectorXf& inout, float threshold)
{
    inout = (threshold < inout.array().abs()).select(inout, 0.0f);
}

但是，比较和 select 方法都没有被 Eigen (as of now) 向量化。

如果速度很重要，您需要编写一些手动 SIMD 代码，或者编写支持 packet 方法的自定义仿函数（这使用 Eigen 的内部功能，因此不能保证稳定!):

template<typename Scalar> struct threshold_op {
  Scalar threshold;
  threshold_op(const Scalar& value) : threshold(value) {}
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const{
    return threshold < std::abs(a) ? a : Scalar(0); }
  template<typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const {
    using namespace Eigen::internal;
    return pand(pcmp_lt(pset1<Packet>(threshold),pabs(a)), a);
  }
};
namespace Eigen { namespace internal {

template<typename Scalar>
struct functor_traits<threshold_op<Scalar> >
{ enum {
    Cost = 3*NumTraits<Scalar>::AddCost,
    PacketAccess = packet_traits<Scalar>::HasAbs };
};

}}

然后可以将其传递给 unaryExpr:

inout = inout.unaryExpr(threshold_op<float>(threshold));

Godbolt-Demo（应该与SSE/AVX/AVX512/NEON/...一起工作）：https://godbolt.org/z/bslATI

这实际上可能是您减速的唯一原因是低于正常数字。在那种情况下，一个简单的

_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);

应该可以解决问题（参见：Why does changing 0.1f to 0 slow down performance by 10x?）

在 Eigen 中将低于某个阈值的数字四舍五入为零

Rounding numbers below a certain threshold to zero in Eigen

c++

matrix

eigen