Numpy массив и для l oop: как его улучшить

Question

Я новичок с Python, и у меня есть вопрос о том, чтобы l oop ускорился.

Пусть "u" будет numpy массивом измерений (N, K) и пусть "kernel_vect" будет numpy массивом измерений (K,) обоих чисел float64. Я хотел бы ускорить следующий код (исключив, например, для l oop)

Kernel_appo = np.zeros((N**2,))
    for k in range(K):
        uk = u[:,k]
        Mat_appo = np.outer(uk,uk)
        Kernel_appo = Kernel_appo  + kernel_vect[k] * routines.vec(Mat_appo)

Есть идеи? Спасибо!

Answer 1

Здесь более быстрая реализация без al oop over k:

# Version 2
Kernel_appo = np.zeros((N**2,))
for n1 in range(N):
    for n2 in range(N):
        Kernel_appo[n1*N+n2] = (u[n1,:] * u[n2,:] * kernel_vect).sum()

Мы можем сделать это быстрее, используя симметрию продукта u:

# Version 3
Kernel_appo = np.zeros((N,N))
for n1 in range(N):
    for n2 in range(n1,N):
        Kernel_appo[n1,n2] = (u[n1,:] * u[n2,:] * kernel_vect).sum()
Kernel_appo = np.triu(Kernel_appo, 1) + np.tril(Kernel_appo.transpose(), 0) # make the matrix symmetric
Kernel_appo = np.ravel(Kernel_appo, order='C')

Вот версия, которая удаляет одну из l oop:

# Version 4
Kernel_appo = np.zeros((N,N))
for n1 in range(N):
    Kernel_appo[n1,n1:N] = ((u[n1,:] * kernel_vect) * u[n1:N,:]).sum(axis=1)
Kernel_appo = np.triu(Kernel_appo, 1) + np.tril(Kernel_appo.transpose(), 0)
Kernel_appo = np.ravel(Kernel_appo, order='C')

У нас все еще есть al oop над N. Однако кажется разумным сохранить его, так как N мало. Его удаление, безусловно, заставит numpy создавать огромные матрицы в памяти, что должно привести к падению производительности (и даже может привести к падению sh, если N и K очень большие).

Обратите внимание, что версия 4, вероятно, не будет Быть таким быстрым, если K намного больше (поскольку временные матрицы numpy не могут поместиться в кэш процессора).

ОБНОВЛЕНИЕ : я просто обнаруживаю, что можно использовать awesome np.einsum в этом случае:

# Version 5
Kernel_appo = np.ravel(np.einsum('ji,ki,i->jk', u, u, kernel_vect, optimize=True), order='C')

Будьте готовы, потому что эта более простая реализация также намного быстрее (потому что numpy может векторизовать код и запустить его параллельно).

Вот результаты производительности с N = 50 и K = 5000 на моей машине:

Initial code: 58.15 ms
Version 2:    19.94 ms
Version 3:    10.11 ms
Version 4:     5.08 ms
Version 5:     0.57 ms

Окончательная реализация теперь примерно в 100 раз быстрее, чем начальный!