圧縮行格納方式(CRS)の疎行列のデータ構造

今回も疎行列のお話です。
前回の記事で登場したcrsとcscについて、具体的にどのようなデータ構造なのかを紹介します。
ちなみにcrsとcscはそれぞれ、
圧縮行格納方式 (Compressed Sparse Row) と、
圧縮列格納方式 (Compressed Sparse Column) の略です。
ほぼ同じ処理を行方向に行うか列方向の違いしかなく、転置を取るとそれぞれ入れ替わるので、 CSRの方を紹介します。

ちなみに、 wikipediaの説明で理解したので、それをみながら記事を書いています。
例として取り上げる行列はこれ。(wikipediaの例と同じ。)
$$
\left(
\begin{matrix}
1 & 2 & 3 & 0 \\
0 & 0 & 0 & 1 \\
2 & 0 & 0 & 2 \\
0 & 0 & 0 & 1 \\
\end{matrix}
\right)
$$

まず、csr形式のデータで作りましょう。
今回はarrayで作って変換するのではなく、お作法にしたがい、lil形式で生成してから変換します。


from scipy import sparse
# 成分が全て0の 4 * 4 の lil形式の疎行列を作成。
M_lil = sparse.lil_matrix((4, 4))
# 各成分を代入。
M_lil[0, 0] = 1
M_lil[0, 1] = 2
M_lil[0, 2] = 3
M_lil[1, 3] = 1
M_lil[2, 0] = 2
M_lil[2, 3] = 2
M_lil[3, 3] = 1
# M_csr形式に変換
M_csr = sparse.csr_matrix(M_lil)
# 確認
print(M_csr)
# 出力
'''
  (0, 0)    1.0
  (0, 1)    2.0
  (0, 2)    3.0
  (1, 3)    1.0
  (2, 0)    2.0
  (2, 3)    2.0
  (3, 3)    1.0
'''

これで、csr形式の変数、M_csrに例の行列が格納されました。
printすると整形されて表示されるのですが、実際のデータ構造はこうはなっていません。
wikipediaの説明と、ドキュメントをみながら確認します。
まず、 実際のデータは、次の3つの属性に格納されています。

data ・・・ CSR format data array of the matrix
indices ・・・ CSR format index array of the matrix
indptr ・・・ CSR format index pointer array of the matrix

具体例を見てから説明します。


print(M_csr.data)
# [1. 2. 3. 1. 2. 2. 1.]
print(M_csr.indices)
# [0 1 2 3 0 3 3]
print(M_csr.indptr)
# [0 3 4 6 7]

まず、data が 疎行列の0では無い要素の値を、左上から行方向(右側へ)に順番に並べたものです。
(csrのrが対応。 cscの場合はここで列方向(下向き)に並べたものになります。)
そして、indices が、それぞれの要素が、何列目の要素なのかを示す配列です。

明らかにわかるように、data と indices の要素の数はその疎行列の0では無い成分の個数です。
あとは、dataの各要素が何行目なのかがわかれば、行列を復元できますが、
それを担っているのが、indptr です。
これだけ、wikipediaの説明と異なっていて非常にわかりにくいですが、次のように解釈できます。


# 行列の最初の行のデータは、indptrの最初の2個のデータで作ったスライスの値
print(M_csr.data[0: 3])
# [1. 2. 3.]
# 次の行のデータは、indptrの一つずらした2個のデータで作ったスライスの値
print(M_csr.data[3: 4])
# [1.]
# 以下繰り返し
print(M_csr.data[4: 6])
# [2. 2.]
print(M_csr.data[6: 7])
# [1.]

明らかにわかる通り、 indptr の要素の個数は行の数より1つ大きくなります。

これで、csr_matrixの中のデータの構造がわかりました。
また、data属性の中に行単位でデータが固まって存在してて、
行単位の取り出しや演算が得意なことにも納得できると思います。

csc_matrixとcsr_matrix

前回の記事に続いて疎行列の話です。
参考:scipyにおける疎行列

今回は疎行列を実装する型(クラス)のうち、csc_matrixcsr_matrixの紹介です。
前回紹介した、lil_matrixには効率的に疎行列を構成できるメリットがありますが、
計算が遅いというデメリットがあります。
ドキュメントに次のように書いてある通り、演算するときはcsrかcsc形式に変換した方が良いです。
arithmetic operations LIL + LIL are slow (consider CSR or CSC)

とりあえず本当に遅いのか時間を測っておきましょう。
今回は 10000*10000の大きめの疎行列を作り、実験します。
(arrayの場合の時間も測定したいので、最初にarrayで作成して変換していますが、本来は最初からlil型で生成してメモリを節約の恩恵を受けた方が良いです。)
また、csrと、cscはそれぞれ行方向、列方向の取り出しに有利という特徴があるので、実験用データによって有利不利無いように対称行列を作りました。


