Optimiere mittlere äußere Produkte

Question

Ich schreibe gerade ein kurzes Programm, um eine Analyse von zufälligen Matrixeigenwertverteilungen durchzuführen, aber die Parameterauswahl, die für meine Analyse benötigt wird, hat sich als extrem langsam herausgestellt. Grundsätzlich sollte ich die untenstehende Funktion durchlaufen, idealerweise für ca. 5000 mal, und schließlich die komplette Liste der Eigenwerte am Ende sammeln.

C = np.zeros((N,N)) 
time_series = np.random.normal(mu,sigma, (N + B*(M-1)) ) 

for k in range(int(M)): 
    C += np.outer(time_series[k*B : (N) + k*B], time_series[k*B : (N) + k*B]) 
C = C/M 

eg_v = np.linalg.eigvalsh(C) 

wo ich brauche N = 1000, B um 10, M = 100. jedoch mit dieser Wahl der Parameter das Programm dauert etwas wie 4-5 Stunden auf meinem ziemlich leistungs Laptop laufen zu lassen.

Abgesehen von Hardware-Einschränkungen, gibt es etwas, was ich in Bezug auf den Code tun kann, um das Ganze zu beschleunigen?

Vielen Dank im Voraus!

Answer 1

Sie die Schleife mit einer vektorisiert Lösung ersetzen kann np.tensordot

So verwenden die folgende -

C = np.zeros((N,N)) 
for k in range(int(M)): 
    C += np.outer(time_series[k*B : (N) + k*B], time_series[k*B : (N) + k*B]) 

durch ersetzt werden könnte -

# Get the starting indices for each iteration 
idx = (np.arange(M)*B)[:,None] + np.arange(N) 

# Get the range of indices across all iterations as a 2D array and index 
# time_series with it to give us "time_series[k*B : (N) + k*B]" equivalent 
time_idx = time_series[idx] 

# Use broadcasting to perform summation accumulation 
C = np.tensordot(time_idx,time_idx,axes=([0],[0])) 

Die tensordot durch eine ersetzt werden könnte einfaches Punktprodukt:

C = time_idx.T.dot(time_idx) 

Runtime Test

Funktionen:

def original_app(time_series,B,N,M): 
    C = np.zeros((N,N)) 
    for k in range(int(M)): 
     C += np.outer(time_series[k*B : (N) + k*B], time_series[k*B : (N) + k*B]) 
    return C 

def vectorized_app(time_series,B,N,M): 
    idx = (np.arange(M)*B)[:,None] + np.arange(N) 
    time_idx = time_series[idx] 
    return np.tensordot(time_idx,time_idx,axes=([0],[0])) 

Eingänge:

In [115]: # Inputs 
    ...: mu = 1.2 
    ...: sigma = 0.5 
    ...: N = 1000 
    ...: M = 100 
    ...: B = 10 
    ...: time_series = np.random.normal(mu,sigma, (N + B*(M-1)) ) 
    ...: 

Timings:

In [116]: out1 = original_app(time_series,B,N,M) 

In [117]: out2 = vectorized_app(time_series,B,N,M) 

In [118]: np.allclose(out1,out2) 
Out[118]: True 

In [119]: %timeit original_app(time_series,B,N,M) 
1 loops, best of 3: 1.56 s per loop 

In [120]: %timeit vectorized_app(time_series,B,N,M) 
10 loops, best of 3: 26.2 ms per loop 

So sehen wir eine 60x Beschleunigung für die in der Frage aufgeführten Eingänge!