sparse_matrix¶

動機¶

sparseとよく聞くけど、概要しか分かっていないのでnumpyの力を借りて学ぶことにした

scipyだった
https://ja.wikipedia.org/wiki/%E7%96%8E%E8%A1%8C%E5%88%97 はあまり詳しく書いてない
なので、http://docs.scipy.org/doc/scipy/reference/sparse.html から各項目について学ぶのがよさそう
データ構造とその利点は学べるけど、その構造を活かした計算については調べてない
概要は、利用されていない・関係を表さない0という値ばかりに対してメモリを使いすぎるのは無駄だからコンパクトにする

疎行列(sparse matrix) の準備¶

In [2]:

import numpy as np

sparse_arr = np.zeros((4, 4))
sparse_arr[2, 0] = 1
sparse_arr[1, 1] = 1
sparse_arr[3, 0] = 1
sparse_arr[3, 3] = 1
sparse_arr

Out[2]:

array([[ 0.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 1.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  1.]])

疎行列用のデータ構造¶

In [3]:

from scipy.sparse import (
    dok_matrix,
    lil_matrix,
    coo_matrix,
    csr_matrix,
    csc_matrix,
    dia_matrix,
    bsr_matrix
)

DOK(Dictionary of Key)¶

ゼロでない値を (行, 列): 値となる辞書に対応させる

constructorからの作成¶

なにも値を持たない状態の疎行列を作成し、値を設定していく

In [4]:

# http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.dok_matrix.html#scipy.sparse.dok_matrix
S = dok_matrix((2, 3), dtype=np.float32)
S

Out[4]:

<2x3 sparse matrix of type '<class 'numpy.float32'>'
    with 0 stored elements in Dictionary Of Keys format>

In [9]:

S[0, 1], S[1, 2] = 1, 1
S

Out[9]:

<2x3 sparse matrix of type '<class 'numpy.float32'>'
    with 2 stored elements in Dictionary Of Keys format>

In [10]:

for k, v in S.items():
    print(k, v)

(0, 1) 1.0
(1, 2) 1.0

In [11]:

S[0, 1] = 0
for k, v in S.items():
    print(k, v)

(1, 2) 1.0

ndarrayからの変換¶

In [15]:

S = dok_matrix(sparse_arr)
S[3, 3] = 0
S.shape, S

Out[15]:

((4, 4),
 <4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements in Dictionary Of Keys format>)

In [13]:

for k, v in S.items():
    print(k, v)

(3, 0) 1.0
(2, 0) 1.0
(1, 1) 1.0

Row-based linked list sparse matrix¶

In [22]:

# http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.lil_matrix.html#scipy.sparse.lil_matrix
S = lil_matrix(sparse_arr)
S.data, S.rows, sparse_arr

Out[22]:

(array([[], [1.0], [1.0], [1.0, 1.0]], dtype=object),
 array([[], [1], [0], [0, 3]], dtype=object),
 array([[ 0.,  0.,  0.,  0.],
        [ 0.,  1.,  0.,  0.],
        [ 1.,  0.,  0.,  0.],
        [ 1.,  0.,  0.,  1.]]))

sparse matrixからの変換¶

In [108]:

S = lil_matrix(dok_matrix(sparse_arr))
S

Out[108]:

<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 4 stored elements in LInked List format>

ndarrayへの変換¶

In [62]:

S.toarray()

Out[62]:

array([[ 0.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 1.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  1.]])

In [63]:

S.todense()

Out[63]:

matrix([[ 0.,  0.,  0.,  0.],
        [ 0.,  1.,  0.,  0.],
        [ 1.,  0.,  0.,  0.],
        [ 1.,  0.,  0.,  1.]])

sparse matrixへの変換¶

In [71]:

S.tolil()

Out[71]:

<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 4 stored elements in LInked List format>

In [64]:

S.tobsr()

Out[64]:

<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 4 stored elements (blocksize = 1x1) in Block Sparse Row format>

In [65]:

S.tocoo()

Out[65]:

<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 4 stored elements in COOrdinate format>

In [66]:

S.tocsc()

Out[66]:

<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 4 stored elements in Compressed Sparse Column format>

In [67]:

S.tocsr()

Out[67]:

<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 4 stored elements in Compressed Sparse Row format>

In [68]:

S.todia()

Out[68]:

<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 7 stored elements (3 diagonals) in DIAgonal format>

In [69]:

S.todok()

Out[69]:

<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 4 stored elements in Dictionary Of Keys format>

A sparse matrix in COOrdinate format.¶

Also known as the ‘ijv’ or ‘triplet’ format.

In [82]:

# http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html#scipy.sparse.coo_matrix
S = coo_matrix(sparse_arr)
S.data, S.row, S.col

Out[82]:

(array([ 1.,  1.,  1.,  1.]),
 array([1, 2, 3, 3], dtype=int32),
 array([1, 0, 0, 3], dtype=int32))

ここまでは単純なのですぐわかった

Compressed Sparse Row matrix¶

In [81]:

# http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html#scipy.sparse.csr_matrix
S = csr_matrix(sparse_arr)
S.data, S.indices, S.indptr

Out[81]:

(array([ 1.,  1.,  1.,  1.]),
 array([1, 0, 0, 3], dtype=int32),
 array([0, 0, 1, 2, 4], dtype=int32))

Compressed Sparse Column matrix¶

In [91]:

# http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html#scipy.sparse.csc_matrix
S = csc_matrix(sparse_arr)
S.data, S.indices, S.indptr

Out[91]:

(array([ 1.,  1.,  1.,  1.]),
 array([2, 3, 1, 3], dtype=int32),
 array([0, 2, 3, 3, 4], dtype=int32))

In [26]:

S.toarray()

Out[26]:

array([[ 0.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 1.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  1.]])

よくわからないので、Wikipediaの例でやってみる¶

In [34]:

a = np.array([
 [1, 2, 3, 0],
 [0, 0, 0, 1],
 [2, 0, 0, 2],
 [0, 0, 0, 1],
    # 以下、追加
 [0, 0, 0, 0],
 [0, 0, 0, 1],
])
a

Out[34]:

array([[1, 2, 3, 0],
       [0, 0, 0, 1],
       [2, 0, 0, 2],
       [0, 0, 0, 1],
       [0, 0, 0, 0],
       [0, 0, 0, 1]])

A = [1 2 3 1 2 2 1] 非零要素リスト
IA = [1 1 1 2 3 3 4] i行
JA = [1 2 3 4 1 4 4] j列

これはCOOなので、以下より求められる

In [35]:

coo_a = coo_matrix(a)
coo_a.data, coo_a.row, coo_a.col

Out[35]:

(array([1, 2, 3, 1, 2, 2, 1, 1]),
 array([0, 0, 0, 1, 2, 2, 3, 5], dtype=int32),
 array([0, 1, 2, 3, 0, 3, 3, 3], dtype=int32))

indexの基準がずれてるけど、ずらせば結果は同じ

配列IAというのは同じ数字がずっと続くので、これを圧縮して何番目の要素からiが始まるかを書きならべる
この場合だと、IAの中で、
1が初めに出現するのは1番目
2が初めに出現するのは4番目
3が初めに出現するのは5番目
4が初めに出現するのは7番目
よってIA' = [1 4 5 7]

実際にCSRにしてみる

In [36]:

csr_a = csr_matrix(coo_a)
csr_a.data, csr_a.indices, csr_a.indptr

Out[36]:

(array([1, 2, 3, 1, 2, 2, 1, 1], dtype=int64),
 array([0, 1, 2, 3, 0, 3, 3, 3], dtype=int32),
 array([0, 3, 4, 6, 7, 7, 8], dtype=int32))

ポインタ部分の結果がおなじになった
indicesは列を表す
列のどの部分までが、どの行を表すかを圧縮しつつ表現できるようになった
- 0行は0から3未満
- 1行は3から4未満
- 2行は4から6未満
- 3行は6から7未満

列もやってみる¶

やることは変わらない
ただし列方向になるので転置して考える

In [37]:

a.T

Out[37]:

array([[1, 0, 2, 0, 0, 0],
       [2, 0, 0, 0, 0, 0],
       [3, 0, 0, 0, 0, 0],
       [0, 1, 2, 1, 0, 1]])

In [38]:

coo_at = coo_matrix(a.T)
coo_at.data, coo_at.row, coo_at.col

Out[38]:

(array([1, 2, 2, 3, 1, 2, 1, 1]),
 array([0, 0, 1, 2, 3, 3, 3, 3], dtype=int32),
 array([0, 2, 0, 0, 1, 2, 3, 5], dtype=int32))

In [39]:

csc_a = csc_matrix(coo_a)
csc_a.data, csc_a.indices, csc_a.indptr

Out[39]:

(array([1, 2, 2, 3, 1, 2, 1, 1], dtype=int64),
 array([0, 2, 0, 0, 1, 2, 3, 5], dtype=int32),
 array([0, 2, 3, 4, 8], dtype=int32))

転置に対して考えるので、indicesは行になる
あとはCSRと同じ

Sparse matrix with DIAgonal storage¶

In [95]:

# http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.dia_matrix.html#scipy.sparse.dia_matrix
S = dia_matrix(sparse_arr)
S.data, S.offsets

Out[95]:

(array([[ 1.,  0.,  0.,  0.],
        [ 1.,  0.,  0.,  0.],
        [ 0.,  1.,  0.,  1.]]), array([-3, -2,  0], dtype=int32))

In [131]:

sparse_arr

Out[131]:

array([[ 0.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 1.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  1.]])

Block Sparse Row matrix¶

https://ja.wikipedia.org/wiki/%E5%8C%BA%E5%88%86%E8%A1%8C%E5%88%97
どのときに役に立つか書いてある資料を見つけられなかった
- 公式の導入に役に立つものもあるし、とりあえず保留

In [97]:

# http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.bsr_matrix.html#scipy.sparse.bsr_matrix
S = bsr_matrix(sparse_arr)
S.data, S.indices, S.indptr, S.blocksize

Out[97]:

(array([[[ 1.]],

        [[ 1.]],

        [[ 1.]],

        [[ 1.]]]),
 array([1, 0, 0, 3], dtype=int32),
 array([0, 0, 1, 2, 4], dtype=int32),
 (1, 1))