Быстрые циклы Numpy

Ответ 1

Здесь код для outer:

def outer(a, b, out=None):    
    a = asarray(a)
    b = asarray(b)
    return multiply(a.ravel()[:, newaxis], b.ravel()[newaxis,:], out)

Таким образом, каждый вызов outer включает в себя несколько вызовов python. В конечном итоге они вызывают скомпилированный код для выполнения умножения. Но каждый из них несет накладные расходы, которые не имеют никакого отношения к размеру ваших массивов.

Таким образом, 200 (200 ** 2?) звонков на outer будут иметь все эти служебные данные, тогда как один вызов outer со всеми 200 строками имеет один служебный набор, за которым следует одна операция с быстрой компиляцией.

cython и numba не компилируются и не обходят код Python в outer. Все, что они могут сделать, это упорядочить итерационный код, который вы написали, - и это не требует много времени.

Не вдаваясь в подробности, MATLAB jit должен иметь возможность заменить "внешний" на более быстрый код - он перезаписывает итерацию. Но мой опыт работы с MATLAB датируется временем до его jit.

Для реальных улучшений скорости с помощью cython и numba вам нужно использовать примитивный код numpy/python до конца. Или еще лучше сосредоточьте свои усилия на медленных внутренних кусках.

Замена outer упрощенной версией сокращает время выполнения примерно в два раза:

def foo1(N):
        size = N
        np.random.seed(1000031212)
        bar = np.random.rand(size, size)
        moo = np.zeros((size,size), dtype = np.float)
        for i in range(0,size):
                for j in range(0,size):
                        val = bar[j]
                        moo += val[:,None]*val   
        return moo

С полным N=200 ваша функция заняла 17 секунд за цикл. Если я заменю внутренние две строки на pass (без вычисления), время падает до 3 мс за цикл. Другими словами, внешний механизм цикла не является большим потребителем времени, по крайней мере, не по сравнению со многими вызовами outer().

Ответ 2

При наличии памяти вы можете использовать np.einsum для выполнения этих тяжелых вычислений векторизованным способом, например:

moo = size*np.einsum('ij,ik->jk',bar,bar)

Можно также использовать np.tensordot -

moo = size*np.tensordot(bar,bar,axes=(0,0))

Или просто np.dot -

moo = size*bar.T.dot(bar)

Ответ 3

Многие учебники и демонстрации Cython, Numba и т.д. показывают, что эти инструменты могут ускорить ваш код автоматически, но на практике это часто бывает не так: вам нужно немного изменить код извлеките лучшую производительность. Если вы уже выполнили некоторую степень векторизации, это обычно означает выписывание ВСЕХ циклов. Причины Операции массива Numpy не являются оптимальными:

• Много временных массивов создаются и зацикливаются;
• Значительные накладные расходы для каждого вызова, если массивы малы;
• Логика короткого замыкания не может быть реализована, поскольку массивы обрабатываются в целом;
• Иногда оптимальный алгоритм не может быть выражен с использованием выражений массива, и вы соглашаетесь на алгоритм с худшей временной сложностью.

Использование Numba или Cython не будет оптимизировать эти проблемы! Вместо этого эти инструменты позволяют писать код цикла, который намного быстрее, чем простой Python.

Кроме того, для Numba вы должны знать разницу между "object mode" и "nopython mode" . Тесные петли из вашего примера должны работать в режиме nopython, чтобы обеспечить значительное ускорение. Однако numpy.outer еще не поддерживается Numba, в результате чего функция компилируется в режиме объекта. Украсьте с помощью jit(nopython=True), чтобы такие случаи выдавали исключение.

Пример для демонстрации ускорения действительно возможен:

import numpy as np
from numba import jit

@jit
def foo_nb(bar):
    size = bar.shape[0]
    moo = np.zeros((size, size))
    for i in range(0,size):
        for j in range(0,size):
            val = bar[j]
            moo += np.outer(val, val)
    return moo

@jit
def foo_nb2(bar):
    size = bar.shape[0]
    moo = np.zeros((size, size))
    for i in range(size):
        for j in range(size):
            for k in range(0,size):
                for l in range(0,size):
                    moo[k,l] += bar[j,k] * bar[j,l]
    return moo

size = 100
bar = np.random.rand(size, size)

np.allclose(foo_nb(bar), foo_nb2(bar))
# True

%timeit foo_nb(bar)
# 1 loop, best of 3: 816 ms per loop
%timeit foo_nb2(bar)
# 10 loops, best of 3: 176 ms per loop

Ответ 4

Пример, который вы нам показываете, - это неэффективный алгоритм, так как вы вычисляете один и тот же внешний продукт несколько раз. Результирующая временная сложность O (n ^ 4). Его можно свести к n ^ 3.

for i in range(0,size):
    val = bar[i]
    moo += size * np.outer(val, val)