NumPy векторизации двойного Python для петли

Question

V (п, р) NumPy массив, как правило, размеры п ~ 10, р ~ 20000

Кода я теперь выгляжу

A = np.zeros(p) 
for i in xrange(n): 
    for j in xrange(i+1): 
     A += F[i,j] * V[i,:] * V[j,:] 

Как бы я переписывать это, чтобы избежать двойного питона для цикла?

Answer 1

Хотя Исаакиевская ответ представляется перспективным, так как она удаляет те две вложенные для петель, вы вынуждены создать промежуточный массив M, который n раз превышает размер исходного V массива. Python для петель не дешев, но доступ к памяти не является бесплатным либо:

n = 10 
p = 20000 
V = np.random.rand(n, p) 
F = np.random.rand(n, n) 

def op_code(V, F): 
    n, p = V.shape 
    A = np.zeros(p) 
    for i in xrange(n): 
     for j in xrange(i+1): 
      A += F[i,j] * V[i,:] * V[j,:] 
    return A 

def isaac_code(V, F): 
    n, p = V.shape 
    F = F.copy() 
    F[np.triu_indices(n, 1)] = 0 
    M = (V.reshape(n, 1, p) * V.reshape(1, n, p)) * F.reshape(n, n, 1) 
    return M.sum((0, 1)) 

Если теперь принять как для испытательной поездки:

In [20]: np.allclose(isaac_code(V, F), op_code(V, F)) 
Out[20]: True 

In [21]: %timeit op_code(V, F) 
100 loops, best of 3: 3.18 ms per loop 

In [22]: %timeit isaac_code(V, F) 
10 loops, best of 3: 24.3 ms per loop 

Так Извлечение для петель обходятся вам ое замедление. Не очень хорошая вещь ... На данный момент вы можете даже подумать, нужна ли функция, занимающая около 3 мс для оценки, дальнейшая оптимизация. В случае, если вы делаете, есть небольшое улучшение, которое можно было с помощью np.einsum:

def einsum_code(V, F): 
    n, p = V.shape 
    F = F.copy() 
    F[np.triu_indices(n, 1)] = 0 
    return np.einsum('ij,ik,jk->k', F, V, V) 

А теперь:

In [23]: np.allclose(einsum_code(V, F), op_code(V, F)) 
Out[23]: True 

In [24]: %timeit einsum_code(V, F) 
100 loops, best of 3: 2.53 ms per loop 

Так что примерно на 20% ускорить, который вводит код, который может очень хорошо не быть таким же читаемым, как ваши петли. Я бы сказал, не стоит ...

Answer 2

Сложная часть об этом заключается в том, что вы хотите только взять сумму элементов с помощью j <= i. Если бы не это, то вы можете сделать следующее:

M = (V.reshape(n, 1, p) * V.reshape(1, n, p)) * F.reshape(n, n, 1) 
A = M.sum(0).sum(0) 

Если F симметрично (если F[i,j] == F[j,i]), то вы можете использовать симметрию M над следующим:

D = M[range(n), range(n)].sum(0) 
A = (M.sum(0).sum(0) - D)/2.0 + D 

Это сказал, что это действительно не отличный кандидат для векторизации, так как у вас есть n << p, и поэтому ваши for -loops не окажут большого влияния на скорость этого вычисления.

Редактировать: Как сказано ниже Билл, вы можете просто убедиться, что элементы F, что вы не хотите использовать устанавливаются в ноле первый, а затем M.sum(0).sum(0) результата будет то, что вы хотите.

Answer 3

Выражение может быть записано в виде

и, таким образом, вы можете подвести его, как это с помощью np.newaxis -construct:

na = np.newaxis 
X = (np.tri(n)*F)[:,:,na]*V[:,na,:]*V[na,:,:] 
X.sum(axis=1).sum(axis=0) 

Здесь строится 3D-массив X[i,j,p], а затем суммируются две первые оси, что приводит к 1D-массиву A[p]. Дополнительно F умножилось на треугольную матрицу, чтобы ограничить суммирование по задаче.