numpy 在SymPy中使用带lambdify的公共子表达式消除的最佳实践

llycmphe  于 2023-04-30  发布在  其他
关注(0)|答案(3)|浏览(144)

我目前正在尝试使用SymPy来生成和数值计算函数及其梯度。为了简单起见,我将使用以下函数作为示例(请记住,真实的函数要长得多):

import sympy as sp
def g(x):
    return sp.cos(x) + sp.cos(x)**2 + sp.cos(x)**3

数值计算这个函数及其导数很容易:

import numpy as np
g_expr = sp.lambdify(x,g(x),modules='numpy')
dg_expr = sp.lambdify(x,sp.diff(g(x)),modules='numpy')

print g_expr(np.linspace(0,1,50))
print dg_expr(np.linspace(0,1,50))

然而,对于我的真实的函数,lambdify很慢,无论是在生成数值函数还是在计算方面。由于我的函数中的许多元素都是相似的,我想在lambdify中使用公共子表达式消除(cse)来加速这个过程。我知道SymPy有一个内置函数来执行cse,

>>> print sp.cse(g(x))
([(x0, cos(x))], [x0**3 + x0**2 + x0])

但是不知道要使用什么语法来在我的lambdify函数中使用这个结果(我仍然希望使用x作为输入参数):

>>> g_expr_fast = sp.lambdify(x,sp.cse(g(x)),modules='numpy')
>>> print g_expr_fast(np.linspace(0,1,50))
Traceback (most recent call last):
  File "test3.py", line 34, in <module>
    print g_expr1(nx1)
  File "<string>", line 1, in <lambda>
NameError: global name 'x0' is not defined

任何关于如何在lambdify中使用cse的帮助都将不胜感激。或者,如果有更好的方法来加速我的梯度计算,我也会很感激听到这些。
如果它是相关的,我使用Python 2。7.3和SymPy 0。7.6.

u1ehiz5o

u1ehiz5o1#

计算速度可以提高:

  • 通过提取公因子并且不重复计算,可以获得一点(在我的示例中是36倍的改进)。
  • 更大的加速(高达100 x-1000 x)可以通过例如使用numba创建扩展模块来获得。

前言

我假设这是“计算函数在渐近一次和使用它在不同的项目以后多次”类型的情况。因此,有一些手动复制粘贴和创建文件包括在内。* 然而 * 它可以自动创建新文件的功能,也编译步骤,但我离开了这一点。
我也遇到了类似的问题,我对不同的方法做了一些基准测试。我使用的函数很长(len(str(expr)) = 45857),cse(expr)将其分解为72个子表达式。在这里复制粘贴太长了,但这里有一些步骤,可以使使用sympy创建的函数的速度提高100 - 1000倍。

基准测试

A)求单浮点数

对每个参数使用 * 一个浮点值 * 来计算函数的时间。使用timeit myfunc(*params)

  • 基线:使用modules="numpy"进行lambda定义:277µs
  • (1)将str(expr)复制粘贴到函数定义:275µs(无差异)
  • (2)cse后表达式的复制粘贴:8.2 µs(提高33倍)
  • (3)cse(optimizations="basic")后表达式的复制粘贴:7.6µs(提高36倍)
  • (4)使用numba将代码编译为func_numba_f():0.25µs(提高1090倍)
  • (5)使用sympy autowrap:0.47 µs(提高589倍)
    B)评估np.数组1000浮点数
  • (1)将str(expr)复制粘贴到函数定义:15100 µs|15.1µs/值
  • (2)cse后表达式的复制粘贴:493微秒|每个值0.49µs(31倍改进)
  • (3)cse(optimizations="basic")后表达式的复制粘贴:413微秒|每个值0.41µs(37倍改进)
  • (4)使用numba将代码编译为func_numba_arr():114µs|0.11µs/值(132倍改善)
  • (5)用NP自动 Package 。矢量化:480µs|0.48µs/值(31倍改善)

(1)str(expr)的复制粘贴

  • 将表达式字符串复制粘贴到新函数中,并返回值.
  • 将函数保存在另一个文件中。

(2)cse后表达式的复制粘贴

  • 想法:通过识别公共部分使代码更短。
  • 首先,复制粘贴公共部分:
repl, redu = cse(K)
for variable, expr in repl:
    print(f"{variable} = {expr}")
  • 然后,复制粘贴返回值:print(redu[0]) .
  • 创建另一个文件,并粘贴到函数定义中

