numpy/scipyで非常に反復的な行列を疎な行列に追加していますか？

Question

私は、単一の（トレーニング）ベクトルと多数の他の（観測）ベクトルの間でJensen-Shannon divergenceを計算するために、NumPy/Scipyに関数を実装しようとしています。観測ベクトルは非常に大きい（500,000x65536）Scipy sparse matrixに格納されます（密行列はメモリに収まらない）。

I Tが訓練ベクトルであり、各観測ベクトルO I_{、のためにT + O Iを計算するために必要なアルゴリズムの一部として。私はNumPyの通常の放送規則を使ってこれを行う方法を見つけることができませんでした。なぜなら、疎な行列はそれらをサポートしていないように思われるからです（Tが密な配列として残っていると、Scipyは、私がTを疎な行列にすると、T + O iは形が矛盾するので失敗します）。 〜}

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) 

しかし、それは唯一の保存だとき、これは、（6ギガバイトの周りに）メモリの膨大な量を占める：

は現在、私は500,000x65536スパース行列に訓練ベクトルをタイリングする著しく非効率的なステップを取っています1500の「本当の」要素。それはまた、構築がかなり遅いです。

stride_tricksを使用して、CSRマトリクスのindptrとデータメンバーが繰り返しデータに余分なメモリを使用しないようにするために、賢明にしようとしました。

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 

しかし、ストライドビューmtraining.dataとmtraining.indicesが間違った形をしている（とthis answerに応じて右の形状にする方法はありません）ので、これは動作しません。 .flatイテレータを使用してフラットに見せようとすると、配列のように見えないため（dtypeメンバを持たないなど）、flatten（）メソッドを使用するとコピーが終了します。

これを行う方法はありますか？

Answer 1

私が考えていなかったあなたの他の選択肢は、あなた自身が疎な形式で合計を実装することです。これにより、配列の周期的な性質をフルに活用することができます。だから、

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

は、あなたもあなたのトレーニングベクトルと観測行列の列のそれぞれの間の同時発生を探す必要はありません：これは、あなたがscipyのダウンロードのスパース行列のこの特異な振る舞いを乱用した場合に、やることは非常に簡単にすることができますそれらを追加するには、そこに正しいポインタを持つすべてのデータを詰め込み、集計する必要があるデータは、データにアクセスすると合計されます。

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 

それは以前のように早くはありませんが、それはで10,000行を行うことができます約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)