用 sympy 缓慢替换符号矩阵答案

【问题标题】：Slow substitution of symbolic matrix with sympy用 sympy 缓慢替换符号矩阵
【发布时间】：2017-08-18 11:41:00
【问题描述】：

我正在使用 sympy 处理大小为 QxQ 的符号雅可比矩阵 J。该矩阵的每个系数包含Q 符号，从f[0] 到f[Q-1]。我想做的是将J 的每个系数中的每个符号替换为已知值g[0] 到g[Q-1]（不再是符号）。我发现最快的方法如下：

for k in range(Q):
    J = J.subs(f[k], g[k])

但是，我发现这个“基本”操作很长！例如，使用这个 MCVE：

import sympy
import numpy as np
import time

Q = 17
f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, f16 = \
    sympy.symbols("f0 f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12 f13 f14 f15 f16")
f = [f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, f16]
e = np.array([0., 1., 0., -1., 0., 1., -1., -1., 1.,
              2., -2., -2., 2., 3., 0., -3., 0.])
u = np.sum(f * e) / np.sum(f)
function = e * np.sum(f) * (1. + u * e + (u * e)**2 - u * u)
F = sympy.Matrix(function)
g = e * (1. + 0.2 * e + (0.2 * e)**2)

start_time = time.time()
J = F.jacobian(f)
print("--- %s seconds ---" % (time.time() - start_time))

start_time = time.time()
for k in range(Q):
    J = J.subs(f[k], g[k])
print("--- %s seconds ---" % (time.time() - start_time))

在我的计算机上替换大约需要 5 秒，而雅可比矩阵的计算只需要 0.6 秒。在另一个（更长的）代码中，Q=37 的替换需要 360 秒（而雅可比计算需要 20 秒）！

此外，当我查看正在运行的进程时，我可以看到 Python 进程有时会在矩阵替换期间停止工作。

有谁知道这可能来自哪里？
有没有办法让这个操作更快？

【问题讨论】：

欢迎使用 SO，不要忘记使用tour :) 关于这个问题，我怀疑雅可比可能不会被简化，导致大量重复的变量/计算。

标签： python numpy sympy substitution differential-equations

【解决方案1】：

您可能想为此尝试Theano。它实现了一个jacobian 函数，它比sympy 并行且更快。

以下版本实现了3.88的加速！现在替换时间不如秒。

import numpy as np
import sympy as sp
import theano as th
import time


def main_sympy():
    start_time = time.time()

    Q = 17
    f = sp.symbols(('f{} ' * Q).format(*range(Q)))

    e = np.array([0., 1., 0., -1., 0., 1., -1., -1., 1.,
                  2., -2., -2., 2., 3., 0., -3., 0.])
    u = np.sum(f * e) / np.sum(f)
    ue = u * e
    phi = e * np.sum(f) * (1. + ue + ue*ue - u*u)
    F = sp.Matrix(phi)
    J = F.jacobian(f)

    g = e * (1. + 0.2*e + (0.2*e)**2)

    for k in range(Q):
        J = J.subs(f[k], g[k])

    print("--- %s seconds ---" % (time.time() - start_time))
    return J


def main_theano():
    start_time = time.time()

    Q = 17
    f = th.tensor.dvector('f')

    e = np.array([0., 1., 0., -1., 0., 1., -1., -1., 1., 2.,
                  -2., -2., 2., 3., 0., -3., 0.])
    u = (f * e).sum() / f.sum()
    ue = u * e
    phi = e * f.sum() * (1. + ue + ue*ue - u*u)
    jacobi = th.gradient.jacobian(phi, f)
    J = th.function([f], jacobi)

    g = e * (1. + 0.2*e + (0.2*e)**2)
    Jg = J(g)  # evaluate expression

    print("--- %s seconds ---" % (time.time() - start_time))
    return Jg


J0 = np.array(main_sympy().tolist(), dtype='float64')
J1 = main_theano()

print(np.allclose(J0, J1))  # compare results

【讨论】：