Lanczos 算法

                     

贡献者: addis

预备知识 氢原子球坐标薛定谔方程数值解,施密特正交归一化(向量与矩阵)

1. 背景介绍

  1在计算含时薛定谔方程时,若已知某时刻 $t$ 的波函数 $\Psi( \boldsymbol{\mathbf{r}} ,t)$,要求 $\Psi( \boldsymbol{\mathbf{r}} ,t+\Delta t)$,通常使用传播算符

\begin{equation} U(\Delta t) = \exp\left(-\frac{ \mathrm{i} }{\hbar} \int_{t}^{t+\Delta t} \,\mathrm{d}{t} \boldsymbol\cdot H(t)\right) \approx \exp\left(- \mathrm{i} H(t+\Delta t/2) \Delta t\right) ~, \end{equation}
其中 $H$ 为哈密顿算符。若 $H$ 不含时,约等号变为等号
\begin{equation} \Psi( \boldsymbol{\mathbf{r}} , t+\Delta t) = \exp\left(- \mathrm{i} H\Delta t/\hbar\right) \Psi( \boldsymbol{\mathbf{r}} , t)~. \end{equation}
通常情况下,用有限个正交归一基底 $\Psi_0,\dots\Psi_{K-1}$ 近似展开 $\Psi( \boldsymbol{\mathbf{r}} , t+\Delta t)$,这时,$H$ 可以表示成矩阵的形式。
\begin{equation} H_{ij} = \left\langle \Psi_i \right\rvert H \left\lvert \Psi_j \right\rangle ~. \end{equation}
这样,幺正算符 $ \exp\left(- \mathrm{i} H\Delta t/\hbar\right) $ 也可以表示成矩阵。根据定义
\begin{equation} \exp\left(- \mathrm{i} H\Delta t/\hbar\right) = 1+ (- \mathrm{i} H\Delta t/\hbar) + \frac{1}{2!} (- \mathrm{i} H\Delta t/\hbar)^2\dots~ \end{equation}
此时若把 $H$ 对角化(求解本征方程),得到本征矢列矩阵 $ \boldsymbol{\mathbf{P}} $(单位正交阵)以及本征值矩阵 $ \boldsymbol{\mathbf{\Lambda}} $(对角矩阵),则可进行基底变换变到 $H$ 矩阵的本征空间求解上式
\begin{equation} \begin{aligned} \exp\left(- \mathrm{i} H\Delta t/\hbar\right) &= \boldsymbol{\mathbf{P}} \times \rm{diag}( \mathrm{e} ^{- \mathrm{i} E_1 t/\hbar}, \mathrm{e} ^{- \mathrm{i} E_2 t/\hbar}\dots \mathrm{e} ^{- \mathrm{i} E_N t/\hbar})\times \boldsymbol{\mathbf{P}} ^{-1}\\ &= \sum_{j=0}^{K-1} P_{ij} \mathrm{e} ^{- \mathrm{i} E_j t/\hbar} P_{jk}^{-1} = \sum_{j=0}^{K-1} P_{ij}P_{kj} \mathrm{e} ^{- \mathrm{i} E_j t/\hbar} ~. \end{aligned} \end{equation}
要说明的是,这种算法的误差(除数值计算误差)来源于两个方面,第一是用 $ \exp\left( \mathrm{i} H\Delta t/\hbar\right) $ 代替 $ \exp\left(-\frac{ \mathrm{i} }{\hbar} \int_{t}^{t+\Delta t} \,\mathrm{d}{t} \cdot H(t)\right) $,但如果不含时就没有该误差,第二是有限数量的基底不具有完备性,这个误差可以随着基底数量增加而减小。

2. Lanczos 算法

   在以上的算法中,选取施密特正交归一化 的 Krilov 基底作为基底。令 $\Psi_0 \equiv \Psi( \boldsymbol{\mathbf{r}} , t)$, $K$ 阶的 Krilov 基底为

\begin{equation} \left\{ \left\lvert \Psi_0 \right\rangle , H \left\lvert \Psi_0 \right\rangle , H^2 \left\lvert \Psi_0 \right\rangle \dots H^{K-1} \left\lvert \Psi_0 \right\rangle \right\} ~. \end{equation}
我们把正交归一化后的基底记为
\begin{equation} \left\{ \left\lvert \Psi_0 \right\rangle , \left\lvert \Psi_1 \right\rangle , \left\lvert \Psi_2 \right\rangle \dots \left\lvert \Psi_{K-1} \right\rangle \right\} ~. \end{equation}
这样,用其展开哈密顿算符
\begin{equation} H_{ij}= \left\langle \Psi_i \middle| H \middle| \Psi_j \right\rangle ~. \end{equation}
在算法上,本来我们可以按部就班地按照以上步骤做,然而由于 Krilov 基底的性质,可以通过一些定理(见词条最后)大大减少计算量。正交归一化步骤如下

   【1】如果 $\Psi( \boldsymbol{\mathbf{r}} ,t)$ 没有归一化,将其进行归一化

\begin{equation} \left\lvert \Psi_0 \right\rangle =\Psi( \boldsymbol{\mathbf{r}} ,t)/\sqrt{ \left\langle{\Psi( \boldsymbol{\mathbf{r}} ,t)}\middle| \Psi( \boldsymbol{\mathbf{r}} ,t) \right\rangle }~. \end{equation}

   【2】把 $\Psi_1$ 进行施密特正交化,令

\begin{equation} \lvert{\tilde\Psi_1}\rangle = (1 - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert )H \left\lvert \Psi_0 \right\rangle ~, \qquad \left\lvert \Psi_1 \right\rangle = \lvert{\tilde\Psi_1}\rangle /\sqrt{ \langle{\tilde\Psi_1}|{\tilde\Psi_1}\rangle }~. \end{equation}

   【3】现在起我们令 $\beta_j = \sqrt{ \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle }$。

定理 1 

   对以上定义的基底

\begin{equation} \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \left\lvert \Psi_{j - 1} \right\rangle \left\langle \Psi_{j - 1} \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) \ H^{j + 1} \left\lvert \Psi_0 \right\rangle ~ \end{equation}
\begin{equation} \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \left\lvert \Psi_{j - 1} \right\rangle \left\langle \Psi_{j - 1} \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) \ H \left\lvert \Psi_j \right\rangle ~. \end{equation}
共线。

   故现在起采用后者进行正交化。

定理 2 

\begin{equation} \left\langle \Psi_{j-1} \right\rvert H \left\lvert \Psi_j \right\rangle = \sqrt{ \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle } = \beta_j~. \end{equation}

   【4】把 $H^2 \left\lvert \Psi_0 \right\rangle $ 进行施密特正交化,令

\begin{equation} \alpha_i = \left\langle \Psi_i \middle| H \middle| \Psi_i \right\rangle ~. \end{equation}
根据以上两个定理,正交化结果也可以写成
\begin{equation} \begin{aligned} \lvert{\tilde\Psi_2}\rangle &= \left(1 - \left\lvert \Psi_1 \right\rangle \left\langle \Psi_1 \right\rvert - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H \left\lvert \Psi_1 \right\rangle \\ &= H \left\lvert \Psi_1 \right\rangle - \alpha_1 \left\lvert \Psi_1 \right\rangle - \beta_1 \left\lvert \Psi_0 \right\rangle ~.\\ \end{aligned} \end{equation}

   【5】再次归一化 $ \left\lvert \Psi_2 \right\rangle = \lvert{\tilde\Psi_2}\rangle / \sqrt{ \langle{\tilde\Psi_2}|{\tilde\Psi_2}\rangle } = \lvert{\tilde\Psi_2}\rangle /\beta_2~.$

定理 3 定理 3

\begin{equation} \begin{aligned} \lvert{\tilde\Psi_{j+1}}\rangle &= \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H \left\lvert \Psi_j \right\rangle \\ &= \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \left\lvert \Psi_{j-1} \right\rangle \left\langle \Psi_{j-1} \right\rvert \right) H \left\lvert \Psi_j \right\rangle \\ &= H \left\lvert \Psi_j \right\rangle - \alpha_j \left\lvert \Psi_j \right\rangle - \beta_j \left\lvert \Psi_{j-1} \right\rangle ~, \end{aligned} \end{equation}

   这就是最简洁的正交化公式。比起最原始的正交化

