pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

Python(NumPy)による行列の計算

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

pianofisica.hatenablog.com



行列の入力

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

 \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の『行列』をあつかうには2次元配列の 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()

 

ベクトルの入力

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

 \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 の対角成分を1とするように、対角行列との積 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)

により確かめられます。このように、行列を対角成分が1の下三角形行列 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 次直交行列 Qm\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 次ユニタリ行列 Un 次ユニタリ行列 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、行列、ベクトル、線形代数