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

【问题标题】：Rounding numbers below a certain threshold to zero in Eigen在 Eigen 中将低于某个阈值的数字四舍五入为零
【发布时间】：2019-02-03 14:22:18
【问题描述】：

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

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

有没有更有效的方法来使用Eigen 库？

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

【问题讨论】：

标签： c++ matrix eigen

【解决方案1】：

写你想要的最短的方法是

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);

应该做的伎俩（cf：Why does changing 0.1f to 0 slow down performance by 10x?）

【讨论】：

感谢您的提示和信息。我知道减速低于正常的数字正在产生，这也是我明确避免小于 1e-15 的数字的部分原因。因为它在正常范围内接壤，它们对我的计算没有用。
我的意思是，您可能可以省去手动将值刷新为零的工作，而只需依靠 CPU 的自动功能。
也谢谢你，非常感谢。对不起，如果我的评论听起来有点简洁。英语不是我的母语。我试图说“你是对的，我也在努力防止不正常的浮动，这是正确的！”。实际上，我现在正在寻找以跨平台方式启用它的方法。
看看你的godbolt-demo，.select() 版本是由GCC自动矢量化的，但不是clang。

【解决方案2】：

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

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

【讨论】：