主页 / coding / 问题 / 79591252
Accepted
Seneral
Asked: 2025-04-25 02:20:57 +0800 CST

除非需要，否则不要扩展子表达式

  • 772

我想求一个复杂函数的导数。不幸的是，当我构造要求导的表达式时，它会将所有子表达式扩展至紧邻的符号。因此，一旦我对最终表达式进行微分，它就会膨胀到难以控制的大小。我曾尝试手动将一些子表达式替换回最终结果中之前的形式，但由于表达式构造过程中的重构，错失了很多简化的机会。

这只是我想要推导的第一阶段，即使尝试替换一些用于构建它的表达式，它已经变得非常臃肿。

from sympy import *

aX, aY, aZ = symbols('aX aY aZ')
rotInc = Matrix(3,1,[aX, aY, aZ])
theta = sqrt((rotInc.T @ rotInc)[0,0])
incQuat = Quaternion.from_axis_angle(rotInc/theta, theta*2)

qX, qY, qZ, qW = symbols('qX qY qZ qW')
baseQuat = Quaternion(qW, qX, qY, qZ)
poseQuat = incQuat * baseQuat

d4 = diff(poseQuat, aX)
d4s = d4.subs({
    incQuat.a: symbols('iW'),
    incQuat.b: symbols('iX'),
    incQuat.c: symbols('iY'),
    incQuat.d: symbols('iZ'),
    theta: symbols('theta')
})

我知道cse（公共子表达式消除），这表明存在某种系统可以保留命名的子表达式。我希望 sympycse在我构造表达式时构建一个类似 returns 的结构，并且只在必要时（例如在微分过程中）用子表达式的组件替换它。这样可以在最终的微分过程中保留大部分这些符号，从而得到更清晰、可立即使用的输出。

是否有这样一种在 sympy 中构造表达式的模式/方法，或者其他可以帮助我保持最终表达式简单的方法？谢谢！

sympy

3 个回答

  1. 这不是一个答案，因为我不知道你想得到什么。但是，我理解你的问题。你知道theta是 的函数aX, aY, aZ， 的分量也是 的函数incQuat。如果我们像这样处理未定义的应用函数会怎么样？

    
from sympy import *

aX, aY, aZ = symbols('a_X a_Y a_Z')
rotInc = Matrix(3,1,[aX, aY, aZ])
theta = Function("theta")(aX, aY, aZ)
incQuat = Quaternion.from_axis_angle(rotInc/theta, theta*2)

# substitution dictionary, will be used later
subs_d1 = {
    theta: sqrt((rotInc.T @ rotInc)[0,0])
}
symb_components = [Function(f"f_{i}")(aX, aY, aZ) for i in range(4)]
# substitution dictionary, will be used later
subs_d2 = {f: c for f, c in zip(symb_components, incQuat.to_Matrix())}
incQuat_symbolic = Quaternion(*symb_components)

qX, qY, qZ, qW = symbols('q_X q_Y q_Z q_W')
baseQuat = Quaternion(qW, qX, qY, qZ)
poseQuat = incQuat_symbolic * baseQuat

d4 = diff(poseQuat, aX)
d4

    所得方程

    这是第一步的表达式。然后，当你准备好进行实际计算时，请代入字典：

    d4.subs(subs_d2).doit().subs(subs_d1).doit()

    最终你会得到非常长的表达式......

    • 1
  2. Best Answer

    我以为这在 SymPy 的某个地方，但也许在准备编译函数（或类似的东西）时用到了它。不过，按照你的想法（在这个问题上可能有点多余）：

    from sympy import symbols, Function, Derivative, cse

def dep(r, x):
    reps = {}
    for v, e in r:
        if e.xreplace(reps).has(x):
            reps[v] = Function(str(v))(x)
    return reps

def differentiate_with_cse(expr, x, backsub=False):
    r, e = cse(expr)
    reps = dep(r, x)
    dr = []
    for v, s in r:
        if v in reps:
            f = reps[v]
            df = Derivative(f, x).doit()
            ds = s.xreplace(reps).diff(x)
            dr.append((df, ds))
    dexpr = e[0].subs(reps).diff(x).subs(dr).expand()
    for k, v in reversed(reps.items()):
        dexpr = dexpr.subs(v, k)
    return dexpr.subs(list(reversed(r))) if backsub else dexpr,dict(r)

...your code up to derivative

