我在为sympy.MatrixSymbol的元素似乎无法与sympy的分化例程很好地交互而苦苦挣扎。

我尝试使用sympy.MatrixSymbol元素而不是“普通” sympy符号的事实是因为我想自动包装一个大函数,这似乎是克服参数限制并启用a输入的唯一方法单阵列。

为了让读者了解对可能解决方案的限制,我将首先概述我的意图。但是,仓促的读者也可以跳到下面的代码块,这些代码块说明了我的问题。


声明某种向量或变量数组。
从前者的元素中构建一些表达;这些表达式将构成所述向量的向量值函数的组成部分。除了此功能外,我还想获得Jacobian w.r.t.向量。
使用自动换行(与cython后端一起)来获得矢量函数及其Jacobian的数值实现。这对前面的步骤施加了一些限制:(a)希望函数的输入作为向量而不是符号列表给出。 (这两者都是因为自动封装函数的输入数量似乎受到限制,并且以后为了减轻与scipy的交互,即避免经常将numpy矢量解包到列表中)。


在旅途中,我遇到了两个问题:


Cython似乎不喜欢某些sympy函数,其中sympy.Max是我严重依赖的函数。自动换行的“帮助程序”组合似乎无法一次处理多个帮助程序。
就我自己使用cython容易理解的abs()或sign()来规避它本身而言,这并不是什么大不了的事情。


(另请参见上面的this question


如前所述,autowrap / cython不能接受超过509个符号形式的参数,至少在我的编译器设置中不这样。 (另请参见here
因为无论如何我都希望给向量而不是列表作为函数的输入,所以我寻找一种方法来获取包装函数以numpy数组作为输入(与DeferredVector + lambdify相比)。看起来很自然的方法是sympy.MatrixSymbol。 (请参阅上面的链接。我不确定是否有其他选择,如果可以的话,欢迎提出建议。)
然后,我的最新问题就从这里开始:我意识到sympy.MatrixSymbol的元素在许多方面都不像“其他” sympy符号那样工作。必须分别为属性分配实和可交换的属性,尽管如此似乎工作良好。但是,当我尝试获得Jacobian时,我的真正麻烦就开始了。 sympy似乎没有开箱即用地衍生元素:




import sympy

X= sympy.MatrixSymbol("X",10,1)

for element in X:
    element._assumptions.update({"real":True, "commutative":True})

X[0].diff(X[0])
Out[2]: Derivative(X[0, 0], X[0, 0])


X[1].diff(X[0])
Out[15]: Derivative(X[1, 0], X[0, 0])


以下块是我想做的一个最小示例,但是这里使用普通符号:
(我认为它捕获了我需要的所有内容,如果我忘记了一些内容,我会在以后添加。)

import sympy
from sympy.utilities.autowrap import autowrap

X = sympy.symbols("X:2", real = True)
expr0 = X[1]*( (X[0] - abs(X[0]) ) /2)**2
expr1 = X[0]*( (X[1] - abs(X[1]) ) /2)**2
F = sympy.Matrix([expr0, expr1])
J = F.jacobian([X[0],X[1]])
J_num = autowrap(J, args = [X[0],X[1]], backend="cython")


这是我(当前)使用sympy.MatrixSymbol的最佳猜测,它当然会失败,因为Derivative中的J-表达式:

X= sympy.MatrixSymbol("X",2,1)
for element in X:
    element._assumptions.update({"real":True, "commutative":True, "complex":False})
expr0 = X[1]*( (X[0] - abs(X[0]) ) /2)**2
expr1 = X[0]*( (X[1] - abs(X[1]) ) /2)**2
F = sympy.Matrix([expr0, expr1])
J = F.jacobian([X[0],X[1]])
J_num = autowrap(J, args = [X], backend="cython")


运行上述命令后,J看起来像这样:

J
Out[50]:
Matrix([
[(1 - Derivative(X[0, 0], X[0, 0])*X[0, 0]/Abs(X[0, 0]))*(-Abs(X[0, 0])/2 + X[0, 0]/2)*X[1, 0], (-Abs(X[0, 0])/2 + X[0, 0]/2)**2],
[(-Abs(X[1, 0])/2 + X[1, 0]/2)**2, (1 - Derivative(X[1, 0], X[1, 0])*X[1, 0]/Abs(X[1, 0]))*(-Abs(X[1, 0])/2 + X[1, 0]/2)*X[0, 0]]])


毫无疑问,自动换行不喜欢:

[...]
wrapped_code_2.c(4): warning C4013: 'Derivative' undefined; assuming extern returning int
[...]
wrapped_code_2.obj : error LNK2001: unresolved external symbol Derivative


如何告诉sympy X[0].diff(X[0])=1X[0].diff(X[1])=0?甚至是abs(X[0]).diff(X[0]) = sign(X[0])

还是可以使用sympy.MatrixSymbol和still get a cythonized function(输入是单个向量而不是符号列表)来解决问题?

对于任何输入都将非常有用,很可能是上述过程的任何步骤的解决方法。谢谢阅读!



编辑:
简短的一句话:我可以提出的一个解决方案是:
使用普通符号构造FJ;然后用一些sympy.MatrixSymbol的元素替换两个表达式中的符号。这似乎可以完成工作,但是替换需要花费大量时间,因为J可以达到〜1000x1000及更高的尺寸。因此,我宁愿避免这种方法。

最佳答案

经过更广泛的研究后,似乎我上面描述的问题已在development / github版本中解决。相应地更新后,所有涉及DerivativeMatrixElement术语都可以正确解析!

请参见here以供参考。

关于python - sympy.MatrixSymbol:区分问题,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/44377903/

10-13 09:14