pybind11 与 numpy 的矩阵乘积

Question

编辑 2（2018 年 6 月 18 日）

我使用了Matrixclass中提出的

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> 的要好得多。图表（纵坐标为秒，横坐标为方阵大小）：

Numpy 在 intel mkl 下相对来说有点小。在大小区域 [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);

intel core i5 4570 上的初始代码运行。运行这次在 intel core i7 6700 HQ 上的所有三个案例都给出了：

两点说明 :

1) 对于相同的 Python 3.6.5 32 位，Python 在我笔记本电脑的核心 i7 上使用 numpy 比在我工作时使用的旧核心 i5 台式机上慢。朴素的 c++ 稍微快一点。非常st运行ge.

2) 在核心 i7 上，c++ intel mkl 与 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 运行domly 生成矩阵的乘积的平均（在 500 个产品上）时间，并得到以下结果：

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

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

所以我猜想这个因素主要是numpy数组"conversions"变成py::array_t<double>和逆向转换造成的。

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

Answer 1

I don't think the factor would have been drastically improved would I have used intel mkl for the c++ part of the product"

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

pybind11 与 numpy 的矩阵乘积

pybind11 vs numpy for a matrix product

c++

python

numpy

python-3.x

pybind11