递归符号计算可以提高性能

Question

在我的研究中，我试图解决科尔莫戈罗夫倒向方程，即对 $$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

然而，我对它们的速度都不满意。我希望至少进行1000次迭代，而在运行simpy方法3天后，我得到了30:/

我希望第二种方法（数值）可以优化，以避免每次都重新计算表达式。不幸的是，我自己看不到这个解决方案。另外，我也试过枫树，但又一次运气不佳

Answer 1

概述

这里有两个关于导数的公式很有趣：

1. Faa di Bruno's formula这是一种快速找到f(g(x))的n阶导数的方法，看起来很像Multinomial theorem
2. General Leibniz rule是一种快速找到f(x)*g(x)的n阶导数的方法，看起来很像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 visualizeSymPy的调用树。它实际上更大，但这里是它的一部分：

我们还可以使用py-spy模块生成火焰图，以查看花费的时间：

据我所知，花在{}上的时间占34%，花在函数{}上的时间占10%，花在{}上的时间占12%，花在{}上的时间占12%

优化#5

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

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

因此，我将提到的代码移植到SymEngine.py，并注意到它的运行速度比power=8的Symphy版本快98倍（而power=30的运行速度则快4320倍）

所需模块可通过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所说，对数图上的直线对应于单项式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的非常粗略的近似值是O(n^5.69)，而space complexity的近似值是O(2^2.62)

这里有{a19}关于决定增长率是多项式还是指数（它涉及绘制半对数和对数图，并查看数据显示为直线的位置）

使用已定义的b和g函数的性能

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

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

time slope = 2.95
memory slope = 1.35

记录数据拟合已知增长率

我想进一步研究匹配观察到的资源消耗（时间和内存），因此我编写了以下Python module，将每个增长率（从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（在上面提到的符号引擎代码中）运行此命令，则系数将如下所示：

[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]

因此，事实证明，第二列与A048473重合，根据OEIS，这是“在n个铭文之后，塞尔皮斯基三角形中三角形的数量（各种尺寸，包括孔）。

这方面的所有代码也可以在this repo中找到

Answer 2

第i条线的多项式系数与第（i-1）条线的系数之间的关系

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

这可以通过初始化一个单独的2d数组并按如下方式计算度来检查：

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

垂直公式：

然后，如果我们用u(i,j,k)表示x^k在c[i][j]中的系数，我们可以尝试找到u(i,j,k)在u(i-1,_,_)中的公式。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

是的，这似乎很难。与这里提到的递归序列相关的最接近的一类被称为Holonomic sequences（也称为D-有限或P-recursive）。序列c[i][j]不是C-finite，因为它具有多项式系数（在一般情况下，甚至多项式系数的递归的渐近性也是an open problem）

但是c[i][j]的递归关系不符合此条件，因为存在导数。如果我们在c[i][j]公式中省略导数，那么它将符合完整序列的条件。以下是我找到这些解决方案的一些地方：

1. "The Concrete Tetrahedron: Symbolic Sums, Recurrence Equations, Generating Functions, Asymptotic Estimates" by Kauers and Paule-第7章完整序列和幂级数
2. Analytic Combinatorics by Flajolet and Sedgewick-附录B.4完整函数

但是{}也是一个多变量递归，所以这是它不符合上述理论的另一个原因

然而，还有一本书叫做Analytic Combinatorics in Several Variables by Robin Pemantle and Mark C. Wilson，它确实处理了几个变量的重复出现

上面提到的所有书籍都需要很多复杂的分析，它们远远超出了我目前所知道的小数学，所以希望对这类数学有更深入理解的人可以尝试一下

具有生成函数相关操作并可以对此类序列进行操作的最高级CA是Maple和gfun package{a11}（目前仅处理单变量情况）