与 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-loop（可能还有很多其他原因），使我的代码变得幼稚，例如它的运行时非常微弱到 Python 的那个。我也在Is there something like numpy.logspace in C++? 看到了建议。但是，没有明显的差异。那么，如何修改我的代码或编写一个与 python 版本相媲美的新代码？

Answer 1

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

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

#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();
}

Answer 2

至少可以对所示代码进行三个明显的优化。

1) 从logspace 返回时，以 C++17 模式编译以保证复制省略。

2)

std::vector<double> logspace;
for (int i = 0; i < k; i++)

使用logspace.reserve() 预分配向量，以避免在填充此向量时进行无用的重复重新分配。

3)

void logspace_print (std::vector<double> logspace)

此处按值传递会创建向量的完整副本，但没有任何用处。更改此函数，使其通过引用获取logspace 参数。

有一种可能的微优化可能会或可能不会产生任何影响：

logspace.push_back(pow(10, i * (b - a) / (k - 1)));

这个公式的“(b-a)/(k-1)”部分是常数，可以展开出循环。不过，我希望编译器自行完成，这是一个相当基本的优化。