第38卷第1期 大学物理 Vol.38 No.1 2019年1月 COLLEGE PHYSICS Jan.2019 多个含不同质量项的费米子矩阵求逆算法 吴良凯，顾鑫，刘坤，冯龙海 (江苏大学物理系，江苏镇江212013) 摘要：在格点量子色动力学的模拟中，Rational Hybrid Monte Carlo(RHMC)算法是一种精确的，能应用到任意多个味道数 的费米子的方法.它的思想是把费米子行列式展开为有理函数的形式.但该方法会带来很多彼此相差一个常对角矩阵的矩阵 的求逆的计算，消耗大量的时间和计算资源，限制了RHMC算法的应用.本文利用移位多项式，针对共轭梯度法得到多个含有 不同质量项的矩阵求逆的一种方法，该方法可以应用到RHMC算法中. 关键词：费米子：逆矩阵；格点量子色动力学：移位多项式 中图分类号：0413.3文献标识码：A 文章编号：1000-0712(2019)01-0029-03 【D0I】10.16854/j.cmki.1000-0712.180181 量子场论是在大学物理课程中的一门重要的课 U代表规范场，S.是规范场作用量 程.它是研究基本粒子之间相互作用的理论，其方法 格点量子色动力学进行的研究涉及到组态的演 在理论物理和凝聚态物理中得到了广泛的应用.量 化和物理量的计算，这些工作都要求进行费米子矩 子色动力学是其一个分支，是研究夸克，胶子之间作 阵的求逆，因此费米子矩阵的求逆在格点量子色动 用的有力武器.由于色相互作用在低能下的强相互 力学的模拟中居于中心的地位.矩阵求逆有很多的 作用性质，微扰论不能得到应用.诺贝尔奖得主Ken- 算法.事实上，所有的求解线性方程组的算法都可以 neth G.Wilson创立了格点量子色动力学.格点量子 用来求矩阵的逆矩阵.这些方法包括消元法，矩阵 色动力学把连续的4维时空离散化，把费米场定义 分解法，追赶法以及迭代法.消元法中还包括选主 在格点上，把规范场定义在连接格点的链上，其实质 元的消元法；矩阵分解法中包括LU分解法，QR分 是截断对动量的积分以得到有限的结果.通过Wik 解法，Cholesky分解法；迭代法包括Jacob迭代法， 转动从闵氏空间转到欧式空间，使得Monte Carlo模 Gauss--Siedel迭代法，共轭梯度法以及BiCGstab迭 拟能顺利进行格点量子色动力学在强子谱，强子矩 代法回.由于费米子矩阵是维数很大的稀疏矩阵， 阵元和相变的研究中取得了巨大成就0.由于在欧 共轭梯度法以及BiCGstab迭代法在费米子矩阵的 式空间中讨论场论，统计物理的很多概念（包括配 求逆中应用较广 分函数)可以推广到格点量子色动力学中。 组态的演化有多种算法，有R算法，中算法， 在量子场论中，系统的拉氏量可以表达为费米 Hybrid Monte Carlo算法，以及RHMC算法（其全称 场的双线性形式和规范场部分之和： 为Rational Hybrid Monte Carlo算法).RHMC算法是 L=-y.(0.-ig4,)w-rF 一种精确的，能应用到任意个味道费米子的算法.它 (1) 4 的思想是把费米行列式展开为有理函数的形式同 在格点量子色动力学中，在每一个格点上都有一个 由于费米场满足grassman代数，无法用普通的 费米场，费米场的双线性形式可以用矩阵表示为 复数场来表示，在模拟中不能直接处理费米场因此 -中，M中，其中M是费米子矩阵，ij标记不同的格 我们先积分掉费米场，得到费米子行列式，由于行列 点.因此系统的配分函数可写成如下的路径积分的 式的计算量惊人，我们用另外的普通的复数场的积 形式： 分来代替费米子行列式这个普通的复数场，我们用 P来表示（在格点上，P成为一个列矩阵），所以费米 Z=Dr Du DUexp[-中，M-Sg] (2〉 子作用量可以用P表示为回 收稿日期：2018-03-26：修回日期：2018-09-14 基金项目：国家自然科学基金(11347029)和江苏大学大学生科研项目资助 作者简介：吴良凯(1970一)，湖北武汉人，江苏大学物理系教师，博士，主要从事大学物理教学与研究工作 (C)1994-2019 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
书 第 38 卷第 1 期 大 学 物 理 Vol．38 No．1 2019 年 1 月 COLLEGE PHYSICS Jan． 2019 收稿日期: 2018－03－26; 修回日期: 2018－09－14 基金项目: 国家自然科学基金( 11347029) 和江苏大学大学生科研项目资助 作者简介: 吴良凯( 1970—) ，湖北武汉人，江苏大学物理系教师，博士，主要从事大学物理教学与研究工作． 多个含不同质量项的费米子矩阵求逆算法 吴良凯，顾 鑫，刘 坤，冯龙海 ( 江苏大学 物理系，江苏 镇江 212013) 摘要: 在格点量子色动力学的模拟中，Ｒational Hybrid Monte Carlo ( ＲHMC) 算法是一种精确的，能应用到任意多个味道数 的费米子的方法．它的思想是把费米子行列式展开为有理函数的形式．但该方法会带来很多彼此相差一个常对角矩阵的矩阵 的求逆的计算，消耗大量的时间和计算资源，限制了 ＲHMC 算法的应用．本文利用移位多项式，针对共轭梯度法得到多个含有 不同质量项的矩阵求逆的一种方法，该方法可以应用到 ＲHMC 算法中． 关键词: 费米子; 逆矩阵; 格点量子色动力学; 移位多项式 中图分类号: O 413．3 文献标识码: A 文章编号: 1000-0712( 2019) 01-0029-03 【DOI】10．16854 /j．cnki．1000-0712．180181 量子场论是在大学物理课程中的一门重要的课 程．它是研究基本粒子之间相互作用的理论，其方法 在理论物理和凝聚态物理中得到了广泛的应用．量 子色动力学是其一个分支，是研究夸克，胶子之间作 用的有力武器．由于色相互作用在低能下的强相互 作用性质，微扰论不能得到应用．诺贝尔奖得主 Kenneth G． Wilson 创立了格点量子色动力学．格点量子 色动力学把连续的 4 维时空离散化，把费米场定义 在格点上，把规范场定义在连接格点的链上，其实质 是截断对动量的积分以得到有限的结果．通过 Wick 转动从闵氏空间转到欧式空间，使得 Monte Carlo 模 拟能顺利进行．格点量子色动力学在强子谱，强子矩 阵元和相变的研究中取得了巨大成就［1］．由于在欧 式空间中讨论场论，统计物理的很多概念( 包括配 分函数) 可以推广到格点量子色动力学中． 在量子场论中，系统的拉氏量可以表达为费米 场的双线性形式和规范场部分之和: L = －ψγμ( μ－igAμ ) ψ－ 1 4 FμνFμν ( 1) 在格点量子色动力学中，在每一个格点上都有一个 费米场，费米场的双线性形式可以用矩阵表示为 －ψiMijψj ，其中 Mij是费米子矩阵，i、j 标记不同的格 点． 因此系统的配分函数可写成如下的路径积分的 形式: Z = ∫DψDψDUexp［－ ψiMijψj － Sg］ ( 2) U 代表规范场，Sg 是规范场作用量． 格点量子色动力学进行的研究涉及到组态的演 化和物理量的计算，这些工作都要求进行费米子矩 阵的求逆，因此费米子矩阵的求逆在格点量子色动 力学的模拟中居于中心的地位．矩阵求逆有很多的 算法．事实上，所有的求解线性方程组的算法都可以 用来求矩阵的逆矩阵．这些方法包括消元法，矩阵 分解法，追赶法以及迭代法．消元法中还包括选主 元的消元法; 矩阵分解法中包括 LU 分解法，QＲ 分 解法，Cholesky 分解法; 迭代法包括 Jacob 迭代法， Gauss－Siedel 迭代法，共轭梯度法以及 BiCGstab 迭 代法［2］． 由于费米子矩阵是维数很大的稀疏矩阵， 共轭梯度法以及 BiCGstab 迭代法在费米子矩阵的 求逆中应用较广． 组态的演化有多种算法，有 Ｒ 算法， 算法， Hybrid Monte Carlo 算法，以及 ＲHMC 算法( 其全称 为 Ｒational Hybrid Monte Carlo 算法) ．ＲHMC 算法是 一种精确的，能应用到任意个味道费米子的算法．它 的思想是把费米行列式展开为有理函数的形式［3］． 由于费米场满足 grassman 代数，无法用普通的 复数场来表示，在模拟中不能直接处理费米场．因此 我们先积分掉费米场，得到费米子行列式，由于行列 式的计算量惊人，我们用另外的普通的复数场的积 分来代替费米子行列式．这个普通的复数场，我们用 φ 来表示( 在格点上，φ 成为一个列矩阵) ，所以费米 子作用量可以用 φ 表示为［4］
30 大学物理 第38卷 S,= (3) 为了定出系数，以及相应的αB%Y.考察 S为费米子作用量，n为味道数，M为费米子矩阵 等式： 在上式中，矩阵的分数次幂不好处理，我们把 （aAr+Barn+yr-i）= (MM)“展开成有理函数的形式，得到如下的表 ag(A+w）r.+βrnty-' (12) 达式： 对应项的系数相等，即可以确定α、%、y以 S,=<e1a+习M+B (4) 及：的演化公式.为了方便应用于下面要考察的共 轭梯度法，假定Bn+y。=1,B%+y=1,则可以确定 但随之而来的问题是在模拟中会出现很多矩阵 移位多项式的参数为 求逆的计算（例如，在式(4）的求和中有很多项，分 母是两个费米子矩阵的乘积和一个常对角矩阵的 a=a (13) 和.在求S的过程中就涉及到这样的一系列矩阵的 求逆)，带来时间和计算资源的大量消耗，限制 B=(R.-0,) (14) RHMC算法的应用. 如何利用其中一次矩阵求逆的结果而较好较快 ”1积 (15) 的得到其他的类似的，仅仅相差一个常数矩阵的矩 (1-Bn）+-1(Bn-wan） 阵的逆矩阵，是值得研究的问题.本文针对共轭梯 2 共轭梯度法 度法方法，参考了文献5]，利用移位多项式得到 一种解决方案. 对于线性系统： Ax=v (16) 1移位多项式 共轭梯度法是：先任取xo,计算r。=v-Axo,让S。= 我们的目的是在己经知道如何求解线性系统 r。然后按照以下方式进行迭代，考察rn如果r。满 Ax=v (5) 足收敛条件，则终止迭代过程，此时所得到的x。为 的条件下，求解线性系统： 所求的解. (A+@)x=v (6) (rnrn） a=- (17) A为矩阵，w为常对角矩阵w=wI(常对角矩阵可 (s,As) 以写成常数和单位矩阵的乘积，在以下的叙述中，用 Xntl=Xa-dnSn (18) 非黑体的ω来表示和单位矩阵相乘的常数)，v为列 Intl=rn+anAs (19) 矢量，x为待求列矢量.相应于Ax=v的系统，假定某 (ra+1,ra+i） bi=(rnra） (20) 个量的迭代演化过程己知，其迭代过程可以写成如 下多项式形式： Sntl=Fmt1+bntiSn (21) r=P (A)ro (7) 3 互相有耦合的两量的迭代 相应于(A+w)x=v的系统，定义移位多项式 (在以下的叙述中，我们利用σ来标记和ω相应的 在共轭梯度法中，有 系数，多项式等量)： rt1=rta As (22) Pg(A+w）=P.(A) (8) Sa=r+baS-1 (23) 求出了系数：，就能了解相应于(A+w)x=v的系统 两个量的迭代耦合在一起，我们希望得到一个量 的相应量的迭代演化，则可以求解系统(A+w)x=. 的迭代演化表达式通过变化，As.-1=(rn-r-）/a。-1 相应于Ax=v的系统，考察如下的迭代过程： 得到 rntl=anArn+Barn+yara-1 (9) 相应地，相应于(A+ω)x=v的系统，相应量的 Tm+1=(1+ anb da-1 r-1 (24) 演化迭代过程为 和式(9) r=(A+o)r+Baro+yr (10) Ent1=Brn+anArn+ynrn-1 (25) 其相应的移位多项式为 对比可知 Pa(A+0)=g0P (A) (11) n=a。 (26) (C)1994-2019 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
30 大 学 物 理 第 38 卷 Sf = ＜φ| ( M+ M) －nf/4 | φ＞ ( 3) Sf 为费米子作用量，nf 为味道数，M 为费米子矩阵． 在上 式 中，矩阵的分数次幂不好处理，我 们 把 ( M+ M) －nf/4 展开成有理函数的形式，得到如下的表 达式: Sf = ＜ φ | α0 + Σn αn M+ M + βn | φ ＞ ( 4) 但随之而来的问题是在模拟中会出现很多矩阵 求逆的计算( 例如，在式( 4) 的求和中有很多项，分 母是两个费米子矩阵的乘积和一个常对角矩阵的 和． 在求 Sf 的过程中就涉及到这样的一系列矩阵的 求逆) ，带来时间和计算资源的大量消耗，限 制 ＲHMC 算法的应用． 如何利用其中一次矩阵求逆的结果而较好较快 的得到其他的类似的，仅仅相差一个常数矩阵的矩 阵的逆矩阵，是值得研究的问题．本文针对 共轭梯 度法 方法，参考了文献［5］，利用移位多项式得到 一种解决方案． 1 移位多项式 我们的目的是在已经知道如何求解线性系统 Ax = v ( 5) 的条件下，求解线性系统: ( A+ω) x = v ( 6) A 为矩阵，ω 为常对角矩阵 ω = ωI ( 常对角矩阵可 以写成常数和单位矩阵的乘积，在以下的叙述中，用 非黑体的 ω 来表示和单位矩阵相乘的常数) ，ν 为列 矢量，x 为待求列矢量．相应于 Ax = v 的系统，假定某 个量的迭代演化过程已知，其迭代过程可以写成如 下多项式形式: rn =Pn( A) r0 ( 7) 相应于( A+ω) x = v 的系统，定义移位多项式 ( 在以下的叙述中，我们利用 σ 来标记和 ω 相应的 系数，多项式等量) : Pσ n ( A+ω) = ξ σ n Pn( A) ( 8) 求出了系数 ξ σ n ，就能了解相应于( A+ω) x = v 的系统 的相应量的迭代演化，则可以求解系统( A+ω) x = v． 相应于 Ax = v 的系统，考察如下的迭代过程: rn+1 =αnArn+βn rn+γn rn－1 ( 9) 相应地，相应于( A+ω) x = v 的系统，相应量的 演化迭代过程为 r σ n+1 =ασ n ( A+ω) r σ n +βσ n r σ n +γσ n r σ n－1 ( 10) 其相应的移位多项式为 Pσ n ( A+ω) = ξ σ n Pn( A) ( 11) 为了定出系数 ξ σ n ，以及相应的 ασ n 、βσ n 、γσ n ．考察 等式: ξ σ n+1( αnArn+βn rn+γn rn－1 ) = ασ n ( A+ω) ξ σ n rn+βσ n ξ σ n rn+γσ n ξ σ n－1 rn－1 ( 12) 对应项的系数相等，即可以确定 ασ n 、βσ n 、γσ n 以 及 ξ σ n 的演化公式．为了方便应用于下面要考察的共 轭梯度法，假定 βn +γn = 1，βσ n +γσ n = 1，则可以确定 移位多项式的参数为 ασ n =αn ξ σ n+1 ξ σ n ( 13) βσ n = ( βn－ωαn ) ξ σ n+1 ξ σ n ( 14) ξ σ n+1 = ξ σ n－1 ξ σ n ξ σ n ( 1－βn ) +ξ σ n－1( βn－ωαn ) ( 15) 2 共轭梯度法 对于线性系统: Ax = v ( 16) 共轭梯度法是: 先任取 x0，计算 r0 = v－Ax0，让 s0 = r0 ． 然后按照以下方式进行迭代，考察 rn ．如果 rn 满 足收敛条件，则终止迭代过程，此时所得到的 xn 为 所求的解． an = － ( rn，rn ) ( sn，Asn ) ( 17) xn+1 = xn－an sn ( 18) rn+1 = rn+anAsn ( 19) bn+1 = ( rn+1，rn+1 ) ( rn，rn ) ( 20) sn+1 = rn+1+bn+1 sn ( 21) 3 互相有耦合的两量的迭代 在共轭梯度法中，有 rn+1 = rn+anAsn ( 22) sn = rn+bn sn－1 ( 23) 两个量的迭代耦合在一起，我们希望得到一个量 的迭代演化表达式．通过变化，Asn－1 = ( rn －rn－1 ) /an－1， 得到 rn+1 = ( 1+ an bn an－1 ) rn+anArn－ an bn an－1 rn－1 ( 24) 和式( 9) rn+1 = βn rn+αnArn+γn rn－1 ( 25) 对比可知 αn = an ( 26)
第1期 吴良凯，等：多个含不同质量项的费米子矩阵求逆算法 31 月=1+86 snt1=r+1+b+1S。 (38) (27) an-1 有了a.、bn,根据表达式(31)可以计算：，再根据式 abn-1-B。 (29)、(30)得到a、bg则有 Y= (28) dn-1 xat1=x8-q9s (39) 利用式(13)、(14)、(15)，我们可以得到、 si+1=ξ+1ra+1+b+1sg (40) B%:的演化表达式.再由式(26)、(27)、(28)，经过 计算可以得到相应的演化表示式： 5结论 dmo. +1 利用移位多项式得到(A+w)x=的演化算法， (29) 该方法省掉大量的矩阵和矢量的乘积运算，能节约 8ag-, 大量的模拟时间，使得RHMC算法在实际上得到应 ba=b (30) Sa-1da-1 用.对于收敛更快，应用矩阵更广的BiCGstab算 8”1aa-1 法因，利用移位多项式也可以得到相应的算法。 041 (31) anbn(1-5)+5-1an-1（1-wan） 参考文献： 4应用移位多项式的共轭梯度法 ]刘川格点量子色动力学].北京：北京大学出版社， 2017:29. 应用式(29)、（30)、(31)，我们可以得到应用 2] 郑继明，朱伟，刘勇，等数值分析].北京：科学出版 移位多项式的共轭梯度法的算法如下：求解系统 社，2016：6097. Ax=v (32) 31 Clark MA,Kennedy A D.Accelerating Staggered Fer- 和 mion Dynamics with the Rational Hybrid Monte Carlo (A+@)x=v (33) (RHMC)Algorithm [J].Physical Review,2007, 先取初始值如下： D75:011502. x0=x6=0,T0=S0=s6=v, ④ Wong K Y,Woloshyn R M.Dynamical Simulations with a-1=1=8=1,b。=b6=0. Highly Improved Staggered Quarks [EB/OL].https:// arxiv.org/abs/0710.0737 有 [5]Jegerlehner B.Krylov Space Solvers for Shifted Linear Sys- (rnr） a.=-(s,As,) (34) tems [EB/OL].http://inspirehep.net/record/427442 [6]Van der Vorst H A.Bi-CGstab:A Fast and Smoothly Xa+1=Xa-dnSn (35) Converging Variant of Bi-CG for the Solution of Nonsym- ruti=rta As (36) metric Linear Systems []SIAM Journal on Scientific （r+1ra+i） and Statistical Computing,1992,13:631. bn+1= (37) (rnrn） Inversion algorithm for fermion matrix with multi-different mass terms WU Liang-kai,GU Xin,LIU Kun,FENG Long-hai (Department of Physics,Jiangsu University,Zhenjiang,Jiangsu 212013,China) Abstract:In the simulation of lattice gauge theory,RHMC algorithm is an accurate algorithm which can be used to any number of fermions.Its central idea is to expand the fermion determinant into rational functions,the im- plementation of this algorithm results in a lot of matrix inversion of fermion matrix,and brings about numerous time and resources consumption.Therefore,its use is limited.This paper uses the shifted polynomial to obtain an ap- proach to solve this problem for the conjugate gradient method.This approach can be applied into RHMC algorithm. Key words:fermions;inverse matrix;lattice gauge theory;shifted polynomial (C)1994-2019 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
第 1 期 吴良凯，等: 多个含不同质量项的费米子矩阵求逆算法 31 βn = 1+ an bn an－1 ( 27) γn = － an bn an－1 = 1－βn ( 28) 利用式( 13) 、( 14) 、( 15) ，我们可以得到 ασ n 、 βσ n 、ξ σ n 的演化表达式．再由式( 26) 、( 27) 、( 28) ，经过 计算可以得到相应的演化表示式: aσ n = an ξ σ n+1 ξ σ n ( 29) bσ n = bn ξ σ n aσ n－1 ξ σ n－1 an－1 ( 30) ξ σ n+1 = ξ σ n ξ σ n－1 an－1 an bn( ξ σ n－1－ξ σ n ) +ξ σ n－1 an－1( 1－ωan ) ( 31) 4 应用移位多项式的共轭梯度法 应用式( 29) 、( 30) 、( 31) ，我们可以得到应用 移位多项式的共轭梯度法的算法如下: 求解系统 Ax = v ( 32) 和 ( A+ω) x = v ( 33) 先取初始值如下: x0 = xσ 0 = 0， r0 = s0 = s σ 0 = v， a－1 = ξ σ －1 = ξ σ 0 = 1， b0 = bσ 0 = 0． 有 an = － ( rn，rn ) ( sn，Asn ) ( 34) xn+1 = xn－an sn ( 35) rn+1 = rn+anAsn ( 36) bn+1 = ( rn+1，rn+1 ) ( rn，rn ) ( 37) sn+1 = rn+1+bn+1 sn ( 38) 有了 an、bn，根据表达式( 31) 可以计算 ξ σ n ，再根据式 ( 29) 、( 30) 得到 aσ n 、bσ n ．则有 xσ n+1 = xσ n －aσ n s σ n ( 39) s σ n+1 = ξ σ n+1 rn+1+bσ n+1 s σ n ( 40) 5 结论 利用移位多项式得到( A+ω) x = v 的演化算法， 该方法省掉大量的矩阵和矢量的乘积运算，能节约 大量的模拟时间，使得 ＲHMC 算法在实际上得到应 用．对 于 收 敛 更 快，应 用 矩 阵 更 广 的 BiCGstab 算 法［6］，利用移位多项式也可以得到相应的算法． 参考文献: ［1］ 刘川．格点量子色动力学［M］．北京: 北京大学出版社， 2017: 2-9． ［2］ 郑继明，朱伟，刘勇，等．数值分析［M］．北京: 科学出版 社，2016: 60-97． ［3］ Clark M A，Kennedy A D． Accelerating Staggered Fermion Dynamics with the Ｒational Hybrid Monte Carlo ( ＲHMC ) Algorithm ［J］． Physical Ｒeview， 2007， D75: 011502． ［4］ Wong K Y，Woloshyn Ｒ M． Dynamical Simulations with Highly Improved Staggered Quarks ［EB/OL］． https: / / arxiv．org /abs/0710．0737 ［5］ Jegerlehner B． Krylov Space Solvers for Shifted Linear Systems ［EB/OL］． http: / /inspirehep．net /record / 427442 ［6］ Van der Vorst H A． Bi －CGstab: A Fast and Smoothly Converging Variant of Bi－CG for the Solution of Nonsymmetric Linear Systems ［J］． SIAM Journal on Scientific and Statistical Computing，1992，13: 631． Inversion algorithm for fermion matrix with multi－different mass terms WU Liang-kai，GU Xin，LIU Kun，FENG Long-hai ( Department of Physics，Jiangsu University，Zhenjiang，Jiangsu 212013，China) Abstract: In the simulation of lattice gauge theory，ＲHMC algorithm is an accurate algorithm which can be used to any number of fermions． Its central idea is to expand the fermion determinant into rational functions，the implementation of this algorithm results in a lot of matrix inversion of fermion matrix，and brings about numerous time and resources consumption． Therefore，its use is limited． This paper uses the shifted polynomial to obtain an approach to solve this problem for the conjugate gradient method． This approach can be applied into ＲHMC algorithm． Key words: fermions; inverse matrix; lattice gauge theory; shifted polynomial