NumPy Matrix Mnożenie Wydajność dla macierzy o znanej strukturze

Question

Mam dwie macierze NxN, które chcę, aby pomnożyć razem: A i B. W NumPy użyłem:

import numpy as np 
C = np.dot(A, B) 

Jednak zdarza mi się, że dla macierzy B tylko rząd n i kolumna n są niezerowe (wynika to bezpośrednio z formuły analitycznej, która wytworzyła matrycę i jest bez wątpienia zawsze taka sytuacja).

nadziei na wykorzystanie tego faktu i zmniejszenia liczby powieleń wymaganych do wytworzenia c, I otrzymuje powyżej:

import numpy as np 
for row in range(0, N): 
    for col in range(0, N): 
     if col != n: 
      C[row, col] = A[row, n]*B[n, col] #Just one scalar multiplication 
     else: 
      C[row, col] = np.dot(A[row], B[:, n]) 

analitycznie, należy zmniejszyć całkowitą złożoność w następujący sposób: W przypadku ogólnym (bez użycia fantazyjnych sztuczek, tylko podstawowe mnożenie macierzy) C = AB, gdzie A i B są oba NxN, powinno być O (N^3). Oznacza to, że wszystkie N wierszy musi pomnożyć wszystkich kolumn N, a każdy z tych produktów kropek zawiera N mnożenia => O (N N N) = O (N^3). #

Wykorzystanie struktury B jako Zrobiłem powyżej jednak powinien pójść jako O (N^2 + N^2) = O (2N^2) = O (N^2). Oznacza to, że wszystkie wiersze N ​​muszą pomnożyć wszystkie kolumny N, jednak dla wszystkich z nich (z wyjątkiem tych obejmujących "B [:, n]") wymagane jest tylko jedno zwielokrotnienie skalarne: tylko jeden element "B [:, m]" jest niezerowe dla m! = n. Gdy n == m, które wystąpią N razy (jeden raz dla każdego wiersza A, który musi pomnożyć kolumnę n z B), musi wystąpić N multiplikacji skalarnych. #

Jednak pierwszy blok kodu (za pomocą np.dot (A, B)) jest znacznie szybszy. Jestem świadomy (poprzez informacje takie jak: Why is matrix multiplication faster with numpy than with ctypes in Python?), że szczegóły implementacji niskiego poziomu np.dot są prawdopodobnie za to odpowiedzialne. Moje pytanie brzmi następująco: Jak mogę wykorzystać strukturę macierzy B, aby poprawić wydajność mnożenia bez poświęcania wydajności implementacji NumPy, bez budowania własnego mnożenia macierzy niskiego poziomu w c?

Ta metoda jest częścią optymalizacji numerycznej wielu zmiennych, a zatem O (N^3) jest trudny, natomiast O (N^2) prawdopodobnie wykona zadanie.

Dziękuję za pomoc. Ponadto, jestem nowy SO, więc proszę wybaczyć wszelkie błędy początkujących.

Answer 1

Gdybym zrozumiał A i B poprawnie, to nie rozumieją for pętle i dlatego są nie tylko pomnożenie przez dwóch niezerowych wektorów:

# say A & B are like this: 
n, N = 3, 5 
A = np.array(np.random.randn(N, N)) 

B = np.zeros_like(A) 
B[ n ] = np.random.randn(N) 
B[:, n] = np.random.randn(N) 

mieć niezerową wiersza i kolumny B:

rowb, colb = B[n,:], np.copy(B[:,n]) 
colb[ n ] = 0 

wielokrotnie A przez oba wektorze:

X = np.outer(A[:,n], rowb) 
X[:,n] += np.dot(A, colb) 

do zweryfikowania sprawdzenie:

X - np.dot(A, B) 

z N=100:

%timeit np.dot(A, B) 
1000 loops, best of 3: 1.39 ms per loop 

%timeit colb = np.copy(B[:,n]); colb[ n ] = 0; X = np.outer(A[:,n], B[n,:]); X[:,n] += np.dot(A, colb) 
10000 loops, best of 3: 98.5 µs per loop 

Answer 2

ja zmierzyłem czas, a używając sparse jest szybsze:

import numpy as np 
from scipy import sparse 

from timeit import timeit 

A = np.random.rand(100,100) 
B = np.zeros(A.shape, dtype=np.float) 

B[3] = np.random.rand(100) 
B[:,3] = np.random.rand(100) 

sparse_B = sparse.csr_matrix(B) 

n = 1000 

t1 = timeit('np.dot(A, B)', 'from __main__ import np, A, B', number=n) 
print 'dense way : {}'.format(t1) 
t2 = timeit('A * sparse_B', 'from __main__ import A, sparse_B',number=n) 
print 'sparse way : {}'.format(t2) 

Wynik:

dense way : 1.15117192268 
sparse way : 0.113152980804 
>>> np.allclose(np.dot(A, B), A * sparse_B) 
True 

Wraz ze wzrostem liczby wierszy B wzrasta również przewaga czasowa mnożenia przy użyciu macierzy rzadkiej.