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

Posted 2023-02-23

【中文标题】与 numpy.logspace() 函数相比，如何用 C++ 编写快速代码？

【英文标题】：How to write a fast code in C++ compared to numpy.logspace() function?

【发布时间】：2019-09-15 12:15:01

【问题描述】：

这是在 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 版本相媲美的新代码？

【问题讨论】：

你可以一直看numpy的源码吗？

当然。为什么不。但是，正如您也可以看到它的来源，有一系列相互关联的函数需要寻找。

【参考方案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();

【讨论】：

Printf 只有在同步时才比 cout 快（默认情况下是这样）。您可以使用 std::ios::sync_with_stdio(false);禁用它。

"打印通常很昂贵。"打印到文件很昂贵。打印到 终端 需要很长时间。这就是为什么我敢打赌 Python 和 C++ 的性能相当的原因。将输出通过管道传输到文件中，以便更好地比较语言本身。现在他可能正在分析终端。

我没有看过我的答案中的印刷。已经太久了。

@Unapiedra，我很享受您的好回答（此外，感谢您推荐 Google 基准测试。听起来很棒）。

@Unapiedra & @Mooing Duck，我理解在运行时打印的真正糟糕的性能（以及效果）（更准确地说，Python 正在使用一种狡猾的技巧来打印。确实，它不会打印所有输出；））。无论如何，我有一个关于提高for-loop 性能的想法。在尝试之前，我想研究一下它的理论：将区间[0, 10000000]划分为（例如）[0, 1000] U [1000, 10000] U ... U [1000000, 10000000]并对其应用for循环是否是一种好方法？

【参考方案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)”部分是常数，可以展开出循环。不过，我希望编译器自行完成，这是一个相当基本的优化。

【讨论】：

我应用了 2)、3) 和你的微优化建议（不幸的是，我还不能打开 c++17-compiler）。但是，差异小于 1 秒（我在 Linux 中使用time 命令来比较运行时间）。

我会给 +1，除非你忽略了终端的性能。