降低浮点精度的隐式蒙特卡洛方法

《Journal of Computational Physics》：Reduced Floating-Point Precision Implicit Monte Carlo

　　摘要：本文提出了利用精度降低的浮点数实现来精确求解热辐射传输问题的算法，采用了隐式蒙特卡洛方法。研究了多种属于算术操作和缩放方法范畴的技术，以评估其提升低精度计算准确性的能力。对比了不同热辐射基准问题在半精度和双精度实现下的计算结果。

引言：热辐射传输为低精度计算提供了一个具有挑战性且有趣的应用领域。由于涉及的非线性物理现象以及物理量跨度极大的特点，使得精确进行低精度计算面临巨大困难。本研究考虑采用精度降低的隐式蒙特卡洛方法来解决辐射传输问题[1]。要在低精度环境下成功应用该算法并得到精确结果，需要运用两类技术：算术操作和缩放方法。算术操作包括使用不同的数值求和算法，以及通过调整运算顺序来处理多个值的乘积。缩放技术可以是静态的，也可以是动态的，其方法包括改变基本单位、基于相对值重新调整距离，以及使用多种经算法选定的能量尺度。所有这些技术的核心目标都是将低精度变量保持在合适的数值范围内，以避免下溢/上溢现象，并最小化计算中的舍入误差。开展此项研究的动机有两个：一是寻找利用低精度浮点数据类型加速科学计算的方法；二是为防止未来计算架构无法支持更高精度浮点数据类型而导致仿真功能丧失做好准备。本研究建立在洛斯阿拉莫斯国家实验室早期关于在低精度浮点数环境中实现隐式蒙特卡洛方法的研究基础之上[2]。具体而言，我们扩展了之前的距离缩放方法，使其能与多种能量尺度配合使用，同时用一种更强大的方法替代了三重乘积运算顺序重排机制，该方法能够处理任意数量的数值和缩放因子。本研究的方法与结果是通过Julia编程语言实现的[3]，相关软件的维护链接为https://github.com/simonbutson/MixedPrecisionIMC.jl，其存档版本则按照MIT许可协议存储在Zenodo平台上，地址为https://zenodo.org/records/20650656[4]。

