Skip to content

Instantly share code, notes, and snippets.

@h3y6e

h3y6e/eigenvectors-from-eigenvalues.ipynb

Created November 19, 2019 08:18
Show Gist options
  • Save h3y6e/1db43d4fdb90a9ecf51636619e44d3ce to your computer and use it in GitHub Desktop.
Save h3y6e/1db43d4fdb90a9ecf51636619e44d3ce to your computer and use it in GitHub Desktop.
Download ZIP
"Eigenvectors from Eigenvalues" を Julia で検証する
Display the source blob
Display the rendered blob
Raw
eigenvectors-from-eigenvalues.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Eigenvectors from Eigenvalues\n",
"\n",
" - [Eigenvectors from Eigenvalues](https://arxiv.org/pdf/1908.03795.pdf)\n",
" - Peter B. Denton, Stephen J. Parke, Terence Tao, Xining Zhang\n",
"\n",
"$n$ 番煎じだと思うが，Julia で検証してみる．\n",
"\n",
"この論文の重要な部分は **Lemma 2** である．\n",
"$$\n",
"|v_{i,j}|^2 \\prod_{k = 1; k \\neq i}^n{(\\lambda_i(A) - \\lambda_k(A))} = \\prod_{k = 1}^{n-1}{(\\lambda_i(A) - \\lambda_k(M_j))}\n",
"$$\n",
"ただし，\n",
" - $A$ : $n \\times n$ のエルミート行列\n",
" - $\\lambda_i(A)$ : $A$ の $i$ 番目の固有値\n",
" - $v_{i,j}$ : $\\lambda_i(A)$ に対する固有ベクトル $v_i$ の $j$ 番目の要素\n",
" - $M_j$ : $A$ から第 $j$ 行と第 $j$ 列を除去して得られた主小行列\n",
"\n",
"この関係式により，固有値(と主小行列固有値)から固有ベクトル(の成分の二乗ノルム)を計算する事が出来る．"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 実行環境"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Julia Version 1.2.0\n",
"Commit c6da87ff4b (2019-08-20 00:03 UTC)\n",
"Platform Info:\n",
" OS: macOS (x86_64-apple-darwin18.6.0)\n",
" CPU: Intel(R) Core(TM) i7-8569U CPU @ 2.80GHz\n",
" WORD_SIZE: 64\n",
" LIBM: libopenlibm\n",
" LLVM: libLLVM-6.0.1 (ORCJIT, skylake)\n"
]
}
],
"source": [
"versioninfo()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 準備"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"printarr (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Test\n",
"using LinearAlgebra\n",
"using RandomMatrices\n",
"printarr(arr) = Base.print_array(IOContext(stdout, :compact => true), arr)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## エルミート行列 $A$ を生成\n",
"[JuliaMath/RandomMatrices.jl](https://github.com/JuliaMath/RandomMatrices.jl) パッケージの `GaussianHermite` を用いてランダムなエルミート行列を生成する．"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"GaussianHermite\\{β\\} represents a Gaussian Hermite ensemble with parameter β.\n",
"\n",
"Wigner\\{β\\} is a synonym.\n",
"\n",
"Example of usage:\n",
"\n",
"\\begin{verbatim}\n",
"β = 2 #β = 1, 2, 4 generates real, complex and quaternionic matrices respectively.\n",
"d = Wigner{β} #same as GaussianHermite{β}\n",
"\n",
"n = 20 #Generate square matrices of this size\n",
"\n",
"S = rand(d, n) #Generate a 20x20 symmetric Wigner matrix\n",
"T = tridrand(d, n) #Generate the symmetric tridiagonal form\n",
"v = eigvalrand(d, n) #Generate a sample of eigenvalues\n",
"\\end{verbatim}\n"
],
"text/markdown": [
"GaussianHermite{β} represents a Gaussian Hermite ensemble with parameter β.\n",
"\n",
"Wigner{β} is a synonym.\n",
"\n",
"Example of usage:\n",
"\n",
"```\n",
"β = 2 #β = 1, 2, 4 generates real, complex and quaternionic matrices respectively.\n",
"d = Wigner{β} #same as GaussianHermite{β}\n",
"\n",
"n = 20 #Generate square matrices of this size\n",
"\n",
"S = rand(d, n) #Generate a 20x20 symmetric Wigner matrix\n",
"T = tridrand(d, n) #Generate the symmetric tridiagonal form\n",
"v = eigvalrand(d, n) #Generate a sample of eigenvalues\n",
"```\n"
],
"text/plain": [
" GaussianHermite{β} represents a Gaussian Hermite ensemble with parameter β.\n",
"\n",
" Wigner{β} is a synonym.\n",
"\n",
" Example of usage:\n",
"\n",
"\u001b[36m β = 2 #β = 1, 2, 4 generates real, complex and quaternionic matrices respectively.\u001b[39m\n",
"\u001b[36m d = Wigner{β} #same as GaussianHermite{β}\u001b[39m\n",
"\u001b[36m \u001b[39m\n",
"\u001b[36m n = 20 #Generate square matrices of this size\u001b[39m\n",
"\u001b[36m \u001b[39m\n",
"\u001b[36m S = rand(d, n) #Generate a 20x20 symmetric Wigner matrix\u001b[39m\n",
"\u001b[36m T = tridrand(d, n) #Generate the symmetric tridiagonal form\u001b[39m\n",
"\u001b[36m v = eigvalrand(d, n) #Generate a sample of eigenvalues\u001b[39m"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@doc GaussianHermite"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.0903755+0.0im 0.416189+0.229176im 0.645621+0.182867im 0.267279+0.259613im 0.138913-0.06654im \n",
" 0.416189-0.229176im 0.296458+0.0im 0.227503+0.36022im 0.289634-0.414161im -0.0693944+0.370095im \n",
" 0.645621-0.182867im 0.227503-0.36022im 0.31481+0.0im -0.10495-0.13137im -0.202772-0.396498im \n",
" 0.267279-0.259613im 0.289634+0.414161im -0.10495+0.13137im 0.129101+0.0im -0.0765937+0.00273656im\n",
" 0.138913+0.06654im -0.0693944-0.370095im -0.202772+0.396498im -0.0765937-0.00273656im -0.585852+0.0im "
]
}
],
"source": [
"N = 5\n",
"A = rand(GaussianHermite(2), N)\n",
"\n",
"printarr(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## $A$ の固有値, 固有ベクトルを求める\n",
"固有値及び固有ベクトルは Standard Library である [LinearAlgebra](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/) の `eigvals` 関数 と `eigvecs` 関数 で求められる．また，`eigen` 関数でも取得することが出来る．"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"固有値(eigvals)\n",
" -1.10838 \n",
" -0.605808 \n",
" -0.0674164\n",
" 0.696097 \n",
" 1.3304 \n",
"固有ベクトル(eigvecs)\n",
" -0.37486-0.0869043im 0.555318+0.125579im 0.203765+0.267967im -0.165965-0.344795im 0.20392-0.475858im\n",
" 0.185334-0.3305im -0.26207-0.282287im -0.103122+0.433395im 0.0493836+0.430492im 0.0775921-0.561967im\n",
" 0.345016+0.288291im -0.25803+0.162003im -0.301826-0.152449im 0.0476126-0.535074im -0.199351-0.512226im\n",
" 0.0488565-0.0582875im -0.272034+0.477345im 0.374051-0.403587im -0.451654+0.280989im 0.226014-0.235649im\n",
" 0.707401+0.0im 0.364418+0.0im 0.520646+0.0im 0.30707-0.0im 0.0376719-0.0im \n",
"\n",
"または，\n",
"\n",
"固有値(eigen)\n",
" -1.10838 \n",
" -0.605808 \n",
" -0.0674164\n",
" 0.696097 \n",
" 1.3304 \n",
"固有ベクトル(eigen)\n",
" -0.37486-0.0869043im 0.555318+0.125579im 0.203765+0.267967im -0.165965-0.344795im 0.20392-0.475858im\n",
" 0.185334-0.3305im -0.26207-0.282287im -0.103122+0.433395im 0.0493836+0.430492im 0.0775921-0.561967im\n",
" 0.345016+0.288291im -0.25803+0.162003im -0.301826-0.152449im 0.0476126-0.535074im -0.199351-0.512226im\n",
" 0.0488565-0.0582875im -0.272034+0.477345im 0.374051-0.403587im -0.451654+0.280989im 0.226014-0.235649im\n",
" 0.707401+0.0im 0.364418+0.0im 0.520646+0.0im 0.30707-0.0im 0.0376719-0.0im "
]
}
],
"source": [
"println(\"固有値(eigvals)\")\n",
"λ = eigvals(A)\n",
"printarr(λ)\n",
"println(\"\\n固有ベクトル(eigvecs)\")\n",
"v = eigvecs(A)\n",
"printarr(v)\n",
"\n",
"println(\"\\n\\nまたは，\\n\\n固有値(eigen)\")\n",
"F = eigen(A)\n",
"printarr(F.values)\n",
"println(\"\\n固有ベクトル(eigen)\")\n",
"printarr(F.vectors)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## $A$ の主小行列 $M_j$\n",
"\n",
"```julia\n",
"A[1:N .!= j, 1:N .!= j]\n",
"```\n",
"または\n",
"```julia\n",
"A[setdiff(1:N, j), setdiff(1:N, j)]\n",
"```\n",
"のように記述することで，$M_j$ を表現することができる．"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"j = 1\n",
" 0.296458+0.0im 0.227503+0.36022im 0.289634-0.414161im -0.0693944+0.370095im \n",
" 0.227503-0.36022im 0.31481+0.0im -0.10495-0.13137im -0.202772-0.396498im \n",
" 0.289634+0.414161im -0.10495+0.13137im 0.129101+0.0im -0.0765937+0.00273656im\n",
" -0.0693944-0.370095im -0.202772+0.396498im -0.0765937-0.00273656im -0.585852+0.0im \n",
"\n",
"j = 2\n",
" 0.0903755+0.0im 0.645621+0.182867im 0.267279+0.259613im 0.138913-0.06654im \n",
" 0.645621-0.182867im 0.31481+0.0im -0.10495-0.13137im -0.202772-0.396498im \n",
" 0.267279-0.259613im -0.10495+0.13137im 0.129101+0.0im -0.0765937+0.00273656im\n",
" 0.138913+0.06654im -0.202772+0.396498im -0.0765937-0.00273656im -0.585852+0.0im \n",
"\n",
"j = 3\n",
" 0.0903755+0.0im 0.416189+0.229176im 0.267279+0.259613im 0.138913-0.06654im \n",
" 0.416189-0.229176im 0.296458+0.0im 0.289634-0.414161im -0.0693944+0.370095im \n",
" 0.267279-0.259613im 0.289634+0.414161im 0.129101+0.0im -0.0765937+0.00273656im\n",
" 0.138913+0.06654im -0.0693944-0.370095im -0.0765937-0.00273656im -0.585852+0.0im \n",
"\n",
"j = 4\n",
" 0.0903755+0.0im 0.416189+0.229176im 0.645621+0.182867im 0.138913-0.06654im \n",
" 0.416189-0.229176im 0.296458+0.0im 0.227503+0.36022im -0.0693944+0.370095im\n",
" 0.645621-0.182867im 0.227503-0.36022im 0.31481+0.0im -0.202772-0.396498im\n",
" 0.138913+0.06654im -0.0693944-0.370095im -0.202772+0.396498im -0.585852+0.0im \n",
"\n",
"j = 5\n",
" 0.0903755+0.0im 0.416189+0.229176im 0.645621+0.182867im 0.267279+0.259613im\n",
" 0.416189-0.229176im 0.296458+0.0im 0.227503+0.36022im 0.289634-0.414161im\n",
" 0.645621-0.182867im 0.227503-0.36022im 0.31481+0.0im -0.10495-0.13137im \n",
" 0.267279-0.259613im 0.289634+0.414161im -0.10495+0.13137im 0.129101+0.0im \n",
"\n"
]
}
],
"source": [
"M = zeros(Complex{Float64}, N-1, N-1, N)\n",
"for j = 1:N\n",
" M[:,:,j] = A[1:N .!= j, 1:N .!= j]\n",
" println(\"j = $j\")\n",
" printarr(M[:,:,j])\n",
" println(\"\\n\")\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Lemma 2 (再掲)\n",
"$$\n",
"|v_{i,j}|^2 \\prod_{k = 1; k \\neq i}^n{(\\lambda_i(A) - \\lambda_k(A))} = \\prod_{k = 1}^{n-1}{(\\lambda_i(A) - \\lambda_k(M_j))}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 左辺"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.340903 -0.221091 0.0677844 -0.166597 1.12213 \n",
" 0.330558 -0.101196 0.118709 -0.213624 1.34738 \n",
" 0.465399 -0.0633124 0.0683906 -0.328318 1.26486 \n",
" 0.0133172 -0.205889 0.181113 -0.321919 0.446352 \n",
" 1.15209 -0.0905785 0.162138 -0.10728 0.00594161"
]
}
],
"source": [
"lhs = abs2.(v) .* [prod(λ[i] - λ[k] for k = 1:N if k != i) for j = 1:N, i = 1:N]\n",
"printarr(lhs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 右辺"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.340903 -0.221091 0.0677844 -0.166597 1.12213 \n",
" 0.330558 -0.101196 0.118709 -0.213624 1.34738 \n",
" 0.465399 -0.0633124 0.0683906 -0.328318 1.26486 \n",
" 0.0133172 -0.205889 0.181113 -0.321919 0.446352 \n",
" 1.15209 -0.0905785 0.162138 -0.10728 0.00594161"
]
}
],
"source": [
"rhs = [prod(λ[i] - eigvals(M[:,:,j])[k] for k = 1:N-1) for j = 1:N, i = 1:N]\n",
"printarr(rhs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 両辺比較"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@test lhs ≈ rhs"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.2.0",
"language": "julia",
"name": "julia-1.2"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.2.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment