Lanczos 算法

                     

贡献者： addis

1. 背景介绍

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

其中 $H$ 为哈密顿算符。若 $H$ 不含时，约等号变为等号 通常情况下，用有限个正交归一基底 $\Psi_0,\dots\Psi_{K-1}$ 近似展开 $\Psi( \boldsymbol{\mathbf{r}} , t+\Delta t)$，这时，$H$ 可以表示成矩阵的形式。 这样，幺正算符 $ \exp\left(- \mathrm{i} H\Delta t/\hbar\right) $ 也可以表示成矩阵。根据定义 此时若把 $ \boldsymbol{\mathbf{H}} $ 对角化（求解本征方程），得到本征矢列矩阵 $ \boldsymbol{\mathbf{P}} $（单位正交阵）以及本征值矩阵 $ \boldsymbol{\mathbf{\Lambda}} $（对角矩阵），则可进行基底变换变到 $H$ 矩阵的本征空间求解上式 要说明的是，这种算法的误差（除数值计算误差）来源于两个方面，第一是用 $ \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 基底为

我们把正交归一化后的基底记为 这样，用其展开哈密顿算符 在算法上，本来我们可以按部就班地按照以上步骤做，然而由于 Krilov 基底的性质，可以通过一些定理（见文章最后）大大减少计算量。正交归一化步骤如下

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

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

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

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

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

根据以上两个定理，正交化结果也可以写成

　　 【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 只用到了一个矩阵-矢量乘法，两个矢量数乘和两个矢量减法。

　　 接下来只要对 $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}} $ 是对称三对角矩阵。所以我们已经顺便求出了所有的矩阵元！

该矩阵具有维数小，易求本征问题的优势。令 那么 注意在 $ \left\lvert \Psi_i \right\rangle $ 基底下，$\Psi( \boldsymbol{\mathbf{r}} , t)$ 就是 $(1, 0, 0, \dots)$。所以在该基底下 这里的 $ \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）

推导：把式 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 $，所以（以下所有系数 只是用来表示线性组合，具体值不重要）

把 $H^{j+1} \left\lvert \Psi_0 \right\rangle $ 施密特正交归一化

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 $ 而

由于上式中 $ \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~,$

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 的证明类似，我们知道

所以对 $m \geqslant j + 2$ 有

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$

易证 $ \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 $ 是实数。而 显然是实数。以上的思路是数学归纳法的逆过程。

---

1. ^ 参考 Renate Pazourek 的论文
2. ^ 为什么只看第一列？因为当前波函数就是第一个 Krylov 基底。