python - 与 numpy.logspace() 函数相比，如何用 C++ 编写更快的代码？

Question

这是 Python 中的代码，可以非常快速地生成对数空间值:

import numpy
print(numpy.logspace(0,1,num=10000000))

我尝试用 C++ 模拟其输出，如下所示:

#include <iostream>
#include <cmath>
#include <vector>
std::vector<double> logspace (const double &a, const double &b, const int &k)
{
    std::vector<double> logspace;
    for (int i = 0; i < k; i++)
    {
        logspace.push_back(pow(10, i * (b - a) / (k - 1)));
    }
    return logspace;
}
void logspace_print (std::vector<double> logspace)
{
    for (auto ls : logspace)
    {
        std::cout << ls << "\n";
    }
    std::cout << "\n";
}
int main ()
{
    std::vector<double> my_ls = logspace(0, 1, 10000000);
    logspace_print(my_ls);
}

放弃浮点运算，使用函数 pow(., .) 和 for 循环(可能还有很多其他原因)，使我的代码作为一个天真的人，例如它的运行时相对于 Python 的人来说非常微弱。我在 Is there something like numpy.logspace in C++? 看到了推荐还。但是，没有明显的区别。那么，我该如何修改我的代码或编写一个与 python 版本相当的新代码呢？

Answer 1

有趣的问题!我的答案在顶部有不同版本的功能。以下只是基准测试代码。使用 google-benchmark 作为库。

• 我的中间结果也可以在这里找到:1 Quick-Bench.com 通常是一个很棒的网站。
• 您没有说明是否要将打印到标准输出作为用例的一部分进行测量。打印通常很昂贵。您避免了 std::endl 的刷新，这很好!此外，printf 可能比 std::cout 更快。另请查看 fmtlib 2 .它快速且易于使用。
• 通常，Numpy 使用的方法是最快的。 (在我的版本中命名为 logspace_v3。)它包括首先运行 linspace，然后就地计算 10 的幂。
• 不过，我强烈地觉得我在这里遗漏了很多东西。使用适当的标志(-march=native -mtune=native 和 fast-math)矢量化应该开始。但我不相信它会。这是一些带矢量化的 Godbolt(第 590 行)3 .
• 禁食的是摆脱 pow 调用。请注意，这会累积浮点错误并导致不准确的结果。
• 次要:通过常量引用传递 double 或整数没有任何好处。

#include <algorithm>
#include <benchmark/benchmark.h>
#include <cmath>
#include <iostream>
#include <numeric>
#include <vector>

#include <gtest/gtest.h>

std::vector<double> logspace(double a, double b, int k) {
  std::vector<double> logspace;
  for (int i = 0; i < k; i++) {
    logspace.push_back(pow(10, i * (b - a) / (k - 1)));
  }
  return logspace;
}

// Pre-allocate the correct size using .reserve()
std::vector<double> logspace_v1(double a, double b, int k) {
  std::vector<double> logspace;
  logspace.reserve(k);
  for (int i = 0; i < k; i++) {
    logspace.push_back(pow(10, i * (b - a) / (k - 1)));
  }
  return logspace;
}

/// Manually extract the constant factor.
std::vector<double> logspace_v2(double a, double b, int k) {
  std::vector<double> logspace;
  logspace.reserve(k);
  const auto exp_scale = (b - a) / (k - 1);
  for (int i = 0; i < k; i++) {
    logspace.push_back(pow(10, i * exp_scale));
  }
  return logspace;
}

/// Copy the impl behavior of numpy.linspace: First linspace then power.
std::vector<double> logspace_v3(double a, double b, int k) {
  /*
  y = linspace(start, stop, num=num, endpoint=endpoint, axis=axis)
  if dtype is None:
      return _nx.power(base, y)
  return _nx.power(base, y).astype(dtype, copy=False)
  */
  const auto exp_scale = (b - a) / (k - 1);
  std::vector<double> logspace;
  logspace.reserve(k);
  for (int i = 0; i < k; i++) {
    logspace.push_back(i * exp_scale);
  }
  std::for_each(logspace.begin(), logspace.end(),
                [](double &x) { x = pow(10, x); });
  return logspace;
}

