C++ vs python numpy 复杂数组的性能

Posted 2023-02-23

【中文标题】C++ vs python numpy 复杂数组的性能

【英文标题】：C++ vs python numpy complex arrays performance

【发布时间】：2020-04-09 19:12:06

【问题描述】：

谁能告诉我为什么这两个程序在运行时间上有巨大的差异？我只是将两个大型复杂数组相乘并比较 python (numpy) 和 c++ 中的时间。我正在使用带有 g++ 的 -O3 标志来编译这个 C++ 代码。我发现只有当我在 C++ 中使用 complex 浮点数时才会出现巨大的差异，它在 numpy 中的速度要快 20 倍以上。

python 代码：

import numpy as np
import time


if __name__ == "__main__":

    # check the data type is the same
    a = np.zeros((1), dtype=np.complex128)
    a[0] = np.complex(3.4e38,3.5e38)
    print(a)
    b = np.zeros((1), dtype=np.complex64)
    b[0] = np.complex(3.4e38,3.5e38)
    print(b)  # imaginary part is infinity

    length = 5000;
    A = np.ones((length), dtype=np.complex64) * np.complex(1,1)
    B = np.ones((length), dtype=np.complex64) * np.complex(1,0)

    num_iterations = 1000000
    time1 = time.time()
    for _ in range(num_iterations):
        A *= B

    time2 = time.time()
    duration = ((time2 - time1)*1e6)/num_iterations
    print(duration)

C++ 代码：

#include <iostream>
#include <complex>
#include <chrono>
using namespace std::chrono;
using namespace std;

int main()


  // check the data type is the same
  complex<double> a = complex<double>(3.4e38, 3.5e38);
  cout << a << endl;
  complex<float> b = complex<float>(3.4e38, 3.5e38);
  cout << b << endl;  // imaginary part is infinity

  const int length = 5000;
  static complex<float> A[length];
  static complex<float> B[length];

  for(int i=0; i < length; i++) 
    A[i] = complex<float>(1,1);
    B[i] = complex<float>(1,0);
  

  int num_iterations = 1000000;
  auto time1 = high_resolution_clock::now();
  for(int k=0; k < num_iterations; k++)
    for(int i=0; i < length; i++)
      A[i] *= B[i];

  auto time2 = high_resolution_clock::now();

  auto duration = duration_cast<microseconds>(time2 - time1);
  cout << "average time:" << duration.count() / num_iterations << endl;

【问题讨论】：

你的实现应该尽可能的相同。 C++ 版本使用float，通常是 32 位。您的 Python 使用 64 位复数。尝试使用std::complex<double>。

如果您将阵列容量设为const，您可以执行以下操作：static std::complex<double> A[length]; 与阵列 B 相同。

您的delete 错误。由于您分配了一个数组，因此语法应该是delete [] A，类似于B。

您需要在计时周期内多次运行测试。 C++ 时钟的分辨率不足以为一次执行计时。此外，时间不够准确也不够一致。我建议至少运行 1E06 次测试。这应该允许时间不准确以及其他开销（例如可能的任务切换）。在互联网上搜索“C++ 基准测试”以获取有关准确分析或基准测试的更多信息。

谢谢Thomas，我将代码更改为响应您的cmets，但仍然存在较大的性能差异。 numpy complex64 数字与 c++ complex 等价，因为两者都是两个 32 位浮点数的组合，如代码所示。我现在正在平均超过 1000 次迭代的时间。

【参考方案1】：

C++ 编译器正在为您做一些额外的检查，以便正确处理 NaN 和其他此类“标准”行为。 如果您添加 -ffast-math 优化标志，您将获得更理智的速度，但更少的“标准”行为。例如complex<float>(inf,0)*complex<float>(inf,0) 不会被评估为 complex<float>(inf,0)。你真的在乎吗？

numpy 正在做有意义的事，不受狭隘地阅读 C++ 标准的阻碍。

例如直到 非常 最近的 g++ 版本，除非使用 -ffast-math，否则以下函数中的后者要快得多。

complex<float> mul1( complex<float> a,complex<float> b)

    return a*b;


complex<float> mul2( complex<float> a,complex<float> b)

    float * fa = reinterpret_cast<float*>(&a);
    const float * fb = reinterpret_cast<float*>(&b);
    float cr = fa[0]*fb[0] - fa[1]*fb[1];
    float ci = fa[0]*fb[1] + fa[1]*fb[0];
    return complex<float>(cr,ci);

您可以在 https://godbolt.org/z/kXPgCh 上尝试此程序集输出以及前一个函数如何默认调用 __mulsc3

附：准备好对 C++ 标准对std::complex<T> 的另一波愤怒了吗？你能猜到std::norm必须默认实现吗？一起玩。点击链接，花十秒钟思考一下。

剧透：它可能是使用 sqrt 然后对其进行平方。