from scipy import sparse
import numpy as np
# 疎行列作成
M_array = np.random.choice(
    [0, 1, 2, 3, 4, 5],
    p=[0.99, 0.002, 0.002, 0.002, 0.002, 0.002],
    size=[10000, 10000]
)
# 自身の転置行列と足し合わせて対称行列にする
M_array = M_array + M_array.T

# それぞれの型に変換する
M_lil = sparse.lil_matrix(M_array)
M_csc = sparse.csc_matrix(M_lil)
M_csr = sparse.csr_matrix(M_lil)

# 和をとるのにかかる時間を計測
%timeit M_array + M_array
# 406 ms ± 2.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit M_lil + M_lil
# 618 ms ± 7.04 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit M_csc + M_csc
# 7.92 ms ± 17.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit M_csr + M_csr
# 7.94 ms ± 51.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

ご覧の通り、 lil 同士の演算は通常のarray形式よりも遅く、一方でcscやcsr同士の演算は圧倒的に速くなりました。

あとは、cscとcsrの違いですが、それぞれ、列方向のスライスと行方向のスライスが効率的に行えます。
(あまり指摘されないようですが、スライスをとるだけであればarrayが一番早い。)

これもそれぞれみておきましょう。


# 列方向のスライス
%timeit M_array[:, 400]
# 225 ns ± 3.06 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit M_lil[:, 400]
# 11.4 ms ± 201 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit M_csc[:, 400]
# 151 µs ± 2.73 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit M_csr[:, 400]
# 2.86 ms ± 127 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

列方向のスライスについては、array が最速ではありますが、 csc が csr に比べてかなり早いことが確認できます。

次は行方向。


# 行方向のスライス
%timeit M_array[400, :]
# 235 ns ± 10.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit M_lil[400, :]
# 36.6 µs ± 1.37 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit M_csc[400, :]
# 2.9 ms ± 112 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit M_csr[400, :]
# 66 µs ± 3.39 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

今度は、csr が csc より速くなりました。
とはいえ、lil が csrよりも速く処理できています。(array最速は変わらず。)

これは、前回の記事で説明した通り、lilがrows属性に、行ごとの値を保有しているからのようです。

csrの方が演算が速く、わざわざ疎行列を作ったうえでなんの演算もせずにスライスで値を取り出すだけという場面もあまり無いので、
csrかcscを選んで使うことが多いと思いますが、行方向のスライスだけならlilが速いことは一応覚えておこうと思います。

scipyにおける疎行列

自然言語処理やレコメンドシステムの開発などをやっていると、サイズは非常に大きいが成分のほとんどは0という行列が頻繁に登場します。
そのような行列を疎行列と言います。
疎行列 – Wikipedia
これをそのまま扱うと、メモリは無駄に消費するし無駄な計算は増えるしで非常に効率の悪いことになります。
そのため、scipyでは疎行列を専門に扱う scipy.sparse というモジュールが用意されています。

ドキュメント: Sparse matrices (scipy.sparse)

リンク先を見たらわかる通り、疎行列を表す型は結構色々な種類があります。(Sparse matrix classes 参照)
それぞれ、他の型との変換が高速とか、行方向/列方向の取り出しが早いとか個別にメリットを持っていますが、
共通しているのはデータ量を大幅に節約できる点です。

今後の記事でいくつか紹介する予定ですが、とりあえずデータ量を削減できるって特徴だけでも確認しておきましょう。

疎行列のメリットを感じるにはかなり小さいのですが、乱数で10*10の疎行列(array型)を作って、
それを lil_matrix に変換して、中のデータを見てみましょう。

まずデータ作成。
(今回はサンプルとして、array型でデータを作って変換していますが、
省メモリの恩恵をきちんと受けるには最初から疎行列でデータを生成するべきである点には注意してください。)


from scipy import sparse
import numpy as np
M_array = np.random.choice(
    [0, 1, 2, 3, 4, 5],
    p=[0.9, 0.02, 0.02, 0.02, 0.02, 0.02],
    size=[10, 10]
)
M_array