系统中的热辐射传输取决于吸收/发射材料与辐射场中的光子之间的能量交换。被吸收的光子会为材料增加能量，而发射的光子则会减少材料的能量。热辐射传输所涉及的物理过程可以通过耦合的辐射传输方程和物质能量平衡方程进行数学建模。辐射传输由玻尔兹曼传输方程描述，该方程用于表示光子特定强度作为空间r、角度Ω、时间t以及频率ν的函数：1c?I(r,Ω,ν,t)?t+Ω·?I(r,Ω,ν,t)+σa(r,ν,T)I(r,Ω,ν,t)=σa(r,ν,T)B(ν,T)+Qr(r,ν,t)4π。该方程与物质能量平衡方程相耦合：?Um(r,t,T)?t=∫0∞∫4πσa(r,ν,T)I(r,Ω,ν,t)dΩdν?∫0∞σa(r,ν,T)B(ν,T)dν+Qm(r,t)，其中物质能量密度Um(r,t,T)是根据热容cv(r, T)定义的：Um(r,T)=∫0Tcv(r,T′)dT′，而平衡态辐射能量密度Ur(r,T)则为：Ur(r,T)=aTr(r,t)4，其中a为辐射常数：a=8π5kB415hc3，T和Tr分别代表物质温度和辐射温度。此外，还可能存在非均匀的辐射和物质能量源Qr(r, ν, T)与Qm(r, t)。这些方程还需满足描述系统初始状态的初始条件，以及可一般表示为以下形式的边界条件：I(r,Ω,ν,0)=I0(r,Ω,ν)，T(r,0)=T0(r)，I(r,Ω,ν,t)=IB(r,Ω,ν,t)，其中r属于?V且Ω·n<0。这些方程中的一个重要假设是物质处于“局部热力学平衡”状态，其发射率可用温度为T的黑体来描述。此外，还假设忽略诸如不同离子-光子相互作用、流体动力学运动、热传导以及相对论效应等复杂物理现象[5]。黑体的普朗克谱表达式为：B(ν,T)=2hν3c21exp(hνkBT)?1。在辐射传输应用和高能量密度物理研究中，常用的能量单位是Jk，其值为1×109 J，时间单位则是shakes，值为1×10?8 s。表1列出了热辐射传输计算中使用的各类参数。解决热辐射传输问题的标准方法是由Fleck和Cummings在1971年提出的隐式蒙特卡洛方法[1]。在该方法中，热辐射传输方程被线性化，从而得到一组离散的辐射传输方程和辐射-物质耦合方程。这些方程可通过对在离散时间步长tn内被材料吸收或发射的粒子（即辐射能量束）进行蒙特卡洛模拟来数值求解。隐式蒙特卡洛方法引入了Fleck因子f，该因子从物理意义上对应于一个时间步长内的辐射吸收及各向同性再发射过程。Fleck因子的存在使得有效散射得以实现，因此可以使用更大的时间步长进行模拟，而不会在蒙特卡洛辐射传输模拟中出现不稳定性。隐式蒙特卡洛方法还利用了物质热容与辐射热容之间的比值，该比值为：β(r,t)=?Ur?Um=dUr/dTdUm/dT=cv,r(r,T)cv,m(r,T)=4aT3cv(r,T)。通过参考[1]中的完整推导，可得到Fleck因子的表达式：f=11+αcσa,pnΔtβn，其中α属于[0, 1]区间，是时间步长差分法的隐含参数，σa,p为普朗克加权透射率：σa,p=∫0∞σa(ν,T)B(ν,T)dν∫0∞B(ν,T)dν=∫0∞σa(ν,T)B(ν,T)dνcUr(T)。隐式蒙特卡洛方法形式下的热辐射传输方程为：1c?I(r,Ω,ν,t)?t+Ω·?I(r,Ω,ν,t)+σanI(r,ν,T,t)=χ(ν)σa,pnfUrn+χ(ν)(1?f)(∫0∞∫4πσan(r,ν,T)I(r,Ω,ν,t)dΩ′dν+Qr(r))，而对应的隐式蒙特卡洛方法形式的物质能量平衡方程为：Umn+1(r,T)?Umn(r,T)Δt=1Δt∫tntn+1∫0∞∫4πfσan(r,ν,T)I(r,Ω,ν,t)dΩ′dνdt?fcσa,pUrn(r,T)+fQm(r,t)，其中χ为发射谱：χ(ν,T)=σa(ν,T)B(ν,T)∫0∞σa(ν,T)B(ν,T)dν。虽然上述内容给出了隐式蒙特卡洛方法的通用形式，但本文的研究将采用灰度近似方法，即用普朗克/灰度透射率σa,p替代任何与频率相关的透射率。这种近似方法能够在使用低精度浮点数表示频率相关透射率时避免许多数值上的难题，同时仍能准确模拟各种热辐射物理现象。下文算法1概述了隐式蒙特卡洛算法的运作原理。该算法的基本流程是模拟粒子在离散时间步长内的发射、吸收和散射过程。每个时间步长内被吸收或沉积的能量用于计算更新后的物质温度，这些温度变化又会影响依赖温度的物理量以及下一个时间步长内用于发射粒子的黑体分布。作为一种蒙特卡洛模拟，单个粒子之间的相互作用是通过概率方式建模的，诸如散射碰撞等事件的发生及其结果均由随机抽样决定。像沉积能量这类物理量则是通过统计一个时间步长内发生的多次能量沉积事件来计算的。由于需要精确汇总大量浮点数值（尤其是在低精度情况下），这就催生了2.1小节中讨论的各种技术。关于隐式蒙特卡洛算法流程的更详细说明可见[1]、[5]文献。浮点数是用于计算机计算的实数R的一个子集。由于浮点数据类型是基于固定位数实现的，因此它们只能表示实数连续体中离散的、有限的一部分数值。虽然可以对浮点数执行各种算术运算，但其运算结果通常无法与实数算术的结果完全一致，计算得到的数值往往需要四舍五入到最接近的可表示浮点数值。在本研究中，我们采用IEEE 754标准所规定的浮点算术规则[6]。一个IEEE 754二进制浮点数可由一个符号位s、e个指数位以及d个有效数字位构成。存储指数和有效数字的实际二进制数值分别记为E和D。为了能够表示小于1和大于1的数值，指数部分会加上2e?1?1的偏移量。有效数字部分则表示为1加上一个分数值，该分数值是通过将有效数字位除以2d得到的。从0到d?1的位用于存储有效数字，从d?1到d?1+e的位用于存储指数，而符号位位于最末的d+e位置。一个二进制数可以用bi的表示法来表示，其中i为对应位的位置，其表达式为：(?1)bd+e×2(bd+e?1bd+e?2...bd?1)2?(2e?1?1)×(1.bd?1bd?2...b0)2，相应的浮点数值表达式为：fpv=(?1)s×2(E?(2e?1?1))×(1+D2d)，其中E=∑i=0e?1bd+i2i+1∈{1,2,...,(2e?1)?1}，D=∑i=1dbd?i2?i∈{1,1+2?d,...,2?2?d}，s=bd+e?1。除了标准可表示数值范围之外，IEEE 754浮点数还有几种特殊值，即Inf和NaN[7]。Inf和-Inf标识符代表无穷大值±∞，通常在尝试计算±1/0或±log(0)这类量时会出现这些值。NaN则用于表示如0/0这类结果不确定的计算情况。此外，该标准还规定了非规格化数的概念，这类数允许表示小于最小规格化浮点数的数值，实现方式是允许有效数字的最高位为0，而非通常情况下的1。非规格化数主要用于避免出现这样一种情况：两个浮点数x和y虽然不相等，但它们的差值x?y却小于最小规格化浮点数的数值，此时若直接计算x?y，结果会下溢为0，从而导致错误地认为x=y，进而给代码缺陷的排查带来困难。IEEE 754标准定义了多种不同的二进制数据类型，包括binary16/half、binary32/single、binary64/double以及binary128/quadruple。下表2列出了半精度、单精度和双精度变量相关参数的对比情况：从表中可以看出，半精度浮点数的表示范围明显小于单精度和双精度浮点数。这一限制促使人们采用第3节中讨论的缩放技术。IEEE 754标准旨在保证算术的交换律成立（生成NaN或Inf等特殊值的计算除外），不过实数的其他代数性质，如结合律和分配律，则不一定成立。表3总结了这些性质在实数域R和浮点算术中的适用性。

