Wie kann ich die Erstellung der Übergangsmatrix in Numpy beschleunigen?

Question

Im Folgenden ist der einfachste Weg den ich kenne, Übergänge in einer Markow-Kette zu zählen und es zu verwenden, einen Übergang Matrix zu besiedeln:

def increment_counts_in_matrix_from_chain(markov_chain, transition_counts_matrix): 
    for i in xrange(1, len(markov_chain)): 
     old_state = markov_chain[i - 1] 
     new_state = markov_chain[i] 
     transition_counts_matrix[old_state, new_state] += 1 

ich habe versucht, bis auf 3 verschiedene Arten zu beschleunigen:

1) eine spärliche Matrix Einzeiler auf diesem Code Matlab basierte Anwendung:

transition_matrix = full(sparse(markov_chain(1:end-1), markov_chain(2:end), 1)) 

die in Numpy/SciPy, sieht wie folgt aus:

def get_sparse_counts_matrix(markov_chain, number_of_states): 
    return coo_matrix(([1]*(len(markov_chain) - 1), (markov_chain[0:-1], markov_chain[1:])), shape=(number_of_states, number_of_states)) 

Und ich habe noch ein paar Python zwickt versucht, wie mit Reißverschluss():

for old_state, new_state in zip(markov_chain[0:-1], markov_chain[1:]): 
    transition_counts_matrix[old_state, new_state] += 1 

And Queues:

old_and_new_states_holder = Queue(maxsize=2) 
old_and_new_states_holder.put(markov_chain[0]) 
for new_state in markov_chain[1:]: 
    old_and_new_states_holder.put(new_state) 
    old_state = old_and_new_states_holder.get() 
    transition_counts_matrix[old_state, new_state] += 1 

Aber keine dieser drei Methoden beschleunigt Dinge. In der Tat war alles außer der Zip() Lösung mindestens 10x langsamer als meine ursprüngliche Lösung.

Gibt es noch andere Lösungen, die sich zu untersuchen lohnt?

Modified Lösung
Die beste Antwort auf die obige Frage war eine Übergangsmatrix aus vielen Ketten für den Aufbau speziell von DSM. Doch für alle, die eine Übergangsmatrix auf der Grundlage einer Liste von Millionen von Markow-Ketten füllen will, ist der schnellste Weg, um dies:

def fast_increment_transition_counts_from_chain(markov_chain, transition_counts_matrix): 
    flat_coords = numpy.ravel_multi_index((markov_chain[:-1], markov_chain[1:]), transition_counts_matrix.shape) 
    transition_counts_matrix.flat += numpy.bincount(flat_coords, minlength=transition_counts_matrix.size) 

def get_fake_transitions(markov_chains): 
    fake_transitions = [] 
    for i in xrange(1,len(markov_chains)): 
     old_chain = markov_chains[i - 1] 
     new_chain = markov_chains[i] 
     end_of_old = old_chain[-1] 
     beginning_of_new = new_chain[0] 
     fake_transitions.append((end_of_old, beginning_of_new)) 
    return fake_transitions 

def decrement_fake_transitions(fake_transitions, counts_matrix): 
    for old_state, new_state in fake_transitions: 
     counts_matrix[old_state, new_state] -= 1 

def fast_get_transition_counts_matrix(markov_chains, number_of_states): 
    """50% faster than original, but must store 2 additional slice copies of all markov chains in memory at once. 
    You might need to break up the chains into manageable chunks that don't exceed your memory. 
    """ 
    transition_counts_matrix = numpy.zeros([number_of_states, number_of_states]) 
    fake_transitions = get_fake_transitions(markov_chains) 
    markov_chains = list(itertools.chain(*markov_chains)) 
    fast_increment_transition_counts_from_chain(markov_chains, transition_counts_matrix) 
    decrement_fake_transitions(fake_transitions, transition_counts_matrix) 
    return transition_counts_matrix 

Answer 1

