6

NumPy / Scipyに関数を実装して、単一の（トレーニング）ベクトルと他の多数の（観測）ベクトルの間のイェンセンシャノンの発散を計算しようとしています。観測ベクトルは、非常に大きな（500,000x65536）Scipyスパース行列に格納されます（密行列はメモリに収まりません）。

アルゴリズムの一部として、観測ベクトルOiごとにT+ O iを計算する必要があります。ここで、Tはトレーニングベクトルです。スパース行列はそれらをサポートしていないように見えるため、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)

ただし、これは、最大1500の「実際の」要素しか格納しない場合、大量のメモリ（約6GB）を消費します。構築もかなり遅いです。

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は間違った形状であるため、これは機能しません（この回答によると、正しい形状にする方法はありません）。.flatイテレータを使用してそれらをフラットに見せようとすると、配列のように見えないため（たとえば、dtypeメンバーがないため）失敗し、flatten（）メソッドを使用するとコピーが作成されます。

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

4

1 に答える 1

3

私が考えもしなかった他のオプションは、配列の周期的な性質を最大限に活用できるように、自分でスパース形式で合計を実装することです。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]])

したがって、トレーニングベクトルと観測行列の各行の間の一致を探してそれらを合計する必要はありません。適切なポインターを含むすべてのデータを詰め込むだけで、合計する必要があるものが合計されます。データにアクセスしたとき。

編集

最初のコードの速度が遅いことを考えると、次のようにメモリを速度と交換できます。

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

以前ほど高速ではありませんが、約1秒で10,000行を実行できます。

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

編集このコードは非常に遅いです

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)
于 2013-03-07T02:11:38.540 に答える