Разница в производительности между numpy и matlab

Ответ 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 показали, что фрагменты кода выполнялись в в основном идентичные количества времени.

Ответ 2

Я бы начал с операций inplace, чтобы не выделять новые массивы каждый раз:

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])

с более простой и быстрой (благодаря Bi Rico):

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 и значительно увеличить скорость без векторизации (и без использования большего объема памяти).

Ответ 3

Разница в производительности между numpy и matlab всегда разочаровывала меня. В конце концов они часто сводятся к основным библиотекам лап. Насколько я знаю, Matlab использует полный атлас лапак по умолчанию, а numpy использует яркий свет. Matlab считает, что люди не заботятся о пространстве и массе, в то время как numpy считает, что люди это делают. Подобный вопрос с хорошим ответом.