Python: scipy.sparseで疎行列計算する

疎行列計算用のモジュール scipy.sparse について解説します。

https://docs.scipy.org/doc/scipy/reference/sparse.html

疎行列クラス

疎行列とは要素のほとんど(90%以上)が0で構成される行列です。
レコメンドやソーシャルグラフなどの分野ではしばしば現れ、総じて巨大な行列(数千×数千〜)のため、そのままではCPUキャッシュやメモリに乗り切らずに非効率な計算になりやすく、扱いが難しいデータです。

こうした疎行列を効率的に扱うために、Pythonではscipy.sparseがよく使われます。

scipy.sparseで提供されている疎行列クラスはいくつかあります。
詳細は公式ドキュメントを参照ですが、代表的な3つのクラス(lil_matrix, csr_matrix, csc_matrix)について紹介します。

https://docs.scipy.org/doc/scipy/reference/sparse.html#sparse-matrix-classes

lil_matrix

lil_matrixは、リストやnumpyのクラス(ndarrayやmatrixなど)と相互変換するために使われることが多いクラスです。行単位のリンクトリストで実装されてます。

lil_matrixのインスタンスをprintで表示すると、非ゼロ要素のインデックスと値だけしか持っていないことがわかります
numpy.ndarrayに変換する時はtoarray()、numpy.matrixに変換する時はtodense()を使います
csr_matrixに変換する時はtocsr()、csc_matrixに変換する時はtocsc()を使います

import numpy as np
from scipy import sparse


# リスト、numpy.ndarray、numpy.matrixのデータ
matrix1 = [[1, 0, 3], [0, 5, 0]]
matrix2 = np.array([[1, 0, 3], [0, 5, 0]])
matrix3 = np.matrix([[1, 0, 3], [0, 5, 0]])

# lil_matrixの生成
lil1 = sparse.lil_matrix(matrix1, dtype=np.float64) # リストから生成(dtypeで要素の型を指定できる)
lil2 = sparse.lil_matrix(matrix2) # numpy.ndarrayから生成
lil3 = sparse.lil_matrix(matrix3) # numpy.matrixから生成

# lil_matrixの表示
print(type(lil1), lil1)
# <class 'scipy.sparse.lil.lil_matrix'>
#  (0, 0)  1.0
#  (0, 2)  3.0
#  (1, 1)  5.0

# 要素にアクセス
print(lil1[0, 0], lil1[0, 1], lil1[1, 1]) # 1.0 0.0 5.0
# スライスで1〜2列目にアクセス
print(lil1[:, 0:2])
#  (0, 0)  1.0
#  (1, 1)  5.0

# lil_matrixからnumpy.ndarrayへの変換
lil1_ndarray = lil1.toarray()
print(type(lil1_ndarray), lil1_ndarray)
# <class 'numpy.ndarray'> [[ 1.  0.  3.] [ 0.  5.  0.]]

# lil_matrixからnumpy.matrixへの変換
lil1_matrix = lil1.todense()
print(type(lil1_matrix), lil1_matrix)
# <class 'numpy.matrixlib.defmatrix.matrix'> [[ 1.  0.  3.] [ 0.  5.  0.]]

# lil_matrixからcsr_matrixへの変換
lil1_csr = lil1.tocsr()
print(type(lil1_csr)) # <class 'scipy.sparse.csr.csr_matrix'>

# lil_matrixからcsc_matrixへの変換
lil1_csc = lil1.tocsc()
print(type(lil1_csc)) # <class 'scipy.sparse.csc.csc_matrix'>

csr_matrix

csr_matrixは、行ごとに圧縮されており、行単位なら高速にアクセスできます。また、csc_matrix同士の加減算やベクトル積・行列積も高速です。

# csr_matrixの生成
csr1 = sparse.csr_matrix(matrix1, dtype=np.float64) # リストから生成(dtypeで要素の型を指定できる)
csr2 = sparse.csr_matrix(matrix2) # numpy.ndarrayから生成
csr3 = sparse.csr_matrix(matrix3) # numpy.matrixから生成
csr4 = sparse.csr_matrix(lil1) # lil_matrixから生成

