首页

共轭梯度法

星期三, 2009-03-04 14:34 — satchel1979

算法基础

      在各种优化算法中，共轭梯度法（Conjugate Gradient）是非常重要的一种。其优点是所需存储量小，具有$ N $步收敛性，稳定性高，而且不需要任何外来参数。下面对此算法作简要地介绍。

      首先介绍最速下降法（Steepest Decent）：对于一个目标函数$ F $，其优化过程首先需要确定优化方向。直观地考虑，此优化方向应该选择为$ F $在当前点处的负梯度方向$ d=-\nabla F $,因为$ F $沿这个方向数值下降最快。确定优化方向d之后，利用一维搜索算法（如解析式、进退法、三次函数法等等） 找出沿d方向上的最小值及其对应点。此后将该点作为新的出发点重复上述过程，直至$ \|d\|\leq \varepsilon $, $ \varepsilon $是预设的收敛判据,此时优化过程结束。

      最速下降法的收敛速度正比于$ F $的条件数（最大本征值与最小本征值的比值），因此如果$ F $的条件数太大，考虑一个二维的图像， 则$ F $等值线是一系列离心率极大的椭圆。在距离极值点较远的地方其梯度方向并不指向极值点所在位置，而是沿着该椭圆的斜长轴方向，在这种情况下执行最速下降， 其结果是搜索方向沿着"之"字形曲折的逼近极值点。显然这种方法的效率比较低下。

共轭梯度法在确定优化方向d上采用了和最速下降法不同的策略。最速下降法的连续两次过程中的搜索方向$ d^{n+1} $和$ d^{n} $之间相互独立。但是在共轭梯度法中$ d^{n+1} $和$ d^{n} $通过矩阵$ \mathbf{A} $相互联系。具体来说，为方便讨论，我们选择一个二次函数$ F $，形如$ F=\frac{1}{2}x^{T}Ax-bx $。则满足方程

                                             $ (d^m)^T\cdot A\cdot d^n = 0 (m\neq n) $    (1)

 的两个矢量$ d^m $和$ d^n $彼此关于$ \mathbf{A} $共轭。可以通过下述方法构建彼此共轭的一系列搜索方向：

                          $ d^{n+1}=-\nabla F^{n+1}+\frac{\|\nabla F^{n+1}\|^2}{\|\nabla F^{n}\|^2}\cdot d^n ; d^0=-\nabla F\mid_{x=x_0} $      (2)

与最速下降法相同，每次确定搜索方向之后，利用一维搜索算法或者解析式（仅对二次型函数有效）得出$ F $在此方向上的极小值。共轭梯度法具有$ N $步收敛性， 也即对于一个自由度为N的体系而言，理论上可以在$ N $步内找到极小值。与此前一样，考虑一个二维的图像，对于一个正定的二次函数$ F $，通过配方可 以将其表示为$ F(x)=\frac{1}{2}(x-\bar{x})A(x-bar{x}) $，因此F的等值面是以$ \bar{x} $为中心的椭圆。设在某等值面上有一点$ x^{(1)} $，则其此处的导数为

                                                  $ \nabla F(x^{(1)})=A(x^{(1)}-\bar{x}) $

 定义$ d^{(2)}=x^{(1)}-\bar{x} $，而$ d^{(1)} $是过$ x^{(1)} $的该等值面上的切向量，则$ d^{(1)} $与$ \nabla F(x^{(1)}) $正交，因此
 
                                                          $ (d^{(1)})^TAd^{(2)}=0 $
 
 也即$ F $等值面上一点处的切向量与该点指向函数极小点的向量关于$ \mathbf{A} $共轭。

       关于这一点的讨论也可以直观的想象：假设$ F $的等值线是一系列离心率极大的椭圆，第一步中两种方法的搜索方向一致，但是第二步中，共轭梯度法所确定的搜索方向直接指向极值点。这是因为考虑共轭性之后，相当于将这一系列的椭圆经过一系列各向异性的伸缩变换，
 转化为一系列圆形。而第二步所确定的$ d^1 $，相当于沿着直径指向该圆形的圆心。因此，共轭梯度法的实现算法如下：
 
 1. 确定起始点$ x_0 $处$ F $的负梯度方向：$ d^0=-\nabla F\mid_{x=x_0} $，若$ \|d^0\|^2 \leq \varepsilon $，则优化结束，转道步骤4；

 2. 设$ i=0 $，沿$ d^i $进行一维搜索，确定$ F $在$ d^i $方向上的最优位置:$ x^{i+1}=x^i+\lambda d^i $，对于正定的二次函数$ F $，$ \lambda $有解析形式:
                             $ \lambda = -\frac{(g^k)^T\cdot d^k}{(d^k)^TAd^k} $

 3. 依据方程（2）构建新的共轭梯度方向$ d^{i+1} $，若$ \|d^{i+1}\|^2 \leq \varepsilon $则优化结束，转到步骤4，否则设$ i=i+1 $，回到步骤2；

 4. 计算$ F\mid_{x=x_i} $并保存$ x_i $。

       注意：使用共轭梯度法要保证初始的搜索方向$ d^0 $必须是目标函数在当前点的负梯度方向。

