Точка Numpy в цикле for

Question

У меня есть следующий код, где я использую numpy.dot, чтобы ускорить мои вычисления.

u = numpy.zeros((l, l)) 
wp = numpy.zeros((l,2)) 
# some code which edits u and wp 
for x in range(N): 
    wavg = numpy.dot(wp[:, 0], wp[:, 1]) 
    wp[:, 0] = 1.0/wavg*numpy.dot(u, numpy.multiply(wp[:, 0], wp[:, 1])) 

Для малых l самая медленная часть - это внешний контур. Теперь я спрашиваю себя, есть ли способ избавиться от этого цикла?

Edit: В математических терминах этот код будет выглядеть следующим образом

Answer 1

Небольшое улучшение путем предварительного расчета на срок

u = np.zeros((l, l)) 
wp = np.zeros((l,2)) 
# some code which edits u and wp 

m = u*wp[:, 1] 
for x in range(N): 
    wavg = np.dot(wp[:, 0], wp[:, 1]) 
    wp[:, 0] = 1.0/wavg*np.dot(m, wp[:, 0]) 

Но мы можем сделать лучше - вместо вычисления среднего каждый раз, мы можем сделать это только на последней итерации:

m = u*wp[:, 1] 
for x in range(N - 1): 
    wp[:, 0] = np.dot(m, wp[:, 0]) 

wavg = np.dot(wp[:, 0], wp[:, 1]) 
wp[:, 0] = 1.0/wavg*np.dot(m, wp[:, 0]) 

Но есть еще одна вещь, которую мы можем сделать - это цикл может быть заменен матрицей экспоненциации:

m = u*wp[:, 1] 
wp[:, 0] = np.linalg.matrix_power(m, N-1).dot(wp[:, 0]) 
wavg = np.dot(wp[:, 0], wp[:, 1]) 
wp[:, 0] = 1.0/wavg*np.dot(m, wp[:, 0]) 

К сожалению, это выглядит медленнее. Однако, если вы можете прекомпотировать np.linalg.matrix_power(m, N-1), то это будет быстрее.

Answer 2

Я до сих пор путают о форме ваших массивов, а также соотношение между кодом и уравнением.

Но только глядя на уравнение, я думаю, что это может быть вычислена как:

In [515]: n,m = 3,4 

In [516]: U = np.ones((n,m)) 

In [517]: w = np.ones((m,)) 

In [518]: f = np.ones((m,)) 

In [519]: np.einsum('ij,j,j->i',U,w,f) 
Out[519]: array([ 4., 4., 4.]) 

На данный момент я заинтересован в получении размеров, чтобы соответствовать; а не на конечных значениях. Расчет достаточно прост, что ему не нужна нотация Эйнштейна, но einsum делает перевод почти механическим.

dot эквивалент

In [520]: np.dot(U, w*f) 
Out[520]: array([ 4., 4., 4.]) 

Так как это развивается с течением времени, f итерации в зависимости от значения предыдущего (и это внешнее значения w(t)), это трудно удалить эту петлю; мы можем просто сделать содержимое быстрее.