行列対角化の計算方法を紹介~LaTeXをMarkdownへ書き替えてみる~

今回の記事では$\LaTeX$ の文章をMarkdown記法への書き換えを試行してみました.試行で使用した文章は,数値計算の手法（ランチョス法）の文章をベースとしています.

# 本記事の目的
本記事は過去に私が執筆した修士論文[1]から計算方法の第3章のみを抜き出して記載しています.計算手法であれば発表時期が古くても今現在記事にするにもニーズはあるのではと考えたこと,投稿論文等では数行等で終わる計算手法は修士論文の中では細かく記載をしているので,ご紹介するメリットがあると考えたのが本記事記載の目的の一つです.そして,Markdown形式（$\LaTeX$）で行列式等を作るのは面倒であり,時間がかかります.数式のソースを載せるのは参考になるのではと考えたのが目的の二つめです.以下、論文調で読みにくい文体ですが,ご容赦ください.

# 第3章　計算方法
この章では,計算手法を述べる.
## 3.1 ランチョス法
銅酸化物高温超伝導体は代表的な強相関電子系である．すなわち，電子のクーロン相互作用と運動エネルギーを正確に考慮していく必要があり，その取り扱いは簡単なものではない．高温超伝導体の理論的研究は様々な数値的手法によってなされてきた．手法の一つとして高温超伝導体のCuO$_2$面の電子状態を数値的に調べる厳密対角化法がある．この手法は有限サイズのクラスターにおいて，$t\hbox{-}J$モデルやハバードモデルによって与えられたハミルトニアンを厳密に対角化する方法であり，多体効果を正確に取り入れられるという大きな利点を持っている．そして，固有状態や固有エネルギーが求まるだけでなく，一粒子励起スペクトルなどの動的な物理量を求められることも特徴の一つである．数値計算によって一粒子励起スペクトルを求めて，角度分解光電子分光(ARPES)実験のスペクトルと比較することで，実験との整合性を常に意識しながら理論を構築することができる．


厳密対角化法はハミルトニアンを三重対角化して，三重対角行列を対角化する手法である．この三重対角化を行うときに，状態をすべて求めていくハウスホルダー法がある．この方法ではコンピュータのメモリの制限より，高々$10^3$次元程度の行列の対角化しか実行できない．よって，大規模疎行列($10^7$次元)を計算するためにはランチョス法(Lancz$\ddot{\rm{o}}$s Method)[2]を用いて対角化を実行する．ランチョス法には，基底状態と低エネルギー励起状態のみが精度良く求められるという制限があるが，本研究で興味があるのは低励起状態なので，ランチョス法を用いて計算することは非常に有効である．以下に厳密対角化の手順を示す．

## 3.1.1 ランチョスステップ

