计算机代数与Python的符号转换优化

天文学与计算机42（2023）100666全长产品sympy 2c：从符号表达式到Python中的快速C/C++函数和联合施密特，a，b，b，B。Moserb，C.S. Lorenzb，A. RéfrégierbaScientific IT Services，ETH Zürich，Binzmühlestrasse 130，CH-8092 Zürich，Switzerlandb粒子物理和天体物理研究所，苏黎世联邦理工学院，Wolfgang-Pauli-Strasse 27，CH-8093 Zürich，瑞士ar t i cl e i nf o文章历史记录：2022年3月23日收到2022年11月5日接受2022年11月23日在线提供保留字：代码生成ODEsolver Python计算机代数a b st ra ct计算机代数系统在科学中起着重要的作用，因为它们促进了新的理论模型。所得到的符号方程通常在编译的编程语言中实现，以便为实际应用提供快速和可移植的代码。我们描述了sympy2c，一个新的Python软件包，旨在弥合之间的差距，符号的定义和数值实现的理论模型。sympy 2c将SymPyPython包中实现的符号方程转换为使用符号转换优化的C/C++代码。生成的函数可以方便地用作Python中的扩展模块。sympy 2c在PyCosmoPython包中用于求解爱因斯坦-玻尔兹曼方程，这是一个描述宇宙中线性扰动演化的大型常微分方程系统。在回顾了sympy2c的功能和用法之后，我们描述了它的实现和优化策略。这包括，特别是，一种新的方法来生成优化的ODE求解器，利用符号雅可比矩阵的稀疏性。我们使用爱因斯坦-玻尔兹曼方程作为测试用例来证明其性能。sympy2c是通用的，并且对计算物理的各个领域都有潜在的用处。sympy2c可在https://cosmology.ethz.ch/research/software-lab/sympy2c.html上公开获得。©2022作者（S）。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。1. 介绍计算机代数系统（CAS），例如Mathematica（Inc. ，2021）和SymPy（Meurer et al. ，2017年），在数学和理论物理等科学学科中发挥着重要作用，因为它们促进了新的或修改的大多数情况下，所得到的方程在编译的编程语言中实现，以便为实际应用提供快速、鲁棒在本文中，我们描述了sympy2c，一个新的Python包，旨在弥合理论模型的符号开发和数值实现为此，sympy 2c将PythonCASSymPy中实现的符号方程转换为快速的C/C++代码，然后可以从Python中作为扩展模块使用sympy2c 的 开 发 始 于 计 算 宇 宙 学 领域， 作 为 Python 包PyCosmo 1的一部分（Moseret al. ，2021; Refregier等人，2017;Tarsitano et al. ，2020年）。即时编译器HOPE的概念（Akeretet al. ，2015年）通讯作者：Scientific IT Services，ETH Zürich，Binzmüh- lestrasse 130，CH-8092 Zürich，Switzerland.电子邮件地址：uwe. id.ethz.ch（美国）Schmitt）。1 https://cosmology.ethz.ch/research/software-lab/PyCosmo.htmlhttps://doi.org/10.1016/j.ascom.2022.100666在sympy2c的发展之前。在其他功能中，PyCosmo为爱因斯坦-玻尔兹曼方程提供了一个快速求解器为了改进PyCosmo的代码结构，并使代码生成器可供更广泛的受众使用，我们将代码创建部分与PyCosmo的其他功能分离，从而创建了单独的Python包sympy2c。sympy 2c扩展了SymPy的基本C/C++代码生成功能，例如通过支持特殊函数，数值积分，插值和ODE的数值解。我们特别优化了sympy2c代码生成器，用于具有稀疏雅可比矩阵的高维刚性常微分方程。一般稀疏线性系统的直接求解器是ODE求解器的核心，将在下面详细描述Sympy2c可在http：//cosmology处公开获得。我的天ch/research/software-lab/sympy2c。HTML.本文的组织结构如下。在第2节中，我们介绍了Python在科学编程中的作用，性能相关的方面以及sympy2c在此背景下的作用。在第3节中，我们概述了sympy2c提供的功能以及如何使用该库。我们在第4节中列出了访问sympy2c包的资源，它的源代码和文档。在第5节中，我们讨论了实现和优化细节2213-1337/©2022作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表天文学与计算期刊主页：www.elsevier.com/locate/ascom联合 施密特湾 Moser，C.S. Lorenz等人天文学与计算机42（2023）1006662=sympy2c。在第6节中，我们展示了sympy2c.第7章总结了我们的结论。2. 使用Python进行Python是一种解释型、高级、通用编程语言，专注于可读性和高效编程。如今，它在许多科学学科中发挥着重要作用。它在科学上的成功因素是它的开放源码许可证，它使用C/C++的可扩展性，以及高质量和易于使用的科学软件包的可用性。numpy和scipy等基本软件包在多学科科学期刊上有介绍（Harris et al. ，2020; Virtanen et al. ，2020年）。用于机器学习的Python库，如TensorFlow（Developers，2022），PyTorch（Paszke et al. ，2019）和scikit-learn（Pedregosaet al. ，2011）被广泛使用（Raschka et al. ，2020年）。另一个重要的软件包是SymPy，一个用Python编写的开源计算机代数系统（CAS）SymPy可以直接作为Python库使用，并且不实现自己的编程语言，但可以被视为嵌入式领域特定语言（Hudak，1997）。将符号编译嵌入到Python中允许在Python中扩展SymPy，并且还可以将SymPy与其他Python库一起使用，而不会因跨越语言障碍而造成困难 这种嵌入到Python中的方法对于sympy2c的开发至关重要，并且也被类似的工具所使用，例如devito（Luporini et al. ，2020年）。Python尤其广泛用于天体物理学。 实例包括astropy（The Astropy Collaboration et al. ，2013; Price-Whelan et al. ，2018年），一个用于天文学的Python包，截至3月，Web of Science2上的引用超过5800次2022年，以及事件视界望远镜的数据处理管道（Akiyama etal. ，2019），LIGO天文台（Abbottet al. ，2016）和时空遗产调查（LSST）（Jurićet al. ，2015）。我们还建议读者参考文章（Faes，2012），该文章概述了Python在天文学和科学中的作用。由于Python是一种解释型和动态类型的编程语言，它提供了很大的灵活性，并支持敏捷开发过程。然而，这也意味着在运行时降低速度和更高的内存消耗，这些特性在许多科学应用中都有影响。为了避免这一点，已经开发了提高执行速度的解决方案，并且可以分类如下：用低级语言（如C/C++或Rust ）实现部分代码，或将Python绑定到现有的C/C++代码。例如，大部分numpy（Harris et al. ，2020）和scipy（Virtanen et al. ，2020）由BLAS（Blackford et al. ，2002）和其他已建立的数字库。 简化Python和低级语言之间绑定的著名工具包括Cython （ Behnel et al. ， 2011 ） ， swig （ Beazley ，1996 ） ， pybind11 （ Jakob et al. ， 2017 ） 、 f2py（Peterson，2009）或PyO3。使用HOPE提供的即时编译（Akeretet al. ，2015）、numba（ Lam et al. ， 2015 ） 或 PyPy （ Rigo andPedroni ，2006）。Python代码自动翻译为C/C++，由Pythran支持（Guelton，2018）。2 https://clarivate.com/products/web-of-science/sympy 2c包属于第三类，但与上述工具相反，它从SymPy表达式而不是现有的Python函数创建C/C++代码。SymPy已经能够将基本表达式转换为C代码，但这种功能是有限的。例如，SymPy代码生成不支持特殊函数，也无法生成SymPy无法计算符号反导数的积分代码软件包pyodesys（Dahlgren，2018）遵循类似的方法来生成C/C++代码，用于评估ODE的右侧及其雅可比矩阵。pyodesys将这些函数委托给现有的ODE求解器，例如pygslodeiv2（Bjorn，2020），它不考虑稀疏性。SymPy的autowrap功能允许将表达式编译到不同的后端，如C或FORTRAN。然而，autowrap不生成没有显式符号解的积分或ODE的代码，因此不足以在PyCosmo中使用（Moser et al. ，2021; Refregier等人，2017; Tar-sitano etal. ，2020年）。sympy2c与上述工具不同通过提供一个非常快速的ODE求解器，通过考虑雅可比矩阵的稀疏性，并通过实现数值积分和样条插值的例程3. 功能在本节中，我们描述并演示了sympy2c的主要功能。有关sympy2cAPI的详细文档，请参阅sympy2c在线文档，网址为https://cosmo-docs.phys.ethz.ch/sympy2c/。更多关于sympy2c内部工作的细节将在第5节中介绍。3.1. 功能创建和使用函数的步骤如下：1. 我们通过提供函数应该求值的符号表达式和函数接受的2. 我们声明一个扩展模块并将此函数添加到模块中。3. 我们触发代码生成过程和已创建代码的编译。4. 我们导入编译好的Python扩展模块。5. 我们现在可以从Python调用声明的C-Function作为此模块中的函数。我们在清单1中演示了sympy2c的使用模式，其中我们创建并使用了一个Python扩展模块，该模块包含一个函数来计算给定高度h和半径r的圆柱体V hπ r2的体积：第4行：与Mathematica相反，符号在Python中不是一等公民，因此我们必须声明r和h。第6行：这声明了一个扩展模块。第7行：我们声明了一个名为“volume_cylinder”的函数，它计算h * pi * r ** 2的值，并接受参数r和h。第8行：我们将此函数添加到扩展模块第10行：这里我们触发代码生成和生成代码的编译。我们还将生成的扩展模块作为imported_module导入。第12行：现在，名为volume_cylinder的函数可以作为imported_module模块的一部分使用。·········联合 施密特湾 Moser，C.S. Lorenz等人天文学与计算机42（2023）1006663∫•∞1从sympy2c导入符号，函数，模块2importnumpy asnp 34r， h=符号（“rh”）56模块decl = Module（）7f=函数（“volume_cylinder”，h* np.pi * r ** 2，r，h）8module_decl.add（f）910进口模块= module_decl.compile_and_load（）11打印（12“半径为1，高度为2的圆柱体的体积为”，13import_module.volume_cylinder（1.0，2.0））14）清单1：sympy2c中的函数声明示例。1从sympy2c导入符号，函数，积分，2模块3importnumpy asnp 45t=符号（“t”）67integral= Integral（exp（-（t ** 2）），t，-oo，oo）8module_decl.add（Function（“gauss”，integral））910进口模块= module_decl.compile_and_load（）11高斯积分=imported_module.gauss（）12打印（13“数值误差为”，14abs（gauss_integral - np.sqrt（np.pi）15）清单2：不定积分的数值计算。3.2. 积分sympy 2c的一个重要特性是创建用于积分数值计算的C/C++代码。sympy2c提供了一个函数Integral，它采用一个表达式，积分变量和积分下限和上限的表达式。当使用时，例如在一个sympy 2c函数中，sympy 2c创建C/C++代码来计算数值近似。当反导数的封闭形式未知或无法通过SymPy计算时，这特别有用。在内部，sympy2c调用完善的QUADPACK（Piessens etal. ，2012）GNU科学库（gsl）（Galassi，2009）中提供的例程。这些例程还支持不定积分的计算。函数Integral从sympy2c返回一个符号函数，它可以用作任何其他表达式，但稍后将被转换为一个例程，用于积分的数值近似。清单2显示了如何计算高斯分布的一个示例3.3. 三次样条插值插值函数可以帮助加速数值计算。sympy 2c提供了一个函数InterpolationFunc-tion 1D，它将创建C/C++代码来从gsl调用样条插值。 这样的插值函数也可以用在被积函数中，或者用在前面在3.2节中介绍的数值积分的极限中，或者用在后面在3.4节中介绍的常微分方程的符号表示的右手边。清单3展示了sympy2c如何支持这一点：第6行：这声明了一个具有给定名称的插值函数。它可以用作任何其他符号函数，e.g.sympy.sin.现在，它作为一个积分∞−∞e-t2dt，sympy2c：第7行：我们添加使用插值的函数f函数并计算cos_approx（t ** 2）。• 第6行：这定义了符号积分<$∞e−t2dt。t是• 第12行：编译和导入后，模块包含积分变量和符号oo由SymPy使用代表。第7行：我们添加函数gauss，它计算给定的积分，不带参数。第10行：我们从Python中调用函数gauss，现在它将执行数值积分。一个函数set_cos_approx_values（这个名字是从我们在第6行中指定的名字派生出来的）来设置插值。这里我们提供x i值0的值。0，0。1，. . . ，1. 0.第15行：我们检查插值对实际结果的影响。····联合 施密特湾 Moser，C.S. Lorenz等人天文学与计算机42（2023）10066642=·==1从sympy2c导入（符号，函数，2插值函数1D，模块）3importnumpy asnp 45t=符号（“t”）6cos_approx=InterpolationFunction1D（“cos_approx”）7module_decl.add（Function（“f”，cos_approx（t** 2），t））89 import_module = module_decl.compile_and_load（）1011Xi= np.linspace（0，1.0，11）12import_module.set_cos_approx_values（Xi，np.cos（Xi））13打印（14“x=0.5时的插值误差为”，15abs（imported_module.f（0.5）-np.cos（0.5 ** 2））16）清单3：使用插值函数。3.4. 常微分方程Sympy2c还生成用于刚性和非刚性常微分方程的数值解的有效代码。显式一阶常微分方程的形式为ystec（t）=f（t，y（t）），t≥t0y（t0）=y0∈ Rn.函数f被称为常微分方程的右边，y0是初始条件的向量。这种常微分方程的数值解由在用户指定的时间点ti处的y（ti）的近似组成。计算这种近似的算法通常需要用户提供的函数来评估f和初始条件的向量。一些算法还需要用户提供函数来评估f的雅可比矩阵。sympy2c使用LSODA3 算法（Petzold，1983）。我们选择LSODA是因为它可以检测ODE的刚度，并自动在BDF4之间切换 方法（Curtiss和Hirschfelder，1952）和选择LSODA的另一个原因是，多年来它已被证明是可靠和有效的，这也是为什么它也是scipy的ODE求解器odeint中的默认方法。LSODA实现自适应步长控制，以满足本地误差要求，并确保稳定性。我们使用这种方法来取代现有的PyCosmo手工制作的BDF求解器，以受益于强大而有效的步长控制以及刚度检测和两个LSODA积分器阶数的自动自适应。为了证明这一特征，我们考虑了罗伯逊问题（Robertson，1966;Hairer和Wanner，1996），这是描述浓度为y1，y2和y3的三种反应物的自催化化学反应动力学的刚性方程的一个常见例子。这个问题通常被用作测试问题，以比较刚性常微分方程的求解器。罗伯逊问题的方程如下：ystec1=−k1y1+k3y2y3ystec2=k1y1−k2y2−k3y2y31 LSODA是LSODE（LivermoreS olver for Ordinary ifferentialEquations）的一个变体，具有自动方法切换。2 Backward分化公式ystec3=k2y2，其中反应系数的值由k10给出。04,k23107和k3104。清单4展示了如何使用sympy2c实现这个ODE：1将numpy作为np2fromsympy_to_cimport模块，OdeFast，符号34 y1，y2，y3，t=符号（“y1 y2y3t”）5 k1，k2，k3= 1 e-4，3e7，1 e46y1点=-k1 * y1+ k3 * y2 * y37 y2点=k1 * y1- k2 * y2- k3 * y2 * y38 y3点=k2 * y2 ** 2 910模块decl = Module（）11state_vars= [y1，y2，y3]12rhs= [y1dot，y2dot，y3dot]13module_decl.add（OdeFast（“robertson“，t，14state_vars，rhs））15导入模块= module_decl.compile_and_load（）1617y_start= np.array（[1.0，0.0，0.0]）18tvec= 0.4 * 10 ** np.arange（0，6）19rtol=1e-620 atol= np.array（[1e-8，1e-8，1e-10]）2122结果，诊断=23imported_module.solve_fast_robertson（24y_start，tvec，rtol=rtol，atol=atol 25）26打印（结果）[[1.00000000e+000.00000000e+000.0000000e+00][9.99640065e-013.33213355e-121.20024984e-15] [9.96047827e-01 3.32015942e-121.31485884e-14] [9.60826498e-01 3.20275499e-121.28035090e-13] [6.70346206e-01 2.23448735e-129.17743119e-13] [1.83165556e-02 6.10551854e-141.66613769e-12]第6-8行第13行：使用sympy2c编译ODE求解器的声明。我们提供了一个名字的常微分方程，时间变量，一个列表的状态变量和列表或右手边的表达式的常微分方程。··联合 施密特湾 Moser，C.S. Lorenz等人天文学与计算机42（2023）1006665≥=M=（1）011O=-=⎛⎞--第17-19行：初始值的声明，用于评估的时间网格和LSODA自适应步长控制的公差设置。第21函数名solve_fast_robertson来自第13行声明中使用的ODE名称，并调用LSODA算法和右侧的编译代码，雅可比矩阵和优化的线性方程求解器，稍后将详细介绍4. 使用sympy2c软件包托管在https://pypi.org/project/sympy2c上，因此可以使用pipinstallsympy2c安装。sympy2c的源代码在https上公开托管：//cosmo-gitlab.phys.ethz.ch/cosmo_public/sympy2c 并根据 GPL v3许可。为了减少软件包的大小并避免潜在的许可冲突，sympy2c将在第 一 次 调 用 时 下 载 并 编 译 外 部 C 代 码 ， 例 如 GNU 科 学 库（Galassi，2009）。 文档可在https：//cosmo-docs上找到。phys.ethz.ch/sympy2c/和https://cosmology.ethz.ch/research/software-lab/sympy2c.html5. 执行在本节中，我们将描述symphy2c实现的主要特征。由于快速ODE求解器是sympy2c的核心，也是实现该软件包的主要驱动力之一，因此我们更详细地介绍了其优化策略。5.1. Python扩展模块如第2节所述，Python可以使用扩展模块（Anon，0000）进行扩展。这种方法的目的是提高程序关键部分的性能，或者使用现有的外部C/C++代码。程序员可以使用Cython（Behnel et al. ，2011年），它是Python语言的超集，增加了对变量、函数参数和返回值的可选类型声明的支持。Cython项目提供了将Cython源代码翻译成C/C++的工具，包括与Python解释器交互所需的代码。此外，使用Cython来实现扩展模块有助于避免由于Python对象的引用计数而导致的错误，并且创建的代码可以在所有主要操作系统上编译，并且可以在所有Python 2.6版本中运行，而无需调整原始的Cython源代码。sympy 2c在内部使用Cython及其工具，使生成的C/C++代码可以从Python访问，并创建独立于不同Python版本之间更改操作系统和编译器。我们选择Cython而不是ctypes或cffi，因为它也可以很好地与C++一起工作，并且使用户不必使用单独的C或C++构建系统来编译共享库。5.2. 快速ODE求解器刚性常微分方程的隐式求解器需要在每个时间步求解一个涉及常微分方程右侧雅可比矩阵的线性方程。这些系统通常是稀疏的，其结构取决于常微分方程各组成部分的相互作用。求解稠密线性系统的一个常用工具是LUP分解（Gene和Golub，1996），它将矩阵A分解为A LUP，其中L是下三角矩阵，U是上三角矩阵，P是置换矩阵。尸体腐烂-对于一个大小为n×n的矩阵，时间复杂度为O（n3）。由于P−1=P T 对于置换矩阵，可以通过首先求解LzPTb然后求解Uxz来有效地执行求解Axb。由于L和U都是下三角矩阵（分别是上三角矩阵），所以这两个步骤都只涉及前向（分别是后向）替换。为了解决这样的系统，LSODA使用LAPACK（Andersonetal. ，1999年）例程getrf计算一般矩阵的LUP分解，dgbtrf计算带状矩阵。LSODA不支持其他稀疏矩阵结构。由于运行时的复杂性，（n3），LUP分解占主导地位的所有运行时间的LSODA较大的系统与非带状雅可比矩阵Sympy2c能够通过用专门的和快速的实现替换getrf例程来获得主要的速度改进，该实现考虑了从ODE的符号表示导出的先验已知的稀疏性结构。在sympy2c的开发过程中，我们将我们的原型与直接稀疏求解器SuperLU（Li，2005）进行了比较。SuperLU比LSODA的默认求解器快，但不如我们的方法快。我们还在代码生成器中使用sympy.simplify和sympy.cse函数来消除公共子表达式，以减少浮点运算的总数。5.2.1. 线性解算器中的循环展开作为对代码生成器的介绍，我们忽略LUP分解中的排列，首先关注LU分解。sympy2c展开LU分解计算中出现的循环。为了说明这个想法，我们假设单位矩阵P和n=4的矩阵1X0 0123400x y然后，M的LU分解的第一步由以下嵌套循环组成：对于（k=0;kn- 1;k ++）对于（inti= k+1;in; i++）{ m[i][k]/=m[k][k];for（intj=k+1;jn;j++）m[i][j]-= m[i][k]* m[k][j]}使用矩阵项的专用变量而不是数组，并在代码生成期间展开循环，可以将前面的代码转换为doublem_0 =1.0;doublem_0_1 =x;...m_1_0 /= m_0_0;m_1_1 -= m_1_0 * m_0_1;m_1_2 -= m_1_0 * m_0_2;m_1_3 -= m_1_0 * m_0_3;m_2_0 /= m_0_0;...m_3_2 -= m_3_1 * m_1_2;m_3_3 -= m_3_1 * m_1_3;由于m0，2m0，3m2，1m3，0m3，10，sympy2c减少计算的数量，m_1_0 /= m_0_0;m_1_1 -= m_1_0 * m_0_1;m_2_0 /= m_0_0;m_2_1 -= m_2_0 * m_0_1;··联合 施密特湾 Moser，C.S. Lorenz等人天文学与计算机42（2023）1006666==O2×2（）下一页（）下一页CDX2B2我我=2m_2_1 /= m_1_1;m_2_2 -= m_2_1 * m_1_2;m_2_3 -= m_2_1 * m_1_3;第一个使用for循环和数组的实现需要14个整数加法、14个整数比较和31个浮点运算，而最后一个专门的实现需要只有11个浮点运算，没有整数加法或比较。这个概念与编译器中使用的稀疏条件常数传播（Cooper，20125.2.2. 置换处理行来实现部分旋转是控制数值误差所必需的（Gene和Golub ，1996），并且不能被忽略。 为了检查我们是否需要交换行，LUP求解器检查每行对角线上的相应元素是否比下面同一列中的元素具有更高的幅度 对角线。如果不是这种情况，则交换行检测交换的数学公式是k>i|姆基|>>|MII|.（二）这对我们的方法提出了挑战，因为mij的值在LUP分解期间是最新的，并且不能提前有效地计算。我们在PyCosmo中使用sympy2c来一次又一次地解决具有不同参数的相同系统。因此，我们实施了以下自适应策略来缓解此问题：1. 尝试使用优化的求解器求解线性系统Ax b，而不使用主元。如果需要排列，则专用求解器实现必要的检查，并且如果需要，则回退到通用LUP在这种情况下，我们记录所需的排列。2. 之后，我们为记录的排列Pi生成并编译APT bPT的更优化的求解器。3. 在下一次运行中，如果可行，求解器将在现有的专用求解器之间切换。如果ODE的参数变化需要先前未出现的排列，则我们回到步骤1中的通用求解器并记录排列。4. 继续执行步骤2。我们求解具有变化参数的爱因斯坦-玻尔兹曼方程的实验已经表明，仅需要非常少的所描述的迭代来捕获所涉及的置换。对于MCMC采样中可能出现的物理上非常不可能的参数组合，情况也是如此（Foreman-Mackey et al. ，2013）。为了减少所需的排列数量，我们放宽了上述检查，可配置的安全因子C≥1，5.2.3. 分裂大型系统可以为展开的LUP求解器创建数百万行代码的C/C++函数。这可能会对编译器的优化器造成挑战，从而导致编译时间长和内存消耗高。我们缓解这些问题的方法是使用Schur补将线性系统Mx b1. 我们用方阵A和D（可能有不同的大小）将矩阵M分成块A、B、C、D，并相应地将x和bMx=（A B）（x1）=（b1）（4）2. 使用分块高斯消去法，我们可以分两步求解这个系统：（A−BD−1C）x1=b1−BD−1b2（5）D x2=b2−Cx1（6）代码生成器可以符号化地计算所涉及的矩阵，并最终创建两个较小的C/C++函数，而不是一个较大的函数。sympy2c可以递归地应用这个思想，在生成的代码中创建更多更小的函数另一个好处是，回退通用LUP求解器现在可在较小的系统上运行。因此，在不考虑出现排列的情况下，求解时间受到的影响显著较小：不是求解具有运行时间（n3）的大小为n的系统，而是将大小为n的系统拆分为两个非耦合系统大小为n的函数将运行时间减少到原来的四分之一。这种方法的缺点是减少了旋转：在极端情况下，将大小为n n的矩阵M分裂为n个非耦合求解器，sympy2c根本不会考虑旋转。在我们的实验中，对于合理的块大小，我们没有遇到精度问题。5.3. 代码生成优化由于SymPy是在Python中实现的，因此sympy2c执行的大量符号操作可能会减慢代码生成过程。这尤其适用于较大的ODE系统。另一个影响包含sympy 2c的代码运行时间的因素是生成的C/C++ 代 码 的 编 译 。 为 了 改 善 这 两 种 情 况 下 的 运 行 时 间 ，sympy2c广泛使用基于磁盘的缓存来避免重复计算，从而显著提高后续执行和符号表达式的较小修改的速度。为了提高分裂方法中所需的矩阵的符号求逆的性能，我们通过递归应用以下等式来替换SymPy的逆函数，其中M的大小为n×n，二次矩阵A的大小为n×n，D的大小为k>i|姆基|> C|M II|.（三）尺寸范围这一变化背后的动机是，我们假设，如果受影响的矩阵条目在不同的数量级上不同我们使用C=5或C=A B−1M=C D=A−1+A−1BRCA−1−A−1BR−RCA−1R（七）10在我们的实验中，并观察到减少行交换计算ODE解决方案之间没有显着差异为了进一步加速回退LUP求解器，我们实现了一个变体，该变体考虑了线性系统的符号表示中的先验已知带状结构。这避免了在已知的零条目上循环，但是以在LUP算法中的矩阵更新期间跟踪频带限制为代价。这个功能在我们的应用程序中很有用，但是可以关闭，因为它在给定系统没有绑定的情况下增加了运行时间R=（D-CA-1B）-1。（八）6. 性能我们使用一个特定的高维问题来测量和讨论快速ODE求解器的性能，该ODE的变体和扩展的基准可以在Moser等人 （2021年）。在上面的例子中，为了简单起见，我们忽略了用于排列行的P矩阵。然而，在实践中，联合 施密特湾 Moser，C.S. Lorenz等人天文学与计算机42（2023）1006667++ ×+≈6.1. 设置图1.一、l_max=20的爱因斯坦-玻尔兹曼方程的雅可比矩阵的稀疏模式，零项显示为白色。在LSODA中强制较小的步长，从而导致更长的计算时间。为了测试sympy 2c的性能，我们考虑在PyCosmo中实现的Einstein-Boltzmann 方 程 （ Ma 和 Bertschinger ， 1995;Dodelson，2003）的数值解（Refregier et al. ，2017; Moser etal. ，2021年）。这是一个一阶线性齐次微分方程系统，描述了宇宙中线性扰动的演化。方程组的整体结构为以下形式：y′（t）=J（t）y（t）（9）其中y是扰动场的向量，素数表示关于时间变量t的导数，并且J是时间相关的雅可比矩阵。一般来说，y也取决于波数k，所以这个方程需要求解一个k值的向量这些方程对于精确的宇宙学模型预测至关重要。辐射场、光子因此，ODE系统的大小取决于参数l_max，我们比较不同k∈ {10−5，0}的运行时间。1，1，10}hMpc−1和l_max∈{10，50，100}，得到大小为n∈ {38，158，308}的方程组。6.2. 时间测量表10.1 ，1. 0和10。0 hMpc−1。显示的列为：• n是系统的大小。• 工作模式：– full表示我们禁用了所有优化，LSODA始终使用回退通用LUP求解器。– 带状表示我们禁用了所有优化，并且LSODA使用具有如前所述的雅可比矩阵的带状结构的优化的回退LUP求解器– 优化是指使用展开循环的优化线性求解器，并避免在-旋转已知的零。导致系统的大小为53（l_max 1）。在实际情况下，扰动场y的最终维数可以是几个• TLUP 是在完整LUP线性求解器中花费的总时间这几百除了方程的刚性性质之外，这种大的维度使得它们的快速数值解具有对于这个方程组的理论分析，我们参考Nadkarni-Ghosh and Refregier（2017）.所考虑的方程的矩阵J的稀疏模式如图1所示。可以看到，矩阵是相当稀疏的（5%的非零条目），并具有复杂的结构，使得专门的求解器，例如。不能适用于带状系统我们在Moser等人（2021）中发表了爱因斯坦-玻尔兹曼方程不同变体的时间测量下面，我们将重点展示和比较所讨论的优化的效果。所提出的时间测量是在苏黎世联邦理工学院的高性能计算集群Euler我们分配了一个配备AMD EPYC 7742处理器的完整节点此外，我们为每个配置测量了五次时间，并始终报告这些运行中最快的一次。如上所述，系统取决于波数k。较高的k值会增加溶液中的振荡，当优化的求解器被禁用时（完全和带状模式），或者当优化的求解器遇到未知的行排列并切换到回退LUP求解器时（优化模式），使用求解器在优化模式下报告值0.00表示在求解ODE系统期间不需要回退求解程序。• Toptim是在优化线性求解器中花费的总时间。Ttotal是求解微分方程所需的总时间。除了在求解线性方程中所花费的时间之外，测量为TLUP，最好，这还包括在实际LSODA算法中花费的时间。S是求解常微分方程的加速比。该系数描述了与基线模式完全相比运行时间的减少。表1-的=··联合 施密特湾 Moser，C.S. Lorenz等人天文学与计算机42（2023）100666810=表1k= 0的时间。1 hMpc−1。n ModeTLUP [s]Toptim[s]Ttotal[s]38满4.53e− 2 5.01e−2 1.00带状1.86e− 2 2.32e−2 2.15优化0.00 2.16e− 3 6.76e−3 7.41158满3.28 3.29 1.001.82 1.83 1.80优化0.00 7.30e− 3 1.80e−2 183.36表5表4不同的n的编译时间。n38 58 158T总计[s] 22.5 58.2 157.8Tparse[s] 1.2 14.8 51.0Tsympy2c[s] 5.3 21.8 64.3Tbuild[s] 14.8 20.3 38.9308完整版2.65e1 2.66e1 1.00k = 10−5 hMpc−1时的计时，使用LUP回退求解器。n ModeTLUP [s]Toptim[s]Ttotal[s]1.29e1 1.30e1 2.05表2优化0.00 1.71e−2 3.71e−2 716.66−138满5.22e−3 5.78e−3 1.00带状3.50e− 3 4.40e−3 1.31优化1.65e− 3 6.24e− 5 2.29e−3 2.5358满3.47e− 1 3.49e−1 1.00k=1的时间。0 hMpc。带状1.93e−1 1.95e