Kernel Polynomial Method简介
Kernel Polynomial Method (KPM)
Kernel polynomial method (KPM) 是一种利用数值 Chebyshev 多项式展开来降低计算复杂度的方法。传统计算谱时常需 \(O(N^3)\) 的对角化方法，而 KPM 则将问题转化为矩阵与一组较小向量的乘法，大幅降低复杂度。同期发展的还有 Lanczos 算法，Lanczos 适合求本征态，KPM 更适合谱性质和关联函数的计算 [1]。
KPM 源自函数的正交多项式展开。定义权重函数 \(w(x)\) 及区间 \([a,b]\) 后，两个可积函数 \(f,g:[a,b]\to\mathbb{R}\) 的内积为：
\[ \begin{equation} \langle f|g \rangle = \int_a^b w(x) f(x) g(x) dx \end{equation} \]
通过 Schmidt 正交化可得正交多项式 \(p_n(x)\)，满足：
\[ \begin{equation} \langle p_n|p_m \rangle = \delta_{nm}/h_n \end{equation} \] 其中 \(h_n = ⟨p_n|p_n⟩\)。任意函数可展开为： \[ \begin{equation} f(x) = \sum_{n=0}^\infty \alpha_n p_n(x), \quad \alpha_n = \frac{1}{h_n}\langle p_n|f \rangle \end{equation} \] 通常意义上任何正交多项式都可以用来展开函数从而使用KPM，但经过历史和经验的检验，第一类和第二类Chebyshev多项式由于其优秀的收敛性以及与傅里叶变换的相似性而被认为是最佳选择。 ### Chebyshev 多项式 Chebyshev 多项式定义在区间 \([a,b]=[-1,1]\) 上。第一类 Chebyshev 多项式 \(T_n(x)\) 的权重函数为 \(w(x)=\frac{1}{\pi\sqrt{1-x^2}}\)，第二类 \(U_n(x)\) 的权重函数为 \(w(x)=\pi\sqrt{1-x^2}\)。可以验证：
\[ \begin{equation} \begin{split} T_n(x) &= \cos\left[n \cos^{-1}x\right] \\ U_n(x) &= \frac{\sin\left[(n+1)\cos^{-1}x\right]}{\sin\left[\cos^{-1}x\right]} \end{split} \end{equation} \]
两类 Chebyshev 多项式都存在递推关系。对于第一类：
\[ \begin{equation} \begin{split} T_0(x) &= 1, \\ T_{-1}(x) &= T_1(x) = x, \\ T_{n+1}(x) &= 2xT_n(x)-T_{n-1}(x) \end{split} \end{equation} \]
对于第二类：
\[ \begin{equation} \begin{split} U_0(x) &= 1, \\ U_{-1}(x) &= 0, \\ U_{n+1}(x) &= 2xU_n(x)-U_{n-1}(x) \end{split} \end{equation} \]
此外还有乘积关系：
\[ \begin{equation} \begin{split} 2T_m(x)T_n(x) &= T_{m+n}(x)+T_{m-n}(x), \\ 2(x^2-1)U_{m-1}(x)U_{n-1}(x) &= T_{m+n}(x)-T_{m-n}(x) \end{split} \end{equation} \]
两种权重函数分别给出两种内积定义：
\[\begin{equation} \begin{split} \langle f|g \rangle_1 &= \int_{-1}^{1} \frac{f(x)g(x)}{\pi\sqrt{1-x^2}} dx, \\ \langle f|g \rangle_2 &= \int_{-1}^{1} \pi\sqrt{1-x^2} f(x)g(x) dx \end{split} \end{equation}\]
对应的正交性关系为：
\[\begin{equation} \begin{split} \langle T_n|T_m \rangle_1 &= \frac{1+\delta_{n,0}}{2}\delta_{n,m} \\ \langle U_n|U_m \rangle_2 &= \frac{\pi^2}{2}\delta_{n,m} \end{split} \end{equation}\]
定义 \(\phi_n(x)=\frac{T_n(x)}{\pi\sqrt{1-x^2}}\)，则在第二种内积下有：
\[\begin{equation} \langle \phi_n|\phi_m \rangle_2 = \frac{1+\delta_{n,0}}{2}\delta_{nm} \end{equation}\]
可以验证任意函数都可以展开为：
\[\begin{equation} \begin{split} f(x) &= \sum_{n=0}^{\infty} \frac{\langle f|\phi_n \rangle_2}{\langle \phi_n|\phi_n \rangle_2} \phi_n(x) \\ &= \frac{1}{\pi\sqrt{1-x^2}}\left(\mu_0 + 2\sum_{n=1}^{\infty} \mu_n T_n(x)\right) \end{split} \end{equation}\]
其中系数\(\mu_n\)叫做矩，定义为：
\[\begin{equation} \mu_n = \langle f|\phi_n \rangle_2 = \int_{-1}^{1} f(x) T_n(x) dx \end{equation}\]
核多项式与 Gibbs 现象
实际计算中通常需要截断展开，直接截断会导致严重的 Gibbs 现象。为此常引入阻尼 kernel 系数 \(\mu_n \to g_n \mu_n\) 修正展开系数：
\[\begin{equation} \begin{split} f(x) &= \frac{1}{\pi\sqrt{1-x^2}}\left(\mu_0 + 2\sum_{n=1}^{\infty} \mu_n T_n(x)\right) \\ &\to \frac{1}{\pi\sqrt{1-x^2}}\left(\mu_0 g_0 + 2\sum_{n=1}^{N-1} \mu_n g_n T_n(x)\right) \end{split} \end{equation}\]
常用的 Jackson 核为：
\[\begin{equation} g_n^{\text{J}} = \frac{1}{N+1}\left[(N-n+1)\cos\frac{\pi n}{N+1} + \sin\frac{\pi n}{N+1}\cot\frac{\pi}{N+1}\right] \end{equation}\]
格林函数计算中常用 Lorentz 核：
\[\begin{equation} g_n^L = \frac{\sinh\left[\lambda\left(1-\frac{n}{N}\right)\right]}{\sinh \lambda} \end{equation}\]
其中参数 \(\lambda \in \mathbb{R}\)，通常取 \(3\sim5\)。
态密度计算
区间缩放
Chebyshev 展开定义在 \([-1,1]\)，对于矩阵运算需保证本征值 \(\{E_k\}\) 也在 \([-1,1]\)。可作线性变换：
\[\begin{equation} \begin{split} \tilde{H} &= \frac{H-b}{a} \\ \tilde{E} &= \frac{E-b}{a} \end{split} \end{equation}\]
其中 \(a, b\) 由哈密顿量最大最小本征值 \(E_{\text{max}}, E_{\text{min}}\) 决定：
\[\begin{equation} a = \frac{E_{\text{max}}-E_{\text{min}}}{2-\epsilon} \\ b = \frac{E_{\text{max}}+E_{\text{min}}}{2} \end{equation}\]
参数\(\epsilon\) 用于避免谱边界刚好处于或超过 \([-1,1]\)，通常取 \(0.01\)。
最大最小本征值可用 Lanczos 算法（如 Matlab 的 eigs）获得：
1 | E_max = eigs(H_sparse, 1, 'largestreal') |
离散余弦傅里叶变换
在对函数 \(f(x)\) 展开作截断并利用核多项式修正后，对于给定点 \(x\)，有如下展开式：
\[ \begin{equation} f(x) = \frac{1}{\pi\sqrt{1-x^2}}\left(\mu_0 g_0 + 2\sum_{n=1}^{N-1} \mu_n g_n T_n(x)\right) \end{equation} \]
利用 Chebyshev 多项式的性质，还可以进一步减少计算时间。定义核多项式修正后的矩系数为
\[ \begin{equation} \tilde{\mu}_n = \mu_n g_n \end{equation} \]
此时展开式可以写为
\[ \begin{equation} f(x) = \frac{1}{\pi\sqrt{1-x^2}}\left(\tilde{\mu}_0 + 2\sum_{n=1}^{N-1} \tilde{\mu}_n T_n(x)\right) \end{equation} \]
如果采样点 \(x\) 取为特殊的 \(x_k\)，满足
\[ \begin{equation} x_k = \cos\frac{\pi(k+1/2)}{\tilde{N}},\quad k=0,\ldots,(\tilde{N}-1) \end{equation} \]
则有
\[ \begin{equation} T_n(x_k) = \cos\left[n\cos^{-1}x_k\right] = \cos\left[\frac{n\pi(k+1/2)}{\tilde{N}}\right] \end{equation} \]
这些点被称为 Chebyshev 节点 [2]。\(\tilde{N}\) 通常取 \(2N\)。在这些点上，\(f(x_k)\) 可以通过对 \(\tilde{\mu}_n\) 作第三类离散余弦变换（DCT-III）得到：
\[ \begin{equation} y[k] = \frac{1}{2}x_0 + \sum_{n=1}^{N-1} x_n \cos\left[\frac{n\pi(k+1/2)}{N}\right] \end{equation} \]
由于第三类余弦变换可以利用快速傅里叶变换（FFT）加速，因此比直接乘积更快。Matlab 示例代码：
1 | omega_dct = cos(pi * ((0:N_points-1)+0.5)/N_points); |
注意：第二行对矩的第一个元素除以了 \(\sqrt{2}\)，这是由于 Matlab 的 dct 定义与上述第三类余弦变换的系数在第一项与剩余项之间相差一个 \(\sqrt{2}\)，若不修改，会导致余弦变换的系数收敛到零的速度变慢，给出错误的态密度。
矩的计算
KPM 展开可以被应用到算符上，例如哈密顿量 \(\hat{H}\)。在利用 KPM 计算时需要关注两种形式的矩：
其一是两个态的期望： \[ \begin{equation} \mu_n = \langle \beta | T_n(\tilde{H}) | \alpha \rangle \end{equation} \]
其二为算符 \(\hat{A}\) 与多项式的迹： \[ \begin{equation} \mu_n = \mathrm{Tr}[A T_n(\tilde{H})] \end{equation} \]
对于第一种情况，可以利用迭代关系直接计算。令 \(|\alpha_n\rangle = T_n(H)|\alpha\rangle\)，有
\[ \begin{equation} \begin{split} |\alpha_0\rangle &= |\alpha\rangle \\ |\alpha_1\rangle &= \tilde{H}|\alpha_0\rangle \\ |\alpha_{n+1}\rangle &= 2\tilde{H}|\alpha_n\rangle - |\alpha_{n-1}\rangle \end{split} \end{equation} \]
随即可算出矩 \(\mu_n = \langle \beta | \alpha_n \rangle\)。整个过程只涉及到矩阵乘法，且可以利用稀疏矩阵方法提升大尺寸运算效率。如果 \(|\beta\rangle = |\alpha\rangle\)，还可以进一步简化。利用乘积关系可验证，对于偶数次矩有
\[ \begin{equation} \begin{split} 2\langle \alpha_n | \alpha_n \rangle - \mu_0 &= 2\langle \alpha | T_n(\tilde{H})^\dagger T_n(\tilde{H}) | \alpha \rangle - \mu_0 \\ &= 2\langle \alpha | T_n(\tilde{H}) T_n(\tilde{H}) | \alpha \rangle - \mu_0 \\ &= \langle \alpha | \left(T_0(\tilde{H}) + T_{2n}(\tilde{H})\right) | \alpha \rangle - \mu_0 \\ &= \langle \alpha | \alpha_0 \rangle + \langle \alpha | \alpha_{2n} \rangle - \mu_0 \\ &= \mu_{2n} \end{split} \end{equation} \]
奇数次矩有
\[ \begin{equation} \begin{split} 2\langle \alpha_{n+1} | \alpha_n \rangle - \mu_1 &= 2\langle \alpha | T_{n+1}(\tilde{H})^\dagger T_n(\tilde{H}) | \alpha \rangle - \mu_1 \\ &= 2\langle \alpha | T_{n+1}(\tilde{H}) T_n(\tilde{H}) | \alpha \rangle - \mu_1 \\ &= \langle \alpha | \left(T_1(\tilde{H}) + T_{2n+1}(\tilde{H})\right) | \alpha \rangle - \mu_1 \\ &= \langle \alpha | \alpha_1 \rangle + \langle \alpha | \alpha_{2n+1} \rangle - \mu_1 \\ &= \mu_{2n+1} \end{split} \end{equation} \]
因此当 \(|\alpha\rangle = |\beta\rangle\) 时，可将矩的计算减半：
\[ \begin{equation} \begin{split} \mu_{2n} &= 2\langle \alpha_n | \alpha_n \rangle - \mu_0 \\ \mu_{2n+1} &= 2\langle \alpha_{n+1} | \alpha_n \rangle - \mu_1 \end{split} \end{equation} \]
对于第二种情况，以态密度为例。假设本征态 \(\tilde{H}|k\rangle = E_k|k\rangle\)，态密度定义为
\[ \begin{equation} \rho(E) = \sum_k \delta(E-E_k) = \mathrm{Tr}\left[\delta(E-\tilde{H})\right] \end{equation} \]
\(\delta\) 函数可用 Chebyshev 多项式展开为矩：
\[ \begin{equation} \begin{split} \mu_n &= \int_{-1}^{1} \rho(E) T_n(E) dE \\ &= \int_{-1}^{1} \sum_k \delta(E-E_k) T_n(E) dE \\ &= \sum_k T_n(E_k) \end{split} \end{equation} \]
利用 \(T_n(H)|k\rangle = T_n(E_k)|k\rangle\)，可得
\[ \begin{equation} \mu_n = \sum_k \langle k | T_n(H) | k \rangle = \mathrm{Tr}[T_n(H)] \end{equation} \]
由于 Chebyshev 多项式存在递推关系，有
\[ \begin{equation} \begin{split} |\alpha_n(k)\rangle &= T_n(H)|k\rangle \\ |\alpha_{n+1}(k)\rangle &= 2H|\alpha_n(k)\rangle - |\alpha_{n-1}(k)\rangle \end{split} \end{equation} \]
此时矩变为
\[ \begin{equation} \mu_n = \sum_k \langle \alpha_0(k) | \alpha_n(k) \rangle = \sum_k \langle k | T_n(H) | k \rangle \end{equation} \]
注意到之前公式中我们可以将态密度转换为了展开系数的计算，而上式表明，展开系数的计算只需要计算两个向量的内积，又因为\(|\alpha_n(k)\rangle\)各阶项实际上是一个矩阵乘法运算，通过上述的步骤，我们可以跳过高复杂度的对角化方法，仅通过矩阵乘法就可以得到态密度的展开系数。
实际上，计算迹时并不需要对所有本征态求和，只需对少量随机归一化向量 \(|r\rangle\) 取平均即可获得对 \(\mu_n\) 的良好估计：
\[ \begin{equation} \mu_n = \mathrm{Tr}[A T_n(H)] \approx \frac{1}{R} \sum_{r=1}^R \langle r | A T_n(H) | r \rangle \end{equation} \]
其中 \(|r\rangle = \frac{1}{\sqrt{D}} \sum_{i=1}^D \xi_{rk} |k\rangle\)，\(\xi_{rk}\) 为复随机的归一化向量，满足
\[ \begin{equation} \begin{split} \mathbb{E}[\xi_{rk}] &= 0 \\ \mathbb{E}[\xi_{rk}\xi_{r'l}] &= 0 \\ \mathbb{E}[\xi^*_{rk}\xi_{r'l}] &= \delta_{rr'}\delta_{kl} \end{split} \end{equation} \]
若 \(\xi_{rk}=e^{i\phi}\)，\(\phi\in[0,2\pi]\) 满足均匀或高斯分布，可以获得相当精确的结果。