Jak wyrazić to skomplikowane wyrażenie za pomocą numpy

Chcę zaimplementować następujące wyrażenie w Pythonie: gdzie x i y są tablicami liczbowymi o rozmiarze n , a k to tablica liczbowa o wielkości n × n . Rozmiar n może wynosić do około 10000, a funkcja jest częścią wewnętrznej pętli, która będzie oceniana wiele razy, więc szybkość jest ważna.

Idealnie chciałbym całkowicie uniknąć pętli for, choć myślę, że to nie koniec świata, jeśli taki istnieje. Problem polega na tym, że mam problem z widzeniem, jak to zrobić bez posiadania kilku zagnieżdżonych pętli, co może spowolnić.

Czy ktoś może zobaczyć, jak wyrazić powyższe równanie za pomocą numpy w sposób, który jest wydajny, a najlepiej także czytelny? Mówiąc bardziej ogólnie, jaki jest najlepszy sposób podejścia do tego rodzaju rzeczy?

Oto rozwiązanie Numba. Na mojej maszynie wersja Numba jest> 1000x szybsza niż wersja python bez dekoratora (dla matrycy 200x200, „k” i wektora o długości 200 „a”). Możesz także użyć dekoratora @autojit, który dodaje około 10 mikrosekund na połączenie, dzięki czemu ten sam kod będzie działał z wieloma typami.

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

Ujawnienie: Jestem jednym z programistów Numba.

Oto początek. Po pierwsze przepraszam za błędy.

$i$$i-1$

Edycja: Nie, górny limit był prawidłowy, jak podano w pytaniu. Zostawiłem go tak, jak jest tutaj, ponieważ inna odpowiedź używa teraz tego samego kodu, ale poprawka jest prosta.

Najpierw zapętlona wersja:

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

Zrobiłem z tego jedną pętlę z numpy plasterkami:

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)

$n=5000$

Potem napisałem cytowaną wersję (bardziej czytelnego) zapętlonego kodu.

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

Na moim laptopie ten jest około 200 razy szybszy niż wersja zapętlona (i 8 razy szybszy niż wersja wektoryzowana w 1 pętli). Jestem pewien, że inni mogą zrobić lepiej.

Grałem z wersją Julii i wydawało mi się (jeśli odpowiednio ją zaplanowałem) porównywalną z kodem Cython.

To, czego chcesz, wydaje się być splotem; Myślę, że najszybszym sposobem na osiągnięcie tego celu byłaby numpy.convolvefunkcja.

Być może będziesz musiał naprawić indeksy zgodnie z własnymi potrzebami, ale myślę, że chcesz spróbować czegoś takiego:

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

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