Pythonの疎行列ライブラリscipy.sparseの基本

By | 2014年9月26日

今回は珍しく(?)初心者フレンドリーに書こうと思う。先日のPyCon JPでの発表で、Pythonで使える疎行列のライブラリ「scipy.sparse」について反響があったので、改めてここで使い方の基本をまとめてみることとする。そもそも、便利で使い方もそれほど難しくないのに、日本語の情報がほとんどないようなので敬遠されていることもあるかもしれない。ここでは、使い方の基本を解説する。

基本: 疎行列とは?

要素のほとんどが0であるような行列を疎行列と呼ぶ。そうではない普通の行列を、区別したいときには密行列と呼ぶ。疎行列では、非ゼロ要素だけを覚えておけば良いので、メモリ、計算時間ともに大幅な節約になる。特に非ゼロ要素の割合が小さい場合は、その節約は大きくなる。

特に機械学習の分野だと、計算過程で大規模な疎行列が必要になることがよくある。

使い方

scipy.sparseには、疎行列を表すクラスがいくつかあるが、一般の行列について演算をしたいのならlil_matrix, csc_matrix, csr_matrixだけを覚えれば良い。

典型的な使い方としては
1) lil_matrixを用意する
2) 値をセットする
3) csr_matrix(またはcsc_matrix)に変換する(どちらを使うべきかの判断は後述)
4) 演算をする
という流れになる。

つまりまとめると、

  • 値をセットするときはlil_matrixを使う
  • 計算するときはcsr_matrixまたはcsc_matrixを使う

lil_matrixでも演算はできるが、計算効率を考えるとそれは避けた方がいい。

コード例は以下のようになる。

from scipy.sparse import lil_matrix, csr_matrix

# 疎行列aを用意する
a=lil_matrix((3,3))

# 非ゼロ要素を設定する
a[0,0]=1.; a[0,2]=2.

# lil_matrixをcsr_matrixに変換する
a=a.tocsr()

# 疎行列bを用意する
b=lil_matrix((3,3))

# 非ゼロ要素を設定する
b[1,1]=3.; b[2,0]=4.; b[2,2]=5.

# lil_matrixをcsr_matrixに変換する
b=b.tocsr()

# aとbの積を計算する
c=a.dot(b)

# aとbの和を計算する
d=a+b

このコードでは、
\[
A=
\begin{pmatrix}
1&0&2\cr
0&0&0\cr
0&0&0
\end{pmatrix},\quad
B=
\begin{pmatrix}
0&0&0\cr
0&3&0\cr
4&0&5
\end{pmatrix}
\]
という行列について
\[
C=AB,\quad D=A+B
\]

という計算をしている。もちろん、これはユースケースを示しているだけで、実際に疎行列が威力を発揮するのは大規模行列のときである。3×3程度の大きさだと疎行列を使うメリットはない。

lil_matrixのコンストラクタでは、引数として(m,n)のように行列のサイズだけをしていすれば、まず全要素が0の行列ができる。そうしてできた行列をaとすると、a[i,j]=xの形で非ゼロ要素だけセットする。その後、a.tocsr()またはa.tocsc()で、csr_matrixまたはcsc_matrixに変換してから、和・積などの計算をする。

CSRかCSCか?

では、計算するときにはcsr_matrixにすべきか、csc_matrixにすべきか、という問題が残るが、それは以下のルールを考慮した上、適宜よい方を選択すればよい。

  • csr_matrixは〜行目を取り出す操作、csc_matrixは〜列目を取り出す操作が高速である
  • 同じ型同士の和・積は高速である。つまりcsr_matrix同士の和・積、csc_matrix同士の和・積はともに高速である
  • csr_matrix, csc_matrixは転置行列を取ると型が移り合う。つまりcsr_matrixの転置をとるとcsc_matrixになり、csr_matrixの転置をとるとcsc_matrixになる。

以上の説明で、「高速になる」というのは、大規模データの場合はすべて「それ以外はやってはいけない」と読み替えても間違いではない。例えば、csr_matrixとcsc_matrixの和・積の計算はやってはいけない。

ではなぜこうなっているのかというと、内部構造の理解が必要になるが、それはPyCon JPでもしゃべったが次回の記事で書こうと思う。

更新履歴:
2016/10/24 最初の例で行列Aの説明が間違っているとの指摘を受けましたので、訂正しました。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です