DFT理论中的CG方法应用

      在DFT理论框架下，对于弱相互作用体系，求解体系基态等价为优化下述泛函：

                     $ E=\langle \phi|\hat{H}|\phi \rangle -\varepsilon_n \langle \phi|\hat{S}|\phi \rangle $   (3)

第二项来自于本征函数正交性的约束。在特定的函数基下看待此问题，则本征函数$ |\phi \rangle $ 相当于矢量x，Hamiltonian算符$ \hat{H} $表征为一个矩阵$ \mathbf{A} $，$ E $等于目标函数$ F $，则上述问题等价为优化一个二次函数。因此可以利用共轭梯度法进行求解。实际上，如果考虑Kohn-Sham方程$ \hat{H}|\phi\rangle = \varepsilon \hat{S} |\phi\rangle $，其广义本征值问题与上述问题等价。也即我们可以利用共轭梯度法求出最为接近实际本征值的近似本征值及其相应的本征矢。利用共轭梯度法可以迭代的求解Kohn-Sham方程， 不同于直接对角化Hamiltonian矩阵或者Car-Parrinello动力学方法（这一类方法统称直接求解法），关于共轭梯度直接法将在下文做简要介绍。

      在具体应用共轭梯度法的时候，需要考虑以下几个因素：确定最速下降方向以及共轭方向、正交化处理以及利用预处理技术提升收敛速度。下面分别介绍如下：

 最速下降方向：

      选定某条能带$ m $，根据方程（3），即为优化泛函$ E=\langle \phi_m|\hat{H}-\varepsilon_m\hat{S}|\phi_m\rangle $。由此定义第$ i $次迭代时相应的剩余矢量为：

                                          $ |R(\phi^i_m)\rangle = -(\hat{H}-\varepsilon_m\hat{S}|\phi^i_m\rangle) $  (4)
 
也即此处的最速下降方向，与$ |\phi^i_m $正交。因此此时的Lagrangian乘子可计算如下：
 
                         $ \varepsilon^i_m = \frac{\langle \phi^i_m|\hat{H}|\phi^i_m\rangle}{\langle \phi^i_m|\hat{S}|\phi^i_m\rangle} $   (5)

这个值也是第$ m $个本征值的当前最佳估计值。

正交化：

     正交化的要求源自下列矛盾：共轭梯度法是一种无约束的优化方法，但是如果要进行一系列能量本征值的求解，那么分属不同本征值的本征矢彼此正交。这个约束条件可以通过对最速下降方向进行正交化处理而转化为无约束优化问题。因为每条能带的本征矢都与其他的本征矢正交，所以假设我们已经将小于$ m $以下的所有能带优化完毕，那么第m条能带应该是满足与其余$ m $-1个本征矢正交这个限制条件的最小值。因此，第$ i $次迭代的最速下降方向应该与$ m $-1个本征矢正交。利用Gram-Schmidt正交化方案实现如下：

                                     $ |\varsigma^i_m\rangle=(1-\sum_{n<m}|\phi^i_n\rangle\langle\phi^i_m|S)|\phi^i_m\rangle $          (6)