任意の$N \times N$ のエルミート行列 $\bf{H}$ (ハミルトニアン行列)に対し，任意のユニタリー行列 $\bf{U}$ が存在し，ハミルトニアン行列を三重対角化できる．三重対角化された行列を $\bf{T}$ とすると，次のようにかける．
$\bf{U^{\dagger}HU} = \bf{T} ，\quad \bf{U^{\dagger}U}=\bf{1}.　\tag{1} $　
ここで，
```math
\begin{eqnarray}
\bf{T}= \left(
\begin{array}{ccccccc}
\it{a}_{\rm{1}} & b_{2} &       &        &         &         &\\
b_{2} & a_{2} & b_{3} &        &         &    0    &      \\
      & b_{3} & a_{3} & b_{4}  &         &         &      \\
      &       &       & \ddots &         &         &      \\
      &       &       &        & a_{N-2} & b_{N-1} &      \\ 
      &       &  0    &        & b_{N-1} & a_{N-1} & b_{N} \\
      &       &       &        &         & b_{N}   & a_{N}  \\
\end{array}
\right).\tag{2}
\end{eqnarray}
```
```math
\begin{equation}
\bf{U}=\{ \bf{u_{\rm{1}}，u_{\rm{2}}，u_{\rm{3}}，\ldots，u_{\it{N}}} \}
，\quad \bf{u_{\it{i}}\cdot u_{\it{j}}}=\delta_{\it{ij}}.\tag{3}
\end{equation}
```
のように書くことができる．また，
```math
\begin{equation}
\bf{HU}=\bf{UT}. \tag{4}
\end{equation}
```
となることから，次の漸化式を得る．
```math
\begin{equation} \left\{ 
\begin{array}{lcl}
\bf{Hu}_{\rm{1}}&  =  & a_{\rm{1}}\bf{u}_{\rm{1}}+\it{b_{\rm{2}}}\bf{u}_{\rm{2}}, \\
{\bf Hu}_{i}&   = & b_{i}\bf{u}_{\it{i}-\rm{1}} +\it{a_{i}}\bf{u}_{\it{i}}+\it{b_{i+\rm{1}}}\bf{u}_{\it{i}+\rm{1}}, \tag{5} \\
\bf{Hu}_{\it{N}}&   = & a_{N}\bf{u}_{\it{N}}+\it{b_{N}}\bf{u}_{\it{N}-\rm{1}}.
\end{array} \right.
\end{equation}
```
ここで漸化式を解くために，任意の初期ベクトルとして ${\bf u}_{1} = \bf{u}$ をとり，三重対角行列 $\bf{T}$ の行列要素を決定していく(初期ベクトル $\bf{u}$ は固有ベクトルと直交してはならない)．

1.  $b_{2} \bf  u_{\rm{2}} =\bf{Hu}-\it{a_{\rm{1}}}\bf{u}$を決める．この時$a_{1}$$=\bf  u_{1}^{\dagger}Hu_{1} $であり，$ b_2 $は規格化で決める．  
2. $b_{i+1}\bf u_{\it{i}\rm{+1}}=\bf Hu_{\it{i}}-$ $b_{i} \bf u_{\it{i} - \rm{1}}$ $-a_{i} \bf  u_{\it{i}}$を決める．この時，$a_{i}=\bf  u_{\it{i}}^{\dagger}(Hu_{\it{i}}$ $-b_i \bf  u_{\it{i}-\rm{1}})=\bf u_{\it{i}}^{\dagger}Hu_{\it{i}}$であり，$b_{i}$ は規格化で決める.

 以下，2.の手順を $N-1$ 回繰り返せば三重対角行列 $\bf{T}$ が原理的には得られる．しかし，実際 $N-1$ 回繰り返しても数値誤差が積み重なるので, $\bf U$ のユニタリー性が崩れ，固有値を精度良く求めることはできない．むしろ，ランチョス法の利点は $N-1$ 回よりもはるかに小さな$100$回程度の反復で三重対角行列の途中までの部分行列の最低の固有値が，元の行列の最低の固有値に収束することにある．

## 3.1.2 固有エネルギーと固有状態
前節より，行列が三重対角化されたのでこれを対角化することによりハミルトニアン行列の固有値が求まる．固有値は，二分法(Bisection method)やQL法より求める．こうして固有値(固有エネルギー)や三重対角行列の固有状態が求まるので，次にハミルトニアンの固有状態を求めることを考える．
反復過程の回数(三重対角行列の次元)を$n$，三重対角行列の $\nu$ 番目の固有値を$\lambda_{\nu}$ そしてその固有状態を$\bf{x}_{\nu}$とする．
```math
\begin{equation}
\bf{U}^{\dagger}\bf{HUx}_{\nu}=\lambda_{\nu}\bf{x}_{\nu}. \tag{6}
\end{equation}
```
```math
\begin{equation}
\bf{HUx}_{\nu}=\lambda_{\nu}\bf{Ux}_{\nu}. \tag{7}
\end{equation}
```
よって，ハミルトニアンの固有状態は次のようになる．
```math
\begin{equation}
\bf{Ux}_{\nu}=[\bf{u}_{\rm{1}}，\bf{u}_{\rm{2}}，\bf{u}_{\rm{3}}，\ldots，\bf{u}_{\rm{n}}]
 \cdot
 \left(
 \begin{array}{c}
 \it{x}_{1} \\
 \vdots  \\
 x_{n}
 \end{array}
\right)_{\nu}. \tag{8}
\end{equation}
```

