为什么我的python numpy代码比c ++快？

Question

有人可以告诉我为什么这个Python Numpy代码:

import numpy as np
import time

k_max = 40000
N = 10000

data = np.zeros((2,N))
coefs = np.zeros((k_max,2),dtype=float)

t1 = time.time()
for k in xrange(1,k_max+1):
    cos_k = np.cos(k*data[0,:])
    sin_k = np.sin(k*data[0,:])
    coefs[k-1,0] = (data[1,-1]-data[1,0]) + np.sum(data[1,:-1]*(cos_k[:-1] - cos_k[1:]))
    coefs[k-1,1] = np.sum(data[1,:-1]*(sin_k[:-1] - sin_k[1:]))
t2 = time.time()

print('Time:')
print(t2-t1)

比这个C++代码更快:

#include <cstdio>
#include <iostream>
#include <cmath>
#include <time.h>

using namespace std;

// consts
const unsigned int k_max = 40000;
const unsigned int N = 10000;

int main()
{
    time_t start, stop;
    double diff;
    // table with data
    double data1[ N ];
    double data2[ N ];
    // table of results
    double coefs1[ k_max ];
    double coefs2[ k_max ];
    // main loop
    time( & start );
    for( unsigned int j = 1; j<N; j++ )
    {
        for( unsigned int i = 0; i<k_max; i++ )
        {
            coefs1[ i ] += data2[ j-1 ]*(cos((i+1)*data1[ j-1 ]) - cos((i+1)*data1[ j ]));
            coefs2[ i ] += data2[ j-1 ]*(sin((i+1)*data1[ j-1 ]) - sin((i+1)*data1[ j ]));
        }
    }
    // end of main loop
    time( & stop );
    // speed result
    diff = difftime( stop, start );
    cout << "Time: " << diff << " seconds";
    return 0;
}

第一个显示:"时间:8秒",而第二个显示:"时间:11秒"

我知道numpy是用C语言编写的,但我仍然认为C++示例会更快.我错过了什么吗？有没有办法改进C++代码(或python代码)？预先感谢您的帮助!

编辑: 我已经按照其中一条评论中的建议更改了C++代码(动态表到静态表).C++代码现在更快,但仍然比Python版本慢得多.

编辑2: 我已经从调试模式更改为发布模式,并将'k'从4000增加到40000.现在numpy稍微快一些(8秒到11秒).

Answer 1

我发现这个问题很有意思,因为每次我遇到类似numpy速度的话题时(与c/c ++相比)总会有回答,例如"它是一个薄的包装器,它的核心是用c编写的,所以它是脂肪",但这个没有解释为什么c应该比带有附加层(甚至是薄层)的c慢.

答案是:正确编译时,你的c ++代码并不比你的python代码慢.

我做了一些基准测试,起初看起来numpy的速度要快得多.但我忘了用gcc优化编译.

我再次计算了所有内容,并将结果与​​代码的纯c版本进行了比较.我使用的是gcc版本4.9.2和python2.7.9(使用相同的gcc从源代码编译).编译我用过的c ++代码g++ -O3 main.cpp -o main,编译我用过的c代码gcc -O3 main.c -lm -o main.在所有例子中,我data用一些数字(0.1,0.4)填充变量,因为它改变了结果.我也改变了np.arrays以使用doubles(dtype=np.float64),因为在c ++示例中有双打.我的纯c版代码(类似):

#include <math.h>
#include <stdio.h>
#include <time.h>

const int k_max = 100000;
const int N = 10000;

int main(void)
{
    clock_t t_start, t_end;
    double data1[N], data2[N], coefs1[k_max], coefs2[k_max], seconds;
    int z;
    for( z = 0; z < N; z++ )
    {
        data1[z] = 0.1;
        data2[z] = 0.4;
    }

    int i, j;
    t_start = clock();
    for( i = 0; i < k_max; i++ )
    {
        for( j = 0; j < N-1; j++ )
        {
            coefs1[i] += data2[j] * (cos((i+1) * data1[j]) - cos((i+1) * data1[j+1]));
            coefs2[i] += data2[j] * (sin((i+1) * data1[j]) - sin((i+1) * data1[j+1]));
        }
    }
    t_end = clock();

    seconds = (double)(t_end - t_start) / CLOCKS_PER_SEC;
    printf("Time: %f s\n", seconds);
    return coefs1[0];
}

对于k_max = 100000, N = 10000以下结果:

• python 70.284362 s
• c ++ 69.133199 s
• c 61.638186 s

Python和c ++基本上是同一时间,但请注意,有一个长度为k_max的python循环,与c/c ++相比应该慢得多.它是.

因为k_max = 1000000, N = 1000我们有:

• python 115.42766 s
• c ++ 70.781380 s

用于k_max = 1000000, N = 100:

• python 52.86826 s
• c ++ 7.050597 s

所以差异随着分数的增加而增加k_max/N,但是即使N比大得多k_max,python也不会更快,例如k_max = 100, N = 100000:

• python 0.651587 s
• c ++ 0.568518 s

显然,c/c ++和python之间的主要速度差异在for循环中.但是我想找出numpy和c中数组上简单操作之间的区别.在代码中使用numpy的优点包括:1.将整个数组乘以数字,2.计算整个数组的sin/cos,3.对数组的所有元素求和,而不是分别对每个单独的项进行这些操作.所以我准备了两个脚本来比较这些操作.

Python脚本:

import numpy as np
from time import time

N = 10000
x_len = 100000

def main():
    x = np.ones(x_len, dtype=np.float64) * 1.2345

    start = time()
    for i in xrange(N):
        y1 = np.cos(x, dtype=np.float64)
    end = time()
    print('cos: {} s'.format(end-start))

    start = time()
    for i in xrange(N):
        y2 = x * 7.9463
    end = time()
    print('multi: {} s'.format(end-start))

    start = time()
    for i in xrange(N):
        res = np.sum(x, dtype=np.float64)
    end = time()
    print('sum: {} s'.format(end-start))

    return y1, y2, res

if __name__ == '__main__':
    main()

# results
# cos: 22.7199969292 s
# multi: 0.841291189194 s
# sum: 1.15971088409 s

C脚本:

#include <math.h>
#include <stdio.h>
#include <time.h>

const int N = 10000;
const int x_len = 100000;

int main()
{
    clock_t t_start, t_end;
    double x[x_len], y1[x_len], y2[x_len], res, time;
    int i, j;
    for( i = 0; i < x_len; i++ )
    {
        x[i] = 1.2345;
    }

    t_start = clock();
    for( j = 0; j < N; j++ )
    {
        for( i = 0; i < x_len; i++ )
        {
            y1[i] = cos(x[i]);
        }
    }
    t_end = clock();
    time = (double)(t_end - t_start) / CLOCKS_PER_SEC;
    printf("cos: %f s\n", time);

    t_start = clock();
    for( j = 0; j < N; j++ )
    {
        for( i = 0; i < x_len; i++ )
        {
            y2[i] = x[i] * 7.9463;
        }
    }
    t_end = clock();
    time = (double)(t_end - t_start) / CLOCKS_PER_SEC;
    printf("multi: %f s\n", time);

    t_start = clock();
    for( j = 0; j < N; j++ )
    {
        res = 0.0;
        for( i = 0; i < x_len; i++ )
        {
            res += x[i];
        }
    }
    t_end = clock();
    time = (double)(t_end - t_start) / CLOCKS_PER_SEC;
    printf("sum: %f s\n", time);

    return y1[0], y2[0], res;
}

// results
// cos: 20.910590 s
// multi: 0.633281 s
// sum: 1.153001 s

Python结果:

• cos:22.7199969292 s
• 多:0.841291189194 s
• 总和:1.15971088409 s

C结果:

• cos:20.910590 s
• 多:0.633281 s
• 总和:1.153001

你可以看到numpy非常快,但总是比纯c慢一点.

Answer 2

事实上，我很惊讶没有人提到像 BLAS LAPACK MKL 这样的线性代数库和所有......

