Производительность умножения матриц MATLAB в 5 раз выше, чем у NumPy

Question

Я установил два одинаковых теста в 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 по умолчанию установлен из pip, использующего OpenBLAS.
• Я работаю на Intel NUC i3.

Код MATLAB выполняется за 1 мс, а Python за 5,8 мс, и я не могу понять, почему, поскольку, похоже, они оба используют BLAS.

---

EDIT

От Анаконды:

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 = [...]

От 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) и получил 2x хуже время, например, около 11 мс.Это действительно странно.

Answer 1

Ваш код MATLAB использует массивы с плавающей запятой, но код NumPy использует целочисленные массивы.Это существенно меняет время.Для сравнения «яблоки с яблоками» между MATLAB и NumPy в коде Python / NumPy также должны использоваться массивы с плавающей запятой.

Однако это не единственная существенная проблема.Существует недостаток NumPy, обсуждаемый в выпуске 7569 (и снова в выпуске 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, это занимает примерно 1/8 времени A @ B.

Answer 2

Не могли бы вы попробовать еще раз с недавно выпущенным NumPy 1.16?Мы переработали matmul, чтобы использовать BLAS для внутренних двух измерений, что должно ускорить код.

Answer 3

Я думаю, что это проблема упорядочения памяти.Matlab zeros(a, b, c) похож на numpy zeros((a, b, c), order='F'), который не является значением по умолчанию.

Конечно, как вы уже определили, @ работает на разных осях, вплоть до mtimesx.Чтобы сделать сравнение справедливым, вы должны убедиться, что ваши массивы находятся в порядке matlab, а затем транспонировать для обработки разницы в семантике

# note: `order` in reshape actually changes the resulting array data,
# not just its memory layout
A = np.arange(a * b * c).reshape([b, c, a], order='F').transpose((2, 0, 1))
B = np.arange(a * c * d).reshape([c, d, a], order='F').transpose((2, 0, 1))