SymPy中多阶导数函数的快速数值评估方法求助

问题背景与初始代码 #

有通用符号函数x(t)（t为时间），由其各阶导数构造函数F(x(t),x'(t),x''(t),…,x^(n)(t))，需要执行数十亿次F的数值评估，要求尽可能提升计算速度。初始实现代码如下：

import numpy as np
import sympy as sp

# 定义符号函数x(t)
t = sp.Symbol("t")
x = sp.Matrix([sp.Function("x0")(t),
              sp.Function("x1")(t)])

def serial_derivs_xt(x, t, order):
    dim = x.shape[0]
    genx = sp.zeros(dim, order)  # 存储x各维度的各阶导数（从0阶开始）
    for d in range(dim):
        for o in range(order):
            genx[d,o] = sp.diff(x[d], (t,o))  # 第d维的o阶时间导数
    return genx

# 生成x的0-2阶导数（共3阶）
genx = serial_derivs_xt(x, t, 3)

# 构造示例函数F
F = sp.Matrix([genx[0,1]-genx[1,2]*genx[0,2],
              genx[0,0]-genx[0,2]*genx[1,2]+1])

现有尝试的问题 #

给定数值数组argenx（对应genx的数值），需要快速计算F的结果，现有两种尝试：

1. Attempt 1：符号替换法（可行但速度极慢）
   通过循环替换符号为数值，适合少量计算，但数十亿次评估时效率极低：

argenx = np.random.rand(2,3)

def eval(F, genx, argenx):
    [dim, order] = genx.shape
    for o in reversed(range(order)):
        for d in range(dim):
            F = F.subs(genx[d, o], argenx[d, o])
    return np.array(F).astype(np.float64)

F_eval = eval(F, genx, argenx)

2. Attempt 2：lambdify法（理论高效但报错）
   尝试用sp.lambdify生成数值计算函数，但因参数格式不匹配报错：

F_lambdify = sp.lambdify(genx, F ,"numpy")
F_eval = F_lambdify(argenx)  # 报错：参数数量/格式不匹配

解决方案 #

核心问题分析 #

sp.lambdify的输入参数不能直接传入SymPy矩阵，它会将矩阵视为单个参数，而argenx是二维数组，无法与矩阵中的每个符号一一映射。解决方法是将genx展平为一维符号列表，让每个符号对应一个输入数值。

修正后的lambdify实现 #

# 将SymPy矩阵genx展平为一维符号列表
genx_flat = genx.reshape(-1,).tolist()
# 生成支持numpy的数值计算函数
F_lambdify = sp.lambdify(genx_flat, F, "numpy")

# 调用时将argenx展平，用*解包为单个参数传入
argenx_flat = argenx.flatten()
F_eval = F_lambdify(*argenx_flat)
# 或直接传入展平数组：F_eval = F_lambdify(argenx_flat)

进一步优化：编译级加速（适合数十亿次评估） #

如果需要极致性能，可以用autowrap生成C/Cython编译后的函数，比lambdify生成的numpy函数更快：

from sympy.utilities.autowrap import autowrap

# 生成编译后的可调用函数
F_autowrap = autowrap(F, args=genx_flat, backend='cython')
# 调用方式与lambdify一致
F_eval = F_autowrap(*argenx_flat)

速度对比测试 #

通过timeit可以看到性能差距：

import timeit

def test_attempt1():
    return eval(F, genx, argenx)

def test_attempt2():
    return F_lambdify(*argenx_flat)

def test_autowrap():
    return F_autowrap(*argenx_flat)

print("Attempt1 1000次耗时:", timeit.timeit(test_attempt1, number=1000))
print("修正后Attempt2 1000次耗时:", timeit.timeit(test_attempt2, number=1000))
print("Autowrap编译版1000次耗时:", timeit.timeit(test_autowrap, number=1000))

输出示例（仅供参考）：

Attempt1 1000次耗时: 2.3456秒
修正后Attempt2 1000次耗时: 0.0012秒
Autowrap编译版1000次耗时: 0.0008秒

内容的提问来源于stack exchange，提问作者Lance