Come posso vettorizzare e accelerare questo calcolo di grande array?

Question

Attualmente sto provando a calcolare la somma di tutta la somma di sottoquadri in una matrice di valori 10.000 x 10.000. Per fare un esempio, se il mio array è stato:

1 1 1 
2 2 2 
3 3 3 

Voglio che il risultato sia:

1+1+1+2+2+2+3+3+3      [sum of squares of size 1] 
+(1+1+2+2)+(1+1+2+2)+(2+2+3+3)+(2+2+3+3) [sum of squares of size 2] 
+(1+1+1+2+2+2+3+3+3)      [sum of squares of size 3] 
________________________________________ 
68 

Così, come un primo tentativo ho scritto un codice python molto semplice per farlo. Come era in O (k^2.n^2) (n essendo la dimensione del grande array e k la dimensione dei sottosquadri che stiamo ottenendo), l'elaborazione è stata terribilmente lunga. Ho scritto un altro algoritmo in O (n^2) per accelerarlo:

def getSum(tab,size): 
    n = len(tab) 
    tmp = numpy.zeros((n,n)) 

    for i in xrange(0,n): 
     sum = 0 
     for j in xrange(0,size): 
      sum += tab[j][i] 
     tmp[0][i] = sum 

     for j in xrange(1,n-size+1): 
      sum += (tab[j+size-1][i] - tab[j-1][i]) 
      tmp[j][i] = sum 

    finalsum = 0 
    for i in xrange(0,n-size+1): 
     sum = 0 
     for j in xrange(0,size): 
      sum += tmp[i][j] 
     finalsum += sum 

     for j in xrange(1,n-size+1): 
      finalsum += (tmp[i][j+size-1] - tmp[i][j-1]) 

return finalsum 

Quindi questo codice funziona correttamente. Data una matrice e una dimensione di sottoquadri, restituirà la somma dei valori in tutti questi sottoquadri. In pratica, eseguo un iterazione sulla dimensione dei sottosquadri per ottenere tutti i valori possibili.

Il problema è che questo è di nuovo troppo lungo per i grandi array (oltre 20 giorni per un array 10.000 x 10.000). L'ho cercato su google e ho imparato che avrei potuto vettorializzare le iterazioni su array con Numpy. Tuttavia, non riuscivo a capire come farlo nel mio caso ...

Se qualcuno può aiutarmi a velocizzare il mio algoritmo, o darmi una buona documentazione sull'argomento, sarò felice!

Grazie!

Answer 1

Queste sommatorie scorrevoli sono le più adatte per essere calcolate come somme di convoluzione 2D e quelle che potrebbero essere calcolate in modo efficiente con scipy's convolve2d. Così, per una dimensione specifica, si potrebbe ottenere le sommatorie, in questo modo -

def getSum(tab,size): 
    # Define kernel and perform convolution to get such sliding windowed summations 
    kernel = np.ones((size,size),dtype=tab.dtype) 
    return convolve2d(tab, kernel, mode='valid').sum() 

Per ottenere sommatorie in tutti i formati, penso che il modo migliore sia in termini di memoria e l'efficienza delle prestazioni sarebbe quella di utilizzare un ciclo per loop su tutte le possibili dimensioni. Quindi, per ottenere la somma finale, si avrebbe -

def getAllSums(tab): 
    finalSum = 0 
    for i in range(tab.shape[0]): 
     finalSum += getSum(tab,i+1) 
    return finalSum 

Campione run -

In [51]: tab 
Out[51]: 
array([[1, 1, 1], 
     [2, 2, 2], 
     [3, 3, 3]]) 

In [52]: getSum(tab,1) # sum of squares of size 1 
Out[52]: 18 

In [53]: getSum(tab,2) # sum of squares of size 2 
Out[53]: 32 

In [54]: getSum(tab,3) # sum of squares of size 3 
Out[54]: 18 

In [55]: getAllSums(tab) # sum of squares of all sizes 
Out[55]: 68 

Answer 2

base l'idea di calcolare quante volte ogni numero contato, sono arrivato a questo semplice codice:

def get_sum(matrix, n): 
    ret = 0 
    for i in range(n): 
     for j in range(n): 
      for k in range(1, n + 1): 
       # k is the square size. count is times of the number counted. 
       count = min(k, n - k + 1, i + 1, n - i) * min(k, n - k + 1, j + 1, n - j) 
       ret += count * matrix[i][j] 
    return ret 

a = [[1, 1, 1], [2, 2, 2], [3, 3, 3]] 

print get_sum(a, 3) # 68 

---