(3)cse(optimizations="basic")后表达式的复制粘贴

  • 与(2)相同,但带有optimizations="basic"
  • 这会创建比(2)稍短的代码

(4)使用numba编译代码

  • 使用numba.pycc.CC来编译代码。因此,创建一个复制粘贴函数,如(3)所示
  • 然后,使用代码创建src_mymodule.py
from numba.pycc import CC

cc = CC("my_numba_module")

@cc.export("func_numba_f", "f8(f8, f8, f8, f8, f8)")
@cc.export("func_numba_arr", "f8[:](f8[:],f8[:],f8[:],f8[:],f8[:])")
def myfunc(x1, x2, x3, x4, x5):
    # your function definition here
    return value

if __name__ == "__main__":
    cc.compile()
  • 函数func_numba_f()中有 * 5个 * 浮点值输入变量和一个浮点值输出变量。f8表示浮点数。
  • func_numba_arr()是处理np的版本。dtype="float64"dtype="float32"的数组,具体取决于编译时使用的内容。
  • 然后通过运行python src_mymodule.py编译一次代码。这将创建my_numba_module.cp38-win_amd64.pyd或类似。它只能与文件名中的 * 相同的python版本和位数 * 一起使用。
  • 然后,在另一个python文件中,导入函数并使用它们,如:
from my_numba_module import func_numba_f, func_numba_arr

out = func_numba_f(4,3,2,1,100)

# or:
args = [np.array([x]*N, dtype='float64') for x in (4,3,2,1,100)]
out_arr = func_numba_arr(*args)

(5)使用sympy autowrap

  • 这很简单。在installing和配置Cython之后,我需要两行代码来使用autowrap创建一个函数。
from sympy.utilities.autowrap import autowrap
func = autowrap(expr, backend='cython')
  • 通过指定temp_dir参数,它将保存所有源文件(。c,.h,.pyx)和a .pyd(win)/ .so(unix)文件,可用于稍后导入函数(假设temp_dirsys.path中):
from wrapper_module_1 import autofunc_c
  • 这是迄今为止最简单的方法,尽管生成的C代码不是高度优化的。如果需要,可以在其他步骤中重命名自动捕获的输出。
  • 该函数仅接受标量,但可以使用np.vectorize进行矢量化
func = np.vectorize(func)
gab6jxml

gab6jxml2#

所以这可能不是最好的方法,但对于我的小例子来说,它是有效的。
下面代码的思想是对每个公共子表达式进行lambdifying,并生成一个可能包含所有参数的新函数。我添加了一些额外的sin和cos项,以添加来自先前子表达式的可能依赖项。

import sympy as sp
import sympy.abc
import numpy as np
import matplotlib.pyplot as pl

def g(x):
    return sp.cos(x) + sp.cos(x)**2 + sp.cos(x)**3 + sp.sin(sp.cos(x)+sp.sin(x))**4 + sp.sin(x) - sp.cos(3*x) + sp.sin(x)**2

repl, redu=sp.cse(g(sp.abc.x))

funs = []
syms = [sp.abc.x]
for i, v in enumerate(repl):
    funs.append(sp.lambdify(syms,v[1],modules='numpy'))
    syms.append(v[0])

glam = sp.lambdify(syms,redu[0],modules='numpy')

x = np.linspace(-1,5,10)
xs=[x]

for f in funs:
    xs.append(f(*xs))

print glam(*xs)
glamlam = sp.lambdify(sp.abc.x,g(sp.abc.x),modules='numpy')
print glamlam(x)
print np.allclose(glamlam(x),glam(*xs))

repl包含:

[(x0, cos(x)), (x1, sin(x)), (x2, x0 + x1)]

和redu包含

[x0**3 + x0**2 + x1**2 + x2 + sin(x2)**4 - cos(3*x)]

所以funs包含所有子表达式lambdified,列表xs包含每个子表达式,这样最终可以正确地馈送glamxs随着每个子表达式的增长而增长,最终可能会成为瓶颈。
您可以对sp.cse(sp.diff(g(sp.abc.x)))的表达式执行相同的方法。

2ul0zpep

2ul0zpep3#

从SymPy 1开始9,lambdify可以使用kwarg cse=True应用公共子表达式消除。
参见: www.example.com

相关问题