Làm cách nào để tăng tốc độ tạo ma trận chuyển tiếp trong Numpy?

Question

Sau đây là cách cơ bản nhất mà tôi biết để đếm chuyển tiếp trong một chuỗi Markov và sử dụng nó để cư một ma trận chuyển tiếp:

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 

Tôi đã cố gắng tăng tốc nó lên trong 3 cách khác nhau:

1) Sử dụng một ma trận thưa thớt một lót dựa trên mã này Matlab:

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

nào trong numPy/scipy, trông như thế này:

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

Và tôi đã cố gắng nhiều hơn một vài tinh chỉnh Python, như sử dụng zip():

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

Và 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 

Nhưng không ai trong số những 3 phương pháp tăng tốc mọi thứ lên. Trong thực tế, tất cả mọi thứ nhưng giải pháp zip() là ít nhất 10X chậm hơn so với giải pháp ban đầu của tôi.

Có giải pháp nào khác đáng xem xét không?

giải pháp thay đổi để xây dựng một ma trận chuyển đổi từ nhiều chuỗi
Câu trả lời tốt nhất cho câu hỏi trên đặc biệt là của DSM. Tuy nhiên, đối với bất kỳ ai muốn điền ma trận chuyển đổi dựa trên danh sách hàng triệu chuỗi markov, cách nhanh nhất là:

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

Làm thế nào về điều này, tận dụng lợi thế của np.bincount? Không siêu mạnh mẽ, nhưng chức năng. [Nhờ @Warren Weckesser cho việc cài đặt.]

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) 

mang đến cho:

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 

[sửa]

Tránh flat (với chi phí không làm việc tại chỗ) có thể tiết kiệm một số thời gian cho các ma trận nhỏ:

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

Answer 2

Đây là phương pháp nhanh hơn. Ý tưởng là đếm số lần xuất hiện của mỗi lần chuyển đổi và sử dụng số đếm trong bản cập nhật vectơ của ma trận. (Tôi giả sử rằng quá trình chuyển đổi tương tự có thể xảy ra nhiều lần trong markov_chain.) Lớp Counter từ thư viện collections được sử dụng để đếm số lần xuất hiện của mỗi lần chuyển đổi.

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

Thời gian Ví dụ, trong 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 

Nó nhanh hơn, nhưng không phải là thứ tự độ lớn nhanh hơn. Để tăng tốc thực, bạn có thể xem xét việc thực hiện điều này trong Cython.

Answer 3

Ok, vài ý tưởng để làm xáo trộn, với một số cải tiến nhẹ (ít chi phí undestanding con người)

Hãy bắt đầu với một véc tơ ngẫu nhiên các số nguyên từ 0 đến 9 có độ dài 3000:

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

Phương pháp của bạn, trên máy của tôi, có hiệu suất thời gian của 11.4 ms.

Điều đầu tiên cho một sự cải thiện chút là tránh để đọc dữ liệu hai lần, lưu trữ nó trong một biến tạm thời:

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

này cung cấp cho bạn một sự cải thiện ~ 10% và giảm thời gian để 10,9 ms.

Một cách tiếp cận involuted hơn sử dụng những bước tiến:

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 

Các bước tiến cho phép bạn đọc các con số liên tiếp của mảng lừa mảng khi nghĩ rằng các hàng bắt đầu theo một cách khác (ok, nó không phải là tốt được mô tả, nhưng nếu bạn dành chút thời gian để đọc về các bước tiến bạn sẽ nhận được) Cách tiếp cận này mất hiệu suất, đi đến 12.2 ms, nhưng đó là hành lang để lừa hệ thống nhiều hơn. phẳng cả ma trận chuyển tiếp và các mảng strided một mảng chiều, bạn có thể tăng tốc độ hiệu suất hơn một chút:

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

này đi xuống 7,75 ms. Đó không phải là thứ tự độ lớn, nhưng dù sao cũng tốt hơn 30% :)

Answer 4

Chỉ để đá và vì tôi đã muốn dùng thử, tôi đã áp dụng Numba cho sự cố của bạn. Trong mã, mà liên quan đến việc chỉ cần thêm một trang trí (mặc dù tôi đã thực hiện một cuộc gọi trực tiếp để tôi có thể kiểm tra JIT biến thể numba rằng cung cấp ở đây):

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

Và sau đó 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 

Các Phương pháp autojit thực hiện một số phương pháp đoán dựa trên đầu vào thời gian chạy và hàm jit có các loại được quyết định. Bạn phải cẩn thận một chút vì tê liệt ở những giai đoạn đầu này không giao tiếp rằng có lỗi với jit nếu bạn nhập nhầm loại đầu vào. Nó sẽ chỉ nhổ ra một câu trả lời không chính xác.

Điều đó nói rằng, mặc dù tăng tốc 35x và 485x mà không cần thay đổi mã và chỉ cần thêm một cuộc gọi đến numba (cũng có thể được gọi là trang trí) là khá ấn tượng trong cuốn sách của tôi. Bạn có thể có thể nhận được kết quả tương tự bằng cách sử dụng cython, nhưng nó sẽ đòi hỏi một chút boilerplate và viết một tập tin setup.py.

Tôi cũng thích giải pháp này vì mã vẫn có thể đọc được và bạn có thể viết nó theo cách bạn nghĩ ban đầu về việc triển khai thuật toán.