## 3.1.3 物理量の計算
基底状態の固有エネルギー，固有状態(波動関数)が求まったので，絶対零度の静的な物理量は求まり，次のように書くことができる．
```math
\begin{equation}
\langle \hat{\mathcal{O}} \rangle =\langle \phi_{0}| \hat{\mathcal{O}} | \phi_{0} \rangle. \tag{9}
\end{equation}
```
ここで，$\hat{\mathcal{O}}$は演算子であり，$| \phi_{0} \rangle$はハミルトニアンの基底状態の固有状態である．
また，ランチョス法を応用して動的な物理量も考えることができる．考えるものは次のような周波数(エネルギー)$\omega$に依存したスペクトル関数である．

```math
\begin{equation}
 I(\omega) = - \frac{1}{\pi} {\rm Im} G(\omega+i \epsilon).　\tag{10}
\end{equation}
```

```math
\begin{equation}
G(\omega)= \langle \phi_{0}| \hat{\mathcal{O}}^{\dagger}
\frac{1}{\omega+E_{0}-H+i \epsilon}
 \hat{\mathcal{O}}| \phi_{0} \rangle.　\tag{11}
\end{equation}
```
ここで，$|\phi_{0} \rangle $はハミルトニアン $H$ の基底状態で固有エネルギー $E_0$ を持つ．$\epsilon$は小さな正の実数である．(11)式のグリーン関数はランチョス法を用いて求められる．
この時，初期ベクトル $u_{1}$ として $\hat{\mathcal{O}}|\phi_0\rangle /\sqrt{\langle \phi_{0}|\hat{\mathcal{O}}^{\dagger}\hat{\mathcal{O}}|\phi_{0} \rangle}$ を採用しランチョス法によって三重対角化を実行する．
```math
\begin{equation}
\langle \phi_0 | \hat{\mathcal{O}}^{\dagger} (z\bf{E} - \bf{H})^{\rm{-1}} \hat{\mathcal{O}} | \phi_0 \rangle = \langle \phi_0 | \hat{\mathcal{O}}^{\dagger} \bf{UU}^{\dagger}(\it{z}\bf{E} - \bf{H})^{\rm{-1}} \bf{UU}^{\dagger}\hat{\mathcal{O}}  | \phi_0 \rangle. \tag{12}
\end{equation}
```
```math
\begin{eqnarray}
\bf{U}^{\dagger}(\it{z}\bf{E}-\bf{H})^{\rm{-1}}\bf{U} 
&=& (\bf{U}^{\rm{-1}}\it{z}\bf{U}-\bf{U}^{\rm{-1}}\bf{H}\bf{U})^{\rm{-1}}  \\
&=& (\it{z}\bf{E}-\bf{T})^{\rm{-1}}. \tag{13}
\end{eqnarray}
```
すなわち，(11)式の値は，(14)式の逆行列の(1，1)成分に対応する．

```math
\begin{equation}
z\bf{E}-\bf{T}= \left(
\begin{array}{ccccc}
\it{z}- a_{\rm{1}} &  -b_{2}  &        &            &        \\
-b_{2}  & z-a_{2}  &        &       0    &        \\
        & - b_{3}  & 　　　 &            &        \\
        &          & \ddots &            &        \\
        &          &        &    　　 　 &        \\ 
        &     0    &        &   z-a_{N-1} & -b_{N} \\
        &          &        &   -b_{N}   & z-a_{N} \\
\end{array}
\right). \tag{14}
\end{equation}
```