>>> deriv, pats = differentiate_with_cse(poseQuat, aX)
>>> deriv.simplify()
(-qX*x7) + qW*x7*i + (-qZ*x7)*j + qY*x7*k
>>> pats
{x0: aX**2, x1: aY**2, x10: aZ*x7, x2: aZ**2, x3: x0 + x1 + x2, x4: sqrt(x3), x5: cos(x4), x6: 1/x3, x7: sin(x4)/(x4*sqrt(x0*x6 + x1*x6 + x2*x6)), x8: aX*x7, x9: aY*x7})
    • 0

  3. 因此，借助@Davide_sd 提供的函数提示，我创建了一个通用方法，可以让我非常轻松地控制子步骤的拆分方式。基本上，我手动推导拆分出的函数，但与 cse 非常相似，将结果保存在字典中以在所有实例之间共享。
    组成计算的基本表达式是输入，永远不会被修改，推导列表以您想要推导的内容为种子（多个表达式也可以），并且它将根据需要使用表达式列表递归地推导它们。
    最后，我仍然可以使用 cse 来 a) 如果需要，将其转换为该格式，以及 b) 分解出更常见的实例。
    它与我的小示例配合得很好，可能会随着我为需要推导的函数添加更多复杂性而更新它。

    from sympy import *

def find_derivatives(expression):
    derivatives = []
    if isinstance(expression, Derivative):
        #print(expression)
        derivatives.append(expression)
    elif isinstance(expression, Basic):
        for a in expression.args:
            derivatives += find_derivatives(a)
    elif isinstance(expression, MatrixBase):
        for i in range(rows):
            for j in range(cols):
                derivatives += find_derivatives(self[i, j])
    return derivatives

def derive_recursively(expression_list, derive_done, derive_todo):
    newly_derived = {}
    for s, e in derive_todo.items():
        print("Handling derivatives in " + str(e))
        derivatives = find_derivatives(e)
        for d in derivatives:
            if d in newly_derived:
                #print("Found derivative " + str(d) + " in done list, already handled!")
                continue
            if d in derive_todo:
                #print("Found derivative " + str(d) + " in todo list, already handling!")
                continue
            if d in expression_list:
                #print("Found derivative " + str(d) + " in past list, already handled!")
                continue
            if d.expr in expression_list:
                expression = expression_list[d.expr]
                print("  Deriving " + str(d.expr) + " w.r.t. " + str(d.variables))
                print("    Expression: " + str(expression))
                derivative = Derivative(expression, *d.variable_count).doit().simplify()
                print("    Derivative: " + str(derivative))
                if derivative == 0:
                    e = e.subs(d, 0)
                    derive_todo[s] = e
                    print("        Replacing main expression with: " + str(e))
                    continue
                newly_derived[d] = derivative
                continue
            print("Did NOT find base expression " + str(d.expr) + " in provided expression list!")
    derive_done |= derive_todo
    if len(newly_derived) == 0:
        return derive_done
    return derive_recursively(expression_list, derive_done, newly_derived)

incRot_c = symbols('aX aY aZ')
incRot_s = Matrix(3,1,incRot_c)
theta_s = Function("theta")(*incRot_c)
theta_e = sqrt((incRot_s.T @ incRot_s)[0,0])
incQuat_c = [ Function(f"i{i}")(*incRot_c) for i in "WXYZ" ]
incQuat_s = Quaternion(*incQuat_c)
incQuat_e = Quaternion.from_axis_angle(incRot_s/theta_s, theta_s*2)
baseQuat_c = symbols('qX qY qZ qW')
baseQuat_s = Quaternion(*baseQuat_c)
poseQuat_c = [ Function(f"p{i}")(*incRot_c, *baseQuat_c) for i in "WXYZ" ]
poseQuat_s = Quaternion(*poseQuat_c)
# Could also do it like this and in expressions just refer poseQuat_s to poseQuat_e, but output is less readable
#poseQuat_s = Function(f"pq")(*incRot_c, *baseQuat_c)
poseQuat_e = incQuat_s * baseQuat_s

expressions = { theta_s: theta_e } | \
    { incQuat_c[i]: incQuat_e.to_Matrix()[i] for i in range(4) } | \
    { poseQuat_c[i]: poseQuat_e.to_Matrix()[i] for i in range(4) }

derivatives = derive_recursively(expressions, {}, { symbols('res'): diff(poseQuat_s, incRot_c[0]) })
print(derivatives)

elements = cse(list(expressions.values()) + list(derivatives.values()))
pprint(elements)
    • 0

相关问题