预处理：

      理论上讲，在$ N $步之内对能带$ m $的优化可以结束。但是可以进行一系列操作提高优化效率。从数学上讲，等于对于矩阵$ \mathbf{A} $进行相似变换，改善其条件数， 使得尽量多的本征值简并。这种操作之所以可以提高效率，是因为如果最速下降方向是偏差（即当前解与精确解之差）乘以一个常数的话，那么沿着最速下降方向移动适当的距离就可以非常精确的到达精确解。

     设当前步骤下最优解与精确解之间的差别为$ |\delta\phi^i_m\rangle $ ，这个量可以用体系的本征值展开为：
 
                                                          $ |\delta\phi^i_m\rangle = \sum_{n}c_{n,m}|\phi_n\rangle $        (7)

因此第$ m $条能带的精确解可以写为$ |\phi_m\rangle = |\phi^i_m\rangle + |\delta\phi^i_m\rangle $，将其带入方程（4），可得

                          $ |R(\phi^i_m)\rangle = -(H-\varepsilon_mS)|\phi_m\rangle + (H-\varepsilon_mS)|\delta\phi_m\rangle $

在比较接近精确解的时候，上述方程的头一项可以忽略，仅考虑第二项即可，将方程（7）带入，则可得：

                     $ |R(\phi^i_m)\rangle = (H-\varepsilon_mS)|\delta\phi_m\rangle = \sum{n}(\varepsilon_n-\varepsilon^i_m)c_{n,m}|\phi_n\rangle $

可以看出，如果$ n $个能级彼此简并的话，则最速下降方向是$ |\delta\phi^i_m $ 的一个常数倍。因此，如前所述，沿着最速下降方向移动适当的距离就可以非常精确的到达精确解。这可以大大提升共轭梯度法的收敛速度。而预处理可以通过乘以一个预处理矩阵$ \mathbf{K} $得以实现。而$ \mathbf{K} $决定于在计算中所采用的基函数。

     以平面波基为例，对于$ G $较高的平面波，动能项为主要项，因此，如果我们要构造一个简并度比较高的变换，那么矩阵$ \mathbf{K} $最方便的选择是一个对角阵，对角元是动能的倒数。但是对于$ G $较低的平面波，动能项不占优势，因此$ \mathbf{K} $应该趋近于1。一般取下面的表达式：

                                       $ \mathbf{K}_{m,n}=\frac{27+18x+12x^2+8x^3}{27+18x+12x^2+8x^3+16x^4}\delta_{m,n} $      (8)

式中的$ x $为平面波动能与$ |\phi^i_m\rangle $动能的比值。

     同时考虑最速下降方向的正交性与预处理，可以由此构造最速下降方向为：

                                            $ |\eta^i_m\rangle = \mathbf{K}|\varsigma^i_m\rangle $

 但是乘以矩阵$ \mathbf{K} $会破坏正交性，因此需要对$ |\eta^i_m\rangle $再额外进行一次Gram-Schmidt正交化：

        $ |\eta'^{i}_m\rangle = |\eta^{i}_m\rangle-\langle\phi^i_m|\eta^i_m\rangle-\sum_{n<m}\langle\phi_n|\eta^i_m\rangle|\phi_n\rangle $        (9)

将$ |\eta'^i_m\rangle $作为最速下降方向。注意（9）式中第三项里的$ |\phi_m\rangle $没有上标，表明$ m $以下的能带均已优化到精确解。

