Python, SciPy（scipy.sparse）で疎行列を生成・変換
SciPy（scipy.sparse）を使うと疎行列（スパース行列）を効率的に扱うことができる。PythonのリストやNumPy配列numpy.ndarrayの密行列（非スパース行列）を疎行列のクラスに変換することも可能。
- 疎行列（スパース行列）と密行列（非スパース行列）
- SciPy（
scipy.sparse）の疎行列の種類- CSR:
scipy.sparse.csr_matrix - CSC:
scipy.sparse.csc_matrix - COO:
scipy.sparse.coo_matrix - LIL:
scipy.sparse.lil_matrix
- CSR:
- SciPyで疎行列クラスを生成・変換
- リスト、
numpy.ndarrayと相互変換 - 値、行、列のリストから生成
- 空の疎行列を生成、要素を追加・削除
- 他の疎行列から生成・変換
- リスト、
- 疎行列生成の処理時間比較
scipy.sparseとnumpy.ndarrayの処理速度・メモリ使用量の比較
以下の説明内容およびサンプルコードはSciPy1.3.1のもの。バージョンが違うと異なる場合があるので注意。
疎行列（スパース行列）と密行列（非スパース行列）
成分のほとんどが0である行列を疎行列（スパース行列）といい、疎行列ではない行列（成分のほとんどが非0）を密行列（非スパース行列）という。全体の何割が0だと疎行列だとかいう決まりはない（たぶん）。
Pythonにおける行列（二次元配列）はリストのリストやNumPy配列numpy.ndarrayで表す。成分のほとんどが0であっても要素数分の値をデータとして持つ。
一方、SciPyのscipy.sparseでは0以外の要素（非ゼロ要素）のみをデータとして持つ。データの格納形式によっていくつかのクラスがある（後述）。
成分のほとんどが0である疎行列はリストやnumpy.ndarrayで表すよりも、scipy.sparseのクラスで表したほうがメモリ使用量が少なく処理速度も高速になる。
あくまでもデータの持ち方の違いであり、scipy.sparseのクラスで成分のほとんどが非ゼロの行列（密行列）を表すことも可能。ただし、密行列の場合はnumpy.ndarrayよりもメモリを消費する。
処理速度やメモリ使用量の実験結果は最後に紹介する。
便宜上、以下の説明では、成分の0の多さとは関係なく、scipy.sparseのクラスで表された行列のことを疎行列と呼び、リストやnumpy.ndarrayで表した行列を密行列と呼ぶ場合がある。疎行列形式で表された行列と密行列形式で表された行列ぐらいのイメージ。
SciPy（scipy.sparse）の疎行列クラスの種類
scipy.sparseにはcsr_matrixやlil_matrixのような様々なクラスがあり、それぞれ異なる方式でデータを格納する。
説明のため以下のようにライブラリをインポートし、元のリストを作成しておく。
from scipy.sparse import csr_matrix, csc_matrix, coo_matrix, lil_matrix
l = [[0, 10, 20],
[30, 0, 40],
[0, 0, 50]]
各要素の値、行、列のリストで表すと以下のようになる。例えば値10が(0, 1)に値40が(1, 2)にあることを意味している。
data = [10, 20, 30, 40, 50]
row = [0, 0, 1, 1, 2]
col = [1, 2, 0, 2, 2]
scipy.sparseではこのような値、および、行・列のインデックスのリストをそのままデータとして保持しているわけではなく、クラスによって様々な方式でデータを格納している。
なお、ここで紹介するものがすべてではない。そのほかのクラスは公式ドキュメントを参照。Wikipediaにも簡単な紹介がある。
- Sparse matrices (scipy.sparse) — SciPy v1.3.0 Reference Guide
- 疎行列 - Wikipedia
- Sparse matrix - Wikipedia
各形式は相互に変換可能。ざっくりとした使い分けの指針は以下の通り。詳細な長所・短所については公式ドキュメントのそれぞれのページを参照。
- データから疎行列を生成 → COO
- 算術演算や行列積計算 → CSR
- 要素の追加・更新 → LIL
- 行の取得 → CSRかLIL
- 列の取得 → CSC
- 要素・部分行列の取得 → LIL
疎行列の生成や変換については本記事の後半、値の取得や演算などについては以下の記事を参照。
CSR: scipy.sparse.csr_matrix
CSRはCompressed Sparse Rowの略。圧縮行格納方式。最もスタンダードな方式。
属性data, indices, indptrにデータが格納されている。
csr = csr_matrix(l)
print(csr)
# (0, 1) 10
# (0, 2) 20
# (1, 0) 30
# (1, 2) 40
# (2, 2) 50
print(type(csr))
# <class 'scipy.sparse.csr.csr_matrix'>
print(csr.data)
# [10 20 30 40 50]
print(csr.indices)
# [1 2 0 2 2]
print(csr.indptr)
# [0 2 4 5]
data, indicesは値、列のインデックスのリスト。indptrは行のインデックスを圧縮したリストとなる。
i行のデータがdata, indicesのindptr[i]:indptr[i+1]に格納されていることを意味している。例の場合、0行目のデータが0:2、1行目のデータが2:4、2行目のデータが4:5。
indptrのサイズ（要素数）は行数 + 1となる。
若干ややこしいが、とりあえずは行のリストが圧縮されていることを覚えておけばよい。
CSR同士だと算術演算や行列積を高速に行える。疎行列の演算を行う場合はCSRにしておけばとりあえず間違いない。データ構造上、行の取得が高速。
CSC: scipy.sparse.csc_matrix
CSCはCompressed Sparse Columnの略。圧縮列格納方式。
属性data, indices, indptrにデータが格納されている。
csc = csc_matrix(l)
print(csc)
# (1, 0) 30
# (0, 1) 10
# (0, 2) 20
# (1, 2) 40
# (2, 2) 50
print(type(csc))
# <class 'scipy.sparse.csc.csc_matrix'>
print(csc.data)
# [30 10 20 40 50]
print(csc.indices)
# [1 0 0 1 2]
print(csc.indptr)
# [0 1 2 5]
data, indicesは値、行のインデックスのリスト。indptrは列のインデックスを圧縮したリストとなる。
CSCはCSRの列バージョン。CSRは行単位でデータが左から右、上から下の順番で格納されていたが、CSCは列単位なので上から下、左から右の順番となる。indptrの圧縮の考え方はCSRと同じ。
indptrのサイズ（要素数）は列数 + 1となる。
CSRと同様に、CSC同士だと算術演算や行列積を高速に行える。ただし行列積はCSR同士のほうが速い（らしい）。データ構造上、列の取得が高速。基本はCSRだが、列を頻繁に取得したい場合はCSC。
COO: scipy.sparse.coo_matrix
COOはCoordinateの略。
属性data, row, colにデータが格納されている。
coo = coo_matrix(l)
print(coo)
# (0, 1) 10
# (0, 2) 20
# (1, 0) 30
# (1, 2) 40
# (2, 2) 50
print(type(coo))
# <class 'scipy.sparse.coo.coo_matrix'>
print(coo.data)
# [10 20 30 40 50]
print(coo.row)
# [0 0 1 1 2]
print(coo.col)
# [1 2 0 2 2]
「Coordinate（座標）」という名前の通り、値、および、行・列のインデックスのリストがそのまま格納されている。
オブジェクトの生成は高速だが算術演算などはサポートされていない。CSRやCSCへの変換も高速なので、COOを生成後、CSRに変換し演算を行うのが基本的な使い方。
LIL: scipy.sparse.lil_matrix
LILはLinked Listの略。
属性data, rowsにデータが格納されている。
lil = lil_matrix(l)
print(lil)
# (0, 1) 10
# (0, 2) 20
# (1, 0) 30
# (1, 2) 40
# (2, 2) 50
print(type(lil))
# <class 'scipy.sparse.lil.lil_matrix'>
print(lil.data)
# [list([10, 20]) list([30, 40]) list([50])]
print(lil.rows)
# [list([1, 2]) list([0, 2]) list([2])]
行ごとに値と列のインデックスのリストがそれぞれ格納されている。
行・列を指定して要素を追加したり、要素の値を更新したりするのを高速に行うことができる。演算は遅い。行ごとに格納されているので行の取得は速いが列の取得は遅い。
SciPyで疎行列を生成・変換
scipy.sparseのcsr_matrixやlil_matrixなどのクラスのオブジェクトを生成および変換する方法について説明する。
リスト、numpy.ndarrayと相互変換
上のサンプルコード中ですでに行っているように、scipy.sparseの各クラスのコンストラクタの第一引数に2次元リストを指定すると、それを元にオブジェクトが生成される。
例はcsr_matrixだが、その他のクラスでも同様。
l = [[0, 10, 20],
[30, 0, 40],
[0, 0, 0]]
csr = csr_matrix(l)
print(csr)
# (0, 1) 10
# (0, 2) 20
# (1, 0) 30
# (1, 2) 40
print(type(csr))
# <class 'scipy.sparse.csr.csr_matrix'>
print(csr.shape)
# (3, 3)
numpy.ndarrayでもOK。
a = np.array(l)
print(a)
# [[ 0 10 20]
# [30 0 40]
# [ 0 0 0]]
print(type(a))
# <class 'numpy.ndarray'>
csr = csr_matrix(a)
print(csr)
# (0, 1) 10
# (0, 2) 20
# (1, 0) 30
# (1, 2) 40
print(type(csr))
# <class 'scipy.sparse.csr.csr_matrix'>
1次元のリストやnumpy.ndarrayは1行の行列とみなされる。3次元以上だとエラー。
print(csr_matrix([0, 1, 2]))
# (0, 1) 1
# (0, 2) 2
print(csr_matrix([0, 1, 2]).shape)
# (1, 3)
# print(csr_matrix([[[0, 1, 2]]]))
# TypeError: expected dimension <= 2 array or matrix
scipy.sparseのクラスからnumpy.ndarrayに変換するにはtoarray()メソッドを使う。これも例はcsr_matrixだが、その他のクラスでも同様。
print(csr.toarray())
# [[ 0 10 20]
# [30 0 40]
# [ 0 0 0]]
print(type(csr.toarray()))
# <class 'numpy.ndarray'>
リストに変換したい場合はさらにnumpy.ndarrayのtolist()を使う。
print(csr.toarray().tolist())
# [[0, 10, 20], [30, 0, 40], [0, 0, 0]]
print(type(csr.toarray().tolist()))
# <class 'list'>
scipy.sparseのクラスにはtodense()というメソッドもある。numpy.ndarrayではなくnumpy.matrixに変換される。
print(csr.todense())
# [[ 0 10 20]
# [30 0 40]
# [ 0 0 0]]
print(type(csr.todense()))
# <class 'numpy.matrix'>
numpy.matrixは行列（2次元配列）に特化したクラス。
値、行、列のリストから生成
値、行、列のリストがそれぞれdata, row, columnとすると、コンストラクタの第一引数に(data, (row, column))と指定することでオブジェクトを生成できる。全体を囲む括弧()も必要なので注意。
この方法が使えるのはcsr_matrix, csc_matrix, coo_matrixのみ。
引数shapeで形状を指定できる。shapeを省略した場合、データが存在する範囲のみの行列が生成される。
data = [10, 20, 30, 40]
row = [0, 0, 1, 1]
col = [1, 2, 0, 2]
print(csr_matrix((data, (row, col))).toarray())
# [[ 0 10 20]
# [30 0 40]]
print(csr_matrix((data, (row, col)), shape=(3, 3)).toarray())
# [[ 0 10 20]
# [30 0 40]
# [ 0 0 0]]
print(csr_matrix((data, (row, col)), shape=(4, 4)).toarray())
# [[ 0 10 20 0]
# [30 0 40 0]
# [ 0 0 0 0]
# [ 0 0 0 0]]
rowとcolumnの組み合わせが重複している場合、csr_matrix, csc_matrixとcoo_matrixで挙動が異なる。coo_matrixは別個の要素としてデータを保持するが、csr_matrix, csc_matrixは重複要素の値が合計される。
data = [10, 20, 30, 40]
row = [0, 0, 1, 1]
col = [1, 2, 2, 2]
print(csr_matrix((data, (row, col))))
# (0, 1) 10
# (0, 2) 20
# (1, 2) 70
print(csc_matrix((data, (row, col))))
# (0, 1) 10
# (0, 2) 20
# (1, 2) 70
print(coo_matrix((data, (row, col))))
# (0, 1) 10
# (0, 2) 20
# (1, 2) 30
# (1, 2) 40
coo_matrixからtoarray()や次に説明するtocsr()などで変換すると重複要素の値が合計される。
print(coo_matrix((data, (row, col))).tocsr())
# (0, 1) 10
# (0, 2) 20
# (1, 2) 70
print(coo_matrix((data, (row, col))).toarray())
# [[ 0 10 20]
# [ 0 0 70]]
csr_matrix, csc_matrixは(data, indices, indptr)と指定してもOK。indicesは圧縮されていない方のインデックスでindptrは圧縮された方のインデックス。
print(csr_matrix(([10, 20, 30, 40], [1, 2, 0, 2], [0, 2, 4, 4]), shape=(3, 3)).toarray())
# [[ 0 10 20]
# [30 0 40]
# [ 0 0 0]]
print(csc_matrix(([30, 10, 20, 40], [1, 0, 0, 1], [0, 1, 2, 4]), shape=(3, 3)).toarray())
# [[ 0 10 20]
# [30 0 40]
# [ 0 0 0]]
data, indices, indptrは元のdata, row, columnから以下の関数で算出できる。もっと効率的に書けるかもしれない。dataや圧縮されていない方のインデックスも並び替えられる場合があるので注意。
def calc_indptr(n, data, row, column, tocsr=True):
if tocsr:
Ai, Aj = row, column
else:
Ai, Aj = column, row
indptr = np.zeros(n + 1, dtype=int)
indices = np.zeros_like(Aj)
data_ = np.zeros_like(data)
for a in Ai:
indptr[a + 1] += 1
indptr = indptr.cumsum()
for i, j, d in zip(Ai, Aj, data):
dest = indptr[i]
indices[dest] = j
data_[dest] = d
indptr[i] += 1
indptr[-1] = 0
return data_, indices, np.roll(indptr, 1)
print(calc_indptr(3, [10, 20, 30, 40], [0, 0, 1, 1], [1, 2, 0, 2]))
# (array([10, 20, 30, 40]), array([1, 2, 0, 2]), array([0, 2, 4, 4]))
print(calc_indptr(3, [10, 20, 30, 40], [0, 0, 1, 1], [1, 2, 0, 2], False))
# (array([30, 10, 20, 40]), array([1, 0, 0, 1]), array([0, 1, 2, 4]))
以下のソースコードを参考にした。
空の疎行列を生成、要素を追加・削除
コンストラクタの第一引数に(行数, 列数)のタプルで形状を指定すると空の疎行列が生成できる。引数dtypeでデータ型を指定可能。
lil = lil_matrix((3, 3), dtype=int)
print(lil.toarray())
# [[0 0 0]
# [0 0 0]
# [0 0 0]]
以下のように[行, 列]で位置（インデックス）を指定して要素に値を代入できる。位置は0始まり。
lil[1, 0] = 10
lil[2, 2] = 30
print(lil)
# (1, 0) 10
# (2, 2) 30
print(lil.toarray())
# [[ 0 0 0]
# [10 0 0]
# [ 0 0 30]]
lil_matrixの場合、0を代入すると要素が削除される。
lil[2, 2] = 0
print(lil)
# (1, 0) 10
print(lil.toarray())
# [[ 0 0 0]
# [10 0 0]
# [ 0 0 0]]
種類の説明の冒頭で書いたように、要素の追加はlil_matrixで行うのが効率的。csr_matrixやcsc_matrixでも可能だが、SparseEfficiencyWarning（非効率だという警告）が出る。
csr = csr_matrix((3, 3), dtype=int)
print(csr.toarray())
# [[0 0 0]
# [0 0 0]
# [0 0 0]]
# csr[1, 0] = 10
# SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
coo_matrixではこのような処理はサポートされていない。エラーが発生する。
coo = coo_matrix((3, 3), dtype=int)
# coo[1, 0] = 10
# TypeError: 'coo_matrix' object does not support item assignment
csr_matrixやcsc_matrix、coo_matrixの要素の値を更新したい場合は次に説明する方法でlil_matrixに変換してから行う。
疎行列の値の取得や更新についての詳細は以下の記事を参照。
他の疎行列から生成・変換
tocsr()やtocoo()、tolil()といったメソッドによって、異なる種類の疎行列クラスに変換できる。
l = [[0, 10, 20],
[30, 0, 0],
[0, 0, 0]]
csr = csr_matrix(l)
print(csr)
# (0, 1) 10
# (0, 2) 20
# (1, 0) 30
print(type(csr))
# <class 'scipy.sparse.csr.csr_matrix'>
lil = csr.tolil()
print(lil)
# (0, 1) 10
# (0, 2) 20
# (1, 0) 30
print(type(lil))
# <class 'scipy.sparse.lil.lil_matrix'>
各コンストラクタに疎行列のオブジェクトを指定しても同様に変換可能。
lil = lil_matrix(csr)
print(lil)
# (0, 1) 10
# (0, 2) 20
# (1, 0) 30
print(type(lil))
# <class 'scipy.sparse.lil.lil_matrix'>
公式ドキュメントを見ると、いずれの場合も引数copy=False（デフォルト）だとデータが共有されると書いてあるが、そのような挙動になっていない。片方の要素の値を更新してももう一方はそのまま。
lil[0, 0] = 100
print(lil.toarray())
# [[100 10 20]
# [ 30 0 0]
# [ 0 0 0]]
print(csr.toarray())
# [[ 0 10 20]
# [30 0 0]
# [ 0 0 0]]
以下のStack Overflowの回答によると、同じ種類の場合のみcopyが効くとのこと。
確かにそうなっている。copy=False（デフォルト）だと片方の要素の値を更新するともう一方の値も更新される。copy=Trueとすると別々に管理される。
lil2 = lil_matrix(lil)
print(lil2.toarray())
# [[100 10 20]
# [ 30 0 0]
# [ 0 0 0]]
lil[0, 0] = 0
print(lil2.toarray())
# [[ 0 10 20]
# [30 0 0]
# [ 0 0 0]]
lil2_copy = lil_matrix(lil, copy=True)
print(lil2_copy.toarray())
# [[ 0 10 20]
# [30 0 0]
# [ 0 0 0]]
lil[0, 0] = 100
print(lil2_copy.toarray())
# [[ 0 10 20]
# [30 0 0]
# [ 0 0 0]]
疎行列生成の処理時間比較
scipy.sparse行列の生成にかかる時間を比較する。
1000行1000列のランダムな位置に値1の要素が1000個あるとする。なお、実際には同じ位置の組み合わせが生成される可能性もあるので要素数は1000以下になる場合もある。
n = 1000
np.random.seed(0)
d = np.ones(n, dtype=int)
i = np.random.randint(0, n, n)
j = np.random.randint(0, n, n)
print(d[:10])
print(i[:10])
print(j[:10])
# [1 1 1 1 1 1 1 1 1 1]
# [684 559 629 192 835 763 707 359 9 723]
# [ 20 683 630 128 484 365 105 706 225 652]
以下の例はJupyter Notebookのマジックコマンド%%timeitを利用しており、Pythonスクリプトとして実行しても計測されない。
また、行列のサイズや非ゼロ要素の割合によって結果は異なる。あくまでも一例。
値、行、列から生成
値、行・列のインデックスのリストがある場合は、(data, (row, column))で指定するのが最速。特にcoo_matrixが速い。
csr_matrixとcsc_matrixはcoo_matrixを生成してから変換したほうが速い。
%%timeit
csr = csr_matrix((d, (i, j)), (n, n))
# 261 µs ± 8.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
csc = csc_matrix((d, (i, j)), (n, n))
# 256 µs ± 4.81 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
coo = coo_matrix((d, (i, j)), (n, n))
# 51.7 µs ± 1.37 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%%timeit
csr = coo_matrix((d, (i, j)), (n, n)).tocsr()
# 182 µs ± 1.98 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%%timeit
csc = coo_matrix((d, (i, j)), (n, n)).tocsc()
# 182 µs ± 1.29 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
空のlil_matrixを生成して要素を一つずつ追加していく方法ははるかに遅い。
%%timeit
lil = lil_matrix((n, n))
for d_, i_, j_ in zip(d, i, j):
lil[i_, j_] = d_
# 4.87 ms ± 53.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
空行列の生成やforループ自体にもある程度の時間がかかるので、非ゼロ要素が少なくてもcoo_matrixなどより速くなることは無いと思われる。
%%timeit
lil = lil_matrix((n, n))
# 409 µs ± 58.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
for d_, i_, j_ in zip(d, i, j):
pass
# 314 µs ± 10.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
一つの要素を表す(値, 行, 列)の組み合わせのリストがある場合も、lil_matrixに一つずつ追加していくよりも、転置して値と行、列に分けてcoo_matrixを生成するほうが高速。
dij = list(zip(d, i, j))
print(dij[:5])
# [(1, 684, 20), (1, 559, 683), (1, 629, 630), (1, 192, 128), (1, 835, 484)]
%%timeit
lil = lil_matrix((n, n))
for d, i, j in dij:
lil[i, j] = d
# 4.78 ms ± 127 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit
d_, i_, j_ = zip(*dij)
coo = coo_matrix((d_, (i_, j_)), (n, n))
# 522 µs ± 26.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
AtCoderなどの競技プログラミングの問題のような(値, 行, 列)の組み合わせが外から一つずつ入力されるような場合でも、一旦リストに格納してcoo_matrixを生成するほうが高速。この例では、まとめてリストに格納してから分割したほうが、それぞれ別々のリストに格納していくよりも速かった。
%%timeit
dij_ = []
for t in dij:
dij_.append(t)
d_, i_, j_ = zip(*dij_)
coo = coo_matrix((d_, (i_, j_)), (n, n))
# 548 µs ± 6.26 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
d_ = []
i_ = []
j_ = []
for d, i, j in dij:
d_.append(d)
i_.append(i)
j_.append(j)
coo = coo_matrix((d_, (i_, j_)), (n, n))
# 636 µs ± 6.56 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
numpy.ndarrayから変換
numpy.ndarrayから変換する場合。
a = coo_matrix((d, (i, j)), (n, n)).toarray()
print(a.shape)
# (1000, 1000)
%%timeit
csr = csr_matrix(a)
# 6.11 ms ± 351 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit
csc = csc_matrix(a)
# 6.16 ms ± 239 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit
coo = coo_matrix(a)
# 6.02 ms ± 433 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit
lil = lil_matrix(a)
# 7.98 ms ± 200 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit
lil = coo_matrix(a).tolil()
# 7.96 ms ± 478 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
lil_matrixは遅いが、他のクラスで生成してから変換しても処理時間は同じ。最初から所望の形式に変換すればよい。
scipy.sparseとnumpy.ndarrayの処理速度・メモリ使用量の比較
最後に、scipy.sparse.csr_matrixとnumpy.ndarrayのメモリ使用量・処理速度を比較する。
処理速度の比較として行列積の処理時間を計測している。処理内容によって結果が変わる可能性もあるのであくまでも目安。
非ゼロ要素の割合にもよるが、100行100列を超えるサイズの疎行列はscipy.sparseを使うメリットがある。
成分のほとんどが0である行列（疎行列）
1000行1000列の単位行列（対角成分が1、他は0）を例とする。
a = np.eye(1000, dtype=np.int64)
print(type(a))
# <class 'numpy.ndarray'>
print(a[:10, :10])
# [[1 0 0 0 0 0 0 0 0 0]
# [0 1 0 0 0 0 0 0 0 0]
# [0 0 1 0 0 0 0 0 0 0]
# [0 0 0 1 0 0 0 0 0 0]
# [0 0 0 0 1 0 0 0 0 0]
# [0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 0 1 0 0 0]
# [0 0 0 0 0 0 0 1 0 0]
# [0 0 0 0 0 0 0 0 1 0]
# [0 0 0 0 0 0 0 0 0 1]]
print(a.shape)
# (1000, 1000)
csr = csr_matrix(a)
print(type(csr))
# <class 'scipy.sparse.csr.csr_matrix'>
csr_matrixのメモリ使用量を取得する関数を以下のように定義する。
def get_size_of_csr(csr):
return csr.data.nbytes + csr.indices.nbytes + csr.indptr.nbytes
当然ながら非ゼロ要素のみをデータとして持つcsr_matrixのほうがメモリ使用量がはるかに少ない。
print(a.nbytes)
# 8000000
print(get_size_of_csr(csr))
# 16004
行列積にかかる処理時間を比較する。scipy.sparseは*で行列積の演算となる。
%%timeit
a @ a
# 2.19 s ± 90.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%%timeit
csr * csr
# 135 µs ± 5.34 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
処理速度もcsr_matrixのほうが速いという結果。
成分のほとんどが0である行列（いわゆる疎行列）はscipy.sparseのクラスのほうがメモリ使用量も処理速度も効率的。
成分のほとんどが0でない行列（密行列）
1000行1000列で要素が1以上100未満の乱数である行列を例とする。すべての要素が非ゼロ。
a_dense = np.random.randint(1, 100, (1000, 1000))
csr_dense = csr_matrix(a_dense)
print(a_dense[:10, :10])
# [[98 46 10 41 53 34 82 91 77 82]
# [30 21 18 89 30 20 60 88 54 5]
# [44 91 92 34 21 75 55 21 60 44]
# [90 30 57 52 14 5 64 62 41 7]
# [59 37 79 21 37 77 45 74 43 61]
# [10 73 92 98 39 72 72 4 82 81]
# [26 48 64 93 81 78 3 6 52 52]
# [28 85 42 19 51 61 23 81 89 79]
# [32 78 24 84 28 67 80 38 85 64]
# [98 11 48 85 17 5 87 30 90 22]]
print(a_dense.dtype)
# int64
print(a_dense.shape)
# (1000, 1000)
一つの要素に対するデータ量はcsr_matrixのほうが多いので、非ゼロ要素が多くなるとメモリ使用量はnumpy.ndarrayより多くなる。
print(a_dense.nbytes)
# 8000000
print(get_size_of_csr(csr_dense))
# 12004004
行列積にかかる処理時間はcsr_matrixのほうが若干速い。
%%timeit
a_dense @ a_dense
# 2.19 s ± 244 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%%timeit
csr_dense * csr_dense
# 1.91 s ± 123 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
密行列の場合、メモリ使用量と処理時間のどちらを重視するかによってscipy.sparseかnumpy.ndarrayのどちらを使うべきかが異なる。
scipy.sparseのクラスを使っても処理時間としてはデメリットがあるわけでもないので、対象の行列が疎行列か密行列か分からない、かつ、メモリの制限がない（厳しくない）ような場合はscipy.sparseを使う選択肢もあり得るだろう。
サイズが小さい行列
10行10列、100行100列の単位行列を例とする。
a_10_10 = np.eye(10, dtype=np.int64)
csr_10_10 = csr_matrix(a_10_10)
print(a_10_10.nbytes)
# 800
print(get_size_of_csr(csr_10_10))
# 164
%%timeit
a_10_10 @ a_10_10
# 2.23 µs ± 421 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%%timeit
csr_10_10 * csr_10_10
# 112 µs ± 455 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
a_100_100 = np.eye(100, dtype=np.int64)
csr_100_100 = csr_matrix(a_100_100)
print(a_100_100.nbytes)
# 80000
print(get_size_of_csr(csr_100_100))
# 1604
%%timeit
a_100_100 @ a_100_100
# 49.3 µs ± 4.32 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%%timeit
csr_100_100 * csr_100_100
# 115 µs ± 790 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
scipy.sparseのクラスのほうがメモリ使用量は少ないが、処理速度は遅い。ただしどちらも絶対値としては小さいので、極限までメモリ使用量や処理時間を切り詰めたいというような場合でもなければどちらを使ってもよいかもしれない。
200行200列の単位行列の例。
a_200_200 = np.eye(200, dtype=np.int64)
csr_200_200 = csr_matrix(a_200_200)
print(a_200_200.nbytes)
# 320000
print(get_size_of_csr(csr_200_200))
# 3204
%%timeit
a_200_200 @ a_200_200
# 6.08 ms ± 126 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit
csr_200_200 * csr_200_200
# 132 µs ± 14.9 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
numpy.ndarrayの処理時間が一気に増大している（単位に注意）。
ざっくりとした目安だが、100行100列を超えるあたりのサイズの疎行列はscipy.sparseを使うメリットがあるといえる。