\begin{equation} \lvert{\tilde\Psi_{j+1}}\rangle = \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H^{j+1} \left\lvert \Psi_0 \right\rangle ~. \end{equation}
定理 3 只用到了一个矩阵-矢量乘法,两个矢量数乘和两个矢量减法。

   接下来只要对 $j = 2,4,5...K - 1$ 不断重复定理 3 中的正交化和归一化,就可以得到所有正交归一基底 $ \left\lvert \Psi_j \right\rangle $。

   现在我们来求该基底下的 $ \boldsymbol{\mathbf{H}} $ 矩阵。根据 $\alpha_i$ 的定义以及定理 2 可以发现它们分别是矩阵的主对角元和副对角元。

定理 4

   在正交归一基底 $ \left\lvert \Psi_j \right\rangle $ 中,矩阵 $ \boldsymbol{\mathbf{H}} $ 是对称三对角矩阵。所以我们已经顺便求出了所有的矩阵元!

\begin{equation} \boldsymbol{\mathbf{H}} = \begin{pmatrix} \alpha_0 & \beta_1 &&& \\ \beta_1 & \alpha_1 & \beta_2 && \\ & \beta_2 & \alpha_2 & \ddots& \\ & & \ddots& \ddots & \beta_{K-1} \\ &&&\beta_{K-1} & \alpha_{K-1} \end{pmatrix}~.\end{equation}
该矩阵具有维数小,易求本征问题的优势。令
\begin{equation} \boldsymbol{\mathbf{H}} = \boldsymbol{\mathbf{P}} \boldsymbol{\mathbf{\Lambda}} \boldsymbol{\mathbf{P}} ^{-1}~, \end{equation}
那么
\begin{equation} \mathrm{e} ^{- \mathrm{i} \boldsymbol{\mathbf{H}} \Delta t} = \boldsymbol{\mathbf{P}} \exp\left(- \mathrm{i} \boldsymbol{\mathbf{\Lambda}} \Delta t\right) \boldsymbol{\mathbf{P}} ^\dagger ~. \end{equation}
注意在 $ \left\lvert \Psi_i \right\rangle $ 基底下,$\Psi( \boldsymbol{\mathbf{r}} , t)$ 就是 $(1, 0, 0, \dots)$。所以在该基底下
\begin{equation} \boldsymbol{\mathbf{a}} \equiv \mathrm{e} ^{- \mathrm{i} \boldsymbol{\mathbf{H}} \Delta t} \Psi( \boldsymbol{\mathbf{r}} , t)= \boldsymbol{\mathbf{P}} \exp\left(- \mathrm{i} \boldsymbol{\mathbf{\Lambda}} \Delta t\right) ( \boldsymbol{\mathbf{P}} ^\dagger )_{\text{第 1 列}} = \boldsymbol{\mathbf{P}} \exp\left(- \mathrm{i} \boldsymbol{\mathbf{\Lambda}} \Delta t\right) ( \boldsymbol{\mathbf{P}} _{\text{第 1 行}}) ^\dagger ~. \end{equation}
这里的 $ \boldsymbol{\mathbf{a}} $ 其实就是 $ \exp\left(- \mathrm{i} \boldsymbol{\mathbf{H}} \,\mathrm{d}{t} \right) $ 的第一列。

3. 误差分析

   要判断 $K$ 是否足够,只需要计算式 21 中的矢量 $ \boldsymbol{\mathbf{a}} $ 的最后几个矩阵元的平方和是否可以忽略不计,如果可以就说明 $K$ 是足够的。

   第一个被截去的系数 $a_{K+1}$ 的模长可以用下式估计(Feist 论文 4.51)

\begin{equation} \left\lvert a_{K+1} \right\rvert \approx \frac{\Delta t^K}{K!} \prod_i \beta_i~. \end{equation}
推导:把式 18 做 $K-1$ 次幂,左下角的矩阵元就是 $\prod_i \beta_i$。此时如果把式 18 拓展到无限大的矩阵,那么做 $K-1$ 次幂后第一列也是相同的。也就是说,$ \exp\left(- \mathrm{i} \boldsymbol{\mathbf{H}} \,\mathrm{d}{t} \right) $ 的泰勒展开的前 $K$ 项的第一列都是精确的。若考察 $K+1$ 项,则无限大矩阵的结果中第一列最后会多出一个元素 $\beta_1\beta_2\dots\beta_K$(可以借助 Mathematica)。接下来的推导就很简单了。

   对于固定的 $K$,该式可以用于薛定谔方程的变步长演化。Feist 用的 Sort Iterative Lanczos 就是使用 $K$ 等于 "a few",以及较小的 $\Delta t$。

