🌾

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

2022/01/29に公開約4,700字

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

  1. パッケージ
  2. ベンチマーク
  3. 行列の宣言
  4. 固有値
  5. 固有ベクトル
  6. 対角化
  7. eigen()の使い方
  8. エルミート行列
  9. そして量子力学へ

パッケージ

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()に行列を渡すと, 固有値と固有ベクトルを同時に求める. 固有値と固有ベクトルの配列はそれぞれvaluesvectorsというキーで個別に取り出すこともできる.

固有値と固有ベクトルを同時に求める例
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による具体的な計算例を挙げながら解説する.

https://zenn.dev/ohno/articles/c48920c9327b16

参考文献

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

https://gist.github.com/ohno/f1daed3883312bc7d1cd01513a7b48ca

Discussion

ログインするとコメントできます