Добавление очень повторяющейся матрицы к разреженной в numpy/scipy?

Question

Я пытаюсь реализовать функцию в NumPy/Scipy, чтобы вычислить Jensen-Shannon divergence между одним (обучающим) вектором и большим числом других (наблюдательных) векторов. Векторы наблюдения хранятся в очень большом (500 000x65536) Scipy sparse matrix (плотная матрица не будет вписываться в память).

В рамках алгоритма мне нужно вычислить T + O _я для каждого вектора наблюдения O _я, где Т представляет собой обучающий вектор. Я не смог найти способ сделать это, используя обычные правила вещания NumPy, поскольку разреженные матрицы, похоже, не поддерживают их (если T остается как плотный массив, Scipy пытается сначала разрезать разреженную матрицу, если я превращаю T в разреженную матрицу, T + O _i терпит неудачу, потому что формы несовместимы).

В настоящее время я беру грубо неэффективный шаг черепицы тренировочную вектора в 500,000x65536 разреженную матрицу:

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) 

Но это занимает огромное количество памяти (около 6GB), когда это только хранение ~ 1500 «реальных» элементов. Это также довольно медленно.

Я попытался получить умный, используя stride_tricks, чтобы заставить indptr и члены данных матрицы CSR не использовать дополнительную память для повторных данных.

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 

Но это не работает, потому что strided просмотров mtraining.data и mtraining.indices являются неправильную форму (и в соответствии с this answer нет никакого способа, чтобы сделать его правильной формы). Пытаясь заставить их выглядеть плоскими, используя итератор .flat, терпит неудачу, потому что он не выглядит достаточно, как массив (например, он не имеет элемента dtype), а использование метода flatten() заканчивает создание копии.

Есть ли способ сделать это?

Answer 1

Ваш другой вариант, который я даже не рассматривал, заключается в том, чтобы реализовать сумму в разреженном формате самостоятельно, чтобы вы могли в полной мере использовать периодический характер вашего массива. Это может быть очень легко сделать, если вы злоупотребляете это своеобразное поведение разреженных матриц SciPy в:

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

Так что вам даже не придется искать совпадения между вашим обучающим вектором и каждым из строк матрицы наблюдений чтобы добавить их: просто скопируйте все данные с помощью правильных указателей там, и то, что нужно получить, суммируется при доступе к данным.

EDIT

Учитывая медленность первого кода, вы можете торговать памятью скорости следующим образом:

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 

Это не так быстро, как и раньше, но он может сделать 10000 строк около 1 с.:

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 Этот код очень, очень медленно

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)