4. 定理 1 证明

   根据施密特正交化的性质,Krilov 基底的前 $j$ 项的与 $ \left\lvert \Psi_0 \right\rangle \dots \left\lvert \Psi_j \right\rangle $ 展开同一空间。所以 $H^j \left\lvert \Psi_0 \right\rangle = c_j \left\lvert \Psi_j \right\rangle +\ldots + c_0 \left\lvert \Psi_0 \right\rangle $,所以(以下所有系数 只是用来表示线性组合,具体值不重要)

\begin{equation} \begin{aligned} H^{j+1} \left\lvert \Psi_0 \right\rangle &= c_j H \left\lvert \Psi_j \right\rangle + H \left(c_{j-1} \left\lvert \Psi_{j-1} \right\rangle \dots + c_0 \left\lvert \Psi_0 \right\rangle \right) \\ & = c_j H \left\lvert \Psi_j \right\rangle + H \left(c'_{j-1} H^{j-1} \left\lvert \Psi_0 \right\rangle + c'_{j-2} H^{j - 2} \left\lvert \Psi_0 \right\rangle \dots + c'_0 \left\lvert \Psi_0 \right\rangle \right) \\ &= c_j H \left\lvert \Psi_j \right\rangle + \left(c'_{j-1} H^j \left\lvert \Psi_0 \right\rangle + c'_{j-2} H^{j-1} \left\lvert \Psi_0 \right\rangle \dots + c'_0 \left\lvert \Psi_0 \right\rangle \right) \\ &= c_j H \left\lvert \Psi_j \right\rangle + \left(c''_{j-1} \left\lvert \Psi_j \right\rangle + c''_{j-2} \left\lvert \Psi_{j-1} \right\rangle + \ldots + c''_0 \left\lvert \Psi_0 \right\rangle \right) ~.\\ \end{aligned} \end{equation}
把 $H^{j+1} \left\lvert \Psi_0 \right\rangle $ 施密特正交归一化
\begin{equation} \begin{aligned} &\phantom{=} \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H^{j+1} \left\lvert \Psi_0 \right\rangle \\ &= c_j \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H \left\lvert \Psi_j \right\rangle \\ &\phantom{=} + \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) \left(c''_{j-1} \left\lvert \Psi_j \right\rangle + c''_{j - 2} \left\lvert \Psi_{j - 1} \right\rangle + \ldots + c''_0 \left\lvert \Psi_0 \right\rangle \right) \\ &= c_j \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H \left\lvert \Psi_j \right\rangle ~. \end{aligned} \end{equation}

5. 定理 2 证明

   要证 $ \left\langle \Psi_{j-1} \right\rvert H \left\lvert \Psi_j \right\rangle = \sqrt{ \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle }$,即证 $ \langle{\Psi_{j - 1}}|{H}|{\tilde \Psi_j}\rangle = \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle $,即证 $ \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle = \langle{\tilde\Psi_j}|{H\Psi_{j-1}}\rangle $ 而