片段总结：
求和算法：蒙特卡洛方法通过多次随机采样的结果求和，得到各种计数值，这些计数值可用于计算不同的感兴趣的物理量。对于由N个数值xi相加得到的通用计数值Sn，可表示为：Sn=∑i=1Nxi。由于实数具有结合律，在实数域中进行计算时，各个数值的相加顺序并不重要。不过，选择不同的数值求和算法会对计算的准确性和性能产生重要影响。
缩放方法：缩放方法的主要目的是对浮点变量进行数值上的重新调整，从而使使用这些变量进行的计算能够落在浮点数据类型的可表示范围内。这对于处理低精度浮点数所具备的较小数值范围来说尤为重要。缩放方法可以是静态的，比如改变距离或时间等物理量的基本单位；也可以是动态的，即在模拟过程中让缩放因子不断变化。
算法复杂度：尽管诸如成对求和以及运算顺序重排算法之类的技术能够提升低精度计算的准确性，但它们会带来算法复杂度的上升。如果使用低精度浮点数的目标是提高计算性能，那么在选择要采用的技术时就必须权衡精度与性能之间的关系。成对求和方法所带来的额外内存需求在很大程度上...
随机游走加速：随机游走加速是一种技术，它用一次较长的随机游走步骤来替代蒙特卡洛模拟中的大量短距离碰撞。在光学厚度较大的介质中，单位长度内会发生大量碰撞，这种情况使用该技术十分有益。虽然每次碰撞对计算结果的贡献都很小，但为了准确模拟问题，必须计算每一次碰撞。而通过一次整体的随机游走，就可以用少得多的计算量来模拟几乎相同的物理现象。
实验结果：本研究共分析了四个测试问题，每个问题都体现了我们所采用的算术操作和缩放技术的不同方面。需要指出的是，在这些低精度实现方案中，都采用了成对求和方法来提高各项计数值的准确性。如果在对低精度情况采用简单的求和方法，过一段时间后，能量沉积的计数值就会增长到如此程度，以至于后续的加法运算都会因为四舍五入而完全失真，进而阻碍物质状态的进一步变化。
结论：我们阐述并展示了多种技术，这些技术可借助隐式蒙特卡洛方法实现热辐射传输问题的高精度低精度模拟。文中介绍了多种数值求和算法形式的算术操作，包括简单求和、成对求和和卡汉求和，同时还介绍了一种用于重新安排乘法运算顺序的算法流程。此外，还介绍了包括静态单位重标定在内的各类缩放技术。马修·克利夫兰：写作——审稿与编辑、监督、资源协调、项目管理、研究工作、资金筹措。亚历克斯·朗：写作——审稿与编辑、监督、项目管理、研究工作、资金筹措、概念设计。托德·帕尔默：写作——审稿与编辑、监督、资源协调、项目管理、研究工作、资金筹措。

利益冲突声明：西蒙·巴特森声明，其研究得到了洛斯阿拉莫斯国家实验室的财务支持及差旅资助。如有其他作者，他们也声明自己不存在任何可能影响本文所述研究的已知利益冲突或个人关系。

西蒙·巴特森|马修·克利夫兰|亚历克斯·朗|托德·帕尔默