Dodawanie macierzy bardzo powtarzalnej do rzadkiej w numpy/scipy?

Question

Próbuję zaimplementować funkcję w NumPy/Scipy, aby obliczyć Jensen-Shannon divergence między pojedynczym (treningowym) wektorem i dużą liczbą innych wektorów (obserwacji). Wektory obserwacyjne są przechowywane w bardzo dużym (500 000 x 65536) Scipy sparse matrix (gęsta macierz nie mieści się w pamięci).

jako część algorytmu potrzebne do obliczenia T + O _i dla każdego wektora obserwacji O _i, gdzie T jest wektorem szkolenia. Nie byłem w stanie znaleźć sposobu, aby to zrobić używając zwykłych reguł rozgłaszania NumPy, ponieważ rzadkie macierze nie wydają się wspierać tych (jeśli T jest pozostawiona jako gęstą tablicę, Scipy próbuje najpierw zrobić gęstą macierz, która działa brak pamięci, jeśli ustawię T na rzadką macierz, T + O _i zawodzi, ponieważ kształty są niespójne).

Obecnie Biorę rażąco nieefektywne krok kafli wektor treningowy do 500,000x65536 rozrzedzony matrycy:

training = sp.csr_matrix(training.astype(np.float32)) 
tindptr = np.arange(0, len(training.indices)*observations.shape[0]+1, len(training.indices), dtype=np.int32) 
tindices = np.tile(training.indices, observations.shape[0]) 
tdata = np.tile(training.data, observations.shape[0]) 
mtraining = sp.csr_matrix((tdata, tindices, tindptr), shape=observations.shape) 

Ale to zajmuje ogromną ilość pamięci (około 6 GB), gdy jest przechowywanie tylko ~ 1500 "prawdziwych" elementów. Jest również dość powolny w konstruowaniu.

Próbowałem uzyskać sprytny przy użyciu stride_tricks, aby indptr w matrycy CSR i członkowie danych nie korzystali z dodatkowej pamięci na powtarzanych danych.

training = sp.csr_matrix(training) 
mtraining = sp.csr_matrix(observations.shape,dtype=np.int32) 
tdata = training.data 
vdata = np.lib.stride_tricks.as_strided(tdata, (mtraining.shape[0], tdata.size), (0, tdata.itemsize)) 
indices = training.indices 
vindices = np.lib.stride_tricks.as_strided(indices, (mtraining.shape[0], indices.size), (0, indices.itemsize)) 
mtraining.indptr = np.arange(0, len(indices)*mtraining.shape[0]+1, len(indices), dtype=np.int32) 
mtraining.data = vdata 
mtraining.indices = vindices 

Ale to nie działa, ponieważ strided odsłony mtraining.data i mtraining.indices są niewłaściwy kształt (i według this answer nie ma sposobu, aby to prawo kształt). Próba sprawienia, by wyglądały płasko, używając iteratora płaskiego, nie działa, ponieważ nie wygląda wystarczająco jak tablica (np. Nie ma elementu typu dtype), a użycie metody spłaszczania() kończy się tworzeniem kopii.

Czy jest jakiś sposób, aby to zrobić?

Answer 1

Inną opcją, której nawet nie brałem pod uwagę, jest samodzielne wprowadzenie sumy w osobnym formacie, abyś mógł w pełni wykorzystać okresową naturę swojej tablicy. Może to być bardzo proste do zrobienia, jeśli nadużycie to osobliwe zachowanie rzadkich macierzy scipy za:

>>> a = sps.csr_matrix([1,2,3,4]) 
>>> a.data 
array([1, 2, 3, 4]) 
>>> a.indices 
array([0, 1, 2, 3]) 
>>> a.indptr 
array([0, 4]) 

>>> b = sps.csr_matrix((np.array([1, 2, 3, 4, 5]), 
...      np.array([0, 1, 2, 3, 0]), 
...      np.array([0, 5])), shape=(1, 4)) 
>>> b 
<1x4 sparse matrix of type '<type 'numpy.int32'>' 
    with 5 stored elements in Compressed Sparse Row format> 
>>> b.todense() 
matrix([[6, 2, 3, 4]]) 