# csr_matrixの表示
print(type(csr1), csr1)
# <class 'scipy.sparse.csr.csr_matrix'>   
#  (0, 0)  1.0
#  (0, 2)  3.0
#  (1, 1)  5.0

# 要素にアクセス
print(csr1[0, 0], csr1[0, 1], csr1[1, 1]) # 1.0 0.0 5.0
# スライスで1〜2列目にアクセス
print(csr1[:, 0:2])
#  (0, 0)  1.0
#  (1, 1)  5.0

# csc_matrixからlil_matrixへの変換
csr1_lil = csr1.tolil()
print(type(csr1_lil), csr1_lil)
# <class 'scipy.sparse.lil.lil_matrix'> 
#  (0, 0)  1.0
#  (0, 2)  3.0
#  (1, 1)  5.0

csc_matrix

csc_matrixは、列ごとに圧縮されており、列単位なら高速にアクセスできます。こちらも、csr_matrix同士の加減算やベクトル積・行列積も高速です。

toarray()でnumpy.ndarray、todense()でnumpy.matrix、tolil()でlil_matrixへそれぞれ変換できます
lil_matrixやcsr_matrixと異なり、csc_matrixクラスは列順で表示されていることがわかります((0, 0)要素の次に(1, 1)要素が表示されている)

# csc_matrixの生成
csc1 = sparse.csc_matrix(matrix1, dtype=np.float64) # リストから生成(dtypeで要素の型を指定できる)
csc2 = sparse.csc_matrix(matrix2) # numpy.ndarrayから生成
csc3 = sparse.csc_matrix(matrix3) # numpy.matrixから生成
csc4 = sparse.csc_matrix(lil1) # lil_matrixから生成

# csc_matrixの表示
print(type(csc1), csc1)
# <class 'scipy.sparse.csc.csc_matrix'>
#  (0, 0)  1.0
#  (1, 1)  5.0
#  (0, 2)  3.0

# 要素にアクセス
print(csc1[0, 0], csc1[0, 1], csc1[1, 1]) # 1.0 0.0 5.0
# スライスで1〜2列目にアクセス
print(csc1[:, 0:2])
#  (0, 0)  1.0
#  (1, 1)  5.0

# csc_matrixからlil_matrixへの変換
csc1_lil = csc1.tolil()
print(type(csc1_lil), csc1_lil)
# <class 'scipy.sparse.lil.lil_matrix'>
#  (0, 0)  1.0
#  (0, 2)  3.0
#  (1, 1)  5.0

計算

scipy.sparseの計算ですが、内積とアダマール積(要素ごとの積)を以外は、numpy.ndarrayと概ね同じです。

lil_matrix、csr_matrix、csc_matrixを混合して計算できます
- 演算子の前(オペランド)の型で評価されます
numpyのndarrayやmatrixと混合して計算することもできます
- その場合、matrixで評価されます
- 性能的な恩恵も受けられにくくなります

加減算

演算子+と-で疎行列同士で加減算ができます。

csr1 = sparse.csr_matrix([[1, 2], [0, 4]])
csr2 = sparse.csr_matrix([[5, 0], [0, 8]])
csc1 = sparse.csc_matrix([[1, 2], [0, 4]])
ndarray1 = np.array([[1, 0], [3, 4]])

# 加算
result_add = csr1 + csr2
print(type(result_add), result_add)
# <class 'scipy.sparse.csr.csr_matrix'>
#  (0, 0)  6
#  (0, 1)  2
#  (1, 1)  12

result_add_ndarray = csr1 + ndarray1 # ndarrayとも加算可
print(type(result_add_ndarray), result_add_ndarray)
# <class 'numpy.matrixlib.defmatrix.matrix'>
# [[2 2]
#  [3 8]]

result_add_csr_csc = csr1 + csc1 # matrix_csrとmatrix_cscで加算可
print(type(result_add_csr_csc), result_add_csr_csc)
# <class 'scipy.sparse.csr.csr_matrix'>   
#  (0, 0)  2
#  (0, 1)  4
#  (1, 1)  8

# 減算
result_sub = csr1 - csr2
print(type(result_sub), result_sub)
# <class 'scipy.sparse.csr.csr_matrix'>
#  (0, 0)  -4
#  (0, 1)  2
#  (1, 1)  -4

