Comment exprimer cette expression compliquée à l'aide de tranches numpy

Je souhaite implémenter l'expression suivante en Python:

$$x_i = \sum_{j=1}^{i-1}k_{i-j,j}a_{i-j}a_j,$$ où et sont des tableaux numpy de taille , et est un tableau numpy de taille . La taille peut aller jusqu'à environ 10000, et la fonction fait partie d'une boucle interne qui sera évaluée plusieurs fois, donc la vitesse est importante.$x$$y$$n$$k$$n\times n$$n$

Idéalement, j'aimerais éviter complètement une boucle for, même si je suppose que ce n'est pas la fin du monde s'il y en a une. Le problème est que j'ai du mal à voir comment le faire sans avoir quelques boucles imbriquées, et cela risque de le rendre plutôt lent.

Quelqu'un peut-il voir comment exprimer l'équation ci-dessus en utilisant numpy d'une manière efficace et de préférence également lisible? Plus généralement, quelle est la meilleure façon d'aborder ce genre de chose?

Voici la solution Numba. Sur ma machine, la version Numba est> 1000x plus rapide que la version python sans le décorateur (pour une matrice 200x200, 'k' et un vecteur de longueur 200 'a'). Vous pouvez également utiliser le décorateur @autojit qui ajoute environ 10 microsecondes par appel afin que le même code fonctionne avec plusieurs types.

from numba import jit, autojit

@jit('f8[:](f8[:,:],f8[:])')
#@autojit
def looped_ver(k, a):
    x = np.empty_like(a)
    for i in range(x.size):
        sm = 0.0
        for j in range(0, i+1):
            sm += k[i-j,j] * a[i-j] * a[j]
        x[i] = sm
    return x

Divulgation: je suis l'un des développeurs de Numba.

Voici un début. Tout d'abord, mes excuses pour toute erreur.

J'ai expérimenté quelques approches différentes. J'étais un peu confus par les limites de la somme - la limite supérieure devrait-elle être , plutôt que i - 1 ?$i$$i-1$

Edit: Non, la limite supérieure était correcte comme indiqué dans la question. Je l'ai laissé tel quel, car une autre réponse utilise maintenant le même code, mais la correction est simple.

D'abord une version en boucle:

def looped_ver(k, a):
    x = np.empty_like(a)
    for i in range(x.size):
        sm = 0
        for j in range(0, i+1):
            sm += k[i-j,j] * a[i-j] * a[j]
        x[i] = sm
    return x

Je l'ai fait une seule boucle avec des tranches numpy:

def vectorized_ver(k, a):
    ktr = zeros_like(k)
    ar = zeros_like(k)
    sz = len(a)
    for i in range(sz):
        ktr[i,:i+1] = k[::-1].diagonal(-sz+i+1)
        a_ = a[:i+1]
        ar[i,:i+1] = a_[::-1] * a_
    return np.sum(ktr * ar, 1)

La version numpy avec une boucle explicite est environ 25 fois plus rapide sur mon ordinateur lorsque .$n=5000$

Ensuite, j'ai écrit une version Cython du code en boucle (plus lisible).

import numpy as np
import cython
cimport numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
def cyth_ver(double [:, ::1] k not None,
              double [:] a not None):
    cdef double[:] x = np.empty_like(a)
    cdef double sm
    cdef int i, j

    for i in range(len(a)):
        sm = 0.0
        for j in range(i+1):
            sm = sm + k[i-j,j] * a[i-j] * a[j]
        x[i] = sm
    return x

Sur mon ordinateur portable, celui-ci est environ 200 fois plus rapide que la version en boucle (et 8 fois plus rapide que la version vectorisée à 1 boucle). Je suis sûr que les autres peuvent faire mieux.

J'ai joué avec une version Julia, et elle semblait (si je la chronométrais correctement) comparable au code Cython.

Ce que vous voulez semble être une convolution; Je pense que le moyen le plus rapide d'y parvenir serait la numpy.convolvefonction.

Vous devrez peut-être fixer les indices en fonction de vos besoins exacts, mais je pense que vous aimeriez essayer quelque chose comme:

import numpy as np
a = [1, 2, 3, 4, 5]
k = [2, 4, 6, 8, 10]

result = np.convolve(a, k*a[::-1])