więc nawet nie trzeba szukać koincydencji między swoim wektora treningowego i każdy z wierszy macierzy obserwacji aby je dodać: po prostu zapakuj wszystkie dane z odpowiednimi wskaźnikami, a to, co wymaga zsumowania, zostanie zsumowane, gdy dane będą dostępne.

EDIT

względu na powolność pierwszego kodu, można handlować pamięci dla prędkości następująco:

def csr_add_sparse_vec(sps_mat, sps_vec) : 
    """Adds a sparse vector to every row of a sparse matrix""" 
    # No checks done, but both arguments should be sparse matrices in CSR 
    # format, both should have the same number of columns, and the vector 
    # should be a vector and have only one row. 

    rows, cols = sps_mat.shape 
    nnz_vec = len(sps_vec.data) 
    nnz_per_row = np.diff(sps_mat.indptr) 
    longest_row = np.max(nnz_per_row) 

    old_data = np.zeros((rows * longest_row,), dtype=sps_mat.data.dtype) 
    old_cols = np.zeros((rows * longest_row,), dtype=sps_mat.indices.dtype) 

    data_idx = np.arange(longest_row) < nnz_per_row[:, None] 
    data_idx = data_idx.reshape(-1) 
    old_data[data_idx] = sps_mat.data 
    old_cols[data_idx] = sps_mat.indices 
    old_data = old_data.reshape(rows, -1) 
    old_cols = old_cols.reshape(rows, -1) 

    new_data = np.zeros((rows, longest_row + nnz_vec,), 
         dtype=sps_mat.data.dtype) 
    new_data[:, :longest_row] = old_data 
    del old_data 
    new_cols = np.zeros((rows, longest_row + nnz_vec,), 
         dtype=sps_mat.indices.dtype) 
    new_cols[:, :longest_row] = old_cols 
    del old_cols 
    new_data[:, longest_row:] = sps_vec.data 
    new_cols[:, longest_row:] = sps_vec.indices 
    new_data = new_data.reshape(-1) 
    new_cols = new_cols.reshape(-1) 
    new_pointer = np.arange(0, (rows + 1) * (longest_row + nnz_vec), 
          longest_row + nnz_vec) 

    ret = sps.csr_matrix((new_data, new_cols, new_pointer), 
         shape=sps_mat.shape) 
    ret.eliminate_zeros() 

    return ret 

To nie jest tak szybki jak wcześniej, ale może to zrobić 10.000 wierszy około 1 s.:

In [2]: a 
Out[2]: 
<10000x65536 sparse matrix of type '<type 'numpy.float64'>' 
    with 15000000 stored elements in Compressed Sparse Row format> 

In [3]: b 
Out[3]: 
<1x65536 sparse matrix of type '<type 'numpy.float64'>' 
    with 1500 stored elements in Compressed Sparse Row format> 

In [4]: csr_add_sparse_vec(a, b) 
Out[4]: 
<10000x65536 sparse matrix of type '<type 'numpy.float64'>' 
    with 30000000 stored elements in Compressed Sparse Row format> 

In [5]: %timeit csr_add_sparse_vec(a, b) 
1 loops, best of 3: 956 ms per loop 

---

EDIT Kod ten jest bardzo powolny

def csr_add_sparse_vec(sps_mat, sps_vec) : 
    """Adds a sparse vector to every row of a sparse matrix""" 
    # No checks done, but both arguments should be sparse matrices in CSR 
    # format, both should have the same number of columns, and the vector 
    # should be a vector and have only one row. 

    rows, cols = sps_mat.shape 

    new_data = sps_mat.data 
    new_pointer = sps_mat.indptr.copy() 
    new_cols = sps_mat.indices 

    aux_idx = np.arange(rows + 1) 

    for value, col in itertools.izip(sps_vec.data, sps_vec.indices) : 
     new_data = np.insert(new_data, new_pointer[1:], [value] * rows) 
     new_cols = np.insert(new_cols, new_pointer[1:], [col] * rows) 
     new_pointer += aux_idx 

    return sps.csr_matrix((new_data, new_cols, new_pointer), 
          shape=sps_mat.shape)