Python, SciPyで疎行列の計算・処理(逆行列、固有値、連結、保存など)
SciPyのscipy.sparse
行列(疎行列)の計算(四則演算、逆行列、固有値など)や各種処理(連結や保存など)を行う方法について、以下の内容を説明する。
- 疎行列の四則演算、行列積
- 疎行列(
scipy.sparse
)同士 - NumPy配列
ndarray
との計算
- 疎行列(
- 疎行列の逆行列、固有値など:
scipy.sparse.linalg
- 逆行列:
inv()
- 固有値・固有ベクトル:
eigs()
- 逆行列:
scipy.sparse
のクラスのメソッド- 平均・合計・最大値・最小値
- 各要素を処理(平方根や三角関数など)
- 非ゼロ要素の位置:
nonzero()
scipy.sparse
の関数- 疎行列の連結:
hstack()
,vstack()
- 上三角行列・下三角行列の取得:
triu()
,tril()
- 疎行列オブジェクトの保存・読み込み:
save_npz()
,load_npz()
- 疎行列の連結:
疎行列の基本については以下の記事を参照。
以下の説明内容およびサンプルコードはSciPy1.3.1
のもの。バージョンが違うと異なる場合があるので注意。
疎行列の四則演算、行列積
疎行列(scipy.sparse)同士
csr_matrix
を例とする。分かりやすいようにtoarray()
でNumPy配列ndarray
に変換して出力する。
+
、-
は要素ごとの足し算・引き算となる。結果も同じくcsr_matrix
。
import numpy as np
from scipy.sparse import csr_matrix, csc_matrix, coo_matrix, lil_matrix
l = [[0, 1, 2],
[3, 0, 4],
[0, 0, 0]]
csr = csr_matrix(l)
print((csr + csr).toarray())
# [[0 2 4]
# [6 0 8]
# [0 0 0]]
print(type((csr + csr)))
# <class 'scipy.sparse.csr.csr_matrix'>
print((csr - csr).toarray())
# [[0 0 0]
# [0 0 0]
# [0 0 0]]
*
は行列積。要素ごとの掛け算はmultiply()
メソッドを使う。dot()
で明示的に行列積を表すことも可能。
print((csr * csr).toarray())
# [[3 0 4]
# [0 3 6]
# [0 0 0]]
print((csr.multiply(csr)).toarray())
# [[ 0 1 4]
# [ 9 0 16]
# [ 0 0 0]]
print((csr.dot(csr)).toarray())
# [[3 0 4]
# [0 3 6]
# [0 0 0]]
/
は要素ごとの割り算。0
で割られた要素はnan
となり、結果はnumpy.matrix
となる。numpy.matrix
は行列(二次元配列)に特化した型。numpy.ndarray
に変換可能。
print(csr / csr)
# [[nan 1. 1.]
# [ 1. nan 1.]
# [nan nan nan]]
print(type(csr / csr))
# <class 'numpy.matrix'>
print(np.array(csr / csr))
# [[nan 1. 1.]
# [ 1. nan 1.]
# [nan nan nan]]
print(type(np.array(csr / csr)))
# <class 'numpy.ndarray'>
すべての要素が非ゼロ要素でも変わらずnumpy.matrix
が返される。
csr_full = csr_matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(csr_full / csr_full)
# [[1. 1. 1.]
# [1. 1. 1.]
# [1. 1. 1.]]
print(type(csr_full / csr_full))
# <class 'numpy.matrix'>
0
の要素が無くなってしまい非効率であるため、スカラー値との足し算・引き算はサポートされていない。すべての要素が非ゼロ要素でも同様。
# print(csr + 10)
# NotImplementedError: adding a nonzero scalar to a sparse matrix is not supported
# print(csr - 10)
# NotImplementedError: subtracting a nonzero scalar from a sparse matrix is not supported
# print(csr_full + 10)
# NotImplementedError: adding a nonzero scalar to a sparse matrix is not supported
# print(csr_full - 10)
# NotImplementedError: subtracting a nonzero scalar from a sparse matrix is not supported
スカラー値との掛け算、割り算はOK。
print((csr * 10).toarray())
# [[ 0 10 20]
# [30 0 40]
# [ 0 0 0]]
print((csr / 10).toarray())
# [[0. 0.1 0.2]
# [0.3 0. 0.4]
# [0. 0. 0. ]]
べき乗は行列積の繰り返しとなる。正の整数のみ指定可能。
print((csr ** 2).toarray())
# [[3 0 4]
# [0 3 6]
# [0 0 0]]
print((csr ** 3).toarray())
# [[ 0 3 6]
# [ 9 0 12]
# [ 0 0 0]]
print((csr * csr * csr).toarray())
# [[ 0 3 6]
# [ 9 0 12]
# [ 0 0 0]]
# print((csr ** -1).toarray())
# alueError: exponent must be >= 0
# print((csr ** 0.5).toarray())
# ValueError: exponent must be an integer
疎行列の演算は結果のクラスに注意。csr_matrix
とcsc_matrix
との演算の結果のクラスは順番によって異なる。また、coo_matrix
, lil_matrix
は自動的にcsr_matrix
に変換される。
csc = csc_matrix(l)
coo = coo_matrix(l)
lil = lil_matrix(l)
print(type(csc + csc))
# <class 'scipy.sparse.csc.csc_matrix'>
print(type(csr + csc))
# <class 'scipy.sparse.csr.csr_matrix'>
print(type(csc + csr))
# <class 'scipy.sparse.csc.csc_matrix'>
print(type(coo + coo))
# <class 'scipy.sparse.csr.csr_matrix'>
print(type(lil + lil))
# <class 'scipy.sparse.csr.csr_matrix'>
なお、処理速度はcsr_matrix
同士、または、csc_matrix
同士が高速。異なるクラスを混在させたままにするよりもどちらかに統一して処理したほうがよい。
NumPy配列ndarrayとの計算
あまり使い所はなさそうだが、scipy.sparse
とnumpy.ndarray
との四則演算も可能。
+
, -
の結果はnumpy.matrix
。どちらが先でもOK。
a = np.array(l)
print(a)
# [[0 1 2]
# [3 0 4]
# [0 0 0]]
print(type(a))
# <class 'numpy.ndarray'>
print(a + csr)
# [[0 2 4]
# [6 0 8]
# [0 0 0]]
print(type(a + csr))
# <class 'numpy.matrix'>
print(type(csr - a))
# <class 'numpy.matrix'>
割り算はndarray / sparse
の順番だとエラーとなる。sparse / ndarray
の順番だとエラーにはならないが警告が出る。0
の要素はnan
となる。結果はnumpy.matrix
。
# print(type(a / csr))
# TypeError: unsupported operand type(s) for /: 'numpy.ndarray' and 'csr_matrix'
print(csr / a)
# [[nan 1. 1.]
# [ 1. nan 1.]
# [nan nan nan]]
#
# /usr/local/lib/python3.7/site-packages/scipy/sparse/base.py:596: RuntimeWarning: invalid value encountered in true_divide
# return np.true_divide(self.todense(), other)
print(type(csr / a))
# <class 'numpy.matrix'>
#
# /usr/local/lib/python3.7/site-packages/scipy/sparse/base.py:596: RuntimeWarning: invalid value encountered in true_divide
# return np.true_divide(self.todense(), other)
*
は行列積となり、結果はnumpy.ndarray
。
print(csr * a)
# [[3 0 4]
# [0 3 6]
# [0 0 0]]
print(type(csr * a))
# <class 'numpy.ndarray'>
print(type(a * csr))
# <class 'numpy.ndarray'>
print(type(csr.dot(a)))
# <class 'numpy.ndarray'>
公式ドキュメントにあるように、NumPy1.7時点で、numpy.ndarray
のメソッドdot()
やnumpy.dot()
の引数にscipy.sparse
行列を指定するとおかしな結果となるので注意。
Warning
As of NumPy 1.7, np.dot is not aware of sparse matrices, therefore using it will result on unexpected results or errors.
Sparse matrices (scipy.sparse) — SciPy v1.3.0 Reference Guide
print(a.dot(csr))
# [[<3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>]
# [<3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>]
# [<3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>]]
要素ごとの掛け算はscipy.sparse
行列のmultiply()
メソッドを使う。結果はcoo_matrix
。こちらもnumpy.multiply()
だとうまくいかない。
print(csr.multiply(a))
# (0, 1) 1
# (0, 2) 4
# (1, 0) 9
# (1, 2) 16
print(type(csr.multiply(a)))
# <class 'scipy.sparse.coo.coo_matrix'>
print(csr.multiply(a).toarray())
# [[ 0 1 4]
# [ 9 0 16]
# [ 0 0 0]]
print(np.multiply(a, csr))
# [[<3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>]
# [<3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>]
# [<3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>
# <3x3 sparse matrix of type '<class 'numpy.int64'>'
# with 4 stored elements in Compressed Sparse Row format>]]
疎行列の逆行列、固有値など: scipy.sparse.linalg
scipy.sparse
行列(疎行列)の逆行列や固有値などを算出するにはサブモジュールscipy.sparse.linalg
の関数を使う。linalg
は線型代数(linear algebra)の意。高度な処理が多く提供されている。
numpy.ndarray
の逆行列や固有値の算出については以下の記事を参照。
逆行列: inv()
逆行列はinv()
。
csc_matrix
以外を引数に指定すると、csc_matrix
を使わないと非効率だという警告SparseEfficiencyWarning
が出る。ここではコメントアウトしているが、エラーではなく警告なので実行することは可能。
import numpy as np
from scipy.sparse import csr_matrix, csc_matrix
from scipy.sparse.linalg import inv, eigs
l = [[2, 5],
[1, 3]]
csr = csr_matrix(l)
csc = csc_matrix(l)
# csr_inv = linalg.inv(csr)
# SparseEfficiencyWarning: splu requires CSC matrix format
# SparseEfficiencyWarning: spsolve is more efficient when sparse b is in the CSC matrix format
csc_inv = inv(csc)
print(csc_inv)
# (0, 0) 3.0
# (1, 0) -1.0
# (0, 1) -5.0
# (1, 1) 2.0
print(type(csc_inv))
# <class 'scipy.sparse.csc.csc_matrix'>
print(csc_inv.toarray())
# [[ 3. -5.]
# [-1. 2.]]
なお、便宜上、例で使っている行列は0
を含まないもの。実際はわざわざscipy.sparse
を使うメリットはない。
固有値・固有ベクトル: eigs()
先にnumpy.linalg.eig()
の結果を示す。固有値と固有ベクトルが返される。
l = [[1, 1, 2],
[0, 2, -1],
[0, 0, 3]]
w, v = np.linalg.eig(l)
print(w)
# [1. 2. 3.]
print(v)
# [[ 1. 0.70710678 0.33333333]
# [ 0. 0.70710678 -0.66666667]
# [ 0. 0. 0.66666667]]
scipy.sparse.linalg
の固有値を算出する関数はeigs()
。
eigs()
で得られるのは1組の固有値と固有ベクトルのみ。
csr_f = csr_matrix(l, dtype=float)
csr_i = csc_matrix(l, dtype=int)
w, v = eigs(csr_f, k=1)
print(w)
# [3.+0.j]
print(v)
# [[ 0.33333333+0.j]
# [-0.66666667+0.j]
# [ 0.66666667+0.j]]
取得する固有値と固有ベクトルの数は引数k
で指定できるが、k < N - 1
(Nは行数)という制約がある。公式ドキュメントにあるように、すべての固有値・固有ベクトルを算出することはできない。
The number of eigenvalues and eigenvectors desired. k must be smaller than N-1. It is not possible to compute all eigenvectors of a matrix.
scipy.sparse.linalg.eigs — SciPy v1.3.0 Reference Guide
k
に大きい値を指定するとエラーとなる。
# w, v = eigs(csr_f, k=2)
# TypeError: Cannot use scipy.linalg.eig for sparse A with k >= N - 1. Use scipy.linalg.eig(A.toarray()) or reduce k.
小さいサイズの行列のすべての固有値と固有ベクトルを求めるような用途にはnumpy.linalg.eig()
やscipy.linalg.eig()
を使うべきで、scipy.sparse.linalg.eigs()
は巨大な疎行列向け、ということだと思われる。
なお、scipy.sparse.linalg.eigs()
の第一引数に指定する行列のデータ型dtype
はfloat
である必要がある。整数int
だとエラー。
# w, v = eigs(csr_i, k=1)
# ValueError: matrix type must be 'f', 'd', 'F', or 'D'
scipy.sparseのクラスのメソッド
csr_matrix
やlil_matrix
などのscipy.sparse
のクラスのメソッドをいくつか紹介する。
共通の基底クラスspmatrix
で定義されているメソッドはどのクラスでも使える。
一方、各クラスで独自に定義されているものは別のクラスでは使えない。以下、適宜注記するが、AttributeError
が出た場合などは公式ドキュメントを参照されたい。
data
属性に各要素の値を単純な配列で格納しているcsr_matrix
, csc_matrix
, coo_matrix
などは共通のメソッドを持っている。
平均・合計・最大値・最小値
numpy.ndarray
と同様に、平均mean()
、合計sum()
、最大値max()
、最小値mean()
を算出できる。
mean()
, sum()
は全クラス共通、max()
, mean()
はlil_matrix
などでは定義されていない。
l = [[0, 1, 2],
[3, 0, 4],
[0, 0, 0]]
csr = csr_matrix(l)
lil = lil_matrix(l)
print(csr.sum())
# 10
print(csr.mean())
# 1.111111111111111
print(csr.max())
# 4
print(csr.min())
# 0
# print(lil.max())
# AttributeError: max not found
各要素を処理(平方根や三角関数など)
要素ごとに平方根sqrt()
、三角関数sin()
, tan()
を適用するメソッドが、csr_matrix
, csc_matrix
, coo_matrix
などで定義されている。lil_matrix
などにはない。
print(csr.sqrt().toarray())
# [[0. 1. 1.41421356]
# [1.73205081 0. 2. ]
# [0. 0. 0. ]]
print(csr.sin().toarray())
# [[ 0. 0.84147098 0.90929743]
# [ 0.14112001 0. -0.7568025 ]
# [ 0. 0. 0. ]]
print(csr.tan().toarray())
# [[ 0. 1.55740772 -2.18503986]
# [-0.14254654 0. 1.15782128]
# [ 0. 0. 0. ]]
# print(lil.sqrt())
# AttributeError: sqrt not found
cos()
がないのは0
が1
になってしまうからだと思われる。
非ゼロ要素のみに対してcos()
を適用したい場合、csr_matrix
, csc_matrix
, coo_matrix
などではdata
属性(numpy.ndarray
)を処理して上書きすればよい。ここでは元の行列を保持するためcopy()
でコピーした行列を処理している。
csr_ = csr.copy()
print(csr_.data)
# [1 2 3 4]
print(type(csr_.data))
# <class 'numpy.ndarray'>
csr_.data = np.cos(csr_.data)
print(csr_.toarray())
# [[ 0. 0.54030231 -0.41614684]
# [-0.9899925 0. -0.65364362]
# [ 0. 0. 0. ]]
cos()
に限らず、要素に対して一括で任意の処理を適用できる。
csr_ = csr.copy()
csr_.data = csr_.data ** 2
print(csr_.toarray())
# [[ 0 1 4]
# [ 9 0 16]
# [ 0 0 0]]
例えばlil_matrix
の場合はdata
属性がリストを要素とするnumpy.ndarray
なので一括処理するのは面倒。
print(lil.data)
# [list([1, 2]) list([3, 4]) list([])]
非ゼロ要素の位置: nonzero()
非ゼロ要素の位置を取得するにはnonzero()
を使う。行・列のインデックスのリストが返される。
print(csr)
# (0, 1) 1
# (0, 2) 2
# (1, 0) 3
# (1, 2) 4
r, c = csr.nonzero()
print(r, c)
# [0 0 1 1] [1 2 0 2]
非ゼロ要素の個数はcount_nonzero()
、行列が保持しているデータ数はnnz
属性で取得できる。
print(csr.count_nonzero())
# 4
print(csr.nnz)
# 4
例えばcsr_matrix
ではある要素の値を0
としたときに値0
の要素としてデータが残るので、count_nonzero()
とnnz
の値が一致しなくなる。
csr[0, 1] = 0
print(csr)
# (0, 1) 0
# (0, 2) 2
# (1, 0) 3
# (1, 2) 4
print(csr.count_nonzero())
# 3
print(csr.nnz)
# 4
r, c = csr.nonzero()
print(r, c)
# [0 1 1] [2 0 2]
一方、lil_matrix
では0
にした要素は削除されるのでcount_nonzero()
とnnz
の値は一致する。種類によって振る舞いが異なるので注意。
print(lil)
# (0, 1) 1
# (0, 2) 2
# (1, 0) 3
# (1, 2) 4
print(lil.nnz)
# 4
lil[0, 1] = 0
print(lil)
# (0, 2) 2
# (1, 0) 3
# (1, 2) 4
print(lil.nnz)
# 3
scipy.sparseの関数
scipy.sparse
内にいくつか関数が定義されている。数は多くないのでscipy.sparse
行列を頻繁に使うなら一度目を通しておくといいだろう。
いくつか紹介する。
疎行列の連結: hstack(), vstack()
scipy.sparse
行列を連結するのがhstack()
(横)、vstack()
(縦)。
返り値のクラスは、デフォルトでは適切なものが選択される。引数format
に文字列で指定可能。
l = [[0, 1, 2],
[3, 0, 4],
[0, 0, 0]]
csr = csr_matrix(l)
lil = lil_matrix(l)
print(hstack([csr, lil]).toarray())
# [[0 1 2 0 1 2]
# [3 0 4 3 0 4]
# [0 0 0 0 0 0]]
print(type(hstack([csr, lil])))
# <class 'scipy.sparse.coo.coo_matrix'>
print(type(hstack([csr, lil], format='csr')))
# <class 'scipy.sparse.csr.csr_matrix'>
vstack()
でも同様。
print(vstack([csr, lil]).toarray())
# [[0 1 2]
# [3 0 4]
# [0 0 0]
# [0 1 2]
# [3 0 4]
# [0 0 0]]
print(type(vstack([csr, lil])))
# <class 'scipy.sparse.coo.coo_matrix'>
print(type(vstack([csr, lil], format='csr')))
# <class 'scipy.sparse.csr.csr_matrix'>
行・列をスライスなどで取り出して連結することも可能。hstack()
では行数、vstack()
では列数が一致していないとエラーになる。
print(vstack([csr, lil[:2]]).toarray())
# [[0 1 2]
# [3 0 4]
# [0 0 0]
# [0 1 2]
# [3 0 4]]
# print(hstack([csr, lil[:2]]))
# ValueError: blocks[0,:] has incompatible row dimensions. Got blocks[0,1].shape[0] == 2, expected 3.
上三角行列・下三角行列の取得: triu(), tril()
上三角行列をtriu()
、下三角行列をtril()
で取得できる。
hstack()
, vstack()
と同様に引数format
を指定可能。
l = [[0, 1, 2],
[3, 0, 4],
[5, 6, 7]]
csr = csr_matrix(l)
print(triu(csr).toarray())
# [[0 1 2]
# [0 0 4]
# [0 0 7]]
print(type(triu(csr)))
# <class 'scipy.sparse.coo.coo_matrix'>
print(type(triu(csr, format='csr')))
# <class 'scipy.sparse.csr.csr_matrix'>
例の通り、デフォルトでは主対角線を含む上側。引数k
で対角線の位置を移動できる。正が上側(右側)、負が下側(左側)。
print(triu(csr, k=1).toarray())
# [[0 1 2]
# [0 0 4]
# [0 0 0]]
print(triu(csr, k=-1).toarray())
# [[0 1 2]
# [3 0 4]
# [0 6 7]]
tril()
でも同様。
print(tril(csr).toarray())
# [[0 0 0]
# [3 0 0]
# [5 6 7]]
print(type(tril(csr)))
# <class 'scipy.sparse.coo.coo_matrix'>
print(type(tril(csr, format='csr')))
# <class 'scipy.sparse.csr.csr_matrix'>
print(tril(csr, k=1).toarray())
# [[0 1 0]
# [3 0 4]
# [5 6 7]]
print(tril(csr, k=-1).toarray())
# [[0 0 0]
# [3 0 0]
# [5 6 0]]
疎行列オブジェクトの保存・読み込み: save_npz(), load_npz()
scipy.sparse
行列のオブジェクトはsave_npz()
, load_npz()
はnpz
形式でそのまま保存・読み込みできる。
csr_matrix
やcsc_matrix
などの形式がそのまま保存される。
l = [[0, 1, 2],
[3, 0, 4],
[0, 0, 0]]
csr = csr_matrix(l)
csc = csc_matrix(l)
save_npz('data/temp/csr.npz', csr)
save_npz('data/temp/csc.npz', csc)
csr_ = load_npz('data/temp/csr.npz')
print(csr_.toarray())
# [[0 1 2]
# [3 0 4]
# [0 0 0]]
print(type(csr))
# <class 'scipy.sparse.csr.csr_matrix'>
csc_ = load_npz('data/temp/csc.npz')
print(csc_.toarray())
# [[0 1 2]
# [3 0 4]
# [0 0 0]]
print(type(csc))
# <class 'scipy.sparse.csc.csc_matrix'>
SciPy1.3.1
時点で対応しているのはcsr_matrix
, csc_matrix
, coo_matrix
, bsr_matrix
, dia_matrix
。lil_matrix
などはエラーになる。
lil = lil_matrix(l)
# save_npz('data/temp/lil.npz', lil)
# NotImplementedError: Save is not implemented for sparse matrix of format lil.
対応していない形式を保存したい場合はtocsr()
などで対応している形式に変換すればOK。
なお、npz
ファイルは複数のnumpy.ndarray
をまとめたもの。scipy.sparse
行列のnpz
ファイルの中身は以下のようになっている。あえて確認する意味はないが参考まで。
npz = np.load('data/temp/csr.npz')
print(npz.files)
# ['indices', 'indptr', 'format', 'shape', 'data']
print(npz['data'])
# [1 2 3 4]
print(npz['indices'])
# [1 2 0 2]
print(npz['indptr'])
# [0 2 4 4]
print(npz['format'])
# b'csr'
print(npz['shape'])
# [3 3]