Python（NumPy）による行列の計算

プログラミング言語 Pythonの数値計算ライブラリNumPyで行列、その計算をあつかう方法を見てみたいと思います。行列といえば連立一次方程式が思い浮かびますが、NumPyで方程式を解く方法は次の過去記事で解説しています。あわせて読んでみてください。

pianofisica.hatenablog.com

行列の入力
ゼロ行列の入力
単位行列の入力
対角行列の入力
行列の線形和・積・冪乗
行列の転置
複素係数の行列
行列の複素共役
行列の随伴共役
ベクトルの入力
ベクトルの線形和
ベクトルの内積
ベクトルのノルム
行列の固有値・固有ベクトル
行列式
行列の固有和（トレース）
逆行列
LU分解
コレスキー分解
QR分解
特異値分解

行列の入力

次のような行列は以下のようにして入力します。

$\quad \displaystyle{A=\begin{pmatrix} 2 & 4 & 5 \\ 1 & 3 & 2 \\ 4 & 2 & 1 \end{pmatrix} \qquad B=\begin{pmatrix} 2 & 5 & 4 \\ 1 & 2 & 3 \\ 4 & 1 & 2 \end{pmatrix}}$

import numpy as np
A = np.array([
    [2, 4, 5],
    [1, 3, 2],
    [4, 2, 1]])
B = np.array([
    [2, 5, 4],
    [1, 2, 3],
    [4, 1, 2]])

注意：NumPyの『行列』をあつかうには２次元配列の numpy.array と行列の numpy.matrix とがあります。ここではサブモジュール linalg でのあつかいを可能にする前者（numpy.array）によって行列を記述することとします。

ゼロ行列の入力

要素がすべてゼロである次のような行列は以下のようにして入力します。