/// Improve on v3 by applying pow directly
std::vector<double> logspace_v4(double a, double b, int k) {
  const auto exp_scale = (b - a) / (k - 1);
  std::vector<double> logspace(k, 0.);
  std::generate(logspace.begin(), logspace.end(),
                [n = -1, exp_scale]() mutable {
                  n++;
                  return pow(10, n * exp_scale);
                });
  return logspace;
}

/// Use generate_n : First linspace then power.
std::vector<double> logspace_v5(double a, double b, int k) {
  const auto exp_scale = (b - a) / (k - 1);
  std::vector<double> logspace(k, 0.);
  std::iota(logspace.begin(), logspace.end(), 0);
  std::for_each(logspace.begin(), logspace.end(),
                [exp_scale](double &x) { x *= exp_scale; });
  std::for_each(logspace.begin(), logspace.end(),
                [](double &x) { x = pow(10, x); });
  return logspace;
}

std::vector<double> logspace_v6(double a, double b, int k) {
  const auto exp_scale = (b - a) / (k - 1);
  const auto factor = pow(10, exp_scale);
  std::vector<double> logspace;
  logspace.reserve(k);

  // val = pow(b, i * exp_scale);
  // = pow(pow(b, exp_scale), i);
  // = pow(f, i); with f := pow(b, exp_scale);
  // next = cur * f;
  // first = pow(b, a);

  double val = pow(10, a);
  for (int i = 0; i < k; i++) {
    logspace.push_back(val);
    val *= factor;
  }
  return logspace;
}

template <std::vector<double> (*F)(double, double, int)>
static void LogspaceBench(benchmark::State &state) {
  for (auto _ : state) {
    benchmark::DoNotOptimize(F(0, 1, state.range(0)));
  }
}
BENCHMARK_TEMPLATE(LogspaceBench, logspace)->Arg(1000);
BENCHMARK_TEMPLATE(LogspaceBench, logspace_v1)->Arg(1000);
BENCHMARK_TEMPLATE(LogspaceBench, logspace_v2)->Arg(1000);
BENCHMARK_TEMPLATE(LogspaceBench, logspace_v3)->Arg(1000)->Arg(10000000);
BENCHMARK_TEMPLATE(LogspaceBench, logspace_v4)->Arg(1000);
BENCHMARK_TEMPLATE(LogspaceBench, logspace_v5)->Arg(1000);
BENCHMARK_TEMPLATE(LogspaceBench, logspace_v6)->Arg(1000)->Arg(10000000);

class LogspaceTest
    : public testing::TestWithParam<
          std::function<std::vector<double>(double, double, int)>> {};

TEST_P(LogspaceTest, IsSame) {
  auto func = GetParam();
  const auto actual = func(0, 1., 1000);
  const auto expected = logspace(0., 1., 1000);
  // TODO: Buggy with (3, 70, 1000) and (0, 1, 1000)
  ASSERT_EQ(expected.size(), actual.size());
  for (int i = 0; i < expected.size(); i++) {
    ASSERT_DOUBLE_EQ(actual[i], expected[i]) << i;
  }
}

INSTANTIATE_TEST_SUITE_P(InstantiationName, LogspaceTest,
                         testing::Values(logspace, logspace_v1, logspace_v2,
                                         logspace_v3, logspace_v4, logspace_v5,
                                         logspace_v6));

int main(int argc, char **argv) {
  ::benchmark::Initialize(&argc, argv);
  ::benchmark::RunSpecifiedBenchmarks();
  ::testing::InitGoogleTest(&argc, argv);
  return RUN_ALL_TESTS();
}