Wie wäre es mit so etwas, unter Ausnutzung der np.bincount? Nicht super-robust, aber funktional. [Dank Weckesser für das Setup @Warren.]

import numpy as np 
from collections import Counter 

def increment_counts_in_matrix_from_chain(markov_chain, transition_counts_matrix): 
    for i in xrange(1, len(markov_chain)): 
     old_state = markov_chain[i - 1] 
     new_state = markov_chain[i] 
     transition_counts_matrix[old_state, new_state] += 1 

def using_counter(chain, counts_matrix): 
    counts = Counter(zip(chain[:-1], chain[1:])) 
    from_, to = zip(*counts.keys()) 
    counts_matrix[from_, to] = counts.values() 

def using_bincount(chain, counts_matrix): 
    flat_coords = np.ravel_multi_index((chain[:-1], chain[1:]), counts_matrix.shape) 
    counts_matrix.flat = np.bincount(flat_coords, minlength=counts_matrix.size) 

def using_bincount_reshape(chain, counts_matrix): 
    flat_coords = np.ravel_multi_index((chain[:-1], chain[1:]), counts_matrix.shape) 
    return np.bincount(flat_coords, minlength=counts_matrix.size).reshape(counts_matrix.shape) 

die gibt:

In [373]: t = np.random.randint(0,50, 500) 
In [374]: m1 = np.zeros((50,50)) 
In [375]: m2 = m1.copy() 
In [376]: m3 = m1.copy() 

In [377]: timeit increment_counts_in_matrix_from_chain(t, m1) 
100 loops, best of 3: 2.79 ms per loop 

In [378]: timeit using_counter(t, m2) 
1000 loops, best of 3: 924 us per loop 

In [379]: timeit using_bincount(t, m3) 
10000 loops, best of 3: 57.1 us per loop 

[Bearbeiten]

flat Vermeidung (auf Kosten der nicht an Ort und Stelle in Betrieb) kann sparen einige Zeit für kleine Matrizen:

In [80]: timeit using_bincount_reshape(t, m3) 
10000 loops, best of 3: 22.3 us per loop 

Answer 2

Hier ist eine schnellere Methode. Die Idee besteht darin, die Anzahl der Vorkommen jedes Übergangs zu zählen und die Zählungen in einer vektorisierten Aktualisierung der Matrix zu verwenden. (Ich nehme an, dass derselbe Übergang mehrfach in markov_chain auftreten kann.) Die Counter-Klasse aus der Bibliothek collections wird verwendet, um die Anzahl der Vorkommen jedes Übergangs zu zählen.

from collections import Counter 

def update_matrix(chain, counts_matrix): 
    counts = Counter(zip(chain[:-1], chain[1:])) 
    from_, to = zip(*counts.keys()) 
    counts_matrix[from_, to] += counts.values() 

Timing-Beispiel in ipython:

In [64]: t = np.random.randint(0,50, 500) 

In [65]: m1 = zeros((50,50)) 

In [66]: m2 = zeros((50,50)) 

In [67]: %timeit increment_counts_in_matrix_from_chain(t, m1) 
1000 loops, best of 3: 895 us per loop 

In [68]: %timeit update_matrix(t, m2) 
1000 loops, best of 3: 504 us per loop 

Es ist schneller, aber nicht um Größenordnungen schneller. Für eine echte Beschleunigung könnten Sie dies in Cython implementieren.

Answer 3

Ok, einige Ideen zu manipulieren, mit einiger leichten Verbesserung (auf Kosten des menschlichen undestanding)

Lassen Sie sich mit einem zufälligen Vektor von ganzen Zahlen zwischen 0 und 9 der Länge 3000 beginnen:

L = 3000 
N = 10 
states = array(randint(N),size=L) 
transitions = np.zeros((N,N)) 

Ihre Methode hat auf meinem Rechner eine Zeiteinheit Leistung von 11,4 ms.

Das erste, was für eine kleine Verbesserung ist, die Daten zu vermeiden, zweimal zu lesen, es in einer temporären Variablen zu speichern:

old = states[0] 
for i in range(1,len(states)): 
    new = states[i] 
    transitions[new,old]+=1 
    old=new 

Dies gibt Ihnen eine ~ 10 Verbesserung% und fällt die Zeit 10,9 ms.

Ein involuted Ansatz verwendet die Schritte:

def rolling(a, window): 
    shape = (a.size - window + 1, window) 
    strides = (a.itemsize, a.itemsize) 
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides) 

state_2 = rolling(states, 2) 
for i in range(len(state_2)): 
    l,m = state_2[i,0],state_2[i,1] 
    transitions[m,l]+=1 

Die Fortschritte können Sie die fortlaufenden Nummern des Arrays lesen das Array austricksen zu glauben, dass die Zeilen in einer anderen Art und Weise zu starten (ok, es nicht gut ist beschrieben, aber wenn Sie etwas Zeit brauchen, um über Fortschritte zu lesen, erhalten Sie es) Dieser Ansatz verliert die Leistung, gehen zu 12,2 ms, aber es ist der Flur, um das System noch mehr zu tricksen. sowohl die Übergangsmatrix und der strided Array eindimensional Arrays Abflachen, können Sie die Leistung etwas mehr beschleunigen:

transitions = np.zeros(N*N) 
state_2 = rolling(states, 2) 
state_flat = np.sum(state_2 * array([1,10]),axis=1) 
for i in state_flat: 
    transitions[i]+=1 
transitions.reshape((N,N)) 

Dies geht bis auf 7,75 ms . Es ist nicht eine Größenordnung, aber es ist sowieso ein 30% besser :)

Answer 4

Nur für Kicks, und weil ich es ausprobieren wollte, bewarb ich mich Numba zu Ihrem Problem. In Code, das umfasst das Hinzufügen nur ein Dekorateur (obwohl ich ein direkter Anruf gemacht habe, so konnte ich die JIT Test-Varianten, die numba hier liefert):

import numpy as np 
import numba 

def increment_counts_in_matrix_from_chain(markov_chain, transition_counts_matrix): 
    for i in xrange(1, len(markov_chain)): 
     old_state = markov_chain[i - 1] 
     new_state = markov_chain[i] 
     transition_counts_matrix[old_state, new_state] += 1 

autojit_func = numba.autojit()(increment_counts_in_matrix_from_chain) 
jit_func = numba.jit(argtypes=[numba.int64[:,::1],numba.double[:,::1]])(increment_counts_in_matrix_from_chain) 

t = np.random.randint(0,50, 500) 
m1 = np.zeros((50,50)) 
m2 = np.zeros((50,50)) 
m3 = np.zeros((50,50)) 

Und dann Timings:

In [10]: %timeit increment_counts_in_matrix_from_chain(t,m1) 
100 loops, best of 3: 2.38 ms per loop 

In [11]: %timeit autojit_func(t,m2)       

10000 loops, best of 3: 67.5 us per loop 

In [12]: %timeit jit_func(t,m3) 
100000 loops, best of 3: 4.93 us per loop 

Die autojit Methode führt einige Raten basierend auf Laufzeit-Eingaben und die jit Funktion hat diktiert. Sie müssen ein wenig vorsichtig sein, da numba in diesen frühen Stadien nicht kommuniziert, dass es einen Fehler mit jit gab, wenn Sie den falschen Typ für eine Eingabe übergeben. Es wird nur eine falsche Antwort ausspucken.

Das sagte aber, eine 35x und 485x Beschleunigung ohne Code-Änderung und nur einen Anruf an numba hinzufügen (kann auch als Dekorateur genannt werden) ist ziemlich beeindruckend in meinem Buch. Sie könnten wahrscheinlich ähnliche Ergebnisse mit Cython erhalten, aber es würde etwas mehr Vortex erfordern und eine setup.py-Datei schreiben.

Ich mag auch diese Lösung, weil der Code lesbar bleibt und Sie es so schreiben können, wie Sie ursprünglich über die Implementierung des Algorithmus dachten.