$\quad \displaystyle{O=\begin{pmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \end{pmatrix}}$

O = np.zeros((3,2))    # 引数はタプル型
print(O)

単位行列の入力

単位行列は以下のようにして入力します。

id3 = np.identity(3)
print(id3)

対角行列の入力

対角行列は以下のようにして入力します。

D = np.diag([1,2,3])
print(D)

行列の線形和・積・冪乗

例えば

$\quad \displaystyle{-2A+2B=-2\begin{pmatrix} 2 & 4 & 5 \\ 1 & 3 & 2 \\ 4 & 2 & 1 \end{pmatrix}+2\begin{pmatrix} 2 & 5 & 4 \\ 1 & 2 & 3 \\ 4 & 1 & 2 \end{pmatrix}=\begin{pmatrix} 0 & 2 & -2 \\ 0 & -2 & 2 \\ 0 & -2 & 2 \end{pmatrix}}$

は

-2*A+2*B

によって計算できます。行列の積（例えば $AB$ ）については

A@B

として計算します。冪乗（例えば $A^3$ ）は

np.linalg.matrix_power(A, 3)

です。実際

np.linalg.matrix_power(A, 3)-A@A@A

です。

行列の転置

B.transpose()

行列の対称部分、反対称部分はそれぞれ

$\quad \displaystyle{B_{Sym}=\frac12(B+B^{T}) \qquad B_{Asym}=\frac12(B-B^{T})}$

によって取り出すことができます。

1/2*(B+B.transpose())

1/2*(B-B.transpose())

複素係数の行列

C = np.array([
    [2, 1+2j],
    [1-2j, 3]])
T2 = np.array([
    [0, -1j],
    [1j, 0]])

T2@C@T2@C.transpose()

行列の複素共役

C.conjugate()

行列の随伴共役

C.conjugate().transpose()

ベクトルの入力

ベクトルは行列の列（または行）が１つしかない特殊な場合です。例えば

$\quad \displaystyle{u=\begin{pmatrix} 3 \\ 2 \\ 1 \end{pmatrix} \qquad v=\begin{pmatrix} -1 \\ 3 \\ -3 \end{pmatrix}}$

は

u = np.array([
    [3],[2],[1]])
v = np.array([
    [-1],[3],[-3]])

と入力します。

ベクトルの線形和

2*u+v

ベクトルの内積

np.vdot(u,v)

ベクトルのノルム

np.linalg.norm(u)

行列の固有値・固有ベクトル

r, S = np.linalg.eig(A)

により、行列 $A$ の固有値が

$\quad \displaystyle{r =[\lambda_1= 7.78728013, \quad \lambda_2=-2.86354908, \quad \lambda_3= 1.07626895]}$

と求まり、対角化行列が

$\quad \displaystyle{S= \begin{pmatrix} -0.7411868 & -0.64930771 & 0.38693013 \\ -0.38466593 & -0.14398058 & -0.75359587 \\ -0.55015838 & 0.74677245 & 0.53139282 \end{pmatrix} }$

と求まります。よって対応する固有ベクトルが

$\quad \displaystyle{u_1=\begin{pmatrix} -0.7411868 \\ -0.38466593 \\ -0.55015838 \end{pmatrix} \qquad u_2=\begin{pmatrix} -0.64930771\\ -0.14398058\\ 0.74677245 \end{pmatrix} \qquad u_3=\begin{pmatrix} 0.38693013\\ -0.75359587\\ 0.53139282 \end{pmatrix}}$

とわかります。固有値ベクトルの定義から $Au_i=\lambda_iu_i$ なので、 $\frac{1}{\lambda_i}Au_i$ を計算してみると

r1 = 7.78728013
u1 = np.array([
    [-0.7411868],[-0.38466593],[-0.55015838]])
1/r1*A@u1

などとして $\frac{1}{\lambda_i}Au_i=u_i$ を数値的に確かめることができます。

行列式

行列式は

np.linalg.det(A)

で計算されます。行列式の性質から、行列式の値は固有値の積に等しいことが知られています。

$\quad \displaystyle{\det(A)=\prod_i\lambda_i}$

これを確かめてみましょう。

r1 = 7.78728013
r2 = -2.86354908
r3 = 1.07626895
np.linalg.det(A)/(r1*r2*r3)

より、公式が数値的に成り立っていることがわかります。

行列の固有和（トレース）

行列の固有和（対角成分の和、トレース）は

np.trace(A)

で計算されます。トレースの性質から、トレースの値は固有値の和に等しいことが知られています。

$\quad \displaystyle{{\rm Tr}(A)=\sum_i\lambda_i}$

確かめてみると

np.trace(A)/(r1+r2+r3)

より、公式が成り立っていることが数値的にわかります。

逆行列

逆行列は

np.linalg.inv(A)

で計算されます。検算してみると

A@np.linalg.inv(A)

np.linalg.inv(A)@A

より、確かに逆行列を与えます。また、上でみた行列 $A$ の対角化行列 $S$ について $S^{-1}AS$ を計算してみると

np.linalg.inv(S)@A@S

より

$\quad \displaystyle{S^{-1}AS={\rm diag}(7.78\dots,\ -2.86\dots,\ 1.07\dots)}$

となって、対角化されることがわかります。

LU分解

行列を下三角形行列 $L$ と上三角形行列 $U$ （と並び替え行列 $P$ ）の積として表すのがLU分解です。

import scipy.linalg as la
P, L, U = la.lu(A)

いまの具体例の場合には

$\quad \displaystyle{P=\begin{pmatrix} 0 & 1 & 0 \\ 0 & 0& 1 \\ 1& 0& 0 \end{pmatrix} \qquad L=\begin{pmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ 0.25 & 0.8333 & 1 \end{pmatrix} \qquad U=\begin{pmatrix} 4 & 2 & 1 \\ 0 & 3 & 4.5 \\ 0 & 0 & -2 \end{pmatrix}}$

と求まります。検算すると

A-(P@L@U)

確かに $A=PLU$ となっていることがわかります。さらに上三角行列 $U$ の対角成分を１とするように、対角行列との積 $U=D\tilde{U}$ で表すと

$\quad \displaystyle{D=\begin{pmatrix} 4 & 0 & 0 \\ 0 & 3 & 0 \\ 0& 0& -2 \end{pmatrix} \qquad \tilde{U}=\begin{pmatrix} 1 & 0.5 & 0.25 \\ 0 & 1 & 1.5 \\ 0 & 0 & 1 \end{pmatrix} }$

がわかります。実際

D = np.diag([4,3,-2])
U0 = np.array([
    [1,0.5,0.25],
    [0,1,1.5],
    [0,0,1]])
U-(D@U0)

により確かめられます。このように、行列を対角成分が１の下三角形行列 $L$ と上三角形行列 $U$ 、そして対角行列 $D$ の積として $LDU$ と表したとき、これを行列のLDU分解と言います。

コレスキー分解

コレスキー分解（Cholesky decomposition）は、正定値（固有値がすべて正値）エルミート行列に対するLU分解の特別の場合で、行列を下三角行列 $L$ とそのエルミート随伴によって得られる上三角形行列 $L^\dagger={L^*}^{\rm T}$ の積によって表したものです。

H = np.array([
    [3,1,1],
    [1,3,-1],
    [1,-1,2]])
Lh = np.linalg.cholesky(H)

実際

Lh@Lh.conjugate().transpose()

によって確かめられます。

QR分解

QR分解は、 $m\times n$ -実行列を、 $m$ 次直交行列 $Q$ と $m\times n$ -上三角行列 $R$ との積として表すことです。

Q, R = np.linalg.qr(B)

より

$\quad \displaystyle{Q=\begin{pmatrix} -0.43\dots & -0.82\dots & -0.36\dots \\ -0.21\dots & -0.29\dots & 0.93\dots \\ -0.87\dots & 0.48\dots & -0.05\dots \end{pmatrix} \qquad R=\begin{pmatrix} -4.58\dots & -3.49\dots & -4.14\dots \\ 0 & -4.22\dots & -3.20\dots \\ 0 & 0 & 1.24\dots \end{pmatrix} }$

と求まります。確かめてみましょう。

Q.transpose()@Q

Q@Q.transpose()

より、行列 $Q$ が直交行列であることがわかります。さらに

Q@R

より、 $B=QR$ が確認できます。

特異値分解

特異値分解は、 $m\times n$ -行列を、 $m$ 次ユニタリ行列 $U$ と $n$ 次ユニタリ行列 $V$ 、そして以下で説明する行列の特異値によって定まる $m\times n$ -行列 $\Sigma$ の積として表すことです。

PythonのNumPyライブラリでの組み込み関数

K = np.array([
    [3,1,1],
    [1,3,-1]])
U, sigma, V = np.linalg.svd(K)

により

$\quad \displaystyle{K=\begin{pmatrix} 3 & 1 & 1 \\ 1 & 3 & -1 \end{pmatrix} }$

の特異値分解 $K=U\Sigma V$ が

$\quad \displaystyle{\begin{aligned} U&=\begin{pmatrix} -0.70\dots & -0.70\dots \\ -0.70 \dots & 0.70 \end{pmatrix} \\ \Sigma&=\begin{pmatrix} 4.00\dots & 0 & 0 \\ 0 & 2.44\dots & 0 \end{pmatrix} \\ V&=\begin{pmatrix} -0.70\dots & -0.70\dots & 0 \\ -0.57\dots & 0.57\dots & -0.57\dots \\ -0.40\dots & 0.40\dots & 0.81\dots \end{pmatrix} \end{aligned}}$

によって与えられると計算されます。確認すると

U@U.conjugate().transpose()

U.conjugate().transpose()@U

V@V.conjugate().transpose()

V.conjugate().transpose()@V

より $U$ と $V$ がユニタリ行列であることがわかります。さらに、特異値ですが、これは正値行列 $KK^\dagger$ の固有値の正のルートになっています。実際、

M = K@K.conjugate().transpose()
np.linalg.eig(M)

より $M=KK^\dagger$ の固有値が

$\quad \displaystyle{\lambda^{(1)}_{KK^\dagger}=16\qquad \lambda^{(2)}_{KK^\dagger}=6}$

と求まりますが、その正のルートは

$\quad \displaystyle{\sigma_1=\sqrt{\lambda^{(1)}_{KK^\dagger}}=4\qquad \sigma_2=\sqrt{\lambda^{(2)}_{KK^\dagger}}=\sqrt{6}=2.449\dots}$

です。すなわち、行列 $\Sigma$ は

$\quad \displaystyle{\Sigma=\begin{pmatrix} \sigma_1 & 0 & 0 \\ 0 & \sigma_2 & 0 \end{pmatrix} }$

という、対角成分に特異値が並んだ形をとります。最後に分解の検算として

Sigma = np.array([
    [sigma[0],0,0],
    [0,sigma[1],0]])
U@Sigma@V

より、 $K=U\Sigma V$ が確認できます。

今回は、Python（NumPy）で行列をあつかう方法についてまとめて解説しました。

キーワード：Python、NumPy、行列、ベクトル、線形代数