\begin{equation} \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle = \langle{\tilde\Psi_j}|{ \left(1- \left\lvert \Psi_{j-1} \right\rangle \left\langle \Psi_{j-1} \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H}|{\Psi_{j-1}}\rangle ~. \end{equation}
由于上式中 $ \langle{{\tilde \Psi }_j}\rvert \left( \left\lvert \Psi_{j - 1} \right\rangle \left\langle \Psi_{j - 1} \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) = 0~,$
\begin{equation} \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle = \langle{\tilde\Psi_j}|{H}|{\Psi_{j-1}}\rangle ~. \end{equation}

6. 定理 3 证明

   要证明定理 3,即证,对 $n \leqslant j - 2$,有 $ \left\langle \Psi_n \right\rvert H \left\lvert \Psi_j \right\rangle = 0$,$H$ 为厄米算符时即证 $ \left\langle \Psi_j \right\rvert H \left\lvert \Psi_n \right\rangle = 0$,即证对 $m \geqslant j + 2$, $ \left\langle \Psi_m \right\rvert H \left\lvert \Psi_j \right\rangle = 0$。在定理 1 的证明类似,我们知道

\begin{equation} \begin{aligned} H \left\lvert \Psi_j \right\rangle &= H \left(c_j H^j \left\lvert \Psi_0 \right\rangle + \ldots + c_0 \left\lvert \Psi_0 \right\rangle \right) = c_j H^{j+1} \left\lvert \Psi_0 \right\rangle + \ldots + c_0 H \left\lvert \Psi_0 \right\rangle \\ &= c'_j \left\lvert \Psi_{j+1} \right\rangle + \ldots + c'_0 \left\lvert \Psi_0 \right\rangle ~,\\ \end{aligned} \end{equation}
所以对 $m \geqslant j + 2$ 有
\begin{equation} \left\langle \Psi_m \right\rvert H \left\lvert \Psi_j \right\rangle = \left\langle \Psi_m \right\rvert \left(c'_j \left\lvert \Psi_{j+1} \right\rangle + \ldots + c'_0 \left\lvert \Psi_0 \right\rangle \right) = 0~. \end{equation}

7. 定理 4 证明

   首先证明 $ \boldsymbol{\mathbf{H}} $ 是三对角矩阵。即证 $m \geqslant j + 2$ 或 $m \leqslant j - 2$ 时 $ \left\langle \Psi_m \right\rvert H \left\lvert \Psi_j \right\rangle = 0$。 $m \geqslant j + 2$ 的情况在证明 3 中已经证明,只需证明另一种情况。对厄米矩阵,$ \langle{H\Psi_m}|{\Psi_j}\rangle = \left\langle \Psi_m \middle| H \middle| \Psi_j \right\rangle = 0$,取复共轭,即 $ \left\langle \Psi_j \middle| H \middle| \Psi_m \right\rangle = 0$。可见左边的角标的确小于等于右边的角标减 2,证毕。

   然后证明 $ \boldsymbol{\mathbf{H}} $ 是实矩阵,考虑到厄米矩阵 $H_{mn} = H_{nm}^* $ 的性质,只需要证明 $H_{j, j+1}$ 是实数。对任意 $j = 2 \dots K-1$

\begin{equation} \begin{aligned} H_{j, j+1} &= \left\langle \Psi_j \middle| H \middle| \Psi_{j+1} \right\rangle = \left\langle \Psi_j \right\rvert H \left(H \left\lvert \Psi_j \right\rangle - \alpha_j \left\lvert \Psi_j \right\rangle - \beta_j \left\lvert \Psi_{j-1} \right\rangle \right) \\ &= \left\langle \Psi_j \middle| H^2 \middle| \Psi_j \right\rangle - \alpha_j \left\langle \Psi_j \right\rvert H \left\lvert \Psi_j \right\rangle - \beta_j \left\langle \Psi_j \right\rvert H \left\lvert \Psi_{j-1} \right\rangle ~. \end{aligned} \end{equation}
易证 $ \left\langle \Psi_j \middle| H^2 \middle| \Psi_j \right\rangle $,$\alpha_j$ 及 $\beta_j$ 是实数,只需要证明 $ \left\langle \Psi_j \middle| H \middle| \Psi_{j-1} \right\rangle $ 是实数即可,即证 $ \left\langle \Psi_{j-1} \middle| H \middle| \Psi_j \right\rangle $ 是实数,重复上述论证,即证 $ \left\langle \Psi_{j-2} \middle| H \middle| \Psi_{j - 1} \right\rangle $ 是实数,……,即证 $ \left\langle \Psi_1 \middle| H \middle| \Psi_0 \right\rangle $ 是实数,即证 $ \left\langle \Psi_0 \middle| H \middle| \Psi_1 \right\rangle $ 是实数,即证 $ \langle{\Psi_0}|{H}|{\tilde\Psi_1}\rangle $ 是实数。而
\begin{equation} \langle{\Psi_0}|{H}|{\tilde\Psi_1}\rangle = \left\langle \Psi_0 \right\rvert H \left(H \left\lvert \Psi_0 \right\rangle - \alpha_0 \left\lvert \Psi_0 \right\rangle \right) = \left\langle \Psi_0 \middle| H^2 \middle| \Psi_0 \right\rangle - \alpha_0 \left\langle \Psi_0 \middle| H \middle| \Psi_0 \right\rangle ~. \end{equation}
显然是实数。以上的思路是数学归纳法的逆过程。


1. ^ 参考 Renate Pazourek 的论文


致读者: 小时百科一直以来坚持所有内容免费,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 10 元,我们一个星期内就能脱离亏损, 并保证在接下来的一整年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利