すなわち，クラメールの定理(Cramer's rule)を用いると，次のようになる．

```math
\begin{equation}
(z\bf{E-T})^{\rm{-1}}\rm{の(1，1)成分}= \frac{\det D_{\rm{1}}}{\det (\it{z}\bf{E-T})}. \tag{15}
\end{equation}
```

ここで，$ \bf{}D_{n} $とは $z\bf{E}-\bf{T}$ の $n$ 行までと$n$列までを除いた部分の行列のことである．また，$\bf{D}_{\rm{n}}$には漸化式として，次の式が成り立つ.

```math
\begin{equation}
\frac{\det \bf{D}_{\rm{n-1}}}{\det \bf{D}_{\rm{n}}} = z - a_{n} - \cfrac{b_{n+1}^{2}}{\cfrac{\det \bf{D}_{\rm{n}}}{\det \bf{D}_{\rm{n+1}}}}. \tag{16}
\end{equation}
```
よって，次のように
```math
\begin{eqnarray}
G(z) & = &\cfrac{\langle \phi_{0}| \hat{\mathcal{O}}^{\dagger}\hat{\mathcal{O}}| \phi_{0} \rangle}{\cfrac{\det (z\bf{E-T})}{\det {\bf D}_{1}}}  \nonumber \\　
& = &\cfrac{\langle \phi_{0}| \hat{\mathcal{O}}^{\dagger}\hat{\mathcal{O}}| \phi_{0} \rangle}{z-a_1
-\cfrac{b_{2}^{2}}{z-a_2 
-\cfrac{b_{3}^{2}}{z-a_3
-\ldots}}}. \tag{17}
\end{eqnarray}
```
連分数展開される．この時適当な$\epsilon$を与えれば，(10)式のスペクトル$I(\omega)$は連続スペクトルになる．

よく知られた関係式($\epsilon \to 0$)
```math
\begin{equation}
\frac{1}{x+i \epsilon}=\rm{P} \left( \frac{1}{x} \right) -i \pi \delta(x).
\end{equation}
```
を用いると，(10)式は次のようにかける．
```math
\begin{equation}
I (\omega) = \sum_{\nu} | \langle \phi_{\nu} | \hat{\mathcal{O}}| \phi_{0} \rangle |^{2} \delta(\omega - (E_{\nu} - E_{0})). \tag{18}
\end{equation}
```
この時，スペクトルは$\delta$関数の集まりであり，$ | \phi_{\nu} \rangle $は$\nu$番目の励起状態であり，$E_{\nu}$はその固有エネルギーで，$\langle \phi_{\nu} | \hat{\mathcal{O}}| \phi_{0} \rangle$はその固有状態の第一成分を表す．
また，$\delta$関数は有限の$\epsilon$を与えることで，ブロードニングされる．
```math
\begin{equation}
\delta(x) \to \frac{1}{\pi}\frac{\epsilon}{x^2+\epsilon^{2}}. \tag{19} 
\end{equation}
```
式(19)のように，$\epsilon$の半値幅でローレンツ型になる．

(以下,第4章に続きますが,本稿では省略します)


# 参考文献
[1] 高橋亮, 修士論文(東北大学)(2000).
[2] C. Lancz$\ddot{\rm{o}}$s, J. Res. Natl .Bur. Std. $\bf{45}$,255(1950).


# おわりに
 今回の記事では過去に$\LaTeX$で作成した文章のMarkdown形式化を試してみました. 実際は, 細かなところで文法が異なり修正が必要なことに気がつきました. 特に文章中に数式を入れる際にBold（\bf）とすると表記が崩れるところは細かく調整しています.また,$\LaTeX$では数式の下にスクロールが出ることがないのに自動で表示されてしまうところは修正できませんでした.解決方法があればコメント頂けると助かります.
以下はMarkdownに変更するにあたり削除対応等した主なコマンドです.
・\cite
・\subsection
・\label
・\begin{enumerate}　\itme \end{}{enumerate}
・\ref


## 参考：(2)式

```
\begin{eqnarray}
\bf{T}= \left(
\begin{array}{ccccccc}
\it{a}_{\rm{1}} & b_{2} &       &        &         &         &\\
b_{2} & a_{2} & b_{3} &        &         &    0    &      \\
      & b_{3} & a_{3} & b_{4}  &         &         &      \\
      &       &       & \ddots &         &         &      \\
      &       &       &        & a_{N-2} & b_{N-1} &      \\ 
      &       &  0    &        & b_{N-1} & a_{N-1} & b_{N} \\
      &       &       &        &         & b_{N}   & a_{N}  \\
\end{array}
\right).\tag{2}
\end{eqnarray}
```
## 参考：(5)式

```
\begin{equation} \left\{ 
\begin{array}{lcl}
\bf{Hu}_{\rm{1}}&  =  & a_{\rm{1}}\bf{u}_{\rm{1}}+\it{b_{\rm{2}}}\bf{u}_{\rm{2}}, \\
{\bf Hu}_{i}&   = & b_{i}\bf{u}_{\it{i}-\rm{1}} +\it{a_{i}}\bf{u}_{\it{i}}+\it{b_{i+\rm{1}}}\bf{u}_{\it{i}+\rm{1}}, \tag{5} \\
\bf{Hu}_{\it{N}}&   = & a_{N}\bf{u}_{\it{N}}+\it{b_{N}}\bf{u}_{\it{N}-\rm{1}}.
\end{array} \right.
\end{equation}
```
## 参考：(10)式

```
\begin{equation}
 I(\omega) = - \frac{1}{\pi} {\rm Im} G(\omega+i \epsilon).　\tag{10}
\end{equation}
```
## 参考：(11)式

```
\begin{equation}
G(\omega)= \langle \phi_{0}| \hat{\mathcal{O}}^{\dagger}
\frac{1}{\omega+E_{0}-H+i \epsilon}
 \hat{\mathcal{O}}| \phi_{0} \rangle.　\tag{11}
\end{equation}
```
## 参考：(15)式

```
\begin{equation}
(z\bf{E-T})^{\rm{-1}}\rm{の(1，1)成分}= \frac{\det D_{\rm{1}}}{\det (\it{z}\bf{E-T})}. \tag{15}
\end{equation}
```

## 参考：(16)式

```
\begin{equation}
\frac{\det \bf{D}_{\rm{n-1}}}{\det \bf{D}_{\rm{n}}} = z - a_{n} - \cfrac{b_{n+1}^{2}}{\cfrac{\det \bf{D}_{\rm{n}}}{\det \bf{D}_{\rm{n+1}}}}. \tag{16}
\end{equation}
```


## 参考：(17)式
```
\begin{eqnarray}
G(z) & = &\cfrac{\langle \phi_{0}| \hat{\mathcal{O}}^{\dagger}\hat{\mathcal{O}}| \phi_{0} \rangle}{\cfrac{\det (z\bf{E-T})}{\det {\bf D}_{1}}}  \nonumber \\　
& = &\cfrac{\langle \phi_{0}| \hat{\mathcal{O}}^{\dagger}\hat{\mathcal{O}}| \phi_{0} \rangle}{z-a_1
-\cfrac{b_{2}^{2}}{z-a_2 
-\cfrac{b_{3}^{2}}{z-a_3
-\ldots}}}. \tag{17}
\end{eqnarray}
```
## 参考：(18)式

```
\begin{equation}
I (\omega) = \sum_{\nu} | \langle \phi_{\nu} | \hat{\mathcal{O}}| \phi_{0} \rangle |^{2} \delta(\omega - (E_{\nu} - E_{0})). \tag{18}
\end{equation}
```


以上