pybind11 vs numpy 用于矩阵乘积

Posted 2023-02-23

【中文标题】pybind11 vs numpy 用于矩阵乘积

【英文标题】：pybind11 vs numpy for a matrix product

【发布时间】：2018-06-15 15:44:45

【问题描述】：

编辑 2（2018 年 6 月 18 日。）

我使用了Matrix 中提出的类

http://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html

使用Matrix 产品，我实现如下：

Matrix product3(const Matrix &s1, const Matrix &s2) // M = M1 x M2

    size_t rowsM1 = s1.rows();
    size_t colsM1 = s1.cols();
    size_t rowsM2 = s2.rows();
    size_t colsM2 = s2.cols();
    assert(colsM1 == rowsM2);
    size_t resDim = rowsM1 * colsM2;
    double * ptr = new double[resDim];
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, rowsM1, colsM2, colsM1, 1.0, s1.data(), rowsM1, s2.data(), colsM2, 0.0, ptr, std::max(rowsM1, colsM2));
    Matrix res(rowsM1, colsM2, ptr);
    return res;

在发布时（在核心 i7 6700 HQ 上）结果如下：

确实比py::array_t<double> 的要好得多。图表（纵坐标为秒，横坐标为方阵大小）：

在 intel mkl 下，Numpy 相对来说有点小。在大小区域 [1500,1600] 中两者都有显着下降，mkl 更陡峭。可以注意到，“numpy time / intel time”因子随着矩阵大小的增加而减少。

这次在核心 i7-7700K 上：

测试python代码是：

import Binder
import numpy as np
import matplotlib.pyplot as plt
import time

rangeMin = 100
rangeMax = 2000
step = 100

X = []
intel = []
numpy = []

for size in range(rangeMin, rangeMax, step):

    X.append(size)

    m1 = np.array(np.random.rand(size,size), copy = False).astype(np.float64)
    M1 = Binder.Matrix(m1)
    m2 = np.array(np.random.rand(size,size), copy = False).astype(np.float64)
    M2 = Binder.Matrix(m2)

    M = Binder.Matrix(size,size)
    N = np.array([size,size])

    #M.print()

    loopSize = 50

    start_time = time.time()
    for x in range(1, loopSize):
        N = m1 @ m2
    time_elapsed = (time.time() - start_time)/loopSize

    print("Size =\t" + repr(size) + "\tnumpy Time =\t" + repr(time_elapsed))
    numpy.append(time_elapsed)

    start_time = time.time()
    for x in range(1, loopSize):
        M = Binder.product3(M1,M2);
    time_elapsed = (time.time() - start_time)/loopSize

    print("Size =\t" + repr(size) + "\tintel Time =\t" + repr(time_elapsed))
    intel.append(time_elapsed)

fig = plt.figure()
ax1 = fig.add_subplot(111)

ax1.scatter(X, numpy, s=10, c='b', marker="s", label='numpy')
ax1.scatter(X, intel, s=10, c='r', marker="o", label='intel')
plt.legend(loc='upper left');
plt.show()

---

编辑 1（2018 年 6 月 16 日。）

我也尝试过，这次是用 intel mkl，将初始代码的 for 循环替换为

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nbRows1, nbCols2, nbCols1, 1.0, ptr1, nbRows1, ptr2, nbCols2, 0.0, ptr, nbRows1);

初始代码在英特尔酷睿 i5 4570 上运行。这次在英特尔酷睿 i7 6700 HQ 上运行所有三个案例：

两句话：

1) 对于相同的 Python 3.6.5 32 位，在我的笔记本电脑的核心 i7 上使用 numpy 的 Python 比在我在工作中使用的旧核心 i5 桌面上要慢。朴素的 c++ 稍微快一点。很奇怪。

2) 在核心 i7 上，c++ intel mkl vs numpy 的因子为 3.41

---

初始问题

我写了这个非常幼稚的 c++ pybind11 代码：

py::array product1(py::array_t<double> m1, py::array_t<double> m2)

    py::buffer_info info1 = m1.request();
    double * ptr1 = static_cast<double *>(info1.ptr);

    py::buffer_info info2 = m2.request();
    double * ptr2 = static_cast<double *>(info2.ptr);

    unsigned int nbRows1 = info1.shape[0];
    unsigned int nbCols1 = info1.shape[1];

    unsigned int nbRows2 = info2.shape[0];
    unsigned int nbCols2 = info2.shape[1];

    assert(nbCols1 == nbRows2);

    int resDim = nbRows1 * nbCols2;

    double * ptr = new double[resDim];

    double localSum = 0.0;
    for (int i = 0 ; i < nbRows1; ++i)
    
        for (int j = 0 ; j < nbCols2; ++j)
        
            for (int l = 0; l < nbCols1; ++l)
            
                localSum += ptr1[nbCols1 * i + l] * ptr2[nbCols2 * l + j];
            
            ptr[nbCols2 * i + j] = localSum;
            localSum = 0.0;
        
    
    py::array_t<double> mRes = py::array_t<double>
                                (
                                    py::buffer_info
                                    (
                                        ptr,
                                        sizeof(double), //itemsize
                                        py::format_descriptor<double>::format(),
                                        2, // ndim
                                        std::vector<size_t>  nbRows1, nbCols2 , // shape
                                        std::vector<size_t> nbRows1 * sizeof(double), sizeof(double) // strides
                                    )
                                );
    delete[] ptr;
    return mRes;

我比较了执行两个固定 500*500 随机生成矩阵的乘积的平均（500 个产品）时间，得到以下结果：

python with numpy :    0.0067s
python with pybind11 : 0.7941s

那个 118 因素让我感到惊讶。当然，我没想到在第一次尝试时就击败了 numpy，但是两次平均时间之间的 100 倍让我感到惊讶。如果我将 intel mkl 用于产品的 c++ 部分或任何其他库，我认为该因素不会得到显着改善。

所以我猜这个因素主要是由numpy数组“转换”成py::array_t<double>s和逆转换来解释的。

我知道 numpy 依赖于 c 代码（很快就会有 c++ 代码），但我真的很想知道这些转换是如何在 numpy 中完成的。我在 github 上浏览了 numpy 的源代码，但没有找到“编组”部分和 c 产品部分。

【参考方案1】：

如果我将 intel mkl 用于产品的 c++ 部分，我认为该因素不会得到显着改善"

绝对会的。 MKL 和类似库中的矩阵矩阵乘积 (GEMM) 是有史以来最高度优化的代码之一。它与您的 ad-hoc 循环完全不同。

【讨论】：

然后让我用mkl检查一下

像往常一样，我遇到了 mkl 和 Visual Studio 的链接问题。十年，一切都没有改变。

链接，C++ 的祸根 ;-) 你也可以试试 OpenBLAS，它的性能与 MKL 相似。

:) 实际上，这本身不是 mkl 问题，请在此处查看我的问题：***.com/questions/50885810/… 这是 Visual Studio python 项目的问题......（至少我认为。）

哦。抱歉，我帮不上忙，我很久以前就离开了 VS 世界:)