Jak mogę wektoryzować i przyspieszyć obliczenia dla tej dużej macierzy?

Question

Obecnie próbuję wykonać obliczyć sumę wszystkich sum podrzędnych w tablicy 10.000 x 10.000 wartości. Jako przykład, jeśli moja tablica była:

1 1 1 
2 2 2 
3 3 3 

Chcę być wynikiem:

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 

Tak, jak za pierwszym razem napisałem bardzo prosty kod Pythona, aby to zrobić. Tak jak było w O (k^2.n^2) (n jest wielkością dużej tablicy i k wielkości podskali, które otrzymujemy), przetwarzanie było strasznie długie. Pisałem w inny algorytm O (n^2), aby ją przyspieszyć:

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 

Więc ten kod działa poprawnie. Biorąc pod uwagę tablicę i wielkość subsquares, zwróci sumę wartości we wszystkich tych subsquares. Zasadniczo iteracja nad wielkością subskry w celu uzyskania wszystkich możliwych wartości.

Problem polega na tym, że to znowu będzie długie na duże tablice (ponad 20 dni dla tablicy 10.000 x 10.000). Przeszukałem go i dowiedziałem się, że mogę wektoryzować iteracje po tablicach z numpy. Jednak nie mogłem wymyślić, jak to zrobić w moim przypadku ...

Jeśli ktoś może mi pomóc przyspieszyć mój algorytm lub dać mi dobrą dokumentację na ten temat, będę zadowolony!

Dziękujemy!

Answer 1

Te suwane podsumowania najlepiej nadają się do obliczenia jako sumy splotów 2D, a te można skutecznie obliczyć za pomocą scipy's convolve2d. Tak więc, dla określonego rozmiaru, można uzyskać sumowanie, jak tak -

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

dostać sumowanie we wszystkich rozmiarach, myślę, że najlepszym sposobem, zarówno jeśli chodzi o pamięci i wydajności wydajności byłoby użyć pętlę pętla nad wszystkimi możliwymi rozmiarami. Dlatego, aby uzyskać ostateczną sumowanie, trzeba -

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

Sample 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

Bazując na pomysł, aby obliczyć, ile razy każdy numer liczy się, doszedłem do tego prostego kodu:

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 

---

rozwiązanie Divakar jest fantastyczny, jednak myślę, że kopalnia może być bardziej skuteczny, przynajmniej w asymp złożoność czasowa totalizmu (O (n^3) w porównaniu z O (n^3logn) Divakara).

---

Dostaję teraz rozwiązanie O (n^2) ...

Zasadniczo możemy, że:

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 

Widać sum(min(k, x) * min(k, y)) można obliczyć w czasie O (1), gdy 1 < = k < = n/2

Więc doszliśmy do tego O (n^2) Kod:

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 

testu:

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 

można przepisać mój O (n^2) kod Pythona do C i wierzę, że spowoduje to bardzo wydajne rozwiązanie ...

Answer 3

Po doskonałej idei @Divakar, chciałbym zaproponować za pomocą integral images do przyspieszenia zwoje. Jeśli macierz jest bardzo duża, musisz ją spętać kilka razy (raz dla każdego rozmiaru jądra). Kilka zwojów (lub oszacowanie sum wewnątrz kwadratu) można bardzo skutecznie obliczyć przy użyciu obrazów całkowych (zwanych również tabelami zsumowanymi obszarami).

Po integralny obraz M obliczana jest suma wszystkich wartości wewnątrz regionu (x0, y0) - (x1, y1) można obliczyć z zaledwie 4 aritmetic obliczenia, niezależnie od rozmiaru okna (zdjęcie z Wikipedii):

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

ten można bardzo łatwo wektorowy w numpy. Obraz całkowy można obliczyć za pomocą cumsum. W związku z na przykład:

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 napawa wiersza i kolumny zer do obsługi pierwszego rzędu (gdzie x0 = 0 lub y0 = 0).

Następnie, biorąc pod uwagę rozmiar okna W, suma każdego okna o rozmiarze W można obliczyć efektywnie iw pełni wektorowy z numpy jak:

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

nocie że vectorized operacja powyżej, oblicza sumę każdego okno, tj. każdy A, B, C i D macierzy. Suma wszystkich oknach jest następnie obliczana jako

total = all_sums.sum() 

Zauważ, że dla N różnych rozmiarach, różnych do zwojów, integralny obraz ma być obliczana tylko raz, a więc kod może być napisany bardzo skutecznie, jak:

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 

wyjście na przykład:

>>> get_all_sums(tab) 
68 

---

Niektóre synchronizacje porównujące skręty z obrazami zintegrowanymi z macierzami o różnych rozmiarach.getAllSums refeers do splotowego metodą Divakar, natomiast get_all_sums całce obrazów metoda oparta opisanych powyżej:

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

1) z R1 10x10 Matrix

>>> %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) Z R2 100x100 Matrix

>>> %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 

Należy pamiętać, że użycie obrazów zintegrowanych jest 300 razy szybsze niż zwoje dla wystarczająco dużych macierzy.