Skip to content

Instantly share code, notes, and snippets.

@ohno

ohno/LinearAlgebra2.ipynb Secret

Last active January 30, 2022 01:40
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ohno/f1daed3883312bc7d1cd01513a7b48ca to your computer and use it in GitHub Desktop.
Save ohno/f1daed3883312bc7d1cd01513a7b48ca to your computer and use it in GitHub Desktop.
Download ZIP
Introduction to Linear Algebra on The Julia Programming Language
Display the source blob
Display the rendered blob
Raw
LinearAlgebra2.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Juliaで線形代数：固有値問題\n",
"\n",
"この記事を読む前に, Juliaでのベクトルや行列の扱いに慣れておこう.\n",
"\n",
"https://zenn.dev/ohno/articles/cb10dc5b3f5bbc\n",
"\n",
"この記事では固有値問題とエルミート行列について, Juliaによる具体的な計算例を挙げながら解説する.\n",
"\n",
"1. [パッケージ](#パッケージ)\n",
"2. [ベンチマーク](#ベンチマーク)\n",
"3. [行列の宣言](#行列の宣言)\n",
"4. [固有値](#固有値)\n",
"5. [固有ベクトル](#固有ベクトル)\n",
"6. [対角化](#対角化)\n",
"7. [`eigen()`の使い方](#eigen()の使い方)\n",
"8. [エルミート行列](#エルミート行列)\n",
"9. [そして量子力学へ](#そして量子力学へ)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## パッケージ\n",
"\n",
"LinearAlgebra.jlはJuliaに最初からインストールされているので, ノートブック上で`using LinearAlgebra`を宣言するだけでよい."
]
},
{
"cell_type": "code",
"execution_count": 99,
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ベンチマーク\n",
"\n",
"[斎藤正彦『線型代数入門』(1966, 東京大学出版会)]( http://www.utp.or.jp/book/b302039.html) 第5章の問題1をベンチマークとする. 対角形の各対角成分は元の行列の固有値である.\n",
"\n",
"<table>\n",
" <tr>\n",
" <th>問題1</th>\n",
" <th>行列</th>\n",
" <th>対角形</th>\n",
" </tr>\n",
" <tr>\n",
" <td>イ)</td>\n",
" <td>\n",
" \\begin{pmatrix}\n",
" -4 & 0 & 6 \\\\\n",
" -3 & 2 & 3 \\\\\n",
" -3 & 0 & 5 \\\\\n",
" \\end{pmatrix}\n",
" </td>\n",
" <td>\n",
" \\begin{pmatrix}\n",
" 2 & 0 & 0 \\\\\n",
" 0 & 2 & 0 \\\\\n",
" 0 & 0 & -1 \\\\\n",
" \\end{pmatrix}\n",
" </td>\n",
" </tr>\n",
" <tr>\n",
" <td>ロ)</td>\n",
" <td>\n",
" \\begin{pmatrix}\n",
" 1 & 0 & -2 \\\\\n",
" 2 & -1 & -2 \\\\\n",
" -2 & 2 & 0 \\\\\n",
" \\end{pmatrix}\n",
" </td>\n",
" <td>\n",
" \\begin{pmatrix}\n",
" 1 & 0 & 0 \\\\\n",
" 0 & -1 & 0 \\\\\n",
" 0 & 0 & 0 \\\\\n",
" \\end{pmatrix}\n",
" </td>\n",
" </tr>\n",
" <tr>\n",
" <td>ハ)</td>\n",
" <td>\n",
" \\begin{pmatrix}\n",
" -3 & -2 & -2 & 1 \\\\\n",
" 2 & 3 & 2 & 0 \\\\\n",
" 3 & 1 & 2 & -1 \\\\\n",
" -4 & -2 & -2 & 2 \\\\\n",
" \\end{pmatrix}\n",
" </td>\n",
" <td>\n",
" \\begin{pmatrix}\n",
" 0 & 0 & 0 & 0 \\\\\n",
" 0 & 2 & 0 & 0 \\\\\n",
" 0 & 0 & 1 & 0 \\\\\n",
" 0 & 0 & 0 & 1 \\\\\n",
" \\end{pmatrix}\n",
" </td>\n",
" </tr>\n",
"</table>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 行列の宣言\n",
"\n",
"行列は次のように宣言する."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Matrix{Int64}:\n",
" -4 0 6\n",
" -3 2 3\n",
" -3 0 5"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [-4 0 6;\n",
" -3 2 3;\n",
" -3 0 5]"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Matrix{Int64}:\n",
" 1 0 -2\n",
" 2 -1 -2\n",
" -2 2 0"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"B = [ 1 0 -2;\n",
" 2 -1 -2;\n",
" -2 2 0]"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Int64}:\n",
" -3 -2 -2 1\n",
" 2 3 2 0\n",
" 3 1 2 -1\n",
" -4 -2 -2 2"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C = [-3 -2 -2 1;\n",
" 2 3 2 0;\n",
" 3 1 2 -1;\n",
" -4 -2 -2 2]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 固有値\n",
"\n",
"事前に`using LinearAlgebra`を宣言しておく必要がある. 行列を`eigvals()`に渡すことで固有を格納した配列が得られる. なお, `eigvals()`は固有値を小さい順に並べるため, 上記の問題とは順序が異なる. また, Juliaに限らず, コンピュータで実数を扱うことは不可能なので, 代わりに[浮動小数点数](https://ja.wikipedia.org/wiki/%E6%B5%AE%E5%8B%95%E5%B0%8F%E6%95%B0%E7%82%B9%E6%95%B0)が採用される. そのため`-1`が`-0.999999999999998`となったり, `0`が`-2.539596429704816e-15`となったりする. 前述のベンチマークに矛盾しない結果が得られている."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" -1.0\n",
" 2.0\n",
" 2.0"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigvals(A)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" -0.999999999999998\n",
" -2.539596429704816e-15\n",
" 1.0000000000000024"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigvals(B)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Vector{Float64}:\n",
" 4.440892098500626e-15\n",
" 0.9999999999999977\n",
" 0.9999999999999997\n",
" 1.9999999999999996"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigvals(C)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 固有ベクトル\n",
"\n",
"固有値は`eigvec()`で求められる. 固有ベクトル(列ベクトル)を並べた行列が出力される."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Matrix{Float64}:\n",
" -0.816497 0.0 -0.707107\n",
" -0.408248 1.0 0.0\n",
" -0.408248 0.0 -0.707107"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigvecs(A)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Matrix{Float64}:\n",
" -0.666667 0.666667 -0.707107\n",
" -0.333333 0.666667 -0.707107\n",
" -0.666667 0.333333 3.66027e-16"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigvecs(B)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" -0.57735 0.235702 0.0116258 -0.408248\n",
" -8.46545e-16 0.471405 -0.712657 0.816497\n",
" 0.57735 -0.707107 0.701031 -1.46869e-15\n",
" -0.57735 0.471405 0.0232516 -0.408248"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigvecs(C)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 対角化\n",
"\n",
"行列$\\pmb{A}$の固有ベクトルを並べた行列$\\pmb{P}$を使うと, 対角形$\\pmb{P}^{-1}\\pmb{A}\\pmb{P}$を求めることができる. 固有値は小さい順に並び替えているため, 対角形も並び替えたものになっていることに注意."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Matrix{Float64}:\n",
" -1.0 0.0 -6.66134e-16\n",
" 0.0 2.0 0.0\n",
" -4.44089e-16 0.0 2.0"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = eigvecs(A)\n",
"inv(P) * A * P"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Matrix{Float64}:\n",
" -1.0 1.33227e-15 -1.77636e-15\n",
" -5.03301e-15 -3.15544e-30 -3.14018e-15\n",
" -5.32907e-15 2.66454e-15 1.0"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = eigvecs(B)\n",
"inv(P) * B * P"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" -2.56395e-16 3.97757e-15 -1.43785e-15 -3.80727e-15\n",
" -4.98164e-15 1.0 -9.74804e-16 -1.97652e-15\n",
" -1.78924e-15 -3.14804e-16 1.0 -1.39507e-15\n",
" 1.23016e-15 -8.76336e-16 -7.77466e-16 2.0"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = eigvecs(C)\n",
"inv(P) * C * P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## eigen()の使い方\n",
"\n",
"`eigen()`に行列を渡すと, 固有値と固有ベクトルを同時に求める. 固有値と固有ベクトルの配列はそれぞれ`values`と`vectors`というキーで個別に取り出すこともできる."
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}\n",
"values:\n",
"3-element Vector{Float64}:\n",
" -1.0\n",
" 2.0\n",
" 2.0\n",
"vectors:\n",
"3×3 Matrix{Float64}:\n",
" -0.816497 0.0 -0.707107\n",
" -0.408248 1.0 0.0\n",
" -0.408248 0.0 -0.707107"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 固有値と固有ベクトルを同時に求める例\n",
"eigen(A)"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" -1.0\n",
" 2.0\n",
" 2.0"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 固有値を取り出す例\n",
"eigen(A).values"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Matrix{Float64}:\n",
" -0.816497 0.0 -0.707107\n",
" -0.408248 1.0 0.0\n",
" -0.408248 0.0 -0.707107"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 固有ベクトルを取り出す例\n",
"eigen(A).vectors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## エルミート行列\n",
"\n",
"詳細は[Wikipedia](https://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%AB%E3%83%9F%E3%83%BC%E3%83%88%E8%A1%8C%E5%88%97)を参照されたい. 次のようなエルミート行列を考える."
]
},
{
"cell_type": "code",
"execution_count": 100,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Matrix{Complex{Int64}}:\n",
" 1+0im 0+2im\n",
" 0-2im 1+0im"
]
},
"execution_count": 100,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"H = [ 1 2im\n",
" -2im 1 ]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> エルミート行列の固有値は実数である.\n",
"\n",
"これは, 量子力学の支配方程式であるSchrödinger方程式が複素数を含む方程式であるにも関わらず, 導かれるエネルギー固有値が必ず実数であることと密接に関係している."
]
},
{
"cell_type": "code",
"execution_count": 101,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{Float64}:\n",
" -1.0\n",
" 3.0"
]
},
"execution_count": 101,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigvals(H)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> 異なる固有値に属する固有ベクトルは直交する.\n",
"\n",
"これも量子力学において重要である."
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0 + 0.0im\n"
]
}
],
"source": [
"P = eigvecs(H)\n",
"\n",
"println(dot(P[:,1], P[:,2]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## そして量子力学へ\n",
"\n",
"量子力学の支配方程式はSchrödinger方程式という微分方程式である. 一見, 線形代数とは関係無いように思われるが, 適当な基底関数で展開することによって, 行列方程式に帰着できる. シリーズ最後の記事では, その方法の一つである線形変分法について, Juliaによる具体的な計算例を挙げながら解説する.\n",
"\n",
"https://zenn.dev/ohno/articles/c48920c9327b16/"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 参考文献\n",
"- [斎藤 正彦『線型代数入門』(1966, 東京大学出版会)]( http://www.utp.or.jp/book/b302039.html) \n",
"- https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"© 2022 Shuhei Ohno\n",
"<br>Source: https://gist.github.com/ohno/f1daed3883312bc7d1cd01513a7b48ca\n",
"<br>License: https://opensource.org/licenses/MIT"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Julia Version 1.6.2\n",
"Commit 1b93d53fc4 (2021-07-14 15:36 UTC)\n",
"Platform Info:\n",
" OS: Windows (x86_64-w64-mingw32)\n",
" CPU: Intel(R) Core(TM) i7-4650U CPU @ 1.70GHz\n",
" WORD_SIZE: 64\n",
" LIBM: libopenlibm\n",
" LLVM: libLLVM-11.0.1 (ORCJIT, haswell)\n"
]
}
],
"source": [
"versioninfo()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.6.2",
"language": "julia",
"name": "julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment