logo标签列表

如何使用SymPy加速长函数的符号导数计算？

pythonoptimizationsympyintegralderivative
10

我正在用Python编写程序，使用自由ICI方法（目前是SICI方法...但最终将变成自由ICI）来解决薛定谔方程。如果您对此不熟悉，那是因为这个主题几乎没有信息可供参考，并且完全没有可用的样本代码。

这个过程涉及到迭代地得出偏微分方程的解。在这个过程中，需要执行许多符号求导操作。问题在于，随着程序运行，需要进行求导的函数不断变得越来越大，以至于到第五次迭代时，计算符号导数需要耗费大量时间。

我需要加快速度，因为我希望能够完成至少30次迭代，并且在退休之前实现这一点。

我已经删除了不必要的重复计算（或者至少我所知道的那些），这有很大帮助。除此之外，我完全不知道如何加速。

下面是包含正在计算导数的函数的代码（inf_integrate函数只是复合Simpson方法，因为它比使用SymPy的integrate方法更快，并且不会由于振荡函数而引发错误）：

from sympy import *


def inf_integrate(fun, n, a, b):
    f = lambdify(r, fun)
    h = (b-a)/n
    XI0 = f(a) + f(b)
    XI1 = 0
    XI2 = 0

    for i in range(1, n):
        X = a + i*h

        if i % 2 == 0:
            XI2 = XI2 + f(X)
        else:
            XI1 = XI1 + f(X)

    XI = h*(XI0 + 2*XI2 + 4*XI1)/3

    return XI


r = symbols('r')

def H(fun):
    return (-1/2)*diff(fun, r, 2) - (1/r)*diff(fun, r) - (1/r)*fun

E1 = symbols('E1')
low = 10**(-5)
high = 40
n = 5000

g = Lambda(r, r)


psi0 = Lambda(r, exp(-1.5*r))

I1 = inf_integrate(4*pi*(r**2)*psi0(r)*H(psi0(r)), n, low, high)
I2 = inf_integrate(4*pi*(r**2)*psi0(r)*psi0(r), n, low, high)

E0 = I1/I2
print(E0)

for x in range(10):

    f1 = Lambda(r, psi0(r))
    f2 = Lambda(r, g(r)*(H(psi0(r)) - E0*psi0(r)))
    Hf1 = Lambda(r, H(f1(r)))
    Hf2 = Lambda(r, H(f2(r)))

    H11 = inf_integrate(4*pi*(r**2)*f1(r)*Hf1(r), n, low, high)
    H12 = inf_integrate(4*pi*(r**2)*f1(r)*Hf2(r), n, low, high)
    H21 = inf_integrate(4*pi*(r**2)*f2(r)*Hf1(r), n, low, high)
    H22 = inf_integrate(4*pi*(r**2)*f2(r)*Hf2(r), n, low, high)

    S11 = inf_integrate(4*pi*(r**2)*f1(r)*f1(r), n, low, high)
    S12 = inf_integrate(4*pi*(r**2)*f1(r)*f2(r), n, low, high)
    S21 = S12
    S22 = inf_integrate(4*pi*(r**2)*f2(r)*f2(r), n, low, high)

    eqn = Lambda(E1, (H11 - E1*S11)*(H22 - E1*S22) - (H12 - E1*S12)*(H21 - E1*S21))

    roots = solve(eqn(E1), E1)

    E0 = roots[0]

    C = -(H11 - E0*S11)/(H12 - E0*S12)

    psi0 = Lambda(r, f1(r) + C*f2(r))

    print(E0)

程序可以正常工作并收敛到预期的结果，但速度太慢。非常感谢任何有助于加快速度的帮助。

好的。谢谢您的清晰解释。我现在开始看到问题所在了：成千上万的术语正在呈组合爆炸式增长。\n这是用于什么目的 - 需要多精确（为了优化目的）？ - EchoCloud
2个回答
9

这里有几件事情可以做：

  1. 如果对代码进行分析，您将注意到大部分时间都花在集成函数 inf_integrate 中，主要是因为使用了手动的Python循环。可以通过将参数转换为矢量化函数并使用SciPy的集成常规（已编译因此速度更快）来改善。

  2. 当您使用嵌套符号表达式时，检查偶尔的显式简化是否有助于控制爆炸性的复杂性，这似乎在这里成立。

  3. 定义的所有 Lamda 函数都不需要。可以简化表达式。我还没有检查这实际上是否影响运行时间，但肯定有助于下一步工作（因为SymEngine还没有 Lambda ）。

  4. 使用SymEngine而不是SymPy。 SymPy（截至目前）仅基于Python，因此速度较慢。 SymEngine正在开发其已编译的核心，并且可以更快。它几乎具有您所需的所有功能。

  5. 每一步都解决一个本质上不会改变其特性的方程：始终是相同的二次方程，只有系数发生变化。通过一般地解决这个问题，您可以节省大量时间，特别是SymPy不必处理复杂系数。

综上所述，我得出以下结论：

from symengine import *
import sympy
from scipy.integrate import trapz
import numpy as np

r, E1 = symbols('r, E1')
H11, H22, H12, H21 = symbols("H11, H22, H12, H21")
S11, S22, S12, S21 = symbols("S11, S22, S12, S21")
low = 1e-5
high = 40
n = 5000

quadratic_expression = (H11-E1*S11)*(H22-E1*S22)-(H12-E1*S12)*(H21-E1*S21)
general_solution = sympify( sympy.solve(quadratic_expression,E1)[0] )
def solve_quadratic(**kwargs):
    return general_solution.subs(kwargs)

sampling_points = np.linspace(low,high,n)
def inf_integrate(fun):
    f = lambdify([r],[fun])
    values = f(sampling_points)
    return trapz(values,sampling_points)

def H(fun):
    return -fun.diff(r,2)/2 - fun.diff(r)/r - fun/r

psi0 = exp(-3*r/2)
I1 = inf_integrate(4*pi*(r**2)*psi0*H(psi0))
I2 = inf_integrate(4*pi*(r**2)*psi0**2)
E0 = I1/I2
print(E0)

for x in range(30):
    f1 = psi0
    f2 = r * (H(psi0)-E0*psi0)
    Hf1 = H(f1)
    Hf2 = H(f2)

    H11 = inf_integrate( 4*pi*(r**2)*f1*Hf1 )
    H12 = inf_integrate( 4*pi*(r**2)*f1*Hf2 )
    H21 = inf_integrate( 4*pi*(r**2)*f2*Hf1 )
    H22 = inf_integrate( 4*pi*(r**2)*f2*Hf2 )

    S11 = inf_integrate( 4*pi*(r**2)*f1**2 )
    S12 = inf_integrate( 4*pi*(r**2)*f1*f2 )
    S21 = S12
    S22 = inf_integrate( 4*pi*(r**2)*f2**2 )

    E0 = solve_quadratic(
            H11=H11, H22=H22, H12=H12, H21=H21,
            S11=S11, S22=S22, S12=S12, S21=S21,
        )
    print(E0)

    C = -( H11 - E0*S11 )/( H12 - E0*S12 )
    psi0 = (f1 + C*f2).simplify()

在我的机器上，这个收敛速度大约几秒钟就可以收敛到-1/2。

谢谢！你的回复非常及时，现在它运行得非常好。你知道我怎么才能让它在E0值上显示更多的数字吗？另外，积分应该是从0到oo，但是使用这些边界会导致一些错误。将上限设置为40大部分相当于“oo”，将下限设置为10e-5大部分相当于（0）。现在这样做很好，但我感觉当我尝试提高精度时，它会引起太多的误差。 - Shrodinger149
你可以始终使用 decimal 来提高精度以显示更多数字。https://docs.python.org/3.7/library/decimal.html 这可能会稍微减慢一些速度，但可能不会严重影响效率。 - EchoCloud
@Shrodinger149：对于数值估算无穷积分，肯定有更适合的方法，特别是自适应方法，因为您可以自由地评估函数。如果您可以提供一些问题的具体细节（函数衰减速度如何；它们是否振荡等），这可能是一个[scicomp.se]的问题。但是，自己深入研究这个主题也是值得的。 - Wrzlprmft
@Shrodinger149：关于精度，您真的确定需要超过双精度浮点精度吗？无论如何，目前限制您精度的是积分的准确性，因此提高数字精度没有任何意义。 - Wrzlprmft
2