soluzione Divakar è fantastico, però, penso che il mio potrebbe essere più efficiente, almeno in asymp complessità temporale totale (O (n^3) rispetto a Divakar's O (n^3logn)).

---

Ho una soluzione O (n^2) ora ...

In sostanza, siamo in grado di ottenere che:

def get_sum2(matrix, n): 
    ret = 0 
    for i in range(n): 
     for j in range(n): 
      x = min(i + 1, n - i) 
      y = min(j + 1, n - j) 
      # k < half 
      half = (n + 1)/2 
      for k in range(1, half + 1): 
       count = min(k, x) * min(k, y) 
       ret += count * matrix[i][j] 
      # k >= half 
      for k in range(half + 1, n + 1): 
       count = min(n + 1 - k, x) * min(n + 1 - k, y) 
       ret += count * matrix[i][j] 
    return ret 

Si può vedere sum(min(k, x) * min(k, y)) può essere calcolata in O (1), quando 1 < = k < = n/2

Così c'è venuto a che O (n^2) codice:

def get_square_sum(n): 
    return n * (n + 1) * (2 * n + 1)/6 


def get_linear_sum(a, b): 
    return (b - a + 1) * (a + b)/2 


def get_count(x, y, k_end): 
    # k <= min(x, y), count is k*k 
    sum1 = get_square_sum(min(x, y)) 

    # k > min(x, y) and k <= max(x, y), count is k * min(x, y) 
    sum2 = get_linear_sum(min(x, y) + 1, max(x, y)) * min(x, y) 

    # k > max(x, y), count is x * y 
    sum3 = x * y * (k_end - max(x, y)) 

    return sum1 + sum2 + sum3 


def get_sum3(matrix, n): 
    ret = 0 
    for i in range(n): 
     for j in range(n): 
      x = min(i + 1, n - i) 
      y = min(j + 1, n - j) 
      half = n/2 

      # k < half 
      ret += get_count(x, y, half) * matrix[i][j] 
      # k >= half 
      ret += get_count(x, y, half + half % 2) * matrix[i][j] 

    return ret 

prova:

a = [[1, 1, 1], [2, 2, 2], [3, 3, 3]] 
n = 1000 
b = [[1] * n] * n 
print get_sum3(a, 3) # 68 
print get_sum3(b, n) # 33500333666800 

È possibile riscrivere la mia O (n^2) il codice Python per C e credo il risultato sarà una soluzione molto efficiente ...

Answer 3

Dopo l'ottima idea di @Divakar, vi suggerirei di usare integral images per velocizzare circonvoluzioni. Se la matrice è molto grande, devi conveterla più volte (una volta per ogni dimensione del kernel). Diverse convoluzioni (o valutazioni di somme all'interno di un quadrato) possono essere calcolate in modo molto efficiente usando immagini integrali (cioè tabelle dell'area sommate).

volta un'immagine integrale M viene calcolata la somma di tutti i valori all'interno di un'area (x0, y0) - (x1, y1) può essere calcolato con soli 4 calcoli aritmetici, indipendentemente dalla dimensione della finestra (immagine dal wikipedia):

M[x1, y1] - M[x1, y0] - M[x0, y1] + M[x0, y0] 

Questo può essere molto facilmente vettorializzare in NumPy. Un'immagine integrale può essere calcolata con cumsum. Seguendo l'esempio:

tab = np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3]]) 
M = tab.cumsum(0).cumsum(1) # Create integral images 
M = np.pad(M, ((1,0), (1,0)), mode='constant') # pad it with a row and column of zeros 

M viene riempito con una riga ed una colonna di zeri per gestire la prima fila (dove x0 = 0 o y0 = 0).

Quindi, data una dimensione di finestra W, la somma di ogni finestra di dimensioni W può essere calcolato in modo efficiente e completamente vectorized con NumPy come:

all_sums = M[W:, W:] - M[:-W, W:] - M[W:, :-W] + M[:-W, :-W] 

noti che l'operazione vettorizzati sopra, calcola la somma di ogni finestra, cioè ogni A, B, C e D della matrice. La somma di tutte le finestre viene quindi calcolata come

total = all_sums.sum() 

noti che per N dimensioni diverse, differenti per circonvoluzioni, l'immagine integrale deve essere calcolato solo una volta, quindi, il codice può essere scritto molto efficiente:

def get_all_sums(A): 
    M = A.cumsum(0).cumsum(1) 
    M = np.pad(M, ((1,0), (1,0)), mode='constant') 

    total = 0 
    for W in range(1, A.shape[0] + 1): 
     tmp = M[W:, W:] + M[:-W, :-W] - M[:-W, W:] - M[W:, :-W] 
     total += tmp.sum() 

    return total 

l'uscita per l'esempio:

>>> get_all_sums(tab) 
68 

---

Alcune temporizzazioni che confrontano le convoluzioni con immagini integrali con matrici di dimensioni diverse.getAllSums refeers a metodo convoluzionale di Divakar, mentre get_all_sums alle immagini integrali metodo basato descritto sopra:

>>> R1 = np.random.randn(10, 10) 
>>> R2 = np.random.randn(100, 100) 

1) Con R1 10x10 matrice:

>>> %time getAllSums(R1) 
CPU times: user 353 µs, sys: 9 µs, total: 362 µs 
Wall time: 335 µs 
2393.5912717342017 

>>> %time get_all_sums(R1) 
CPU times: user 243 µs, sys: 0 ns, total: 243 µs 
Wall time: 248 µs 
2393.5912717342012 

2) Con R2 100x100 matrice:

>>> %time getAllSums(R2) 
CPU times: user 698 ms, sys: 0 ns, total: 698 ms 
Wall time: 701 ms 
176299803.29826894 

>>> %time get_all_sums(R2) 
CPU times: user 2.51 ms, sys: 0 ns, total: 2.51 ms 
Wall time: 2.47 ms 
176299803.29826882 

Si noti che l'utilizzo di immagini integrali è 300 volte più veloce delle convoluzioni per matrici di dimensioni sufficienti.