Разница в производительности между NumPy и Matlab
Я вычисляю backpropagation алгоритм разреженного автоэнкодера. Я реализовал это в Python, используя numpy И в matlab, Код почти такой же, но производительность сильно отличается. Время, которое требуется matlab для выполнения задачи, составляет 0,252454 секунды, а нулевое - 0,973672151566, что почти в четыре раза больше. Я назову этот код несколько раз позже в проблеме минимизации, так что эта разница приводит к задержке в несколько минут между реализациями. Это нормальное поведение? Как я мог улучшить производительность в NumPy?
Numpy реализация:
Sparse.rho - параметр настройки, sparse.nodes - количество узлов в скрытом слое (25), sparse.input (64) - количество узлов во входном слое, theta1 и theta2 - матрицы весов для первого и второй слой соответственно с размерами 25x64 и 64x25, m равен 10000, rhoest имеет размерность (25,), x имеет размерность 10000x64, a3 10000x64 и a2 10000x25.
UPDATEЯ внес изменения в код, следуя некоторым идеям ответов. Производительность сейчас ничтожная: 0,65 против Matlab: 0,25.
partial_j1 = np.zeros(sparse.theta1.shape)
partial_j2 = np.zeros(sparse.theta2.shape)
partial_b1 = np.zeros(sparse.b1.shape)
partial_b2 = np.zeros(sparse.b2.shape)
t = time.time()
delta3t = (-(x-a3)*a3*(1-a3)).T
for i in range(m):
delta3 = delta3t[:,i:(i+1)]
sum1 = np.dot(sparse.theta2.T,delta3)
delta2 = ( sum1 + sum2 ) * a2[i:(i+1),:].T* (1 - a2[i:(i+1),:].T)
partial_j1 += np.dot(delta2, a1[i:(i+1),:])
partial_j2 += np.dot(delta3, a2[i:(i+1),:])
partial_b1 += delta2
partial_b2 += delta3
print "Backprop time:", time.time() -t
Реализация Matlab:
tic
for i = 1:m
delta3 = -(data(i,:)-a3(i,:)).*a3(i,:).*(1 - a3(i,:));
delta3 = delta3.';
sum1 = W2.'*delta3;
sum2 = beta*(-sparsityParam./rhoest + (1 - sparsityParam) ./ (1.0 - rhoest) );
delta2 = ( sum1 + sum2 ) .* a2(i,:).' .* (1 - a2(i,:).');
W1grad = W1grad + delta2* a1(i,:);
W2grad = W2grad + delta3* a2(i,:);
b1grad = b1grad + delta2;
b2grad = b2grad + delta3;
end
toc
3 ответа
Было бы неправильно говорить "Matlab всегда быстрее, чем NumPy" или наоборот. Часто их производительность сопоставима. При использовании NumPy, чтобы получить хорошую производительность, вы должны иметь в виду, что скорость NumPy зависит от вызова базовых функций, написанных на C/C++/Fortran. Он хорошо работает, когда вы применяете эти функции к целым массивам. В целом, вы получаете худшую производительность, когда вызываете функцию NumPy для меньших массивов или скаляров в цикле Python.
Вы спрашиваете, что не так с циклом Python? Каждая итерация цикла Python - это вызов next метод. Каждое использование [] индексация это вызов __getitem__ метод. каждый += это вызов __iadd__, Каждый точечный поиск атрибута (например, как в np.dot) включает вызовы функций. Эти вызовы функций в целом создают значительное препятствие для скорости. Эти хуки дают Python выразительную силу - индексация для строк означает нечто иное, чем, например, индексация для dicts. Тот же синтаксис, разные значения. Волшебство достигается, давая объектам разные __getitem__ методы.
Но эта выразительная сила достигается за счет скорости. Поэтому, когда вам не нужна вся эта динамическая выразительность, чтобы повысить производительность, попробуйте ограничиться вызовами функций NumPy для целых массивов.
Итак, удалите цикл for; по возможности используйте "векторизованные" уравнения. Например, вместо
for i in range(m):
delta3 = -(x[i,:]-a3[i,:])*a3[i,:]* (1 - a3[i,:])
вы можете вычислить delta3 для каждого i все сразу:
delta3 = -(x-a3)*a3*(1-a3)
Тогда как в for-loopdelta3 является вектором, используя векторизованное уравнение delta3 это матрица
Некоторые из вычислений в for-loop не зависит от i и поэтому должен быть поднят за пределы петли. Например, sum2 выглядит как константа:
sum2 = sparse.beta*(-float(sparse.rho)/rhoest + float(1.0 - sparse.rho) / (1.0 - rhoest) )
Вот работающий пример с альтернативной реализацией (alt) вашего кода (orig).
Мой тест timeit показывает 6,8-кратное улучшение скорости:
In [52]: %timeit orig()
1 loops, best of 3: 495 ms per loop
In [53]: %timeit alt()
10 loops, best of 3: 72.6 ms per loop
import numpy as np
class Bunch(object):
""" http://code.activestate.com/recipes/52308 """
def __init__(self, **kwds):
self.__dict__.update(kwds)
m, n, p = 10 ** 4, 64, 25
sparse = Bunch(
theta1=np.random.random((p, n)),
theta2=np.random.random((n, p)),
b1=np.random.random((p, 1)),
b2=np.random.random((n, 1)),
)
x = np.random.random((m, n))
a3 = np.random.random((m, n))
a2 = np.random.random((m, p))
a1 = np.random.random((m, n))
sum2 = np.random.random((p, ))
sum2 = sum2[:, np.newaxis]
def orig():
partial_j1 = np.zeros(sparse.theta1.shape)
partial_j2 = np.zeros(sparse.theta2.shape)
partial_b1 = np.zeros(sparse.b1.shape)
partial_b2 = np.zeros(sparse.b2.shape)
delta3t = (-(x - a3) * a3 * (1 - a3)).T
for i in range(m):
delta3 = delta3t[:, i:(i + 1)]
sum1 = np.dot(sparse.theta2.T, delta3)
delta2 = (sum1 + sum2) * a2[i:(i + 1), :].T * (1 - a2[i:(i + 1), :].T)
partial_j1 += np.dot(delta2, a1[i:(i + 1), :])
partial_j2 += np.dot(delta3, a2[i:(i + 1), :])
partial_b1 += delta2
partial_b2 += delta3
# delta3: (64, 1)
# sum1: (25, 1)
# delta2: (25, 1)
# a1[i:(i+1),:]: (1, 64)
# partial_j1: (25, 64)
# partial_j2: (64, 25)
# partial_b1: (25, 1)
# partial_b2: (64, 1)
# a2[i:(i+1),:]: (1, 25)
return partial_j1, partial_j2, partial_b1, partial_b2
def alt():
delta3 = (-(x - a3) * a3 * (1 - a3)).T
sum1 = np.dot(sparse.theta2.T, delta3)
delta2 = (sum1 + sum2) * a2.T * (1 - a2.T)
# delta3: (64, 10000)
# sum1: (25, 10000)
# delta2: (25, 10000)
# a1: (10000, 64)
# a2: (10000, 25)
partial_j1 = np.dot(delta2, a1)
partial_j2 = np.dot(delta3, a2)
partial_b1 = delta2.sum(axis=1)
partial_b2 = delta3.sum(axis=1)
return partial_j1, partial_j2, partial_b1, partial_b2
answer = orig()
result = alt()
for a, r in zip(answer, result):
try:
assert np.allclose(np.squeeze(a), r)
except AssertionError:
print(a.shape)
print(r.shape)
raise
Совет: обратите внимание, что я оставил в комментариях форму всех промежуточных массивов. Знание формы массивов помогло мне понять, что делает ваш код. Форма массивов может помочь вам выбрать правильные функции NumPy. Или, по крайней мере, уделение внимания формам поможет вам понять, является ли операция разумной. Например, когда вы вычисляете
np.dot(A, B)
а также A.shape = (n, m) а также B.shape = (m, p), затем np.dot(A, B) будет массив формы (n, p),
Это может помочь построить массивы в C_CONTIGUOUS-порядке (по крайней мере, если использовать np.dot). При этом скорость может увеличиться в 3 раза:
Ниже, x такой же как xf Кроме этого x является C_CONTIGUOUS и xf это F_CONTIGUOUS - и те же отношения для y а также yf,
import numpy as np
m, n, p = 10 ** 4, 64, 25
x = np.random.random((n, m))
xf = np.asarray(x, order='F')
y = np.random.random((m, n))
yf = np.asarray(y, order='F')
assert np.allclose(x, xf)
assert np.allclose(y, yf)
assert np.allclose(np.dot(x, y), np.dot(xf, y))
assert np.allclose(np.dot(x, y), np.dot(xf, yf))
%timeit тесты показывают разницу в скорости:
In [50]: %timeit np.dot(x, y)
100 loops, best of 3: 12.9 ms per loop
In [51]: %timeit np.dot(xf, y)
10 loops, best of 3: 27.7 ms per loop
In [56]: %timeit np.dot(x, yf)
10 loops, best of 3: 21.8 ms per loop
In [53]: %timeit np.dot(xf, yf)
10 loops, best of 3: 33.3 ms per loop
Что касается бенчмаркинга в Python:
Это может вводить в заблуждение использовать разницу в парах time.time() вызовы для измерения скорости кода в Python. Вам нужно повторить измерение много раз. Лучше отключить автоматический сборщик мусора. Также важно измерять большие промежутки времени (например, повторения не менее 10 секунд), чтобы избежать ошибок из-за плохого разрешения в таймере часов и уменьшить значимость time.time колл накладные расходы. Вместо того, чтобы писать весь этот код самостоятельно, Python предоставляет вам модуль timeit. По сути, я использую это для определения времени фрагментов кода, за исключением того, что для удобства я вызываю его через терминал IPython.
Я не уверен, влияет ли это на ваши тесты, но помните, что это может иметь значение. В вопросе, на который я ссылаюсь, согласно time.time два фрагмента кода отличались в 1,7 раза при использовании тестов timeit показал, что фрагменты кода выполнялись в практически одинаковом количестве времени.
Я бы начал с операций на месте, чтобы каждый раз не выделять новые массивы:
partial_j1 += np.dot(delta2, a1[i,:].reshape(1,a1.shape[1]))
partial_j2 += np.dot(delta3, a2[i,:].reshape(1,a2.shape[1]))
partial_b1 += delta2
partial_b2 += delta3
Вы можете заменить это выражение:
a1[i,:].reshape(1,a1.shape[1])
с более простым и быстрым (спасибо Би Рико):
a1[i:i+1]
Также эта строка:
sum2 = sparse.beta*(-float(sparse.rho)/rhoest + float(1.0 - sparse.rho) / (1.0 - rhoest))
Кажется, что он одинаков в каждом цикле, вам не нужно пересчитывать его.
И, возможно, небольшая оптимизация, вы можете заменить все вхождения x[i,:] с x[i],
Наконец, если вы можете позволить себе выделить m В разы больше памяти, вы можете следовать предложению unutbu и векторизовать цикл:
for m in range(m):
delta3 = -(x[i]-a3[i])*a3[i]* (1 - a3[i])
с:
delta3 = -(x-a3)*a3*(1-a3)
И вы всегда можете использовать Numba и значительно увеличить скорость без векторизации (и без использования большего количества памяти).
Разница в производительности между NumPy и Matlab всегда меня расстраивала. Они часто в конечном итоге сводятся к лежащим в основе лэпак библиотекам. Насколько я знаю, Matlab по умолчанию использует полный атлас Lapack, а NumPy использует свет Lapack. Matlab считает, что люди не заботятся о пространстве и объеме, в то время как NumPy считает людей. Подобный вопрос с хорошим ответом.