Python软件包：解析延拓的Padé近似法和最大熵法

≥=¯≡ → −计算机物理通信282（2023）108519ana_js：Python软件包，用于解析延拓约瑟夫·考夫曼，卡斯滕·赫尔德固体物理研究所，维也纳工业大学，Wiedener Hauptstrasse 8-10，1040 Vienna，AustriaAr t i cl e i nf o a b st r a ct文章历史：2021年5月24日收到2022年8月26日收到修订版，2022年2022年9月6日在线发布关键词：解析延拓Padé最大熵本文给出了用Padé近似法或最大熵法求解费米子和玻色子多体格林函数解析延拓的Python软件包。详细描述了超参数的确定和实现。该代码在GitHub上公开，其中还提供了文档和学习资源节目概要项目名称：ana_jackCPC Library链接到程序文件：https://doi.org/10.17632/vphzykvvf2.1开发人员https://github.com/josefkaufmann/ana_cont编程语言：Python外部依赖性： Python（3. 6），numpy，scipy，matplotlib，h5py，PyQt5，Cython补充材料：测试用例文件、教程和说明问题的性质：从松原频率/虚时间到实频率的相关函数的解析延拓。求解方法：Padé插值法、最大熵法包括限制和不寻常功能的其他注释：最重要的功能可以通过图形用户界面访问。为了获得更大的灵活性，建议使用代码作为一个库，并编写特定于问题的脚本。2022作者出版社：Elsevier B.V.这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。1. 介绍复分析的美妙之处之一是，如果我们知道一个复变量的解析函数，即全纯函数，在一个（开）子域上，我们知道它在任何连通域上。量子场论[ 1 ]的多体格林次这是可能的，因为解析延拓。但是，它最终需要一个分析的延续回到真正的（物理的）频率和时间。以时间传播e-iHt为例，其中哈密顿量H从时间0到t（普朗克常数h 1 ）。复时间e−Hτ的解析延拓，这里也称为威克旋转[1]t iτ，允许我们在与玻尔兹曼算子e−Hβ（β1/T：逆温度;玻尔兹曼常数kB1）相同的基础上处理时间传播不仅由于统一的时间传播，（半）解析计算通常更容易在虚时间或虚频率中公式化数值方法，如plain-vanilla精确对角化或求解虚松原频率的parquet方程[3，26]，避免了实频率轴上的极点和相关的数值不稳定性或离散化误差。一些数值方法✩ 本文的评论由Blum Volker教授安排。本文及其相关的计算机程序可通过ScienceDirect上的计算机物理通信主页（http://www.sciencedirect.com/science/journal/00104655）获得。*通讯作者。电子邮件地址：josefkaufma@gmail.com（J.Kaufmann）。https://doi.org/10.1016/j.cpc.2022.1085190010-4655/2022作者。 由Elsevier B.V.发布，这是CC BY许可证下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。目录可在ScienceDirect计算机物理通信www.elsevier.com/locate/cpcKaufmann和K. 举行计算机物理通信282（2023）1085192..（· ·· ）∫¯≡..∫=+∫===+=+→+例如量子Monte-Carlo模拟[12，43]甚至不能正确地实时工作。在所有这些情况下，一个人需要，在非常结束时，如果计算物理相关的动力学，则返回到实轴的解析延拓为此，两种最先进的方法是Padé近似[32]和最大熵（MaxEnt）方法[16]。 最近的进展参见。Bergeron和Tremblay [2]，Levy等人[25]，Kraberger等人[20]。另外，替代方案，如稀疏建模在文献中讨论了中间表示中的格林本文介绍了用Padé近似[ 32 ]和最大熵（MaxEnt）方法求格林函数解析延拓的教学方法我们讨论了噪声和参数的选择如何影响结果。最后，我们介绍了开源Python包ana_python，讨论了它的实现和使用。纲要手稿的结构如下：我们首先在第2节中对多体格林函数及其解析性质进行了最小限度的介绍然后，在第三节中，我们回顾了Padé插值方法。然后我们展示了概率方法，使用贝叶斯在第5节中，我们将讨论技术方面的问题-MaxEnt的观点，这使我们能够建立稳定的工作流程来确定超参数。然后，在第6节中，我们介绍了我们的实际实现。详细地，我们描述了数值方面，以及如何使用它的自我能源和可伸缩性。 安装程序在第7节中说明。最后，我们将在第8节中结束我们的工作。2. 多体格林相互作用电子系统的许多性质都被编码在延迟的单粒子格林函数中[ 1 ]。G R（t）=−i©（t）{c（t），c<$（0）}。（一）该定义使用了二次量化，其中运算符或计数器（t）是 一 个 选 择 时 间 ， 而 计 数 器（t）则将其消除。 为了简洁起见，我们省略了进一步的指标，如动量或网站在这里，因为对于解析延拓只有一般的分析性质是有关联的Heaviside函数©（t）对于所有时间t<0都是0，并且对于t>0都是1角括号表示在巨正则系综上的平均。在频域中观察延迟格林函数通常更有指导意义GR（ω）=∞dt eiωtG−∞R（t），（二）其中频率ω等于通过设置h1得到的能量。然而，多体系统的计算方法通常采用松原形式[28]，其中相关函数在虚时间或频率轴上计算。虚数时间τ中的费米子单粒子格林G（τ）=−Tτc（τ）c（0）（3）用时间序算子Tτ，通过傅里叶变换，得到松原格林βG（ iωn）=0dτeiωnτG（τ），（4）其中iωn表示费米子的费米子松原频率ωn（2n1）π/β可以表明，格林函数的所有物理信息G（ z）=∞d−∞ωr A（ωr）.z−ωr（五）通过评估Eq.（5）在松原频率z处i ωn，我们得到松原格林函数（4），并且在z处 ω我们得到延迟格林函数，如图1所示。一般来说，Eq。（5）允许计算整个复平面中的格林函数。 还请注意，为了以后的参考，格林（5）当A（ω）有一个紧凑的支持。Sokhotski-Plemelj定理（通常称为Weierstrass公式）f（x）=f（x）lim <$→0+ dx xi<$= −iπ f（0）+Pdx（六）−∞ −∞等式（5）在zω处i <$+，<$0+，并取其虚部，这只留下方程右侧的第一项。（6）、我们得到了重要的关系式A（ω）= −1Im G R（ω）= −1Im G（ω + i 0+）。Kaufmann和K. 举行计算机物理通信282（2023）1085193（七）π πKaufmann和K. 举行计算机物理通信282（2023）1085194==−||+Fig. 1. 简单单极点格林函数G（z）1/（z） ωp）在复频率z的平面上。 背景密度/轮廓仅显示Im G（z）. 极点位置ωp用红色星号标记。 左上角的插图显示了延迟格林函数（沿着复平面上的蓝线计算），其中实线是虚部，虚线是实部。很明显，滞后格林函数的虚部在ω处有一个峰值ωp左下角的插图显示了松原格林函数（在复平面上的绿色十字处它的实部和虚部分别用点和叉画出3. Padé插值在第2节中，我们已经看到松原格林（5）.由于谱密度包含有价值的信息，因此从松原格林函数中提取谱密度是有意义的。如果格林函数是解析给出的考虑到在大多数计算中函数仅用数值计算，则需要不同的方法逻辑上的下一步是通过解析函数来近似格林函数的数值，该解析函数可以通过如上所述的代入在实轴上进行计算。这是在Padé逼近技术中实现的[32]，其中给定的一组数据点由有理函数插值。 Vidberg和Serene [ 42 ]已经提出了一种针对格林函数设计的高效算法虽然经常被称为拟合，但Padé近似实际上是给定数据的插值。如果数据太多点（松原频率），Padé插值的有理函数将有太多的极点，并创建非物理特征，特别是当数据不是很精确时。因此，在实践中，需要减少松原频率的数量Vidberg和Serene [42]表明，使用Padé逼近进行解析延拓效果非常好。 然而，作者也讨论了该方法的一个主要缺点。也就是说，这是决不是明确的先验，松原频率应选择用于构建插值。相反，他们建议比较几个Padé近似，其中不同的松原频率集被使用。当虚轴上的数据点受到随机不确定性或噪声的影响时，Padé逼近的另一个问题出现了。然后，仍然不可能通过在上半平面中没有极点的有理函数来插值这些点。然而，噪声有可能对插值造成严重破坏我们在图2中通过一个简单的四极点测试函数来说明这种行为。我们用一组松原频率来评估它并通过添加来自具有不同数量级的四个不同标准偏差的从精确数据（蓝色符号）构建的Padé近似完美地再现了原始函数，它的极点与解析函数的真实极点（黑盒子）重合添加10−9级的噪声（橙色符号）只是稍微移动了一些极点，使所得的谱函数实际上保持不变。 10−4级的噪声（绿色符号）在实际的数值计算中并不罕见，但已经有了更剧烈的影响。 对于这种噪声水平，在图2（底部）的上半平面上出现了几乎相互抵消的极点和零点对，这是插值所必需的。噪声测试数据。对于这样的噪声水平，谱函数已经有些恶化，见图12。2（顶部）。4. 概率方法-最大熵4.1. 解析延拓问题通常采用的是量子蒙特卡罗模拟，其结果本质上包含一些不确定性。 虽然我们上面已经表明，即使在有噪声的数据的情况下，Padé近似也可以产生良好的结果，但很难评估其质量，Kaufmann和K. 举行计算机物理通信282（2023）1085195+×+×=.（G− K A）N¸2K图2. 具有四个极点（黑色正方形）的测试函数的Padé近似。用于构造近似式的松原频率用红点标记。重构的极点和零点由和符号，分别与颜色编码的噪音水平，见图例。可以看出，高噪音水平，振幅10−4导致实轴上方的极点（）和零点（）对，它们（几乎）彼此抵消（组合成类似于拉斯克的符号）。(For有关图中颜色的解释，请参阅本文的网页版本真实世界的情况下，确切的解决方案是未知的。插值将带有随机噪声的数据视为精确数据。因此，谱函数与真实的基本谱函数的偏差具有任意性由于这个原因，大多数解析延拓方法并不试图构造一个精确插值的解析函数。惟有一个本文试图回答这样一个问题：假设N个正态分布的数据Gn已被测量，并具有相应的不确定度σn。给定一个线性关系，如Eq。（5）在谱函数A（ω）和Gn之间，∞Gndω−∞K（ iωn，ω）A（ω）或者在运算符符号G=KA中，（八）拟合数据的最可能的谱函数A（ω）是什么？如果谱函数为A（ω），则可以立即写下测量G的概率：Pr [G|A]= （2π）−N/2exp−i2（九）.iσii2σi最大似然法[参见， 例如， Press et al. [32]]是基于这样的假设，即这在数值上等同于概率A是真实的谱函数，在G被测量之后（假设我们没有关于A的先验信息）：Pr [A|G]= Pr [G|A]。（十）然后，最可能的谱函数是使方程的指数最大化的函数。或者换句话说，最小化价值函数χ2χ2[A]= .n=1[Gn−dωK（iνn，ω）A（ω）]2n.（十一）这是一个最小二乘问题[32]。为了进一步分析，让我们将积分离散化dω K（iνn，ω）A（ω）<$.K（iνn，ωk）A（ωk）<$ωk.（十σ2Kaufmann和K. 举行计算机物理通信282（2023）1085196二）仅为演示起见，假设测量不确定度σn对于所有n都是常数，并且ω的离散化是均匀的，<$ωk=<$ω。在这种情况下，最小二乘问题被简化为以下线性方程组Kaufmann和K. 举行计算机物理通信282（2023）1085197.=n=[]（| ]=[]（|[] []D（ω）（|]M图3. 解析延拓的核的奇异值分解[Eq.（15）]。 左图：右奇异向量矩阵Vm（ω）（m：奇异值的个数）。 右图：奇 异 值 的 对数图。0 = Gn −Knk Ak<$ω。（十三）K因此，内核决定了问题的4.2. 核的性质将谱函数与松原频率中的格林函数相关联的核K（iν，ω）1iνn−ω.（十四）这里，ω是实轴上的频率，iνn是（费米子或玻色子）松原频率。核的性质最好通过首先离散实频率轴来理解，这样核就变成了矩阵Kjk K（ iνj，ωk）。然而，它不是二次矩阵，因为实频率的数量必须足够大以解决谱函数的所有特征，而松原频率的数量只需要足够大以达到渐近区域。这些准则是相互独立的，因此将核选择为二次矩阵是不合理的。数值上的必要性也常常限制松原频率的数量。奇异值分解也可以很容易地对一般的非二次矩阵进行Kjk =.U jmm V km，（15）其中，U和V是列正交矩阵，并且m是奇异值的向量。在图3中，矩阵VT的虚部与 奇异值一起示出。V的列可以理解为谱函数的基函数[35]。在此基础上，谱函数的第m个分量对松原数据的贡献然后由第m个奇异值加权。值由于只有很少的奇异值是显着的大小，最小化方程中的χ2的问题（11）是病态的，有一个退化解4.3. 正则项或先验概率如果问题有许多完全不同的解决方案，这可能意味着我们没有足够的信息来实际解决问题。然而，这也可能意味着我们没有以正确的方式使用我们关于问题的所有知识。 我们的推理必须利用所有可用的信息，这可以通过贝叶斯定理来实现PrA GPr[A] Pr[G|A]，（16）Pr[G]其中，我们希望找到谱函数A（假设），如果数据G已经被测量，则该谱函数A（假设）使假设为真的条件概率或后验概率PrA G最大化。似然函数Pr G A 已经在Eq中定义了。 （9）.此外，Pr A 和 Pr G 我们现在分别得到了假设和数据的先验概率。当量（16）与Eq.相比更丰富。（十）我们希望贝叶斯统计比最大似然法得到更好的结果。 然而，我们首先必须对谱函数的先验概率PrA进行建模。这种建模不是单一的[32]，而是一种非常成功的艺术状态熵先验的选择是Pr[A]=αS[A]（17）由熵S[A]=<$dω<$A（ω）−D（ω）−A（ω）logA（ω）<$（18Kaufmann和K. 举行计算机物理通信282（2023）1085198）Kaufmann和K. 举行计算机物理通信282（2023）1085199（| ]=22[][]σ[]−=−（|]=（|]Al（|]2nALn相对于默认模型D（ω）。这导致最大熵方法（MaxEnt）。熵的概念由Shannon [34]引入信息论，Jaynes [17]建立了最大熵原理作为统计力学中的“一种推理方法”。最大熵方法非常成功地用于光学数据的反卷积[8，15]。在进一步发展和应用于图像处理[13]之后，它也被应用于解析延拓问题[38，37，16]。合并等式（16）Eq.（9），Eq.（17），Eq.（18）我们必须最大化概率哪里PrA G1Pr[G]e-Qα[A]，（19）Qα[A]=1×2[A]−αS[A]（20）具有尚未指定的标度超参数α>0。数据的先验概率PrG（也称为证据）是固定的，因为我们正在使用一个测量数据集。因此，它不依赖于谱函数，并且在我们的考虑中可以在Pr [ A]的归一化中吸收|G]。而不是最小二乘问题的方程。（11）我们现在面临泛函Qα A的最小化。熵，这是作为谱函数的先验概率，现在作为一个正则化项病态最小二乘问题。通过观察Q相对于A的海森矩阵，这变得显而易见：Qα=.Knl<$ωlKnm<$ωm+αδlm<$ωl，（21）在此通过选择实频率轴的离散化并通过在这些离散点。的第一项Eq。（21）包含一个乘积KT K，因此它的条件甚至比K的条件更差：它的大多数特征值实际上是零。然而，第二项的Eq。（21）增加一个正对角矩阵，按α缩放。这具有增加特征值的效果，即，Q在其最小值处的曲率，因此更容易实际定位最小值。让我们注意到，除了熵，我们还可以选择不同的正则化项。Otsuki等人。[31]介绍了所谓的稀疏建模技术，如果解决方案在内核的奇异空间中是稀疏的，则优先考虑它们。 隐式正则化效应可以在随机采样方法中实现[33，29，30，10]。5. 最大熵方法5.1. 超参数α的选择我们已经看到，当α>0时，泛函Qα A可以最小化，现在我们必须确定如何实际选择一个好的α值。这种选择的重要性如图4所示。我们选择不对称双峰光谱，如黑色曲线所示在右侧面板中。由等式（14）将其变换为松原频率，并且噪声具有10 − 4的幅度 添加了。假设一个默认模型，我们计算最小化方程的谱函数A（ω）（20）对于大范围的不同α值。α的影响在图4左上角的密度图中显示。在非常高的α值时，我们恢复了常数默认值模型中，没有注意到数据的干扰。随着α的减小，谱权集中在中间，然后我们观察到两个不同的峰的形成。当α值较小时，峰分裂，最后我们得到六个非常尖锐的峰，而不是两个宽峰一个这些峰的数量和位置主要由数据的噪声决定。在我们的测试案例中，我们可以产生噪音，使用不同的随机种子，这将导致在低α极限下看起来完全不同的频谱由于解析延拓的结果在很大程度上受α的控制，因此合理地自动选择α是非常重要的。已经提出了各种方法来解决这个问题。在所谓的历史MaxEnt中，α的选择方式是χ2-偏差（等式2）。（11））近似等于虚轴上的数据点数量[14]。然而，它被发现低于数据，Gull [14]和Skilling [40]都认为它有点特殊。经典的MaxEnt将贝叶斯推理方案扩展到α和D[40，16]：e−Qα [A]Pr [A，α|G，D]-Pr [A|G，D，α] Pr [α]，（二十二）α其中尺度不变Jeffreys随后，依赖于方程中的谱函数。（22）的解，得到PrαG， D的显式表达式。在图4的左下图中，后验概率PrαG，D 显示为橙色曲线。经典MaxEnt采用最可能值处的谱函数α是问题的最佳解。该值通常通过设置<$PrαG， D/<$α 0来获得，经过几次近似后，可以将其解析地转换为方程Ng2αS。Ng被称为在图4的左下角，我们用一条蓝线表示Ng/2αS;它与1的交点是经典的 MaxEnt解，用一条蓝色的垂直线标记相应的光谱如右图所示与之密切相关的是BryanA（ω）=α A（α，ω）Pr [α|G、D]（23）结果通常非常接近经典的MaxEnt，因为α的概率是尖锐的峰值，参见图4中的蓝色和橙色虚线。可以清楚地看到，经典和Bryan的MaxEnt的解决方案可能会导致非物理的额外峰值-这已经是一个缺点[14]《明史》卷Kaufmann和K. 举行计算机物理通信282（2023）10851910→=→∞[][=]= − + −= − = ∈ []图4. 确定α的方法。模型频谱（右图中的黑色曲线）被转换为松原频率，并添加幅度为10 − 4的正态分布噪声。左上图：通过针对α值范围最小化MaxEnt泛函Q α获得的光谱函数A（颜色代码）。左下图：说明了确定α最佳值的四种不同方法。在这里，红色和蓝色垂直线分别标志着历史和经典的MaxEnt。颜色在右侧面板的图例中）。通过“chi2kink”方法获得的α值用绿色垂直线标记。请注意，Bryan的方法在因 为 它是α范围内的加权平均值（在这种情况下，作为权重的概率以橙色绘制）。右下图：与真实（初始）光谱相比，通过四种不同方法获得的然而，还有另一种方法来决定最佳α。它基于所谓的L曲线标准[23，24]。在MaxEnt的背景下，据我们所知，它首先由Bergeron和Tremblay提出并解释[2]，并在下文中由Kraberger等人应用。[20]. 它比经典或Bryan的MaxEnt更具启发性历史MaxEnt）。为了理解这种方法，我们将χ2绘制为图1左下角的绿色曲线第四章在α的极限它会变成一个恒定的高值，而对于α0，则会接近另一个较低值的常数。根据Bergeron和Tremblay的推理[2]，对这种行为有一个明确的解释。在高α极限下，（大）χ2偏差的权重可以忽略不计，默认模型D使Qα A最小化。我们看到的是χ2AD。另一方面，在低α极限下，熵的贡献减少到强制谱的正性;χ2慢慢接近其全局最小值。低α区域（其中χ2相对较小）和随后的χ2急剧增加的区域被称为噪声拟合和信息拟合区域[2]。α的最佳值位于信息拟合和噪声拟合之间的过渡区域，其中logχ2（ logα）具有扭结。 因此，我们称这种确定α的方法为这个最优α由曲率χrr（α）的最大值表示[2，20]。一个数值上更稳定和更灵活的方法是拟合一个函数Bφ（x;a，b，c，d）=a+1+e-d（x-c）（24）log10× 2（log10α）的几个值。方 程 中 φ 的曲率 （24）在 x 处 最大Cln（2103）/dC1. 317/d。由于这通常导致强下拟合，我们建议使用更一般的xCf/d或α10 c−f/d，f的优选值二二5 .图4所对应于F2的值由绿色垂直线标记。总而言之，历史上的MaxEnt似乎只与经典MaxEnt相匹配，因为后者往往严重过度匹配。Bergeron和Tremblay [2]注意到，α的后验概率公式是一个近似值，只有在良好的违约情况下才有效。模型使用。为了使用经典MaxEnt获得良好的结果，因此，通过手动选择的因子来增加标准差σn以将概率峰值移动到更大的α值已成为常见的做法。由于重新缩放因子不是先验已知的（尽管注意到对实际噪声幅度的依赖性），该过程引入了附加程度的任意性。将最佳α确定为噪声拟合开始和信息拟合结束的点（即使相对于历史MaxEnt，α值也更大。尽管缺点是不太“贝叶斯”，它不依赖于额外的近似和噪声幅度的数据必须知道只有一个常数前因子。唯一剩下的任意性来源是上面提到的参数f，它控制着α的实际值。然而，它仅限于一个相当窄范围，而经典MaxEnt中的误差缩放因子可以变化几个数量级因此，至少在当前解析延拓的背景下，我们确信确定α的最佳方法是将其提取为噪声拟合和信息拟合之间的边界5.2. 消除非物理特征：预模糊即使我们现在能够找到α的在图的右侧面板4我们可以看到所有解在ω=0附近都有一个扩大的曲率，并且在降低α时出现边峰。Kaufmann和K. 举行计算机物理通信282（2023）10851911dω=∫≈=====αb2B在MaxEnt在图像处理中的原始应用中已经观察到类似的趋势，其中它表现为巨大的 强度分布中的峰值。为了解决这个问题，Skilling [41]表明熵不需要被评估直接从谱函数，但是从以线性方式与谱函数相关的A= ▲h。（二十五）在标准MaxEnt中，▲只是一个单位。虽然最初要求▲是正交矩阵，但Skilling [41]随后假设高斯卷积具有很好的模糊效果，称之为在后一种情况下，▲是一个矩阵，其中包含一个高斯函数gb（ω）=exp[−ω2/ 2b2]/ 2πb，即▲ij=gb（ωi−ωrj）。在解析延拓的背景下，preblur是由Kraberger et al.[20]，尽管Sandvik [33]已经使用了类似的平滑算法。Eq.的泛函（20），我们必须最小化，因此成为Q[h]=1×2[gh]− α S[h]。（二十六）在α之后，我们现在有第二个超参数b，它控制高斯g的宽度。现在让我们仔细看看等式的χ2项。（26）.它将拟合“KA“与数据G的差相加 解释fit函数现在是K A∞K i（νn，ω） A（ω）−∞∞dω K（ iνn，ω）−∞∞dωrg−∞b（ω-ωr） h（ωr）=∞∞rK irg−勒赫dωdω（νn，ω）b（ω ω）（ω）−∞ −∞我知道（二十七）这意味着隐藏的谱函数以与谱函数类似的方式与松原格林函数相关，唯一的区别是“blurrred“ker n elK。通过使函数Q最小化并稍微修改k，其然后产生隐藏谱函数h。谱函数A本身随后通过h与与上面的高斯gb（ω）虽然Q的最小化本身不会因为引入预模糊而变得更加困难，但我们现在面临的问题是模糊宽度b不是先验已知的，类似于熵缩放参数α。然而，我们知道可以通过分析χ2（α）来确定α的一个好值，现在值得分析作为α和b的函数的χ2偏差。这是在左上面板中完成的，图5.在b的大范围内，χ2仅取决于α。仅在尖锐边界（b处）之后 在我们的例子中是1），χ2偏差随着b急剧增加，这可以在图5的左下图中最好地看到。换句话说，将b增加到某个值不会改变拟合的质量（但可能会改变结果）。只是后来突然不可能得到一个很好的配合。一个启发式的解释这就是，将核与高斯卷积只允许频谱中具有最小宽度2b的峰值。一旦超过这一观察得到了证实，我们发现， 在图5中，最大允许b接近1，实际光谱（黑色曲线）由两个具有标准差的高斯曲线组成1.图5左上角的观察结果可用于提出以下确定超参数的算法：确定不同b值的最佳α，从0开始递增。接受b的最后一个值，其中χ2在最佳α小于，例如，1.5乘以b的最优α时的χ20. α的最佳值在图中以红色虚线标记，图5的左上角，由此确定的（α，b）用红色星号标记。由于α的最优值几乎不改变，b 的合理范围，也可以使用简化的算法：确定b的最佳α0，然后增加b，始终保持相同的α。如果χ2超过，例如， 1.5乘以它在b0处的值，接受在此之前b的最后一个值作为最佳值。 相应的（α，b）在图中标出5绿色的星星正如我们在右上图中所看到的，通过这两个（α，b）的解析延拓得到的相应的红色和绿色曲线实际上是不可区分的。在右下图中，绘制了这些重建光谱的差异，并显示，与实际结果的偏差现在令人满意地小。让我们顺便注意到，上述围绕ω=0的非物理摆动在与我们的示例情况类似的谱函数中最为突出。在谱函数在ω=0处显示间隙或峰值的情况下，通常不需要预模糊5.3. 非对角元素上述考虑是针对具有有限范数的正定谱函数进行的，通常为1：dω A（ω）1> 0。这个性质由相关函数的对角元素来填充。然而在某些情况下，需要进行解析延拓非对角元素，具有零范数 dω A（ω）0，作为费米子反对易关系的直接结果一消失范数与正定函数不一致，因此谱函数对于非对角线元素显然，熵的定义，方程。 （18），必须适应这种新的情况[20]。 在这一点上，让我们提到核K（ iωn，ω）对于对角元素和非对角元素是相同的Kaufmann和K. 举行计算机物理通信282（2023）10851912=D（ω）2+− −2.[]=，图5. 分析前模糊参数b。左上图：MaxEnt fit的χ 2偏差（颜色代码和等值线）作为超参数b和α的函数。的红色虚线表示每个b处的最佳α，绿色虚线表示b处的最佳α0比较。左下图：沿红色虚线和绿色虚线的χ 2 左上角的虚线。右上图：真实的光谱函数（黑色），以及左上图中橙色、红色和绿色斜线标记的点处的MaxEnt重建，以各自的颜色绘制右下图显示了这些重建与真实光谱的差异 橙色的A（ω）（对应于b=0）大到10−2。对于非对角分量，让我们将谱写为：A（ω）=A+（ω）−A−（ω），（28）其中A+和A−都是正定的，并且具有相同的范数。那么我们可以把熵写成1S[A+，A−]=<$d ω<$A+（ω）−D（ω）−A+（ω）logA+（ω）A−（ω）D（ω）A−（ω）logA−（ω）（29）D（ω）不把A+、A-和A看作完全独立的，我们可以使用等式。（28）消去A−，写出Qα [A，A+]= 1 × 2[A]− α S[A，A+]。（三十）由于我们正在寻找关于A和A+的Q的最小值，我们可以使用<$Qα[ A， A+]= −α<$ S[ A， A+]=0（31）A+也可以去掉A+。在等式中代入A−=A−A+（29）并取导数，A+和A−可以由A表示为：A±[A]=ΔA2+4D2±A将其插入Eq中。（29）我们得出了所谓的正负熵（三十二）SAdωA（ω）log A+[A]（ω）。（三十三）D（ω）作为默认模型的合理选择，Kraberger等人。[20]建议以以下方式包含频谱的对角元素：Dllr =AllAlrlr + ε。（三Kaufmann和K. 举行计算机物理通信282（2023）10851913十四）在我们的实现中，没有必要添加一个小的数字ε。此外，已经开发了将矩阵值格林函数作为整体而不是分量继续的方法1 通常，会引入两个默认模型D+和D−，并假设它们在以后的时间点上相等。为了简洁起见，我们省略了这一步。Kaufmann和K. 举行计算机物理通信282（2023）10851914表1用于初始化AnalyticContinuationProblem的关键字参数。关键字值im_axisnumpy array，shape（Nmats，）re_axisnumpy数组，shape（Nreal，）im_datanumpy数组、形状（N个垫子、）、几何形状或复杂形状kernel_mode'freq_fermionic'，'freq_bosonic'，'time_fermionic'，'time_bosonic'其中之一表2使用MaxEnt时AnalyticContinuationProblem.solve（）的可能关键字参数。关键字值方法optimizerone ofstdev错误条：numpy数组，shape（Nmats，）协方差矩阵：numpy数组，形状（Nmats，Nmats，）alpha_start默认值为109alpha_end重复，默认值为10−3alpha_div默认值为10fit_position调整，默认2.5interactiveboolean，defaultFalsepreblurboolean，defaultFalse模糊宽度调整6. 用于分析延续的Python包：ana_js在前面的章节中，我们已经介绍了两种解析延拓方法，Padé逼近和最大熵方法。虽然这两种方法都是成熟的，并且之前已经实现过[例如Bergeron和Tremblay [2]，Levy等人[25]，Kraberger等人[20]]，但我们在这里引入另一种解析延拓代码。它是作为Python包实现的。这种方法的优点是其增加的可伸缩性。做解析延拓的方法太多了，一些命令行参数。可以在https://github.com/josefkaufmann/ana_cont上以开放源代码的方式获得ana_sample包。让我们在这里也提供一个简要的概述其他解析连续码。▲Maxent（https://github.com/amstremblay/OmegaMaxEnt）[2]是一个带有Python脚本的C++代码，它是第一个使用L形准则（见5.1节）进行最大熵方法的代码。TRIQS/maxent（https://github.com/TRIQS/maxent）[21]是一个用Python编写的TRIQS应用程序，使用最大熵方法并允许矩阵值光谱。w2dynamics（https://github.com/w2dynamics/w2dynamics）[43]是一个连续时间量子蒙特卡罗代码，它还包含一个用FORTRAN 编写的普通最大熵代码及其Python包装器。SOM（https://github.com/krivenko/som）[22]是一个基于TRIQS的随机优化代码，用于C++编写的解析延拓。Spektra（https://www.spektra.app）[10，11]是一个在线工具，用于分析延拓，它使用平均频谱方法，通过这些拟合虚时间或频率数据的程度来加权所有实频谱。• SpM（https://github.com/SpM-lab/SpM）[45]使用稀疏建模进行解析延拓，并用C++编写6.1. 封装结构让我们简要概述一下一揽子计划的结构。但是，详细的API文档托管在josefkaufmann.github.io/ana_cont上。分析的中心类是AnalyticContinuationProblem，它包含特定于问题的信息：实频网格、虚频或时间网格、虚轴上的数据、逆温度和核的类型，见表1。 它有一个方法solve，它有几个关键字参数。最重要的是方法。目前，只实现了两个方法，pade和maxent_svd。调用solve将创建求解器对象的实例。在method这个类包含了解决前面定义的解析延拓问题的所有必要函数。表2中列出了可能的关键字参数。我们手头上的延拓问题的内核是由Kernel类提供的。它处理内核的实频离散化、预模糊和旋转到协方差矩阵的特征基（如果提供）。此外，还有一类格林函数，其主要目的是通过应用Kramers-Kronig关系式Eq.，从给定的谱中构造出全复值格林（5）.我们要强调的是，这个包不包含一个“主程序”，可以通过指定一些参数访问所有功能。虽然这听起来很令人沮丧，但我们相信，实际上这减少了用户和维护人员所需的工作量。在用户端，一个工作脚本可以由10行Python代码组成;示例脚本可以在 GitHub repository.提供坐标轴数据的方式完全自由：它们可以从最喜欢的文件格式中读取，也可以在线计算。计算出的实轴数据可以绘制或进一步处理，而无需绕道file-IO。因此，目前的解析连续码可以很容易地嵌入到复杂的后处理程序。同样对于开发商来说这个原则也是一个很大的·····Kaufmann和K. 举行计算机物理通信282（2023）10851915≡[]..Fi F（ωi）. 频率空间上的积分通过下式转换为和：以ωi为中心。dω F（ω）iFii，其中i是频率间隔伊苏岛MKLKL这是因为它们不必考虑复杂的用户界面和流控制。由于解析延拓的知识和经验对于使用任何解析延拓程序都是绝对必要的，我们看不到这种方法的主要缺点然而，对于临时用户，我们已经在脚本中投射了一些标准过程，用于使用图形用户界面（GUI）进行分析。因此，MaxEnt和Padé对玻色子和费米子格林函数的大多数标准解析延拓即使在没有Python知识的情况下也可以完成。GUI的教程可以在GitHub上的存储库文档中找到：https：//josefkaufmann。github.io/ana_cont6.2. Numerics解析延拓是，相比大型量子蒙特卡罗模拟或贝特-萨尔彼得方程反演，一个数值上相当便宜的任务。因此，代码的性能不必是主要目标。ana_blog包完全用Python编写，部分地牺牲了性能以换取可扩展性和可读性，同时仍然保留了后期优化的可能性。然而，为了使代码合理快速，已经采取了一些措施6.2.1.频域离散化对于数值计算，必须将y轴重新定义为：。ω→ωi。对于任意函数F（ω），我们定义6.2.2.奇异值分解MaxEnt最重要的步骤之一是在谱函数空间中最小化泛函Qα A。 如果将实频率轴离散为1000个区间，那么解的空间就有1000维。这对于确定性求解器来说太大了，需要蒙特卡罗方法[33]。然而，Jobell和Gubernatis [16]执行内核的奇异值分解，从而实现了维度的大幅减少，例如在典型情况下从1000减少到20。 然后，优化问题可以通过确定性方法来解决，如牛顿求根或Levenberg-Margaret算法。MaxentSolverSVD遵循核的奇异值分解策略，如等式中所定义。 （15）.截断奇异值的向量，使得仅大于某个阈值（在我们的例子中为10-10）的值将被保留。在奇异空间中，频谱通过以下参数化：Aj = D j exp.V jmum（35）如Jobell和Gubernatis所提出的[16]。6.2.3.最小化问题上述定义被插入到Eq。（20），经过几个中间步骤，平稳条件<$Qα/<$ Am=0导致fm（u）<$αum+<$m。EKUKM。.KklAl（u）l−Gk=0，（36）这个问题必须为你解决。当根函数的导数解析已知时，求根算法工作得更好。因此，我们对f对u求导，得到雅可比矩阵J：fmUkm Ek.Kkl lAl（u）Vli（37）表达式Eq.（36）Eq.（37）可以通过以下定义以更有效的形式进行数值评估：wlexp. .Vlmum（38a）WmlEkUkmmUknnVlnlDl（38b）公里WmilWmlVli（38c）BmUkmmEkGk（38天）KMKaufmann和K. 举行计算机物理通信282（2023）10851916.