Python, NumPyで行列の演算（逆行列、行列式、固有値など）
Pythonで行列の演算を行うにはNumPyを使うと便利。Python組み込みのリストlistでも二次元配列（リストのリスト）を表現できるが、NumPyを使うと行列の積や逆行列、行列式、固有値などを簡単に算出できる。
NumPyはpipコマンドなどでインストールできる。
$ pip install numpy
NumPyには汎用的な多次元配列のnp.ndarrayクラスと行列（二次元配列）に特化したnp.matrixクラスがある。ただし、np.ndarrayでも行列演算は可能で、NumPy 1.26時点ではnp.matrixの使用は推奨されていない。
It is no longer recommended to use this class, even for linear algebra. Instead use regular arrays. The class may be removed in the future. numpy.matrix — NumPy v1.26 Manual
本記事でもnp.matrixには触れず、np.ndarrayを使った行列演算について説明する。
SciPyを使うと要素のほとんどが0である疎行列を効率的に扱うことができる。
本記事のサンプルコードのNumPyのバージョンは以下の通り。バージョンによって仕様が異なる可能性があるので注意。
import numpy as np
print(np.__version__)
# 1.26.2
行列を表すNumPy配列ndarrayの作成: np.array()
NumPy配列ndarrayはnp.array()で作成できる。
第一引数に二次元リスト（リストのリスト）を指定すると、二次元配列が生成される。
a = np.array([[0, 1, 2], [3, 4, 5]])
print(a)
# [[0 1 2]
# [3 4 5]]
print(type(a))
# <class 'numpy.ndarray'>
等差数列を生成する関数np.arange()などでもndarrayを生成できる。さらにreshape()メソッドで形状を変更することも可能。
print(np.arange(6).reshape(2, 3))
# [[0 1 2]
# [3 4 5]]
そのほか、任意の形状を同じ値で初期化するnp.zeros(), np.ones(), np.full()などもある。CSVファイルを読み込むこともできる。
行列の要素の取得と変更（代入）
二次元のNumPy配列ndarrayの要素には[行番号, 列番号]でアクセスできる。行番号・列番号は0始まり。値の取得および変更（代入）が可能。
a = np.array([[0, 1, 2], [3, 4, 5]])
print(a)
# [[0 1 2]
# [3 4 5]]
print(a[0, 1])
# 1
a[0, 1] = 100
print(a)
# [[ 0 100 2]
# [ 3 4 5]]
スライスなどを使って範囲を選択して部分配列を取得することもできる。
print(a[:, 1:])
# [[100 2]
# [ 4 5]]
行列の要素ごとの計算: 算術演算子
+（加算）や-（減算）、*(乗算)、/（除算）、**（べき乗）などの算術演算子をndarray同士に使うと、要素ごとの計算となる。
a1 = np.arange(6).reshape(2, 3)
print(a1)
# [[0 1 2]
# [3 4 5]]
a2 = np.arange(6, 18, 2).reshape(2, 3)
print(a2)
# [[ 6 8 10]
# [12 14 16]]
print(a1 + a2)
# [[ 6 9 12]
# [15 18 21]]
print(a1 - a2)
# [[ -6 -7 -8]
# [ -9 -10 -11]]
print(a1 * a2)
# [[ 0 8 20]
# [36 56 80]]
print(a1 / a2)
# [[0. 0.125 0.2 ]
# [0.25 0.28571429 0.3125 ]]
print(a1**a2)
# [[ 0 1 1024]
# [ 531441 268435456 152587890625]]
スカラー値との計算も可能。
print(a1 * 100)
# [[ 0 100 200]
# [300 400 500]]
異なる形状でも、ブロードキャストと呼ばれる仕組みによって可能な限り自動的に形状が変換されて計算される。
- 関連記事: NumPyのブロードキャスト（形状の自動変換）
print(a1 * [10, 100, 1000])
# [[ 0 100 2000]
# [ 30 400 5000]]
負の値のべき乗の場合、元の配列のデータ型dtypeが整数intだとエラーになる。浮動小数点数floatに変換すればよい。
# print(a2 ** -1)
# ValueError: Integers to negative integer powers are not allowed.
print(a2.astype(float) ** -1)
# [[0.16666667 0.125 0.1 ]
# [0.08333333 0.07142857 0.0625 ]]
*演算子に対応するnp.multiply()など、各演算子に対応する関数も提供されている。
print(np.multiply(a1, a2))
# [[ 0 8 20]
# [36 56 80]]
行列の積（掛け算）: @演算子, np.matmul(), np.dot()
行列の積（掛け算）を算出するには、@演算子、np.matmul(), np.dot()を使う。dot()はndarrayのメソッドとしても提供されている。
- numpy.matmul — NumPy v1.26 Manual
- numpy.dot — NumPy v1.26 Manual
- numpy.ndarray.dot — NumPy v1.26 Manual
@演算子はPython 3.5, NumPy 1.10以降で利用可能で、a @ bはnp.matmul(a, b)と等価。
a1 = np.arange(4).reshape((2, 2))
print(a1)
# [[0 1]
# [2 3]]
a2 = np.arange(6).reshape((2, 3))
print(a2)
# [[0 1 2]
# [3 4 5]]
print(a1 @ a2)
# [[ 3 4 5]
# [ 9 14 19]]
print(np.matmul(a1, a2))
# [[ 3 4 5]
# [ 9 14 19]]
print(np.dot(a1, a2))
# [[ 3 4 5]
# [ 9 14 19]]
print(a1.dot(a2))
# [[ 3 4 5]
# [ 9 14 19]]
三次元以上の多次元配列に対してはnp.matmul()とnp.dot()で処理が異なるがここでは深追いしない。詳細は上記の公式ドキュメントを参照。
なお、np.dot()は名前の通り一般的なドット積（dot product）のための関数であり、行列の積（二次元配列同士に対する処理）にはnp.matmul()または@演算子の使用が好ましいとされている。
If both
aandbare 2-D arrays, it is matrix multiplication, but usingmatmulora @ bis preferred. numpy.dot — NumPy v1.26 Manual
一次元配列と二次元配列の場合
一行のみの行列を行ベクトル、一列のみの行列を列ベクトルと呼ぶことがあるが、ndarrayの一次元配列には行・列の区別はない。一行のみまたは一列のみであることを明確に示すには二次元配列で表す。
a_row = np.arange(3).reshape(1, 3)
print(a_row)
# [[0 1 2]]
print(a_row.shape)
# (1, 3)
a_col = np.arange(3).reshape(3, 1)
print(a_col)
# [[0]
# [1]
# [2]]
print(a_col.shape)
# (3, 1)
二次元配列同士の積は当然ながら通常の行列の積として計算され、結果は二次元配列となる。
a_2d = np.arange(9).reshape(3, 3)
print(a_2d)
# [[0 1 2]
# [3 4 5]
# [6 7 8]]
print(a_row @ a_2d)
# [[15 18 21]]
print(a_2d @ a_col)
# [[ 5]
# [14]
# [23]]
@演算子、np.matmul(), np.dot()で一次元配列と二次元配列の積を算出する場合、第一引数（左オペランド）の一次元配列は行ベクトル、第二引数（右オペランド）の一次元配列は列ベクトルとして扱われた上で、結果は二次元配列ではなく一次元配列となる。
a_1d = np.array([0, 1, 2])
print(a_1d @ a_2d)
# [15 18 21]
print(np.dot(a_1d, a_2d))
# [15 18 21]
print(a_2d @ a_1d)
# [ 5 14 23]
print(np.dot(a_2d, a_1d))
# [ 5 14 23]
一次元配列同士の場合は、次の節で述べるようにベクトルの内積を返す。
ベクトルの内積: np.inner()
行列（二次元配列）というテーマからずれるが、ベクトル（一次元配列）の内積を算出するにはnp.inner()を使う。返り値はスカラー値。
a1 = np.array([0, 1, 2])
a2 = np.array([3, 4, 5])
print(np.inner(a1, a2))
# 14
@演算子、np.matmul(), np.dot()でも引数が両方とも一次元配列だと内積が返される。
print(a1 @ a2)
# 14
print(np.dot(a1, a2))
# 14
逆行列: np.linalg.inv()
逆行列を算出するにはnp.linalg.inv()を使う。元の行列と逆行列の積は単位行列になる。
a = np.array([[2, 5], [1, 3]])
print(a)
# [[2 5]
# [1 3]]
print(np.linalg.inv(a))
# [[ 3. -5.]
# [-1. 2.]]
print(a @ np.linalg.inv(a))
# [[1. 0.]
# [0. 1.]]
擬似逆行列: np.linalg.pinv()
特異行列（逆行列が存在しない行列）に対してはnp.linalg.inv()はエラーとなる。
a_singular = np.array([[0, 0], [1, 3]])
print(a_singular)
# [[0 0]
# [1 3]]
# print(np.linalg.inv(a_singular))
# LinAlgError: Singular matrix
np.linalg.pinv()を使うとムーア-ペンローズの擬似逆行列（Moore-Penrose pseudo-inverse matrix）を算出できる。
print(np.linalg.pinv(a_singular))
# [[0. 0.1]
# [0. 0.3]]
特異行列とその擬似逆行列の行列積は単位行列にならないが、擬似逆行列の擬似逆行列は元の行列と等しいといった性質がある。詳しい定義や性質などは上記のWikipediaのページを参照。
print(a_singular @ np.linalg.pinv(a_singular))
# [[0. 0.]
# [0. 1.]]
print(np.linalg.pinv(np.linalg.pinv(a_singular)))
# [[0. 0.]
# [1. 3.]]
なお、逆行列が存在する行列（正則行列、非特異行列）に対しては、np.linalg.inv()もnp.linalg.pinv()もどちらも逆行列を返す。
a = np.array([[2, 5], [1, 3]])
print(a)
# [[2 5]
# [1 3]]
print(np.linalg.inv(a))
# [[ 3. -5.]
# [-1. 2.]]
print(np.linalg.pinv(a))
# [[ 3. -5.]
# [-1. 2.]]
行列式: np.linalg.det()
行列式を算出するにはnp.linalg.det()を使う。
a = np.array([[0, 1], [2, 3]])
print(a)
# [[0 1]
# [2 3]]
print(np.linalg.det(a))
# -2.0
固有値と固有ベクトル: np.linalg.eig()
固有値と固有ベクトルを算出するにはnp.linalg.eig()を使う。
n次正方行列に対して、n個の固有値が要素となる一次元配列と、各列が各々の固有値に対応する固有ベクトルとなる二次元配列が返される。
i個目（[i]）の固有値と、i列目（[:, i]）の固有ベクトルが対応している。固有ベクトルは長さ1に規格化した値。
a = np.array([[8, 1], [4, 5]])
print(a)
# [[8 1]
# [4 5]]
eigenvalues, eigenvectors = np.linalg.eig(a)
print(eigenvalues)
# [9. 4.]
print(eigenvectors)
# [[ 0.70710678 -0.24253563]
# [ 0.70710678 0.9701425 ]]
print(a @ eigenvectors[:, 0] == eigenvalues[0] * eigenvectors[:, 0])
# [ True True]
print(a @ eigenvectors[:, 1] == eigenvalues[1] * eigenvectors[:, 1])
# [ True True]
最大固有値を求めたい場合は、最大値のインデックスを返すnp.argmax()を使う。
print(eigenvalues[np.argmax(eigenvalues)])
# 9.0
print(eigenvectors[:, np.argmax(eigenvalues)])
# [0.70710678 0.70710678]
以下のような関数を定義すると、対応する固有値と固有ベクトルのペアのタプルのリストを取得できる。固有ベクトルは長さ1で規格化された値だとよく分からないので、各列の0ではない最も絶対値が小さい値で割っている。出力にはf文字列を使用。
def get_eigenpairs(a):
values, vectors = np.linalg.eig(a)
pairs = []
for i, val in enumerate(values):
vec = vectors[:, i] / min(x for x in np.abs(vectors[:, i]) if x != 0)
pairs.append((val, vec))
return pairs
for val, vec in get_eigenpairs(a):
print(f'value: {val}, vector: {vec}')
# value: 9.0, vector: [1. 1.]
# value: 4.0, vector: [-1. 4.]
3 x 3の例。
a_3d = np.array([[1, 1, 2], [0, 2, -1], [0, 0, 3]])
print(a_3d)
# [[ 1 1 2]
# [ 0 2 -1]
# [ 0 0 3]]
for val, vec in get_eigenpairs(a_3d):
print(f'value: {val}, vector: {vec}')
# value: 1.0, vector: [1. 0. 0.]
# value: 2.0, vector: [1. 1. 0.]
# value: 3.0, vector: [ 1. -2. 2.]
複素数も扱うことができる。虚数単位はjで表される。
a = np.array([[3, 2], [-2, 3]])
print(a)
# [[ 3 2]
# [-2 3]]
for val, vec in get_eigenpairs(a):
print(f'value: {val:.1}, vector: {vec}')
# value: (3+2j), vector: [0.-1.j 1.+0.j]
# value: (3-2j), vector: [0.+1.j 1.-0.j]
Pythonでの複素数の扱いについては以下の記事を参照。
転置行列、対角行列、対称行列、三角行列など
転置行列や三角行列、対角行列、対称行列などについては以下の記事を参照。