递归符号计算 - 提高性能

如何解决递归符号计算 - 提高性能

在我的研究中，我试图解决 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,coef):
    return expand(simplify(beta*coef.diff(x)+gamma*coef.diff(x,2)))
def new_coef_second(gamma,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,coef_minus2):
    return lambda x: gamma(x)*coef_minus2(x)
def new_coef_last(gamma,coef_minus2):
    return expand(simplify(gamma*coef_minus2 ))
def new_coef_second_to_last(gamma,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),0)
    expression = expression.subs(Derivative(b,0)
    expression = expression.subs(Derivative(g,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,j-1],j]))
        elif j == (i+1)*2-2:
            charar[i,j] = set_to_zero(new_coef_second_to_last(g,j-2],j-1]))
        elif j == (i+1)*2-1:
            charar[i,j] = set_to_zero(new_coef_last(g,j-2]))
        else:
            charar[i,j] = set_to_zero(new_coef(g,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,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,coef):
    return lambda x: beta(x)*grad(coef)(x)+gamma(x)*grad(grad(coef))(x)
def new_coef_second(gamma,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,coef_minus2):
    return lambda x: gamma(x)*coef_minus2(x)
def new_coef_second_to_last(gamma,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,charar[j])
        elif j == 1:
            new_funct = new_coef_second(g,charar[j-1],charar[j])
        elif j == (i+1)*2-2:
            new_funct = new_coef_second_to_last(g,charar[j-2],charar[j-1])
        elif j == (i+1)*2-1:
            new_funct = new_coef_last(g,charar[j-2])
        else:
            new_funct = new_coef(g,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 这是一种快速找到 f(g(x)) 的 n 阶导数的方法，看起来很像 Multinomial theorem
General Leibniz rule 是一种快速找到 f(x)*g(x) 的 n 阶导数的方法，看起来很像 Binomial theorem

这两个都在 pull request #13892 中讨论过，使用一般莱布尼茨规则加速了 n 阶导数。

我想看看表达式的系数增长有多快

在您的代码中，计算 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% 的时间花在 getit 的函数 assumptions.py 上，12% 的时间花在了 subs(..) 上在 expand(..) 中，12% 的时间花在 power=8

优化 #5

显然，当 pull request #13892 被合并到 SymPy 时，它也引入了性能回归。

在 comments regarding that regression,Ondrej Certik recommends 之一中，使用 SymEngine 来提高大量使用导数的代码的性能。

因此，我将提到的代码移植到 SymEngine.py 并注意到它的运行速度比 power=30（和 43204320 SymPy 版本快 98 倍强>pip3 install --user symengine} 倍快）

所需的模块可以通过 #!/usr/bin/python3 from symengine import * import pprint x=var("x") b=Function("b")(x) g=Function("g")(x) b3,3) def set_to_zero(e): e = e.subs({b3: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,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() 安装。

c[i][j]

优化后的性能 #5

我认为查看 power 中的术语数量以确定表达式的增长速度会非常有趣。这肯定有助于估计当前代码的复杂性。

但出于实际目的，我已经绘制了上面 SymEngine 代码的当前时间和内存消耗，并设法得到了以下图表：

performance_modif6

时间和内存似乎都随着输入（原始代码中的 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",32))) 参数）呈多项式增长。

可以在此处查看相同的图表，但作为 log-log plot：

performance_modif6_loglog

就像 wiki page 所说的，对数图上的一条直线对应于单项式。 This offers a way to recover 单项式的指数。

因此，如果我们考虑两个点 N=16 和 N=32，它们之间的对数图看起来像一条直线

time slope = 5.69
memory slope = 2.62

输出：

O(n^5.69)

time complexity 的非常粗略的近似值是 O(2^2.62)，space complexity 的近似值是 b。

这里有 more details 关于决定增长率是多项式还是指数的（它涉及绘制半对数和对数对数图，并查看数据在何处显示为直线）。

具有定义的 `g` 和 `b` 函数的性能

在第一个原始代码块中，函数 g 和 b=1+x 是未定义的函数。这意味着 SymPy 和 SymEngine 对它们一无所知。

第二个原始代码块定义了 g=1+x+x**2 和 b 。如果我们用已知的 g 和 time slope = 2.95 memory slope = 1.35 再次运行所有这些代码，代码运行得更快，并且运行时间曲线和内存使用曲线比未知函数要好

pip3 install --user matchgrowth

拟合已知增长率的记录数据

我想多看看匹配观察到的资源消耗（时间和内存），所以我写了以下 Python module 来拟合每个增长率（从 catalog of common such growth rates）到记录的数据，然后向用户展示绘图。

可以通过match-growth.py --infile ./tests/modif7_bench.txt --outfile time.png --col1 N --col2 time --top 1

安装

当这样运行时：

power=8

它生成资源使用情况的图表，以及与之匹配的最接近的增长率。在这种情况下，它会发现多项式增长最接近：

其他注意事项

如果您为 [0,0] [1,5,4,1,17,40,31,9,53,292,487,330,106,161,1912,6091,7677,4693,1520,270,25,485,11956,68719,147522,150706,83088,26573,5075,575,36,1457,73192,735499,2568381,4118677,3528928,1772038,550620,108948,13776,1085,49,4373,443524,7649215,42276402,102638002,130209104,96143469,44255170,13270378,2658264,358890,32340,1876,64,1] 运行此程序（在上面提到的 symengine 代码中），系数将如下所示：

{{1}}

事实证明，第 2 列与 A048473 重合，根据 OEIS 是 “n 个铭文后谢尔宾斯基三角形中的三角形（所有大小，包括洞）的数量”.

this repo 中也提供了所有这些代码。

第 i 行的多项式系数与第 (i-1) 行的系数之间的关系

在上一篇文章中计算了 c[i][j]。可以检查 deg(c[i][j])=j+1 .

这可以通过初始化一个单独的二维数组来检查，并像这样计算度数：

deg[i][j] = degree(poly(parse_expr(str(charar[i][j]))))

垂直公式：

那么如果我们用u(i,j,k)表示x^k中c[i][j]的系数，我们可以尝试根据u(i,k)找到u(i-1,_,_)的公式。 u(i,_) 的公式将与 u(i+1,_)（以及所有后续行）的公式相同，因此有一些缓存的机会。

水平公式：

同样有趣的是，当我们修正 i 时，发现 u(i,_) 的公式与 u(i,j+1,_) 的公式看起来一样，除了 k 的最后 3 个值。但我不确定这是否可以利用。

上面提到的缓存步骤可能有助于跳过不必要的计算。

查看更多about this here。

关于分析、闭式解和渐近的一些说明

我正在努力分析得出这一点

是的，这似乎很难。与这里提到的最接近的一类递归序列称为 Holonomic sequences（也称为 D-finite 或 P-recursive）。序列 c[i][j] 不是 C-finite，因为它有多项式系数（在一般情况下，即使是多项式系数的递推渐近线也是 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 和 gfun package gfun repo（目前仅处理单变量情况).