时间戳

首页 > 解决方案 > 正确有效地返回一个 Eigen::Matrix （带有惰性求值）

c++ - 正确有效地返回一个 Eigen::Matrix （带有惰性求值）

问题描述

正如docs 中所解释的auto，在使用带有 eigen 的类型推导时应该注意。例如，以下函数对我来说是段错误，我可以接受：

struct MultiplyAdd{
  Eigen::Matrix<double, DIMS, DIMS> mult;
  Eigen::Matrix<double, DIMS, 1> add;
  auto counterExample(Eigen::Matrix<double, DIMS, -1> const& data) const {
    return (mult.sinh()*data).colwise() + add;
  }
};

据我了解上面的文档，我可以将返回类型拼写为实际矩阵（而不是惰性求值表达式）或使用eval().

#include <benchmark/benchmark.h>
#include <Eigen/Dense>

constexpr int DIMS = 10;

struct MultiplyAdd{
  Eigen::Matrix<double, DIMS, DIMS> mult;
  Eigen::Matrix<double, DIMS, 1> add;
  Eigen::Matrix<double, DIMS, -1> transform(Eigen::Matrix<double, DIMS, -1> const& data) const {
    return ((mult*data).colwise() + add).eval();
  }
};


struct classic_worker{
  static auto dowork(Eigen::Matrix<double, DIMS, -1> const& data, MultiplyAdd const& trafo) {
    const auto flags = ((2.*trafo.transform(data)).row(2).array() > 1.7f).eval();
    return flags;
  }
};

struct allinline{
  static auto dowork(Eigen::Matrix<double, DIMS, -1> const& data, MultiplyAdd const& trafo) {
    const auto flags = ((2.*((trafo.mult * data).colwise() + trafo.add)).row(2).array() > 1.7f).eval();
    return flags;
  }
};

template <typename STRUCT>
static void quicktest(benchmark::State &state) {
  Eigen::Matrix<double, DIMS, -1> data = Eigen::MatrixXd::Random(DIMS, state.range(0));

  Eigen::Matrix<double, DIMS, DIMS> m = Eigen::MatrixXd::Random(DIMS, DIMS);
  Eigen::Matrix<double, DIMS, 1> a = Eigen::MatrixXd::Random(DIMS, 1);

  MultiplyAdd trafo{m,a};

  for (auto _ : state) {
    benchmark::DoNotOptimize(STRUCT::dowork(data, trafo));
  }
}

// clang-format off
BENCHMARK_TEMPLATE(quicktest, classic_worker               )->UseRealTime()->DenseRange(20,320, 50);
BENCHMARK_TEMPLATE(quicktest, allinline                    )->UseRealTime()->DenseRange(20,320, 50);
// clang-format on

BENCHMARK_MAIN();

在这个较长的示例中，我有两个我的实现MultiplyAdd，并将它们与 google-benchmark 进行比较。在一种情况下transform，我使用可以放入中央标题的方法进行操作，在另一种情况下，我在“调用站点”上进行所有低级数学运算。在示例中，我仅使用计算的第 2 行。它出现

./iter_5 --benchmark_filter='.*(classic|all).*320.*' --benchmark_repetitions=10 --benchmark_report_aggregates_only
...
-----------------------------------------------------------------------------------------
Benchmark                                               Time             CPU   Iterations
-----------------------------------------------------------------------------------------
quicktest<classic_worker>/320/real_time_mean         3974 ns         3974 ns           10
quicktest<classic_worker>/320/real_time_median       3944 ns         3944 ns           10
quicktest<classic_worker>/320/real_time_stddev        121 ns          121 ns           10
quicktest<allinline>/320/real_time_mean              3308 ns         3308 ns           10
quicktest<allinline>/320/real_time_median            3285 ns         3285 ns           10
quicktest<allinline>/320/real_time_stddev            68.2 ns         68.2 ns           10

all-at-call-site 版本的速度提高了约 20%（包含所有微基准测试的警告）。

鉴于上述链接，我的天真解释是我正在避免惰性评估并不必要地计算所有行。（速度差异似乎与我在 3 到 100 之间的选择无关DIMS，所以我怀疑我的解释是错误的。）

