python - 使用 scipy.linalg.blas.sgemm 时，numpy 点返回大数组的无效值

Question

我正在尝试计算 A • A^T:

# These are my dummy values for testing
A = np.ones((150000,265),dtype=np.float32, order='F')
A_T = np.ones((265, 150000),dtype=np.float32, order='F')

out = scipy.linalg.blas.sgemm(alpha=1.0, a=A, b=A_T)

两分钟后:

In [7]: out
Out[7]: 
array([[ 265.,  265.,  265., ...,    0.,    0.,    0.],
       [ 265.,  265.,  265., ...,    0.,    0.,    0.],
       [ 265.,  265.,  265., ...,    0.,    0.,    0.],
       ..., 
       [ 265.,  265.,  265., ...,    0.,    0.,    0.],
       [ 265.,  265.,  265., ...,    0.,    0.,    0.],
       [ 265.,  265.,  265., ...,    0.,    0.,    0.]])

In [10]: out.shape
Out[10]: (150000, 150000)

注意到零了吗？我迷路了......我试过使用 64 位 float 并获得相同的输出。从 35468 开始，数组为零。

In [39]: out[0,35468]
Out[39]: 0.0

In [9]: scipy.__version__
Out[9]: '0.12.1'

更新/编辑:

我相当确定，np.dot 正在调用 *gemm 方法本身。

In [1]: A = np.ones((150000,265), dtype=np.float32, order='F')

In [2]: A_T = np.ones((265, 150000),dtype=np.float32, order='F')

In [3]: out = A.dot(A_T)

In [4]: out.shape
Out[4]: (150000, 150000)

In [5]: out
Out[5]: 
array([[ 265.,  265.,  265., ...,  265.,  265.,  265.],
   [ 265.,  265.,  265., ...,  265.,  265.,  265.],
   [ 265.,  265.,  265., ...,  265.,  265.,  265.],
   ..., 
   [   0.,    0.,    0., ...,    0.,    0.,    0.],
   [   0.,    0.,    0., ...,    0.,    0.,    0.],
   [   0.,    0.,    0., ...,    0.,    0.,    0.]], dtype=float32)

Answer 1

在我的机器上，您的示例因 Memory Error 而死，对于 np.dot(A, A.T) 也是如此。

在减小矩阵大小时，它会产生正确的结果。

要进一步调试它，请尝试直接使用 Fortran BLAS 调用，而不使用 python。如果通过，考虑filing a bug in the scipy tracker .

作为旁注，从 scipy 0.13 开始，您可以使用 syrk。