result_sub_ndarray = csr1 - ndarray1 # ndarrayとも減算可
print(type(result_sub_ndarray), result_sub_ndarray)
# <class 'numpy.matrixlib.defmatrix.matrix'>
# [[ 0  2]
#  [-3  0]]

内積

numpy.ndarrayと異なり、scipy.sparseでは演算子*が内積の計算となります(numpyでは要素ごとの積)。また、dot()でも内積を計算できます。
numpyのクラスがオペランドの場合、要素ごとの積になることに注意してください。

# 内積
result_inner_product = csr1 * csr2
print(result_inner_product)
#  (0, 1)  16
#  (0, 0)  5
#  (1, 1)  32

# dot()でも可能
result_dot = csr1.dot(csr2)

# ndarrayがオペランドの場合、要素ごとの積
result_ndarray_csr = np.array([[1, 0], [3, 4]]) * csr1
print(type(result_ndarray_csr), result_ndarray_csr)
# <class 'numpy.ndarray'> 
# [[1 2]
#  [3 22]]

# csr_matrixがオペランドの場合、内積
result_csr_ndarray = csr1 * np.array([[1, 0], [3, 4]])
print(type(result_csr_ndarray), result_csr_ndarray)
# <class 'numpy.ndarray'> 
# [[7 8]
#  [12 16]]

アダマール積(要素ごとの積)

要素ごとの積が必要な場合は、multiply()で計算できます。

result_multiply = csr1.multiply(csr2)
print(result_multiply)
#  (0, 0)  5
#  (1, 1)  32

固有値・固有ベクトル

固有値・固有ベクトルを求める場合、scipy.sparse.linalgをインポートして、eigsで計算します。
引数kは求める固有値・固有ベクトルの数(デフォルトは6)ですが、行列の次元数-1未満である必要があります。全ての固有値・固有ベクトルを求めることができません。

import scipy.sparse.linalg


csr3 = sparse.csr_matrix([[1, 0, 3], [0, 5, 0], [7, 0, 9]], dtype=np.float64)

# 固有値・固有ベクトル
eigen_values, eigen_vectors = sparse.linalg.eigs(csr3, k=1)
print(eigen_values)
# [ 11.08276253+0.j]
print(eigen_vectors)
# [[ -2.85181797e-01+0.j]
#  [  7.54125462e-17+0.j]
#  [ -9.58473444e-01+0.j]]

特異値分解

特異値分解は、scipy.sparse.linalg.svdsで行います。
こちらも引数kで求める特異値の数(デフォルト6)で指定しますが、行列の低い方の次元数(2×3行列の場合は2)未満にする必要があります。

csr4 = sparse.csr_matrix([[1, 0, 3], [0, 5, 0]], dtype=np.float64)

# 特異値分解
u, s, v_t = sparse.linalg.svds(csr4, k=1)
print(u)
# [[  4.04373124e-17]
#  [  1.00000000e+00]]
print(s)
# [ 5.]
print(v_t)
# [[  8.08746247e-18   1.00000000e+00   2.42623874e-17]]

性能

最後にnumpyとscipy.sparseの各クラスで性能を比較してみたいと思います。

numpyとscipy.sparseの性能比較

疎行列を使った計算で一番多いのは内積かと思います。
そこで、1,000×1,000(1,000,000要素)の疎行列(非ゼロ要素が10%)を作り、それぞれで内積を計算してその時間を計測・比較します。
また、疎行列でない場合にndarrayよりも性能劣化するのかどうかを検証するために、同じサイズの密行列(非ゼロ要素が100%)も作って内積を行います。

以下の結果となりました。カッコ内は試行回数で、10回以上の場合は上位3件の平均値です。

疎行列の場合、ndarrayと比較すると、lil_matrixでもcsr_matrixでも高速
- 特にcsr_matrixは、ndarrayの約3,000倍、lil_matrixの約5倍、高速
密行列の場合でも、ndarrayと比較して性能が低下する、ということはなかった
- 疎行列になる"かも"しれないケースにおいても、scipy.sparseのクラスに変換・利用する、というのも1つの手段ということが示唆されます

クラス	疎行列	密行列
ndarray	2.70s (1)	2.66s (1)
lil_matrix	4.57ms (100)	2.49s (1)
csr_matrix	863µs (1000)	2.30s (1)

