🌾

Juliaで線形代数：固有値問題

2022/01/29に公開

Julia

この記事では固有値問題とエルミート行列について, Juliaによる具体的な計算例を挙げながら解説する.

パッケージ
ベンチマーク
行列の宣言
固有値
固有ベクトル
対角化
eigen()の使い方
エルミート行列
そして量子力学へ

パッケージ

LinearAlgebra.jlはJuliaに最初からインストールされているので, ノートブック上でusing LinearAlgebraを宣言するだけでよい.

パッケージの読み込み

using LinearAlgebra

ベンチマーク

斎藤正彦『線型代数入門』(1966, 東京大学出版会) 第5章の問題1をベンチマークとする. 対角形の各対角成分は元の行列の固有値である.

問題1	行列	対角形
イ)	$\begin{pmatrix} -4 & 0 & 6 \\ -3 & 2 & 3 \\ -3 & 0 & 5 \\ \end{pmatrix}$	$\begin{pmatrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & -1 \\ \end{pmatrix}$
ロ)	$\begin{pmatrix} 1 & 0 & -2 \\ 2 & -1 & -2 \\ -2 & 2 & 0 \\ \end{pmatrix}$	$\begin{pmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}$
ハ)	$\begin{pmatrix} -3 & -2 & -2 & 1 \\ 2 & 3 & 2 & 0 \\ 3 & 1 & 2 & -1 \\ -4 & -2 & -2 & 2 \\ \end{pmatrix}$	$\begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{pmatrix}$

行列の宣言

行列は次のように宣言する.

行列の宣言

A = [-4 0 6;
     -3 2 3;
     -3 0 5]

行列の宣言

B = [ 1  0 -2;
      2 -1 -2;
     -2  2  0]

行列の宣言

C = [-3 -2 -2  1;
      2  3  2  0;
      3  1  2 -1;
     -4 -2 -2  2]

固有値

事前にusing LinearAlgebraを宣言しておく必要がある. 行列をeigvals()に渡すことで固有を格納した配列が得られる. なお, eigvals()は固有値を小さい順に並べるため, 上記の問題とは順序が異なる. また, Juliaに限らず, コンピュータで実数を扱うことは不可能なので, 代わりに浮動小数点数が採用される. そのため-1が-0.999999999999998となったり, 0が-2.539596429704816e-15となったりするが, 前述のベンチマークに矛盾しない結果が得られている.

固有値

eigvals(A)

出力

-1.0
2.0
2.0

固有値

eigvals(B)

出力

-0.999999999999998
-2.539596429704816e-15
1.0000000000000024

固有値

eigvals(C)

出力

4.440892098500626e-15
0.9999999999999977
0.9999999999999997
1.9999999999999996

固有ベクトル

固有値はeigvec()で求められる. 固有ベクトル(列ベクトル)を並べた行列が出力される.

固有ベクトル

eigvecs(A)

出力

-0.816497  0.0  -0.707107
-0.408248  1.0   0.0
-0.408248  0.0  -0.707107

固有ベクトル

eigvecs(B)

出力

-0.666667  0.666667  -0.707107
-0.333333  0.666667  -0.707107
-0.666667  0.333333   3.66027e-16

固有ベクトル

eigvecs(C)

出力

-0.57735       0.235702   0.0116258  -0.408248
-8.46545e-16   0.471405  -0.712657    0.816497
0.57735      -0.707107   0.701031   -1.46869e-15
-0.57735       0.471405   0.0232516  -0.408248

対角化

行列 $\pmb{A}$ の固有ベクトルを並べた行列 $\pmb{P}$ を使うと, 対角形 $\pmb{P}^{-1}\pmb{A}\pmb{P}$ を求めることができる. 固有値は小さい順に並び替えているため, 対角形も並び替えたものになっていることに注意.

対角化

P = eigvecs(A)
inv(P) * A * P

出力

-1.0          0.0  -6.66134e-16
0.0          2.0   0.0
-4.44089e-16  0.0   2.0

対角化

P = eigvecs(B)
inv(P) * B * P

出力

-1.0           1.33227e-15  -1.77636e-15
-5.03301e-15  -3.15544e-30  -3.14018e-15
-5.32907e-15   2.66454e-15   1.0

対角化

P = eigvecs(C)
inv(P) * C * P

出力

-2.56395e-16   3.97757e-15  -1.43785e-15  -3.80727e-15
-4.98164e-15   1.0          -9.74804e-16  -1.97652e-15
-1.78924e-15  -3.14804e-16   1.0          -1.39507e-15
1.23016e-15  -8.76336e-16  -7.77466e-16   2.0

eigen()の使い方

eigen()に行列を渡すと, 固有値と固有ベクトルを同時に求める. 固有値と固有ベクトルの配列はそれぞれvaluesとvectorsというキーで個別に取り出すこともできる.

固有値と固有ベクトルを同時に求める例

eigen(A)

出力

values:
    3-element Vector{Float64}:
     -1.0
      2.0
      2.0
vectors:
    3×3 Matrix{Float64}:
     -0.816497  0.0  -0.707107
     -0.408248  1.0   0.0
     -0.408248  0.0  -0.707107

固有値を取り出す例

eigen(A).values

出力

-1.0
2.0
2.0

固有ベクトルを取り出す例

eigen(A).vectors

出力

-0.816497  0.0  -0.707107
-0.408248  1.0   0.0
-0.408248  0.0  -0.707107

エルミート行列

詳細はWikipediaを参照されたい. 次のようなエルミート行列を考える.

エルミート行列の例

H = [  1   2im
     -2im   1  ]

エルミート行列の固有値は実数である.

これは, 量子力学の支配方程式であるSchrödinger方程式が複素数を含む方程式であるにも関わらず, 導かれるエネルギー固有値が必ず実数であることと密接に関係している.

固有値が実であることの確認

eigvals(H)

出力

-1.0
3.0

異なる固有値に属する固有ベクトルは直交する.

これも量子力学において重要である.

固有ベクトルの直交性

P = eigvecs(H)

println(dot(P[:,1], P[:,2]))

出力

0.0 + 0.0im

そして量子力学へ

量子力学の支配方程式はSchrödinger方程式という微分方程式である. 一見, 線形代数とは関係無いように思われるが, 適当な基底関数で展開することによって, 行列方程式に帰着できる. シリーズ最後の記事では, その方法の一つである線形変分法について, Juliaによる具体的な計算例を挙げながら解説する.

参考文献

この記事の元になったノートブックはGistにアップロードしてある.

パッケージ

ベンチマーク

行列の宣言

固有値

固有ベクトル

対角化

eigen()の使い方

エルミート行列

そして量子力学へ

参考文献

Discussion