Python, SciPyで疎行列の計算・処理(逆行列、固有値、連結、保存など)

Posted: | Tags: 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_matrixcsc_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.sparsenumpy.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()の第一引数に指定する行列のデータ型dtypefloatである必要がある。整数intだとエラー。

# w, v = eigs(csr_i, k=1)
# ValueError: matrix type must be 'f', 'd', 'F', or 'D'

scipy.sparseのクラスのメソッド

csr_matrixlil_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()がないのは01になってしまうからだと思われる。

非ゼロ要素のみに対して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_matrixcsc_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_matrixlil_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]

関連カテゴリー

関連記事