# 1000×1000の疎行列の作成
dimensions = 1000
sparse_matrix = np.zeros((dimensions, dimensions), dtype=np.int64)

non_zero_elements = 10000 # 非ゼロ要素は全体の1%
count = 0
rows = np.random.randint(0, dimensions, 20000)
columns = np.random.randint(0, dimensions, 20000)
values = np.random.randint(0, 10, 20000)
for i in range(0, 20000):
    if sparse_matrix[rows[i], columns[i]] == 0:
        sparse_matrix[rows[i], columns[i]] = values[i]
        count += 1
        if count == non_zero_elements:
            break

sparse_lil = sparse.lil_matrix(sparse_matrix)
sparse_csr = sparse.csr_matrix(sparse_matrix)

# 疎行列の内積速度の計測
%timeit np.dot(sparse_matrix, sparse_matrix)
%timeit sparse_lil * sparse_lil
%timeit sparse_csr * sparse_csr

# 1000×1000の密行列の作成
dense_matrix = np.random.randint(1, 10, (dimensions, dimensions))
dense_lil = sparse.lil_matrix(dense_matrix)
dense_csr = sparse.csr_matrix(dense_matrix)

# 密行列の実行速度の計測
%timeit np.dot(dense_matrix, dense_matrix)
%timeit dense_lil * dense_lil
%timeit dense_csr * dense_csr

csr_matrixとcsc_matrixの性能比較

最後に、csr_matrixとcsc_matrixについて、行単位と列単位でアクセスした場合の速度について計測・比較してみます。

scipy.sparseの疎行列クラスで行・列にアクセスするには、スライスを使う方法とメソッドを使う方法の2つがあります。それぞれの方法で計測します。

スライスを使う方法
- i行目にアクセスするなら、matrix[i]
- i列目にアクセスするなら、matrix[:, i]
メソッドを使う方法
- i行目にアクセスするなら、matrix.getrow(i)
- i列目にアクセスするなら、matrix.getcol(i)

結果は以下のとおりです。

総じて、スライスよりも、メソッドを使う方が高速
- 行単位・列単位でアクセスしたい場合は、極力メソッドを使う方が良さそうです
行アクセスでは、csr_matrix > lil_matrix > csc_matrixの順で、予想通りでした
列アクセス(メソッドの場合)では、csc_matrix > csc_matrix >> lil_matrixの順でした
- lil_matrixは遥かに低速で、何らか計算を行う際にはlil_matrixから別のクラスに変換した方が良さそうです
- スライスの場合は、csc_matrixよりもcsr_matrixの方が高速になっており、予想に反した結果となりました

クラス	行アクセス(スライス)	列アクセス(スライス)	行アクセス(getrow())	列アクセス(getcol())
lil_matrix	498ms (1)	7.40s (1)	469ms (1)	2.41s (1)
csr_matrix	330ms (1)	386ms (1)	304ms (1)	367ms (1)
csc_matrix	504ms (1)	442ms (1)	425ms (1)	305ms (1)

sparse_csc = sparse.csc_matrix(sparse_matrix)

# 行単位(スライス)で平均値を算出
%timeit [sparse_lil[i].mean() for i in range(dimensions)]
%timeit [sparse_csr[i].mean() for i in range(dimensions)]
%timeit [sparse_csc[i].mean() for i in range(dimensions)]

# 列単位(スライス)で平均値を算出
%timeit [sparse_lil[:, i].mean() for i in range(dimensions)]
%timeit [sparse_csr[:, i].mean() for i in range(dimensions)]
%timeit [sparse_csc[:, i].mean() for i in range(dimensions)]

# 行単位(getrow())で平均値を算出
%timeit [sparse_lil.getrow(i).mean() for i in range(dimensions)]
%timeit [sparse_csr.getrow(i).mean() for i in range(dimensions)]
%timeit [sparse_csc.getrow(i).mean() for i in range(dimensions)]

# 列単位(getcol())で平均値を算出
%timeit [sparse_lil.getcol(i).mean() for i in range(dimensions)]
%timeit [sparse_csr.getcol(i).mean() for i in range(dimensions)]
%timeit [sparse_csc.getcol(i).mean() for i in range(dimensions)]