Wrzlprmft的最初回答很好。我已经做了一些整理，并用SymPy的积分替换了笨重的集成函数。这在我的原始代码上没有起作用，但在Wrzlprmft的更正/添加后完美地运行。程序有点慢（仍然比我的原始代码快数个数量级），但不再有限制精度的错误。下面是最终的代码：

from symengine import *
from sympy import *
import sympy

r, E1 = symbols('r, E1')
H11, H22, H12, H21 = symbols("H11, H22, H12, H21")
S11, S22, S12, S21 = symbols("S11, S22, S12, S21")
low = 0
high = oo
n = 100000

quadratic_expression = (H11-E1*S11)*(H22-E1*S22)-(H12-E1*S12)*(H21-E1*S21)
general_solution = sympify(sympy.solve(quadratic_expression, E1)[0])


def solve_quadratic(**kwargs):
    return general_solution.subs(kwargs)


def H(fun):
    return -fun.diff(r, 2)/2 - fun.diff(r)/r - fun/r


psi0 = exp(-3*r/2)
I1 = N(integrate(4*pi*(r**2)*psi0*H(psi0), (r, low, high)))
I2 = N(integrate(4*pi*(r**2)*psi0**2, (r, low, high)))
E0 = I1/I2
print(E0)

for x in range(100):
    f1 = psi0
    f2 = r * (H(psi0)-E0*psi0)
    Hf1 = H(f1)
    Hf2 = H(f2)

    H11 = integrate(4*pi*(r**2)*f1*Hf1, (r, low, high))
    H12 = integrate(4*pi*(r**2)*f1*Hf2, (r, low, high))
    H21 = integrate(4*pi*(r**2)*f2*Hf1, (r, low, high))
    H22 = integrate(4*pi*(r**2)*f2*Hf2, (r, low, high))

    S11 = integrate(4*pi*(r**2)*f1**2, (r, low, high))
    S12 = integrate(4*pi*(r**2)*f1*f2, (r, low, high))
    S21 = S12
    S22 = integrate(4*pi*(r**2)*f2**2, (r, low, high))

    E0 = solve_quadratic(
            H11=H11, H22=H22, H12=H12, H21=H21,
            S11=S11, S22=S22, S12=S12, S21=S21,
        )
    print(E0)

    C = -(H11 - E0*S11)/(H12 - E0*S12)
    psi0 = (f1 + C*f2).simplify()
你可以通过将第一行更改为 from symengine import *; import sympy; integrate = lambda *args: sympy.N(sympy.integrate(*args)) 来加快速度。这解决了以下问题：1）显式携带像 1+sqrt(pi) 这样的项浪费了大量时间，这可能对精度没有太大贡献，而且会累积。2）from sympy import * 完全覆盖了 from symengine import *，使你失去了 SymEngine 的速度优势。 - Wrzlprmft
已接受并点赞（我的点赞并不会改变什么，因为我声望不够高）。非常感谢您的所有帮助。对于“integrate = lambda *args: sympy.N(sympy.integrate(*args))”这部分代码，它应该放在哪里？还有，args是什么？ - Shrodinger149
“where is that supposed to go?” - 你可以直接放在导入语句后面。“Also, what is args?” - 它是一个包含所有参数的元组。当你编写一个简单的包装器并希望将所有位置参数原封不动地传递给另一个函数时，名称 args 是标准约定。 - Wrzlprmft
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接