# 出力例
array([[0, 3, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 4, 0, 0, 0, 0],
       [0, 0, 0, 0, 5, 0, 3, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 2, 0, 0],
       [0, 5, 0, 0, 0, 3, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 2],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

これを型変換します。方法は、strとintの変換と似たような感じで、lil_matrix(D)のようにするとできます。


M_lil = sparse.lil_matrix(M_array)
M_lil
# 出力
<10x10 sparse matrix of type ''
	with 9 stored elements in LInked List format>

これで変換できました。

printすると、0ではない成分の座標と値を保有していることが確認できます。
(ただし、実際のデータ構造は結構違います。print用に整形しているようです。)


print(M_lil)
# 以下出力
  (0, 1)	3
  (1, 3)	1
  (1, 5)	4
  (2, 4)	5
  (2, 6)	3
  (3, 7)	2
  (4, 1)	5
  (4, 5)	3
  (8, 9)	2

0行1列目が3, 1行3列目が1,という風に読んでいくと、確かに0以外の成分が記録されているのがわかります。

なお、実際のデータ構造は型ごとに違うのですが、lil_matrixの場合は次ように、
dataとrowsとうい二つの属性を使ってデータを保有しています。


print(M_lil.data)
'''
[list([3]) list([1, 4]) list([5, 3]) list([2]) list([5, 3]) list([])
 list([]) list([]) list([2]) list([])]
'''

print(M_lil.rows)
'''
[list([1]) list([3, 5]) list([4, 6]) list([7]) list([1, 5]) list([])
 list([]) list([]) list([9]) list([])]
'''

data方は、各行の0ではない要素の値、rowsにそれぞれの要素が何番目に入っているのかの情報を保有しています。
合計18個の数字で、10*10=100個の要素を持つ行列を表現できています。

あと、lil_matrixの重要なメリットして、スライスで値を取り出せるし、代入もできるという点があります。


print(M_lil[0, 1])
# 3

(できて当然に思えますが、他の疎行列のデータ型によってはこれができないものがある。)
そのため、 lil_matrix で疎行列を作って、それを(その他のメリットを持つ)他の型に変換して使うというやり方をよくやります。
(Intended Usageにそう書いてあるので従っています。)

最後に、arrayに戻す方法は toarray()です。


M_lil.toarray()
# 出力
array([[0, 3, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 4, 0, 0, 0, 0],
       [0, 0, 0, 0, 5, 0, 3, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 2, 0, 0],
       [0, 5, 0, 0, 0, 3, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 2],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=int64)

単語ごとにその単語が含まれるテキストの数を数える

前回の記事では単語の出現回数を求めましたが、
今回は単語が出現するテキストの数を算出してみます。

参考:テキストデータ中の単語の出現回数を数える

とても似ていますが、要は1つのテキストに同じ単語が何回出ても1回としてしか数えないというものです。
(tf-idfのidfの計算につかわ割れるやつですね。)

コードはとても似たものになります。


from sklearn.feature_extraction.text import CountVectorizer
from sklearn.datasets import fetch_20newsgroups
import numpy as np

# データの読み込み
remove = ('headers', 'footers', 'quotes')
twenty_news = fetch_20newsgroups(
                                subset='all',
                                remove=remove,
                            )
X = twenty_news.data
# 文章の数
print(len(X))  # 18846

# BoW化
tf_vectorizer = CountVectorizer(
                                min_df=50,
                                token_pattern='(?u)\\b\\w+\\b',  # 1文字の単語も対象とする
                               )
tf = tf_vectorizer.fit_transform(X)
# モデルが集計対象とした(min_df回以上出現した)単語の数
print(len(tf_vectorizer.get_feature_names()))  # 4476
print(tf.shape)  # (18846, 4476)

# 単語ごとにその単語が1回以上出現したドキュメント数を求める。
document_frequency = np.array((tf>0).sum(axis=0))[0]
print(document_frequency.shape)  # (4476,)

# 出現回数上位100位までの単語と出現回数を表示
for i in document_frequency.argsort()[:-100:-1]:
    print(document_frequency[i], "\t", tf_vectorizer.get_feature_names()[i])

# 以下出力
'''
15749 	 the
14108 	 to
13948 	 a
13015 	 i
12991 	 and
12809 	 of
11842 	 in
11685 	 is
11029 	 it
10974 	 that
10406 	 for
8722 	 have
8665 	 this
8596 	 on
8447 	 you
8140 	 be
-- (以下略) --
'''

ポイントはここ

np.array((tf>0).sum(axis=0))[0]

(tf>0) すると、行列の各要素が正の数ならTrue,0ならFalse になります。
それをsumすることで、Trueが1として計算されてTrueの数がもとまるというカラクリです。

テキストデータ中の単語の出現回数を数える

テキストが1個だけなら、scikit-lerarnでBoW作って終わりなので、テキストデータは複数あるものとします。
(とは言っても、やることは各テキストをBoWにして足すだけです。)
結果を全部列挙したらすごい量になるのと、数えるだけだと面白くないので、
アウトプットは出現回数が多い単語のランキングにしましょう。

今回もサンプルデータは20newsgroupsを使わせていただきます。
また、この記事では出現回数数えて、上位の単語を列挙するところををゴールとし、
機械学習にかけたりしないのでstop_wordsや、細かな前処理は省きます。

コードは以下のようになりました。


from sklearn.feature_extraction.text import CountVectorizer
from sklearn.datasets import fetch_20newsgroups
import numpy as np

# データの読み込み
remove = ('headers', 'footers', 'quotes')
twenty_news = fetch_20newsgroups(
                                subset='all',
                                remove=remove,
                            )
X = twenty_news.data
# 文章の数
print(len(X))  # 18846

# BoW化(今回は出現頻度が高い単語に関心があるので、min_dfは大きめに設定。
tf_vectorizer = CountVectorizer(
                                min_df=50,
                                token_pattern='(?u)\\b\\w+\\b',  # 1文字の単語も対象とする
                               )
tf = tf_vectorizer.fit_transform(X)
# モデルが集計対象とした(min_df回以上出現した)単語の数
print(len(tf_vectorizer.get_feature_names()))  # 4476
print(tf.shape)  # (18846, 4476)

# 単語ごとに各ドキュメントの出現回数を足し合わせる
term_frequency = np.array(tf.sum(axis=0))[0]
print(term_frequency.shape)  # (4476,)

# 出現回数上位100位までの単語と出現回数を表示
for i in term_frequency.argsort()[:-100:-1]:
    print(term_frequency[i], "\t", tf_vectorizer.get_feature_names()[i])

# 以下出力
'''
173592 	 the
86907 	 to
77192 	 of
73928 	 a
70083 	 and
59683 	 i
50898 	 in
48899 	 is
45942 	 that
38660 	 it
32302 	 for
30492 	 you
24449 	 s
23617 	 on
23519 	 this
22030 	 be
-- (以下略) --
'''

不要語の除去をやっていないので、当然のように冠詞や前置詞など超汎用的な単語が並びました。

最後の出力には以前の記事で紹介したargsortを使っています。
参考:numpyのarrayを並び替えた結果をインデックスで取得する
[:-100:-1]というスライスで、最後の100項を逆順にとっているのもポイント。

単語ごとに各ドキュメントの出現回数を足し合わせるところでは、少しだけ工夫をしました。
まず、BoWが格納されている変数tfですが、 csr_matrix というデータ型になっています。

これは疎行列というメモリを節約できて、うまく使えば計算時間も短縮できる便利なものなのですが、
一見わかりにくいので toarray()でarray型に変化してしまってる例をよく見かけます。

今回のコードでは、np.array(tf.sum(axis=0))[0]という風に、
先にsumした後に、arrayに変換して、1次元に変換しました。
こうすると、toarray()してから和を取るよりもかなり早く処理できます。

jupyterなら %timeit で手軽に時間を測れるので比較しましょう。


%timeit tf.toarray().sum(axis=0)
# 355 ms ± 5.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit np.array(tf.sum(axis=0))[0]
# 1.47 ms ± 75.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

240倍も違いますね。

列の和だけなら、 csr_matrix より csc_matrix の方が計算が早いと思うのですが、
変換にかかる時間もありますし、十分早いのでこれで良いでしょう。

scikit-learnのソースコードも確認しましたが、
csr_matrix で結果を返すようになっていて、
csc_matrix を使うオプションなどはなさそうです。

pythonで、2進法/8進法/16進法で数値を定義する

通常、pythonのコード中で数値を定義したい時は
a=123
のように10進法で数値を表します。

ただ、pythonにおいては、10進法以外で数値を定義する方法も用意されています。
特に、2,8,16進法については、それぞれ数値の前に、
0b, 0o, 0x をつけることで定義できます。


>>> 0b1010101
85
>>> 0o251
169
>>> 0xe3f8
58360

2進法、8進法の方は実用的に使ったことがないのですが、
16進法はユニコード表を読むときなどに使ったことがあります。

numpyのarrayを並び替えた結果をインデックスで取得する

データを大きい順や小さい順に並び替えることはよくある操作ですが、
その結果として、”n番目の値は何だったのか?”よりも、”何番目の要素がn番目だったのか?”を知りたいケースは地味にあります。
[241,631,362,222,2,44] のようなデータがあって、 最大値は631だってことよりも、
最大値のindexは1だ(0から始まることに注意)ってことを得たい場合ですね。

そのような時、ソートした結果をインデックスのリストで得られると便利です。
それは numpyのargsortで得られます。
numpy.argsort

通常の並べ替え結果を返してくれる sort とそれぞれ実行してみましょう。

それぞれドキュメントを見るといくつか引数が用意されていて、多次元の場合の挙動や、ソートアルゴリズムなどを指定できます。
しかし、自分にとってはその中に昇順/降順の指定が”ない”ことの方が重要です。
デフォルトで昇順にソートするので、降順がいい時はどは別途指定します。

それでは、sortとargsortを昇順降順で試します。


import numpy as np
data = [241, 631, 362, 222, 2, 44]

print(list(np.sort(data)))
# [2, 44, 222, 241, 362, 631]
print(list(np.argsort(data)))
# [4, 5, 3, 0, 2, 1]

# 降順にしたい時は[::-1]を使う
print(list(np.sort(data))[::-1])
# [631, 362, 241, 222, 44, 2]
print(list(np.argsort(data))[::-1])
# [1, 2, 0, 3, 5, 4]

うまく動いていますね。

このブログでもすでに次の記事のコード中で利用しています。
参考:pythonでトピックモデル(LDA)

magic function 自体の説明を見る

jupyter notebook の便利な機能に magic functionがあります。
% や %% で始まるあれですね。
このblog記事中でも、処理時間を測ったり、トレジャーデータに接続したりと使ってます。

ただ、個々の関数の使い方の説明はわかっても、そもそも magic functionとは何者か、
という点の理解をこれまでおろそかにしてました。

それで、何を読めばいいのかちょっと調べていたのですが、
jupyter notebookで、%magicを実行すると、magic functionの説明が出てくるようです。

結構な長文が表示されたので、この記事中に貼り付けるようなことはしないのですが、
ザーッと眺めた限りでも自分が知らなかった情報が多く含まれているようです。

よく使うものについてはこのブログでも今後取り上げていきたいですが、
取り急ぎ、%magicを実行して出てきた文章を読んでいただくと色々参考になるのではないでしょうか。

matplotlibのhist()の戻り値

matplotlibでヒストグラムを書く時、次のように、hist()を使います。(リンク先は公式ドキュメント)


import matplotlib.pyplot as plt
import numpy as np
data = np.random.randn(100) * 10
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.hist(data)
plt.show()

ついこの間まで知らなかったのですが、ドキュメントによると、ax.hist(data) は戻り値を返しています。
戻されるのは、各区間の度数と、区間の区切り位置、描写に使われているmatplotlibのオブジェクトの3つです。
これらのうち、度数と区切り位置を取れるのは可視化とその他の集計を整合性を保ちながら行うのに便利そうです。

とりあえず、戻り値を受け取れるようにしてもう一回やってみましょう。
(図は今回省略します。)


import matplotlib.pyplot as plt
import numpy as np
data = np.random.randn(100) * 10
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
n, bins, patches = ax.hist(data)
plt.show()
print(n)
print(bins)
print(patches)

# 以下出力 (図は略)
[ 1.  5.  6. 18. 14. 16. 21. 13.  4.  2.]
[-29.53118806 -24.17059351 -18.80999896 -13.44940441  -8.08880985
  -2.7282153    2.63237925   7.9929738   13.35356835  18.7141629
  24.07475745]

ビンの数は今回何も指定していないので、デフォルトの10個です。
そのため、度数の配列nには10個の要素が含まれ、
区切り位置binsは右端左端が両方入っているので、11この要素を持つ配列になっています。

pythonの関数中でグローバル変数に代入する

前の記事のフィボナッチ数列の実装の中で、関数が呼び出された回数を数えるために少し仕掛けをいれました。
参考:pythonの関数をメモ化する

で、その中で globalと言うのを使っているのでその解説です。

ドキュメントはこちら。
global文

pythonでは、関数の中でグローバル変数を”参照”することができます。
しかし、そのままではグローバル変数に”代入”はできません。

やってみましょう。
まずは、”参照”ができることの確認。


a = 1


def sample1():
    # グローバル変数aの値を参照できるためaの値1を戻す。
    return a


sample1()  # 1

次に、代入を試してみます。ちょっとわかりにくいですが、グローバル変数bを用意し、
関数中でbに代入していますがこのbが実はローカルの変数で、
関数の外でもともと宣言されていたbに影響していないことがわかります。


b = 2


def sample2():
    # グローバル変数bへの代入はできないので、このbはローカル変数
    b = 3
    print(b)


sample2()  # 3
# 変数bは変更されていない
print(b)  # 2

これが、global文を使うと、グローバル変数への代入が可能になります。


c = 10


def sample3():
    global c
    c = 20
    print(c)


sample3()  # 20
# 関数中の変更が反映されている
print(c)  # 20