c - 使用 cblas_sgemm 执行复杂的矩阵操作以进行乘法

Question

我有 100 个 3x3x3 矩阵，我想将它们与另一个大小为 3x5x5 的大矩阵相乘(类似于用多个滤波器对一个图像进行卷积，但不完全一样)。

为了解释起见，这就是我的大矩阵的样子:

>>> x = np.arange(75).reshape(3, 5, 5)
>>> x
array([[[ 0,  1,  2,  3,  4],
        [ 5,  6,  7,  8,  9],
        [10, 11, 12, 13, 14],
        [15, 16, 17, 18, 19],
        [20, 21, 22, 23, 24]],

       [[25, 26, 27, 28, 29],
        [30, 31, 32, 33, 34],
        [35, 36, 37, 38, 39],
        [40, 41, 42, 43, 44],
        [45, 46, 47, 48, 49]],

       [[50, 51, 52, 53, 54],
        [55, 56, 57, 58, 59],
        [60, 61, 62, 63, 64],
        [65, 66, 67, 68, 69],
        [70, 71, 72, 73, 74]]])

在内存中，我假设大矩阵中的所有子矩阵都存储在连续的位置(如果我错了，请纠正我)。我想要做的是，从这个 3x5x5 矩阵中，我想从大矩阵的每个子矩阵中提取 3 列 5x3 列，然后水平连接它们以获得5x9 矩阵(如果这部分不清楚，我很抱歉，如果需要，我可以更详细地解释)。如果我使用 numpy，我会这样做:

>>> k = np.hstack(np.vstack(x)[:, 0:3].reshape(3, 5, 3))
>>> k
array([[ 0,  1,  2, 25, 26, 27, 50, 51, 52],
       [ 5,  6,  7, 30, 31, 32, 55, 56, 57],
       [10, 11, 12, 35, 36, 37, 60, 61, 62],
       [15, 16, 17, 40, 41, 42, 65, 66, 67],
       [20, 21, 22, 45, 46, 47, 70, 71, 72]])

但是，我没有使用 python，所以我无法访问我需要的 numpy 函数，以便将数据 block reshape 为我想要执行乘法的形式...我只能直接调用C 语言中的 cblas_sgemm 函数(来自 BLAS 库)，其中 k 对应于输入 B。

这是我对 cblas_sgemm 的调用:

    cblas_sgemm(    CblasRowMajor, CblasNoTrans, CblasTrans, 
                    100, 5, 9, 
                    1.0,
                    A, 9,
                    B, 9, // this is actually wrong, since I don't know how to specify the right parameter
                    0.0,
                    result, 5);

基本上，ldb 属性是这里的违规者，因为我的数据没有按照我需要的方式被阻止。我尝试了不同的方法，但我无法让 cblas_sgemm 理解我希望它如何读取和理解我的数据。

简而言之，我不知道如何告诉 cblas_sgemm 阅读 x 喜欢 k。有没有一种方法可以在将数据发送到 C 之前在 python 中巧妙地 reshape 数据，以便 cblas_sgemm 可以按照我想要的方式工作吗？

我将通过设置CblasTrans转置k，因此在乘法过程中，B是9x5。我的矩阵 A 的形状为 100x9。希望有帮助。

如有任何帮助，我们将不胜感激。谢谢!

Answer 1

In short, I don't know how to tell cblas_sgemm to read x like k.

你不能。您必须复印一份。

考虑k :

In [20]: k
Out[20]: 
array([[ 0,  1,  2, 25, 26, 27, 50, 51, 52],
       [ 5,  6,  7, 30, 31, 32, 55, 56, 57],
       [10, 11, 12, 35, 36, 37, 60, 61, 62],
       [15, 16, 17, 40, 41, 42, 65, 66, 67],
       [20, 21, 22, 45, 46, 47, 70, 71, 72]])

在二维数组中，内存中元素的间距在每个轴上必须相同。从怎么知道x创建后内存中的连续元素为0, 1, 2, 3, 4, ... ，但是你的第一行 k包含0, 1, 2, 25, 26, .... 。 1 之间没有空格和2 (即内存地址增加了数组一个元素的大小)，但是2之间内存中有一个很大的跳跃和25 。因此，您必须复制一份才能创建 k .

<小时/>话虽如此，还有一种替代方法可以使用一些 reshape (无需复制)和 numpy 的 einsum 来有效地实现您想要的最终结果。功能。

这是一个例子。首先定义x和A :

In [52]: x = np.arange(75).reshape(3, 5, 5)

In [53]: A = np.arange(90).reshape(10, 9)

这是我对你想要实现的目标的理解； A.dot(k.T)是期望的结果:

In [54]: k = np.hstack(np.vstack(x)[:, 0:3].reshape(3, 5, 3))

In [55]: A.dot(k.T)
Out[55]: 
array([[ 1392,  1572,  1752,  1932,  2112],
       [ 3498,  4083,  4668,  5253,  5838],
       [ 5604,  6594,  7584,  8574,  9564],
       [ 7710,  9105, 10500, 11895, 13290],
       [ 9816, 11616, 13416, 15216, 17016],
       [11922, 14127, 16332, 18537, 20742],
       [14028, 16638, 19248, 21858, 24468],
       [16134, 19149, 22164, 25179, 28194],
       [18240, 21660, 25080, 28500, 31920],
       [20346, 24171, 27996, 31821, 35646]])

以下是如何通过切片 x 获得相同结果的方法和 reshape A :

In [56]: x2 = x[:,:,:3]

In [57]: A2 = A.reshape(-1, 3, 3)

In [58]: einsum('ijk,jlk', A2, x2)
Out[58]: 
array([[ 1392,  1572,  1752,  1932,  2112],
       [ 3498,  4083,  4668,  5253,  5838],
       [ 5604,  6594,  7584,  8574,  9564],
       [ 7710,  9105, 10500, 11895, 13290],
       [ 9816, 11616, 13416, 15216, 17016],
       [11922, 14127, 16332, 18537, 20742],
       [14028, 16638, 19248, 21858, 24468],
       [16134, 19149, 22164, 25179, 28194],
       [18240, 21660, 25080, 28500, 31920],
       [20346, 24171, 27996, 31821, 35646]])