Numpy 使用复杂的线性代数库！ 本质上，Numpy 大多数时候都不是基于纯 c/cpp/fortran 代码构建的……它实际上是基于复杂的库构建的，这些库利用最高性能的算法和思想来优化代码。这些复杂的库很难与经典线性代数计算的简单实现相匹配。最简单的第一个改进示例是阻塞技巧。

我从 ETH 的 CSE 实验室获取了下图，他们在其中比较了不同实现的矩阵向量乘法。y 轴表示计算强度（以 GFLOP 为单位）；长话短说，就是计算完成的速度。x 轴是矩阵的维度。

C 和 C++ 是快速语言，但实际上，如果您想模仿这些库的速度，您可能需要更深入地使用 Fortran 或内在函数指令（这可能是您可以在 C++ 中执行的最接近汇编代码的指令） 。

考虑问题Benchmarking (python vs. c++ using BLAS) and (numpy)，其中 @Jfs 提供了非常好的答案，我们观察到：“在我的机器上，C++ 和 numpy 之间没有区别。”

还有一些参考：

为什么简单的 C++ 矩阵乘法比 BLAS 慢 100 倍？

Answer 3

在我的计算机上，您的（当前）Python 代码运行时间为 14.82 秒（是的，我的计算机很慢）。

我将你的 C++ 代码重写为我认为合理的东西（基本上，我几乎忽略了你的 C++ 代码，只是将你的 Python 重写为 C++。这给了我这个：

#include <cstdio>
#include <iostream>
#include <cmath>
#include <chrono>
#include <vector>
#include <assert.h>

const unsigned int k_max = 40000;
const unsigned int N = 10000;

template <class T>
class matrix2 {
    std::vector<T> data;
    size_t cols;
    size_t rows;
public:
    matrix2(size_t y, size_t x) : cols(x), rows(y), data(x*y) {}
    T &operator()(size_t y, size_t x) {
        assert(x <= cols);
        assert(y <= rows);
        return data[y*cols + x];
    }

    T operator()(size_t y, size_t x) const {
        assert(x <= cols);
        assert(y <= rows);
        return data[y*cols + x];
    }
};

int main() {
    matrix2<double> data(N, 2);
    matrix2<double> coeffs(k_max, 2);

    using namespace std::chrono;

    auto start = high_resolution_clock::now();

    for (int k = 0; k < k_max; k++) {
        for (int j = 0; j < N - 1; j++) {
            coeffs(k, 0) += data(j, 1) * (cos((k + 1)*data(j, 0)) - cos((k + 1)*data(j+1, 0)));
            coeffs(k, 1) += data(j, 1) * (sin((k + 1)*data(j, 0)) - sin((k + 1)*data(j+1, 0)));
        }
    }

    auto end = high_resolution_clock::now();
    std::cout << duration_cast<milliseconds>(end - start).count() << " ms\n";
}

这运行了大约 14.4 秒，所以它比 Python 版本略有改进——但考虑到 Python 主要是一些 C 代码的非常薄的包装器，因此只获得轻微的改进几乎是我们应该期望的。

下一个明显的步骤是使用多个内核。要在 C++ 中做到这一点，我们可以添加以下行：

#pragma omp parallel for

...在外for循环之前：

#pragma omp parallel for
for (int k = 0; k < k_max; k++) {
    for (int j = 0; j < N - 1; j++) {
        coeffs(k, 0) += data(j, 1) * (cos((k + 1)*data(j, 0)) - cos((k + 1)*data(j+1, 0)));
        coeffs(k, 1) += data(j, 1) * (sin((k + 1)*data(j, 0)) - sin((k + 1)*data(j+1, 0)));
    }
}

随着-openmp加入到编译器的命令行（虽然确切的标志当然取决于你使用的编译器），该跑了约4.8秒。如果您有 4 个以上的内核，您可能会期待比这更大的改进（相反，如果您的内核少于 4 个，则预计会有较小的改进——但如今，4 个以上的内核更为常见）。