递归符号计算 - 提高性能
Recursive symbolic calculations - improve the performance
在我的研究中，我试图解决 Kolmogorov 逆向方程，即对
$$Af = b(x)f'(x)+\sigma(x)f''(x)$$
对于特定的 b(x) 和 \sigma(x)，我试图了解在计算更高的 Af 幂时表达式的系数增长的速度。我很难从分析上得出这一点，因此试图从经验上看趋势。
首先，我用过sympy:
from sympy import *
import matplotlib.pyplot as plt
import re
import math
import numpy as np
import time
np.set_printoptions(suppress=True)
x = Symbol('x')
b = Function('b')(x)
g = Function('g')(x)
def new_coef(gamma, beta, coef_minus2, coef_minus1, coef):
return expand(simplify(gamma*coef_minus2 + beta*coef_minus1 + 2*gamma*coef_minus1.diff(x)\
+beta*coef.diff(x)+gamma*coef.diff(x,2)))
def new_coef_first(gamma, beta, coef):
return expand(simplify(beta*coef.diff(x)+gamma*coef.diff(x,2)))
def new_coef_second(gamma, beta, coef_minus1, coef):
return expand(simplify(beta*coef_minus1 + 2*gamma*coef_minus1.diff(x)\
+beta*coef.diff(x)+gamma*coef.diff(x,2)))
def new_coef_last(gamma, beta, coef_minus2):
return lambda x: gamma(x)*coef_minus2(x)
def new_coef_last(gamma, beta, coef_minus2):
return expand(simplify(gamma*coef_minus2 ))
def new_coef_second_to_last(gamma, beta, coef_minus2, coef_minus1):
return expand(simplify(gamma*coef_minus2 + beta*coef_minus1 + 2*gamma*coef_minus1.diff(x)))
def set_to_zero(expression):
expression = expression.subs(Derivative(b, x, x, x), 0)
expression = expression.subs(Derivative(b, x, x), 0)
expression = expression.subs(Derivative(g, x, x, x, x), 0)
expression = expression.subs(Derivative(g, x, x, x), 0)
return expression
def sum_of_coef(expression):
sum_of_coef = 0
for i in str(expression).split(' + '):
if i[0:1] == '(':
i = i[1:]
integers = re.findall(r'\b\d+\b', i)
if len(integers) > 0:
length_int = len(integers[0])
if i[0:length_int] == integers[0]:
sum_of_coef += int(integers[0])
else:
sum_of_coef += 1
else:
sum_of_coef += 1
return sum_of_coef
power = 6
charar = np.zeros((power, power*2), dtype=Symbol)
coef_sum_array = np.zeros((power, power*2))
charar[0,0] = b
charar[0,1] = g
coef_sum_array[0,0] = 1
coef_sum_array[0,1] = 1
for i in range(1, power):
#print(i)
for j in range(0, (i+1)*2):
#print(j, ':')
#start_time = time.time()
if j == 0:
charar[i,j] = set_to_zero(new_coef_first(g, b, charar[i-1, j]))
elif j == 1:
charar[i,j] = set_to_zero(new_coef_second(g, b, charar[i-1, j-1], charar[i-1, j]))
elif j == (i+1)*2-2:
charar[i,j] = set_to_zero(new_coef_second_to_last(g, b, charar[i-1, j-2], charar[i-1, j-1]))
elif j == (i+1)*2-1:
charar[i,j] = set_to_zero(new_coef_last(g, b, charar[i-1, j-2]))
else:
charar[i,j] = set_to_zero(new_coef(g, b, charar[i-1, j-2], charar[i-1, j-1], charar[i-1, j]))
#print("--- %s seconds for expression---" % (time.time() - start_time))
#start_time = time.time()
coef_sum_array[i,j] = sum_of_coef(charar[i,j])
#print("--- %s seconds for coeffiecients---" % (time.time() - start_time))
coef_sum_array
然后，研究自动微分并使用 autograd：
import autograd.numpy as np
from autograd import grad
import time
np.set_printoptions(suppress=True)
b = lambda x: 1 + x
g = lambda x: 1 + x + x**2
def new_coef(gamma, beta, coef_minus2, coef_minus1, coef):
return lambda x: gamma(x)*coef_minus2(x) + beta(x)*coef_minus1(x) + 2*gamma(x)*grad(coef_minus1)(x)\
+beta(x)*grad(coef)(x)+gamma(x)*grad(grad(coef))(x)
def new_coef_first(gamma, beta, coef):
return lambda x: beta(x)*grad(coef)(x)+gamma(x)*grad(grad(coef))(x)
def new_coef_second(gamma, beta, coef_minus1, coef):
return lambda x: beta(x)*coef_minus1(x) + 2*gamma(x)*grad(coef_minus1)(x)\
+beta(x)*grad(coef)(x)+gamma(x)*grad(grad(coef))(x)
def new_coef_last(gamma, beta, coef_minus2):
return lambda x: gamma(x)*coef_minus2(x)
def new_coef_second_to_last(gamma, beta, coef_minus2, coef_minus1):
return lambda x: gamma(x)*coef_minus2(x) + beta(x)*coef_minus1(x) + 2*gamma(x)*grad(coef_minus1)(x)
power = 6
coef_sum_array = np.zeros((power, power*2))
coef_sum_array[0,0] = b(1.0)
coef_sum_array[0,1] = g(1.0)
charar = [b, g]
for i in range(1, power):
print(i)
charar_new = []
for j in range(0, (i+1)*2):
if j == 0:
new_funct = new_coef_first(g, b, charar[j])
elif j == 1:
new_funct = new_coef_second(g, b, charar[j-1], charar[j])
elif j == (i+1)*2-2:
new_funct = new_coef_second_to_last(g, b, charar[j-2], charar[j-1])
elif j == (i+1)*2-1:
new_funct = new_coef_last(g, b, charar[j-2])
else:
new_funct = new_coef(g, b, charar[j-2], charar[j-1], charar[j])
coef_sum_array[i,j] = new_funct(1.0)
charar_new.append(new_funct)
charar = charar_new
coef_sum_array
但是，我对他们中任何一个的速度都不满意。我想进行至少一千次迭代，而在 运行 simpy 方法 3 天后，我得到了 30 :/
我希望可以优化第二种方法（数值）以避免每次都重新计算表达式。不幸的是，我自己看不到该解决方案。另外，我试过 Maple，但还是没有成功。
概览
所以，这里有两个关于导数的公式很有趣：
- Faa di Bruno's formula which is a way to quickly find the n-th derivative of
f(g(x)) , and looks a lot like the Multinomial theorem
- General Leibniz rule which is a way to quickly find the n-th derivative of
f(x)*g(x) and looks a lot like the Binomial theorem
这两个都在 pull request #13892 中讨论过 n 阶导数使用一般莱布尼兹规则加速。
I'm trying to see how fast the coefficients of the expression are growing
在您的代码中，计算 c[i][j] 的一般公式是这样的：
c[i][j] = g * c[i-1][j-2] + b * c[i-1][j-1] + 2 * g * c'[i-1][j-1] + g * c''[i-1][j]
（其中 c'[i][j] 和 c''[i][j] 是 c[i][j] 的一阶和二阶导数）
因此，根据上面提到的莱布尼茨法则，我直觉地认为，计算出的系数应该与Pascal's triangle有关（或者至少它们应该有某种组合关系）。
优化 #1
在原始代码中，函数sum_to_coef(f)将表达式f序列化为一个字符串，然后丢弃所有看起来不像数字的东西，然后对剩余数字求和。
我们可以通过遍历表达式树并收集我们需要的东西来避免这里的序列化
def sum_of_coef(f):
s = 0
if f.func == Add:
for sum_term in f.args:
res = sum_term if sum_term.is_Number else 1
if len(sum_term.args) == 0:
s += res
continue
first = sum_term.args[0]
if first.is_Number == True:
res = first
else:
res = 1
s += res
elif f.func == Mul:
first = f.args[0]
if first.is_Number == True:
s = first
else:
s = 1
elif f.func == Pow:
s = 1
return s
优化#2
在函数set_to_zero(expr)中，b的所有二阶和三阶导数，以及g的三阶和四阶导数都被零替换。
我们可以像这样将所有这些替换折叠成一个语句：
b3,b2=b.diff(x,3),b.diff(x,2)
g4,g3=g.diff(x,4),g.diff(x,3)
def set_to_zero(expression):
expression = expression.subs({b3:0,b2:0,g4:0,g3:0})
return expression
优化#3
在原始代码中，对于每个单元格 c[i][j] 我们调用 simplify。事实证明这对性能有很大影响，但实际上我们可以跳过这个调用，因为幸运的是我们的表达式只是导数或未知函数的乘积之和。
所以行
charar[i,j] = set_to_zero(expand(simplify(expr)))
变成
charar[i,j] = set_to_zero(expand(expr))
优化#4
也尝试了以下方法，但效果很小。
对于两个连续的 j 值，我们计算 c'[i-1][j-1] 两次。
j-1 c[i-1][j-3] c[i-1][j-2] c[i-1][j-1]
j c[i-1][j-2] c[i-1][j-1] c[i-1][j]
如果您查看 else 分支中的循环公式，您会发现 c'[i-1][j-1] 已经计算完毕。可以缓存，但是这个优化
在代码的 SymPy 版本中影响不大。
这里还需要指出的是，possible to visualize SymPy 的调用树参与了这些导数的计算。它实际上更大，但这是其中的一部分：
我们还可以使用 py-spy 模块生成火焰图，看看时间花在了哪里：
据我所知，34% 的时间花在 _eval_derivative_n_times 上，10% 的时间花在 assumptions.py 的函数 getit 上，12% 的在 subs(..) 中花费的时间，在 expand(..)
中花费的时间的 12%
优化#5
显然，当 pull request #13892 被合并到 SymPy 中时，它也引入了性能回归。
在 comments regarding that regression, Ondrej Certik recommends 之一中，使用 SymEngine 来提高大量使用导数的代码的性能。
所以我将提到的代码移植到 SymEngine.py 并注意到它 运行s 98 倍 比 [=43 的 SymPy 版本快=]（4320 倍 power=30）
可以通过pip3 install --user symengine安装所需的模块。
#!/usr/bin/python3
from symengine import *
import pprint
x=var("x")
b=Function("b")(x)
g=Function("g")(x)
b3,b2=b.diff(x,3),b.diff(x,2)
g4,g3=g.diff(x,4),g.diff(x,3)
def set_to_zero(e):
e = e.subs({b3:0,b2:0,g4:0,g3:0})
return e
def sum_of_coef(f):
s = 0
if f.func == Add:
for sum_term in f.args:
res = 1
if len(sum_term.args) == 0:
s += res
continue
first = sum_term.args[0]
if first.is_Number == True:
res = first
else:
res = 1
s += res
elif f.func == Mul:
first = f.args[0]
if first.is_Number == True:
s = first
else:
s = 1
elif f.func == Pow:
s = 1
return s
def main():
power = 8
charar = [[0] * (power*2) for x in range(power)]
coef_sum_array = [[0] * (power*2) for x in range(power)]
charar[0][0] = b
charar[0][1] = g
init_printing()
for i in range(1, power):
jmax = (i+1)*2
for j in range(0, jmax):
c2,c1,c0 = charar[i-1][j-2],charar[i-1][j-1],charar[i-1][j]
#print(c2,c1,c0)
if j == 0:
expr = b*c0.diff(x) + g*c0.diff(x,2)
elif j == 1:
expr = b*c1 + 2*g*c1.diff(x) + b*c0.diff(x) + g*c0.diff(x,2)
elif j == jmax-2:
expr = g*c2 + b*c1 + 2*g*c1.diff(x)
elif j == jmax-1:
expr = g*c2
else:
expr = g*c2 + b*c1 + 2*g*c1.diff(x) + b*c0.diff(x) + g*c0.diff(x,2)
charar[i][j] = set_to_zero(expand(expr))
coef_sum_array[i][j] = sum_of_coef(charar[i][j])
pprint.pprint(Matrix(coef_sum_array))
main()
优化后的性能 #5
我认为查看 c[i][j] 中的术语数量以确定表达式增长的速度会非常有趣。这肯定有助于估计当前代码的复杂性。
但出于实际目的，我绘制了上面 SymEngine 代码的当前时间和内存消耗，并设法得到以下图表：
时间和内存似乎都随着输入（原始代码中的 power 参数）呈多项式增长。
相同的图表，但作为 log-log plot 可以在这里查看：
就像 wiki page says, a straight line on a log-log plot corresponds to a monomial. This offers a way to recover 单项式的指数。
所以如果我们考虑两个点 N=16 和 N=32，它们之间的对数图看起来像一条直线
import pandas as pd
df=pd.read_csv("modif6_bench.txt", sep=',',header=0)
def find_slope(col1,col2,i1,i2):
xData = df[col1].to_numpy()
yData = df[col2].to_numpy()
x0,x1 = xData[i1],xData[i2]
y0,y1 = yData[i1],yData[i2]
m = log(y1/y0)/log(x1/x0)
return m
print("time slope = {0:0.2f}".format(find_slope("N","time",16,32)))
print("memory slope = {0:0.2f}".format(find_slope("N","memory",16,32)))
输出：
time slope = 5.69
memory slope = 2.62
time complexity would be O(n^5.69) and an approximation of space complexity 的非常粗略的近似值是 O(2^2.62)。
这里有more details关于判断增长率是多项式还是指数的（它涉及绘制半对数和双对数图，并查看数据显示为直线的位置） .
已定义 b 和 g 函数的性能
在第一个原始代码块中，函数b和g是未定义的函数。这意味着 SymPy 和 SymEngine 对它们一无所知。
第二个原始代码块定义了 b=1+x 和 g=1+x+x**2 。如果我们再次使用已知的 b 和 g 运行 所有这些，代码 运行s 会快得多，并且 运行ning 时间曲线和内存使用曲线比未知功能更好
time slope = 2.95
memory slope = 1.35
记录的数据符合已知的增长率
我想更多地了解匹配观察到的资源消耗（时间和内存），所以我将以下 Python module that fits each growth rate (from a catalog of common such growth rates) 写入记录的数据，然后向用户显示该图。
可以通过pip3 install --user matchgrowth
安装
当运行这样时：
match-growth.py --infile ./tests/modif7_bench.txt --outfile time.png --col1 N --col2 time --top 1
它生成资源使用情况的图表，以及与之匹配的最接近的增长率。在这种情况下，它发现多项式增长最接近：
其他注意事项
如果你运行这个power=8（在上面提到的symengine代码中）系数将如下所示：
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 5, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 17, 40, 31, 9, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 53, 292, 487, 330, 106, 16, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 161, 1912, 6091, 7677, 4693, 1520, 270, 25, 1, 0, 0, 0, 0, 0, 0]
[1, 485, 11956, 68719, 147522, 150706, 83088, 26573, 5075, 575, 36, 1, 0, 0, 0, 0]
[1, 1457, 73192, 735499, 2568381, 4118677, 3528928, 1772038, 550620, 108948, 13776, 1085, 49, 1, 0, 0]
[1, 4373, 443524, 7649215, 42276402, 102638002, 130209104, 96143469, 44255170, 13270378, 2658264, 358890, 32340, 1876, 64, 1]
事实证明，第 2 列与 A048473 重合，根据 OEIS 是 “Sierpiński 三角形中的三角形（所有大小，包括孔）的数量在 n 个铭文之后".
所有相关代码也可在 this repo.
中找到
第 i 行多项式系数与第 (i-1) 行多项式系数之间的关系
在前面postc[i][j]是计算出来的。可以检查 deg(c[i][j])=j+1 .
这可以通过初始化一个单独的二维数组来检查，并像这样计算度数：
deg[i][j] = degree(poly(parse_expr(str(charar[i][j]))))
垂直公式：
然后如果我们用 u(i,j,k) 表示 c[i][j] 中 x^k 的系数，我们可以尝试根据 u(i-1,_,_) 找到 u(i,j,k) 的公式. u(i,j,_) 的公式与 u(i+1,j,_)（以及所有后续行）的公式相同，因此有一些缓存机会。
水平公式：
同样有趣的是，当我们修复 i 时，发现 u(i,j,_) 的公式看起来与 u(i,j+1,_) 的公式相同，除了 k 的最后 3 个值。但我不确定这是否可以利用。
上述缓存步骤可能有助于跳过不必要的计算。
查看更多about this here.
关于解析、闭式解和渐近的一些注释
I'm struggling to derive this analytically
是的，这似乎很难。与此处提到的递归序列最接近的 class 称为 Holonomic sequences (also called D-finite or P-recursive). The sequence c[i][j] is not C-finite because it has polynomial coefficients (in the general case even the asymptotics of recurrences with polynomial coefficients is an open problem).
但是，c[i][j] 的递归关系由于导数的原因不符合此要求。如果我们要在 c[i][j] 的公式中省略导数，那么它就符合完整序列的条件。以下是我为这些问题找到解决方案的一些地方：
- "The Concrete Tetrahedron: Symbolic Sums, Recurrence Equations, Generating Functions, Asymptotic Estimates" by Kauers and Paule - 第 7 章完整序列和幂级数
- Analytic Combinatorics by Flajolet and Sedgewick - 附录 B.4 完整函数
但是 c[i][j] 也是一个多变量递归，所以这是它不符合上述理论的另一个原因。
然而，还有另一本书叫做 Analytic Combinatorics in Several Variables by Robin Pemantle and Mark C. Wilson，它确实处理了几个变量递归。
上面提到的所有书籍都需要大量复杂的分析，而且它们远远超出了我目前所知道的小数学，因此希望对这种数学有更深入了解的人可以尝试一下。
具有与生成函数相关的操作并且可以对此类序列进行操作的最先进的 CAS 是 Maple and the gfun package gfun repo（目前仅处理单变量情况）。
在我的研究中，我试图解决 Kolmogorov 逆向方程，即对 $$Af = b(x)f'(x)+\sigma(x)f''(x)$$
对于特定的 b(x) 和 \sigma(x)，我试图了解在计算更高的 Af 幂时表达式的系数增长的速度。我很难从分析上得出这一点，因此试图从经验上看趋势。
首先，我用过sympy:
from sympy import *
import matplotlib.pyplot as plt
import re
import math
import numpy as np
import time
np.set_printoptions(suppress=True)
x = Symbol('x')
b = Function('b')(x)
g = Function('g')(x)
def new_coef(gamma, beta, coef_minus2, coef_minus1, coef):
return expand(simplify(gamma*coef_minus2 + beta*coef_minus1 + 2*gamma*coef_minus1.diff(x)\
+beta*coef.diff(x)+gamma*coef.diff(x,2)))
def new_coef_first(gamma, beta, coef):
return expand(simplify(beta*coef.diff(x)+gamma*coef.diff(x,2)))
def new_coef_second(gamma, beta, coef_minus1, coef):
return expand(simplify(beta*coef_minus1 + 2*gamma*coef_minus1.diff(x)\
+beta*coef.diff(x)+gamma*coef.diff(x,2)))
def new_coef_last(gamma, beta, coef_minus2):
return lambda x: gamma(x)*coef_minus2(x)
def new_coef_last(gamma, beta, coef_minus2):
return expand(simplify(gamma*coef_minus2 ))
def new_coef_second_to_last(gamma, beta, coef_minus2, coef_minus1):
return expand(simplify(gamma*coef_minus2 + beta*coef_minus1 + 2*gamma*coef_minus1.diff(x)))
def set_to_zero(expression):
expression = expression.subs(Derivative(b, x, x, x), 0)
expression = expression.subs(Derivative(b, x, x), 0)
expression = expression.subs(Derivative(g, x, x, x, x), 0)
expression = expression.subs(Derivative(g, x, x, x), 0)
return expression
def sum_of_coef(expression):
sum_of_coef = 0
for i in str(expression).split(' + '):
if i[0:1] == '(':
i = i[1:]
integers = re.findall(r'\b\d+\b', i)
if len(integers) > 0:
length_int = len(integers[0])
if i[0:length_int] == integers[0]:
sum_of_coef += int(integers[0])
else:
sum_of_coef += 1
else:
sum_of_coef += 1
return sum_of_coef
power = 6
charar = np.zeros((power, power*2), dtype=Symbol)
coef_sum_array = np.zeros((power, power*2))
charar[0,0] = b
charar[0,1] = g
coef_sum_array[0,0] = 1
coef_sum_array[0,1] = 1
for i in range(1, power):
#print(i)
for j in range(0, (i+1)*2):
#print(j, ':')
#start_time = time.time()
if j == 0:
charar[i,j] = set_to_zero(new_coef_first(g, b, charar[i-1, j]))
elif j == 1:
charar[i,j] = set_to_zero(new_coef_second(g, b, charar[i-1, j-1], charar[i-1, j]))
elif j == (i+1)*2-2:
charar[i,j] = set_to_zero(new_coef_second_to_last(g, b, charar[i-1, j-2], charar[i-1, j-1]))
elif j == (i+1)*2-1:
charar[i,j] = set_to_zero(new_coef_last(g, b, charar[i-1, j-2]))
else:
charar[i,j] = set_to_zero(new_coef(g, b, charar[i-1, j-2], charar[i-1, j-1], charar[i-1, j]))
#print("--- %s seconds for expression---" % (time.time() - start_time))
#start_time = time.time()
coef_sum_array[i,j] = sum_of_coef(charar[i,j])
#print("--- %s seconds for coeffiecients---" % (time.time() - start_time))
coef_sum_array
然后，研究自动微分并使用 autograd：
import autograd.numpy as np
from autograd import grad
import time
np.set_printoptions(suppress=True)
b = lambda x: 1 + x
g = lambda x: 1 + x + x**2
def new_coef(gamma, beta, coef_minus2, coef_minus1, coef):
return lambda x: gamma(x)*coef_minus2(x) + beta(x)*coef_minus1(x) + 2*gamma(x)*grad(coef_minus1)(x)\
+beta(x)*grad(coef)(x)+gamma(x)*grad(grad(coef))(x)
def new_coef_first(gamma, beta, coef):
return lambda x: beta(x)*grad(coef)(x)+gamma(x)*grad(grad(coef))(x)
def new_coef_second(gamma, beta, coef_minus1, coef):
return lambda x: beta(x)*coef_minus1(x) + 2*gamma(x)*grad(coef_minus1)(x)\
+beta(x)*grad(coef)(x)+gamma(x)*grad(grad(coef))(x)
def new_coef_last(gamma, beta, coef_minus2):
return lambda x: gamma(x)*coef_minus2(x)
def new_coef_second_to_last(gamma, beta, coef_minus2, coef_minus1):
return lambda x: gamma(x)*coef_minus2(x) + beta(x)*coef_minus1(x) + 2*gamma(x)*grad(coef_minus1)(x)
power = 6
coef_sum_array = np.zeros((power, power*2))
coef_sum_array[0,0] = b(1.0)
coef_sum_array[0,1] = g(1.0)
charar = [b, g]
for i in range(1, power):
print(i)
charar_new = []
for j in range(0, (i+1)*2):
if j == 0:
new_funct = new_coef_first(g, b, charar[j])
elif j == 1:
new_funct = new_coef_second(g, b, charar[j-1], charar[j])
elif j == (i+1)*2-2:
new_funct = new_coef_second_to_last(g, b, charar[j-2], charar[j-1])
elif j == (i+1)*2-1:
new_funct = new_coef_last(g, b, charar[j-2])
else:
new_funct = new_coef(g, b, charar[j-2], charar[j-1], charar[j])
coef_sum_array[i,j] = new_funct(1.0)
charar_new.append(new_funct)
charar = charar_new
coef_sum_array
但是，我对他们中任何一个的速度都不满意。我想进行至少一千次迭代，而在 运行 simpy 方法 3 天后，我得到了 30 :/
我希望可以优化第二种方法（数值）以避免每次都重新计算表达式。不幸的是，我自己看不到该解决方案。另外，我试过 Maple，但还是没有成功。
概览
所以，这里有两个关于导数的公式很有趣：
- Faa di Bruno's formula which is a way to quickly find the n-th derivative of
f(g(x)), and looks a lot like the Multinomial theorem - General Leibniz rule which is a way to quickly find the n-th derivative of
f(x)*g(x)and looks a lot like the Binomial theorem
这两个都在 pull request #13892 中讨论过 n 阶导数使用一般莱布尼兹规则加速。
I'm trying to see how fast the coefficients of the expression are growing
在您的代码中，计算 c[i][j] 的一般公式是这样的：
c[i][j] = g * c[i-1][j-2] + b * c[i-1][j-1] + 2 * g * c'[i-1][j-1] + g * c''[i-1][j]
（其中 c'[i][j] 和 c''[i][j] 是 c[i][j] 的一阶和二阶导数）
因此，根据上面提到的莱布尼茨法则，我直觉地认为，计算出的系数应该与Pascal's triangle有关（或者至少它们应该有某种组合关系）。
优化 #1
在原始代码中，函数sum_to_coef(f)将表达式f序列化为一个字符串，然后丢弃所有看起来不像数字的东西，然后对剩余数字求和。
我们可以通过遍历表达式树并收集我们需要的东西来避免这里的序列化
def sum_of_coef(f):
s = 0
if f.func == Add:
for sum_term in f.args:
res = sum_term if sum_term.is_Number else 1
if len(sum_term.args) == 0:
s += res
continue
first = sum_term.args[0]
if first.is_Number == True:
res = first
else:
res = 1
s += res
elif f.func == Mul:
first = f.args[0]
if first.is_Number == True:
s = first
else:
s = 1
elif f.func == Pow:
s = 1
return s
优化#2
在函数set_to_zero(expr)中，b的所有二阶和三阶导数，以及g的三阶和四阶导数都被零替换。
我们可以像这样将所有这些替换折叠成一个语句：
b3,b2=b.diff(x,3),b.diff(x,2)
g4,g3=g.diff(x,4),g.diff(x,3)
def set_to_zero(expression):
expression = expression.subs({b3:0,b2:0,g4:0,g3:0})
return expression
优化#3
在原始代码中，对于每个单元格 c[i][j] 我们调用 simplify。事实证明这对性能有很大影响，但实际上我们可以跳过这个调用，因为幸运的是我们的表达式只是导数或未知函数的乘积之和。
所以行
charar[i,j] = set_to_zero(expand(simplify(expr)))
变成
charar[i,j] = set_to_zero(expand(expr))
优化#4
也尝试了以下方法，但效果很小。
对于两个连续的 j 值，我们计算 c'[i-1][j-1] 两次。
j-1 c[i-1][j-3] c[i-1][j-2] c[i-1][j-1]
j c[i-1][j-2] c[i-1][j-1] c[i-1][j]
如果您查看 else 分支中的循环公式，您会发现 c'[i-1][j-1] 已经计算完毕。可以缓存，但是这个优化
在代码的 SymPy 版本中影响不大。
这里还需要指出的是，possible to visualize SymPy 的调用树参与了这些导数的计算。它实际上更大，但这是其中的一部分：
我们还可以使用 py-spy 模块生成火焰图，看看时间花在了哪里：
据我所知，34% 的时间花在 _eval_derivative_n_times 上，10% 的时间花在 assumptions.py 的函数 getit 上，12% 的在 subs(..) 中花费的时间，在 expand(..)
优化#5
显然，当 pull request #13892 被合并到 SymPy 中时，它也引入了性能回归。
在 comments regarding that regression, Ondrej Certik recommends 之一中，使用 SymEngine 来提高大量使用导数的代码的性能。
所以我将提到的代码移植到 SymEngine.py 并注意到它 运行s 98 倍 比 [=43 的 SymPy 版本快=]（4320 倍 power=30）
可以通过pip3 install --user symengine安装所需的模块。
#!/usr/bin/python3
from symengine import *
import pprint
x=var("x")
b=Function("b")(x)
g=Function("g")(x)
b3,b2=b.diff(x,3),b.diff(x,2)
g4,g3=g.diff(x,4),g.diff(x,3)
def set_to_zero(e):
e = e.subs({b3:0,b2:0,g4:0,g3:0})
return e
def sum_of_coef(f):
s = 0
if f.func == Add:
for sum_term in f.args:
res = 1
if len(sum_term.args) == 0:
s += res
continue
first = sum_term.args[0]
if first.is_Number == True:
res = first
else:
res = 1
s += res
elif f.func == Mul:
first = f.args[0]
if first.is_Number == True:
s = first
else:
s = 1
elif f.func == Pow:
s = 1
return s
def main():
power = 8
charar = [[0] * (power*2) for x in range(power)]
coef_sum_array = [[0] * (power*2) for x in range(power)]
charar[0][0] = b
charar[0][1] = g
init_printing()
for i in range(1, power):
jmax = (i+1)*2
for j in range(0, jmax):
c2,c1,c0 = charar[i-1][j-2],charar[i-1][j-1],charar[i-1][j]
#print(c2,c1,c0)
if j == 0:
expr = b*c0.diff(x) + g*c0.diff(x,2)
elif j == 1:
expr = b*c1 + 2*g*c1.diff(x) + b*c0.diff(x) + g*c0.diff(x,2)
elif j == jmax-2:
expr = g*c2 + b*c1 + 2*g*c1.diff(x)
elif j == jmax-1:
expr = g*c2
else:
expr = g*c2 + b*c1 + 2*g*c1.diff(x) + b*c0.diff(x) + g*c0.diff(x,2)
charar[i][j] = set_to_zero(expand(expr))
coef_sum_array[i][j] = sum_of_coef(charar[i][j])
pprint.pprint(Matrix(coef_sum_array))
main()
优化后的性能 #5
我认为查看 c[i][j] 中的术语数量以确定表达式增长的速度会非常有趣。这肯定有助于估计当前代码的复杂性。
但出于实际目的，我绘制了上面 SymEngine 代码的当前时间和内存消耗，并设法得到以下图表：
时间和内存似乎都随着输入（原始代码中的 power 参数）呈多项式增长。
相同的图表，但作为 log-log plot 可以在这里查看：
就像 wiki page says, a straight line on a log-log plot corresponds to a monomial. This offers a way to recover 单项式的指数。
所以如果我们考虑两个点 N=16 和 N=32，它们之间的对数图看起来像一条直线
import pandas as pd
df=pd.read_csv("modif6_bench.txt", sep=',',header=0)
def find_slope(col1,col2,i1,i2):
xData = df[col1].to_numpy()
yData = df[col2].to_numpy()
x0,x1 = xData[i1],xData[i2]
y0,y1 = yData[i1],yData[i2]
m = log(y1/y0)/log(x1/x0)
return m
print("time slope = {0:0.2f}".format(find_slope("N","time",16,32)))
print("memory slope = {0:0.2f}".format(find_slope("N","memory",16,32)))
输出：
time slope = 5.69
memory slope = 2.62
time complexity would be O(n^5.69) and an approximation of space complexity 的非常粗略的近似值是 O(2^2.62)。
这里有more details关于判断增长率是多项式还是指数的（它涉及绘制半对数和双对数图，并查看数据显示为直线的位置） .
已定义 b 和 g 函数的性能
在第一个原始代码块中，函数b和g是未定义的函数。这意味着 SymPy 和 SymEngine 对它们一无所知。
第二个原始代码块定义了 b=1+x 和 g=1+x+x**2 。如果我们再次使用已知的 b 和 g 运行 所有这些，代码 运行s 会快得多，并且 运行ning 时间曲线和内存使用曲线比未知功能更好
time slope = 2.95
memory slope = 1.35
记录的数据符合已知的增长率
我想更多地了解匹配观察到的资源消耗（时间和内存），所以我将以下 Python module that fits each growth rate (from a catalog of common such growth rates) 写入记录的数据，然后向用户显示该图。
可以通过pip3 install --user matchgrowth
当运行这样时：
match-growth.py --infile ./tests/modif7_bench.txt --outfile time.png --col1 N --col2 time --top 1
它生成资源使用情况的图表，以及与之匹配的最接近的增长率。在这种情况下，它发现多项式增长最接近：
其他注意事项
如果你运行这个power=8（在上面提到的symengine代码中）系数将如下所示：
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 5, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 17, 40, 31, 9, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 53, 292, 487, 330, 106, 16, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 161, 1912, 6091, 7677, 4693, 1520, 270, 25, 1, 0, 0, 0, 0, 0, 0]
[1, 485, 11956, 68719, 147522, 150706, 83088, 26573, 5075, 575, 36, 1, 0, 0, 0, 0]
[1, 1457, 73192, 735499, 2568381, 4118677, 3528928, 1772038, 550620, 108948, 13776, 1085, 49, 1, 0, 0]
[1, 4373, 443524, 7649215, 42276402, 102638002, 130209104, 96143469, 44255170, 13270378, 2658264, 358890, 32340, 1876, 64, 1]
事实证明，第 2 列与 A048473 重合，根据 OEIS 是 “Sierpiński 三角形中的三角形（所有大小，包括孔）的数量在 n 个铭文之后".
所有相关代码也可在 this repo.
中找到第 i 行多项式系数与第 (i-1) 行多项式系数之间的关系
在前面postc[i][j]是计算出来的。可以检查 deg(c[i][j])=j+1 .
这可以通过初始化一个单独的二维数组来检查，并像这样计算度数：
deg[i][j] = degree(poly(parse_expr(str(charar[i][j]))))
垂直公式：
然后如果我们用 u(i,j,k) 表示 c[i][j] 中 x^k 的系数，我们可以尝试根据 u(i-1,_,_) 找到 u(i,j,k) 的公式. u(i,j,_) 的公式与 u(i+1,j,_)（以及所有后续行）的公式相同，因此有一些缓存机会。
水平公式：
同样有趣的是，当我们修复 i 时，发现 u(i,j,_) 的公式看起来与 u(i,j+1,_) 的公式相同，除了 k 的最后 3 个值。但我不确定这是否可以利用。
上述缓存步骤可能有助于跳过不必要的计算。
查看更多about this here.
关于解析、闭式解和渐近的一些注释
I'm struggling to derive this analytically
是的，这似乎很难。与此处提到的递归序列最接近的 class 称为 Holonomic sequences (also called D-finite or P-recursive). The sequence c[i][j] is not C-finite because it has polynomial coefficients (in the general case even the asymptotics of recurrences with polynomial coefficients is an open problem).
但是，c[i][j] 的递归关系由于导数的原因不符合此要求。如果我们要在 c[i][j] 的公式中省略导数，那么它就符合完整序列的条件。以下是我为这些问题找到解决方案的一些地方：
- "The Concrete Tetrahedron: Symbolic Sums, Recurrence Equations, Generating Functions, Asymptotic Estimates" by Kauers and Paule - 第 7 章完整序列和幂级数
- Analytic Combinatorics by Flajolet and Sedgewick - 附录 B.4 完整函数
但是 c[i][j] 也是一个多变量递归，所以这是它不符合上述理论的另一个原因。
然而，还有另一本书叫做 Analytic Combinatorics in Several Variables by Robin Pemantle and Mark C. Wilson，它确实处理了几个变量递归。
上面提到的所有书籍都需要大量复杂的分析，而且它们远远超出了我目前所知道的小数学，因此希望对这种数学有更深入了解的人可以尝试一下。
具有与生成函数相关的操作并且可以对此类序列进行操作的最先进的 CAS 是 Maple and the gfun package gfun repo（目前仅处理单变量情况）。