我在 MATLAB 和 Python 中设置了两个关于矩阵乘法与广播的相同测试。对于 Python,我使用了 NumPy,对于 MATLAB,我使用了 mtimesx使用 BLAS 的库。
MATLAB
close all; clear;
N = 1000 + 100; % a few initial runs to be trimmed off at the end
a = 100;
b = 30;
c = 40;
d = 50;
A = rand(b, c, a);
B = rand(c, d, a);
C = zeros(b, d, a);
times = zeros(1, N);
for ii = 1:N
tic
C = mtimesx(A,B);
times(ii) = toc;
end
times = times(101:end) * 1e3;
plot(times);
grid on;
title(median(times));
python
import timeit
import numpy as np
import matplotlib.pyplot as plt
N = 1000 + 100 # a few initial runs to be trimmed off at the end
a = 100
b = 30
c = 40
d = 50
A = np.arange(a * b * c).reshape([a, b, c])
B = np.arange(a * c * d).reshape([a, c, d])
C = np.empty(a * b * d).reshape([a, b, d])
times = np.empty(N)
for i in range(N):
start = timeit.default_timer()
C = A @ B
times[i] = timeit.default_timer() - start
times = times[101:] * 1e3
plt.plot(times, linewidth=0.5)
plt.grid()
plt.title(np.median(times))
plt.show()
- 我的 Python 是从使用 OpenBLAS 的
pip
安装的默认 Python。 - 我在英特尔 NUC i3 上运行。
MATLAB 代码在 1 毫秒内运行,而 Python 在 5.8 毫秒内运行,我不明白为什么,因为它们似乎都在使用 BLAS。
编辑
来自 python :
In [7]: np.__config__.show()
mkl_info:
libraries = ['mkl_rt']
library_dirs = [...]
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = [...]
blas_mkl_info:
libraries = ['mkl_rt']
library_dirs = [...]
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = [...]
blas_opt_info:
libraries = ['mkl_rt']
library_dirs = [...]
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = [...]
lapack_mkl_info:
libraries = ['mkl_rt']
library_dirs = [...]
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = [...]
lapack_opt_info:
libraries = ['mkl_rt']
library_dirs = [...]
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = [...]
使用 pip 从 numpy
In [2]: np.__config__.show()
blas_mkl_info:
NOT AVAILABLE
blis_info:
NOT AVAILABLE
openblas_info:
library_dirs = [...]
libraries = ['openblas']
language = f77
define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
library_dirs = [...]
libraries = ['openblas']
language = f77
define_macros = [('HAVE_CBLAS', None)]
lapack_mkl_info:
NOT AVAILABLE
openblas_lapack_info:
library_dirs = [...]
libraries = ['openblas']
language = f77
define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
library_dirs = [...]
libraries = ['openblas']
language = f77
define_macros = [('HAVE_CBLAS', None)]
编辑 2
我尝试将 C = A @ B
替换为 np.matmul(A, B, out=C)
并得到了 2 倍的更差时间,例如大约 11 毫秒。这真是奇怪。
最佳答案
您的 MATLAB 代码使用 float 组,但 NumPy 代码使用整数数组。这在时间上有很大的不同。对于 MATLAB 和 NumPy 之间的“同类”比较,Python/NumPy 代码也必须使用 float 组。
然而,这并不是唯一重要的问题。在 issue 7569 中讨论了 NumPy 的缺陷(并再次在 issue 8957 中)在 NumPy github 站点中。 “堆叠”数组的矩阵乘法不使用快速 BLAS 例程来执行乘法。这意味着二维以上数组的乘法可能比预期慢得多。
二维数组的乘法确实使用快速例程,因此您可以通过在循环中乘以各个二维数组来解决此问题。令人惊讶的是,尽管有 Python 循环的开销,但在许多情况下,它比应用于完整堆栈数组的 @
、matmul
或 einsum
更快。
这是 NumPy 问题中显示的一个函数的变体,它在 Python 循环中执行矩阵乘法:
def xmul(A, B):
"""
Multiply stacked matrices A (with shape (s, m, n)) by stacked
matrices B (with shape (s, n, p)) to produce an array with
shape (s, m, p).
Mathematically equivalent to A @ B, but faster in many cases.
The arguments are not validated. The code assumes that A and B
are numpy arrays with the same data type and with shapes described
above.
"""
out = np.empty((a.shape[0], a.shape[1], b.shape[2]), dtype=a.dtype)
for j in range(a.shape[0]):
np.matmul(a[j], b[j], out=out[j])
return out
我的 NumPy 安装也使用 MKL(它是 Anaconda 发行版的一部分)。下面是 A @ B
和 xmul(A, B)
的时间比较,使用浮点值数组:
In [204]: A = np.random.rand(100, 30, 40)
In [205]: B = np.random.rand(100, 40, 50)
In [206]: %timeit A @ B
4.76 ms ± 6.37 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [207]: %timeit xmul(A, B)
582 µs ± 35.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
尽管 xmul
使用 Python 循环,但它花费的时间大约是 A @ B
的 1/8。
关于python - MATLAB 矩阵乘法性能比 NumPy 快 5 倍,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/52857659/