唯物是真 @Scaled_Wurm

プログラミング(主にPython2.7)とか機械学習とか

Python で疎行列(SciPy)

疎行列とは

疎行列は成分のほとんどがゼロである行列のことです。
\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1\\ \end{bmatrix}
たとえば、文書に登場する単語の頻度を数えたりするとこういった行列になります。
他にも疎なグラフの隣接行列は疎行列になります。
こういった行列は通常の行列(密行列)を使うよりも、少ないメモリで処理できたり、高速に処理できたりすることがあります。

例えば簡単な例として上の行列を考えます。
この時すべての成分を保存するには4 \times 4 = 16個分の数値のメモリが必要になります。
しかし上の行列では非ゼロの要素は3個だけです。
この3個について、以下のように行と列の位置と値だけを覚えておけば3 \times 3 = 9となりすべての成分を記憶しておくよりも効率的なことがわかると思います。
\begin{bmatrix} (3, 1, 1), (2, 2, 1), (4, 4, 1) \end{bmatrix}

疎行列の種類

上では単純な例を示しましたが、実際の疎行列はそれぞれ特徴の異なる様々な表現手法が用いられていて、SciPyにも以下のような種類の疎行列が実装されています。

  1. csc_matrix: Compressed Sparse Column format
  2. csr_matrix: Compressed Sparse Row format
  3. bsr_matrix: Block Sparse Row format
  4. lil_matrix: List of Lists format
  5. dok_matrix: Dictionary of Keys format
  6. coo_matrix: COOrdinate format (aka IJV, triplet format)
  7. dia_matrix: DIAgonal format
http://docs.scipy.org/doc/scipy/reference/sparse.html

行列を作成するときにはlil_matrixが推奨で、他にもdok_matrixやcoo_matrixも効率的らしいです。
また乗算などの操作をするにはcsc_matrixやcsr_matrixが適しているらしいです(適しているというか疎行列の種類と操作によってはそもそも操作ができなかったりします)

それぞれの行列についての詳細はScipyのドキュメントWikipediaの記事(英語)を見たほうがわかりやすいかもしれません。

疎行列の作成

lil_matrix

以下の公式ドキュメントの例のようにインデックスを指定して、行列の値を入れられるのでわかりやすいです。

>>> from scipy.sparse import lil_matrix
>>> from scipy.sparse.linalg import spsolve
>>> from numpy.linalg import solve, norm
>>> from numpy.random import rand
>>> A = lil_matrix((1000, 1000))
>>> A[0, :100] = rand(100)
>>> A[1, 100:200] = A[0, :100]
>>> A.setdiag(rand(1000))
http://docs.scipy.org/doc/scipy/reference/sparse.html

またlil_matrixからcsr_matrixなどへの変換はtocsr関数などを使えば変換出来ます。

>>> A = A.tocsr()
http://docs.scipy.org/doc/scipy/reference/sparse.html

coo_matrix

coo_matrixを使うときには以下のように値Vとインデックスの配列I, Jから作ります。

>>> from scipy import sparse
>>> from numpy import array
>>> I = array([0,3,1,0])
>>> J = array([0,3,1,2])
>>> V = array([4,5,7,9])
>>> A = sparse.coo_matrix((V,(I,J)),shape=(4,4))
http://docs.scipy.org/doc/scipy/reference/sparse.html

疎行列の変換

疎行列から密行列に変換するにはtodense関数を使います。

dense = sparse.todense()

逆に密行列から疎行列に変換するには、基本的に以下のように密行列を疎行列の引数に与えればうまくいきます。

sparse = csr_matrix(dense)

疎行列の保存

numpy.save関数とnumpy.load関数を使えばちゃんと疎行列として扱ってくれます。

numpy.save('array', sparse)
sparse2 = numpy.load('array.npy')

なぜかnumpy.savetxt関数はうまくいかないです。
scipy.ioにある関数を使えばMATLABの形式でも保存できるらしいです。

疎行列用のアルゴリズム

以下のようにSciPyでは疎行列用のアルゴリズムが実装されていて、特異値分解や最短路の探索などが実装されています。

特異値分解(SVD)

使えるアルゴリズムの例として、SVDの応用で自然言語処理でたまに出てくるLSI(LSA)のコードを載せておきます。

3行4列の行列を、3行2列へと次元の削減をしています。

import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import svds

sparse= lil_matrix((3, 4))
sparse[0, 0] = 1
sparse[0, 3] = 1
sparse[1, 0] = 1
sparse[2, 0] = 1
sparse[2, 3] = 1
print svds(sparse, 2)[0]

ついでに紹介しておくとPython特異値分解をやるライブラリにはいろいろあって、以下のブログ記事(英語)で性能の比較がされています

まとめ

疎なデータの時は密行列ではなく、疎行列を使ってメモリを節約しよう。