Ответ 1
Было бы неправильно говорить: "Matlab всегда быстрее, чем NumPy", или наоборот наоборот. Часто их производительность сопоставима. При использовании NumPy, чтобы получить хорошие вы должны иметь в виду, что скорость NumPy происходит от вызова базовые функции, написанные на C/С++/Fortran. Он хорошо работает, когда вы применяете эти функции для целых массивов. В общем, вы получаете более низкую производительность, когда вы вызываете эту функцию NumPy на меньших массивах или скалярах в цикле Python.
Что плохого в цикле Python вы спрашиваете? Каждая итерация через цикл Python
вызов метода next. Каждое использование индексации [] - это вызов
__getitem__. Каждый += является вызовом __iadd__. Каждый пунктирный атрибут
поиск (например, как np.dot) включает вызовы функций. Эти вызовы функций
что значительно затрудняет скорость. Эти крючки дают Python
экспрессивная мощность - индексация для строк означает нечто иное, чем индексирование
например, для dicts. Тот же синтаксис, разные значения. Магия достигается путем предоставления объектам разных методов __getitem__.
Но эта выразительная сила стоит дорого. Поэтому, когда вам не нужны все что динамическая выразительность, чтобы повысить производительность, постарайтесь ограничить себя Функция NumPy вызывает целые массивы.
Итак, удалите for-loop; используйте "векторизованные" уравнения, когда это возможно. Например, вместо
for i in range(m):
delta3 = -(x[i,:]-a3[i,:])*a3[i,:]* (1 - a3[i,:])
вы можете вычислить delta3 для каждого i одновременно:
delta3 = -(x-a3)*a3*(1-a3)
В то время как в for-loop delta3 есть вектор, использование векторизованного уравнения delta3 является матрицей.
Некоторые из вычислений в for-loop не зависят от i и поэтому должны быть подняты вне цикла. Например, sum2 выглядит как константа:
sum2 = sparse.beta*(-float(sparse.rho)/rhoest + float(1.0 - sparse.rho) / (1.0 - rhoest) )
Вот пример с альтернативной реализацией (alt) вашего кода (orig).
Мой тест timeit показывает 6.8x улучшение скорости:
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 module. Я по сути использую это ко времени фрагментов кода, за исключением того, что я вызываю его через , который я связал с, согласно time.time, два фрагмента кода отличались в 1,7 раза, а тесты с использованием timeit показали, что фрагменты кода выполнялись в в основном идентичные количества времени.