编写返回特征矩阵类型的函数的正确模式是什么（在文档中我只看到提到函数参数）？

编辑：我有一个新的例子，在观察上更接近我的真实案例

#include <benchmark/benchmark.h>
#include <Eigen/Dense>

struct Trafo{
  Eigen::Quaternion<double> mult;
  Eigen::Matrix<double, 3, 1> add;

  inline Eigen::Matrix<double, 3, -1> transform(Eigen::Matrix<double, 3, -1> const& data) const noexcept {
    return ((mult.toRotationMatrix()*data).colwise() + add).eval();
  }
  inline auto counterExample(Eigen::Matrix<double, 3, -1> const& data) const noexcept {
    return (mult.toRotationMatrix()*data).colwise() + add;
  }
};


struct classic_worker{
  static auto dowork(Eigen::Matrix<double, 3, -1> const& data, Trafo const& trafo) {
    const auto flags = ((2.*trafo.transform(data)).row(2).array() > 1.7f).eval();
    return flags;
  }
};

struct counterExample_worker{
  static auto dowork(Eigen::Matrix<double, 3, -1> const& data, Trafo const& trafo) {
    const auto flags = ((2.*trafo.counterExample(data)).row(2).array() > 1.7f).eval();
    return flags;
  }
};

struct allinline{
  static auto dowork(Eigen::Matrix<double, 3, -1> const& data, Trafo const& trafo) {
    const auto flags = ((2.*((trafo.mult.toRotationMatrix() * data).colwise() + trafo.add)).row(2).array() > 1.7f).eval();
    return flags;
  }
};

template <typename STRUCT>
static void quicktest(benchmark::State &state) {
  Eigen::Matrix<double, 3, -1> data = Eigen::MatrixXd::Random(3, state.range(0));

  Eigen::Matrix<double, 4, 1> random = Eigen::MatrixXd::Random(4, 1);
  Eigen::Quaternion<double> m{random(0), random(1), random(2), random(3)};
  Eigen::Matrix<double, 3, 1> a = Eigen::MatrixXd::Random(3, 1);

  Trafo trafo{m,a};

  for (auto _ : state) {
    benchmark::DoNotOptimize(STRUCT::dowork(data, trafo));
  }
}


// clang-format off
BENCHMARK_TEMPLATE(quicktest, counterExample_worker        )->UseRealTime()->DenseRange(20,320, 50);
BENCHMARK_TEMPLATE(quicktest, classic_worker               )->UseRealTime()->DenseRange(20,320, 50);
BENCHMARK_TEMPLATE(quicktest, allinline                    )->UseRealTime()->DenseRange(20,320, 50);
// clang-format on

BENCHMARK_MAIN();

有三种方法可以进行我的计算：

用 g++8 编译并且-O3 -march=native（并且都在同一个翻译单元中，所以编译器可以做它想要的所有内联）我看到下面的时序表

quicktest<counterExample_worker>/320/real_time_mean        52612 ns        52612 ns           10
quicktest<counterExample_worker>/320/real_time_median      55218 ns        55217 ns           10
quicktest<counterExample_worker>/320/real_time_stddev       8501 ns         8501 ns           10
quicktest<classic_worker>/320/real_time_mean                 622 ns          622 ns           10
quicktest<classic_worker>/320/real_time_median               619 ns          619 ns           10
quicktest<classic_worker>/320/real_time_stddev              6.89 ns         6.89 ns           10
quicktest<allinline>/320/real_time_mean                      428 ns          428 ns           10
quicktest<allinline>/320/real_time_median                    426 ns          426 ns           10
quicktest<allinline>/320/real_time_stddev                   5.68 ns         5.67 ns           10

看来，eval在这里不打电话是个坏主意。然而，在它自己的功能中进行转换也需要将所有内容写在一行中。所以问题是：有没有办法在不降低运行时减速的情况下将这里的转换放入单独的函数中？我确实注意到（不幸的是）此示例中的结果在很大程度上取决于编译器和编译标志（andg++-8 -msse3之间的差异消失了，似乎提供了更快的 with和 with ...）allinlineclassic_workerclang++-10allinline-msse3-march=native

标签： c++eigengoogle-benchmark

解决方案

推荐阅读