共轭方向:

     已经确定最速下降方向之后，可以依照经典的共轭梯度方法构造共轭方向：

                           $ |\varphi^i_m\rangle=|\eta'^i_m\rangle+\gamma^i_m|\varphi^(i-1)_m\rangle;\gamma^i_m=\frac{\langle\eta'^i_m|\varsigma^i_m\rangle}{\langle\eta'^{i-1}_m|\varsigma^{i-1}_m\rangle} $         (10)

应当注意，共轭方向中的$ \gamma^i_m $不仅仅是一种表达式。比如Dyutiman Das采用了$ \gamma^i_m $的另外一种形式：

                  $ \gamma^i_m=\frac{(\langle\eta'^i_m|-\langle\eta'^{i-1}_m|)|\eta'^i_m\rangle}{(\langle\eta'^i_m|-\langle\eta'^{i-1}_m|)|\varphi^i_m\rangle} $

这些表达式在没有外约束的条件下是彼此等价的，但是因为本征矢正交性的限制，由不同的$ \gamma^i_m $构造出来的$ |\varphi^i_m\rangle $并不相同，很难说哪一种效率更高，需要在具体的问题中通过测试决定。

     最后，当前的共轭方向$ |\varphi^i_m\rangle $还需要与当前第$ m $个能带的本征矢正交，并归一化。因此，最后的共轭方向的表达式为：

                         $ |\varphi''^i_m\rangle = |\varphi^i_m\rangle - \langle\phi^i_m|\varphi^i_m\rangle|\phi^i_m\rangle;\quad |\varphi'^i_m\rangle = \frac{|\varphi''^i_m}{\langle\varphi''^i_m|\varphi''^i_m\rangle} $          (11)

而以$ |\varphi'^i_m\rangle $为最终的共轭梯度方向作为优化方向。

一维搜索:
确定了优化方向，需要沿优化方向求出目标函数的最优解，而相应的本征矢$ |\phi^{i+1}_m\rangle $相当于当前本征矢$ |\phi^i_m\rangle $和优化方向$ |\varphi'^i_m\rangle $的一个线性叠加。考虑到$ |\phi^{i+1}_m\rangle $和$ |\phi^i_m\rangle $的模方应该相等，因此可写为$ |\phi^{i+1}_m\rangle=cos\theta|\phi^{i+1}_m\rangle+sin\theta|\varphi'^i_m\rangle $。因此优化参数$ \theta $即可。前面说过，对于二次正定的函数，优化步长有解析的形式。假设我们采用经验赝势方法，则可写出$ \theta_{min} $的解析式为：

               $ tan2\theta_{min}=\frac{2\langle\varphi'^i_m|H|\phi^i_m\rangle}{\langle\phi^i_m|H|\phi^i_m\rangle-\langle\varphi'^i_m|H|\varphi'^i_m\rangle} $     (12)

如果采用严格的第一原理计算，那么$ \theta_{min} $虽然仍然有解析形式，但是需要考虑实空间内交换关联能以及Hatree项的积分。另外一种方法则需要算出$ \theta=0 $时的函数值以及 一阶导数值，以及另一个$ \theta $（通常取$ \frac{\pi}{300} $）时的函数值，具体的步骤为：

     首先将能量$ E $写为关于$ \theta $的三角级数：$ E(\theta)=E_0+\sum_{n}[A_{n}cos(2n\theta)+B_{n}sin(2n\theta)] $。Payne和Josnnopoulods指出，这个级数中$ n $> 1的项均可以省略。因此上式简化为$ E(\theta)=E_0+A_1cos(2\theta)+B_1sin(2\theta) $。因此，如果要确定$ \theta_{min} $，我们需要首先求出$ E_0 $、$ A_1 $以及$ B_1 $三个参数的值。可以通过下述四个方程求解：

                $ E_0=\frac{E[\frac{\pi}{300}]-\frac{1}{2}\frac{\partial E}{\partial\theta}|_{\theta=0}-E(0)cos[\frac{2\pi}{300}]}{1-cos[\frac{2\pi}{300}]} $
                $ A_1=\frac{E(0)-E[\frac{\pi}{300}]+\frac{1}{2}\frac{\partial E}{\partial\theta}|_{\theta=0}}{1-cos[\frac{2\pi}{300}]} $
                $ B_1=\frac{1}{2}\frac{\partial E}{\partial\theta}|_{\theta=0} $
                $ \frac{\partial E}{\partial\theta}|_{\theta=0}=\langle\varphi'^i_m|H|\phi^i_m\rangle+\langle\phi^i_m|H|\varphi'^i_m\rangle $

则极值点$ \theta_s=\frac{1}{2}\arctan(\frac{B_1}{A_1}) $，取区间$ [0,\frac{\pi}{2}] $中的$ \theta_s $即为所求的$ \theta_0 $。这两种方法的计算量相差无几。

共轭梯度法求解算法：

      综上所述，利用共轭梯度法求解KS方程的本征值具体步骤如下：

 1．预设收敛判据$ \tau $以及$ \lambda $，初始化$ N $个本征矢$ {\phi} $（如利用随机数作为系数），设j=0，$ {\phi^j}={\phi} $，依照原子电荷分布的叠加作为初始电荷密度$ \rho^j_{in} $；

 2. 选择最低的能带$ m $，设$ i=0 $，根据方程（5）求出在上述$ \rho^j_{in} $下的期待值$ \varepsilon^i_m $，并由方程（4）计算剩余矢量$ |R(\phi^i_m)\rangle $；

 3. 依次按照方程（6）~（9）对$ |R(\phi^i_m)\rangle $进行操作，并通过（10）和（11）构造归一化的共轭梯度方向$ |\varphi'^i_m $；

 4．进行一维搜索，利用方程（12）或者上述的两种方法计算优化的本征矢$ \phi^{i+1}_m $，计算出$ \varepsilon^{i+1}_m $以及$ \||R(\phi^{i+1}_m)\rangle\| $，
 若$ \||R(\phi^{i+1}_m)\rangle\|\leq \tau $，优化结束，设$ m=m+1 $，移到下一条能带，否则设$ i=i+1 $，回到步骤2；

 5．重复上述过程，直至收敛或者$ m\geq N $，计算总能$ E $以及能量变化值$ \Delta E^j $，若$ \Delta E^j \leq \lambda $，全部计算结束，转到步骤6，否则设$ j=j+1 $，利用$ {\phi^j} $更新Hamiltonian矩阵以及电荷密度$ \rho^j_{in} $，回到步骤1；

 6．计算并保存电荷密度、本征矢、总能等各种信息。

 子空间转动

      对最速下降方向以及共轭方向进行约束的方法不只上面一种。如果正交化不仅仅对于这些能带进行，而是将$ n $的取值遍历所有能带指标，则利用类似的算法可以得到同样的本征能级， 但是不能保证所得的本征矢正确，实际上，这些本征矢是KS本征矢的线性叠加。因此，对于金属而言，需要在CG过程结束之后再进行一步"子空间转动"（subspace rotaion）的步骤，即以占据数不为零的所有能带对应的本征矢展开子空间，计算Hamiltonian矩阵$ \mathbf{H} $以及交叠矩阵$ \mathbf{S} $，再进行矩阵的直接对角化。所得的本征矢$ \{B\} $乘以此前得到的本征矢$ \{\phi^j\} $，即得最后的结果。

      可以看到，在优化全部$ N $个本征矢的过程中，电荷密度分布$ \rho^j_{in} $是固定的，直到所有$ N $个本征矢优化结束之后才更新$ \rho^j_{in} $。因此，上面介绍的方法是迭代算法。 我们也可以利用共轭梯度法对体系进行直接求解。过程和上述算法大体一致。主要的区别在于对第$ m $条能带优化结束之后，首先要更新$ \rho_{in} $，然后再移到下一条能带。在实际的计算中，特别是平面波基方法中，最开始的几轮优化对于$ \rho_{in} $的更新会导致系统严重偏离基态（因为本征矢的初始化采用随机数），因此这种做法存在可行性方面的困难。一种解决办法是，保持初始$ \rho_{in} $不变，对所有能带做优化，从第二步开始随时更新$ \rho_{in} $，为提高效率，共轭梯度法的收敛判据$ \tau $以及$ \lambda $随着优化不断进行而逐渐减小，而不采用固定值。