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

14

Chcę zaimplementować następujące wyrażenie w Pythonie: gdzie i są tablicami liczbowymi o rozmiarze , a to tablica liczbowa o wielkości . Rozmiar 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.

x_{ja} = \sum_{jot = 1}^{ja - 1} k_{ja - jot, jot} {za}_{ja - jot} {za}_{jot},

$x_i = \sum_{j=1}^{i-1}k_{i-j,j}a_{i-j}a_j,$

x

$x$

y

$y$

n

$n$

k

$k$

n \times n

$n\times n$

n

$n$

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?

python numpy

— Nataniel
źródło

Kilka dni temu miałem podobne pytanie. Zapytałem o to podczas nakładania stosu. Sprawdź ten post . Używam scipy.weave zamiast cython. Czy ktoś wie, czy to powoduje (znaczącą) różnicę w wydajności?

— 2013

17

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.

— Travis Oliphant
źródło

Dzięki, to wygląda całkiem prosto. Nawet nie wiedziałem o Numbie! Cython, PyPy, Numba ... to zagmatwany świat.

— Nathaniel,

3

Travis, bardzo fajnie, czy masz coś przeciwko dodaniu ujawnienia na dole odpowiedzi, że jesteś jednym z programistów Numba?

— Aron Ahmadia

1

n = 200

$n=200$

@NatWilson - jeśli zadasz to pytanie na scicomp, chętnie spróbuję go rozwiązać :)

— Aron Ahmadia

4

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.

— Nat Wilson
źródło

x_{0}

$x_0$

i - 1

$i-1$

O, rozumiem. Zebrałem to z pierwotnego podsumowania, ale nie byłem pewien, czy taki był zamiar.

— Nat Wilson,

1

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])

— Thomas Baruchel
źródło