将二维多项式系数附加到 Python 中的项的更快方法?

Faster way to attach 2d polynomial coefficients to terms in Python?

所以我试图通过将各自的系数 (k_ij) 附加到各自的单项式 (x**i*y**j,其中 xy 是符号变量)。我的目标是 最小化计算时间 因为我的多项式非常大,并且在我的程序中使用计算时间的主要事情是生成这个符号多项式的一组行(我将需要多次调用此步骤)。考虑到我的程序的其余部分相当 lengthy/complex,我很惊讶地意识到这个步骤在程序中花费了多少时间。我永远不会猜到简单地创建单项式并附加系数以创建多项式需要多少时间。

我的多项式不一定要保存为多项式,它只需要是所有项的总和即可。这样,我的意思是我不需要这样的输出...

Out: Poly(7*x + 2*y, x, y, domain='ZZ')

...不过,我也不反对这种格式,只要它能最大限度地减少计算时间。我想要但不一定需要的只是一个看起来像这样的输出(在这种情况下,可以通过简单地说 z = 7*x + 2*y if x and y 来实现已经用符号定义了):

Out: 7*x + 2*y

所以,我有一个系数矩阵(可以很容易地重新排序以适应用于将其附加到多项式的方法),其中包含所需多项式的所有系数,并且我有我的符号变量。这是我第一次尝试的重新创建(计算时间最长):

import sympy
import time
# let's time it, see how long it takes
from time import clock as tc

# time it!
t0 = tc()

order = 33
order_x = order
order_y = order
deg_x = order - 1
deg_y = order - 1

nnn = order_x*order_y

# here is a random coefficient matrix
coefficient_matrix = numpy.random.rand(nnn)

# define symbolic variables
x = sympy.Symbol('x')
y = sympy.Symbol('y')

# now let's populate the surface polynomial
z = 0
for i2 in range(order_y):
    for i3 in range(order_x):
        z = z + coefficient_matrix[i3 + order_x*i2]*(x**(deg_x - i3))*(y**(deg_y - i2))
        # note that this returns high order to low order terms

t1 = tc()
print(t1-t0)

那很慢(37 秒),大概是因为 for 循环的性质。然后我尝试使用 numpy.polynomial.polynomial.polyvander2d 生成一个伪范德蒙矩阵并将其乘以系数矩阵,但这对计算时间几乎没有影响。我尝试的下一个方法(如下所示)能够大大缩短计算时间:

import sympy
import time
from time import clock as tc
import numpy
from numpy.polynomial.polynomial import polyval2d as P2

# time it!
t0 = tc()

order = 33
order_x = order
order_y = order

# here is a random coefficient matrix
coefficient_matrix = numpy.random.rand(order, order)

# define symbolic variables
x = sympy.Symbol('x')
y = sympy.Symbol('y')

# create polynomial
z = P2(x, y, coefficient_matrix)

# make the polynomial a logical sequence of monomials
z = sympy.expand(z)

t1 = tc()
print(z)
print(t1-t0)

这个方法用了8.5秒(这着实让我吃惊,我原以为会短很多),但是没有z = sympy.expand(z)这行大约用了3秒。我有那行的原因是稍后我需要修改函数并提取新系数,所以我希望它以扩展形式供以后使用(以便它以上面列出的格式出现;如果我没有包括这个行,它以某种分解格式出现)。

有没有办法让 Python 更快地将项与系数匹配,并且 return 它们作为单项式序列?

Constructing a sy.Poly 会花费更少的时间:

using_loop          : 43.56
using_P2            : 12.64
using_poly          :  0.03

from timeit import default_timer as tc
import numpy as np
import sympy as sy
from numpy.polynomial.polynomial import polyval2d as P2


def using_poly(coefficient_matrix, S=sy.S):
    order = coefficient_matrix.shape[0]
    x = sy.Symbol('x')
    y = sy.Symbol('y')
    dct = {i:S(val) for i, val in np.ndenumerate(coefficient_matrix)}
    z = sy.Poly(dct, x, y)
    return z

def using_loop(coefficient_matrix):
    order = coefficient_matrix.shape[0]
    coefficient_matrix = coefficient_matrix.T.ravel()[::-1]
    order_x = order
    order_y = order
    deg_x = order - 1
    deg_y = order - 1
    x = sy.Symbol('x')
    y = sy.Symbol('y')
    z = 0
    for i2 in range(order_y):
        for i3 in range(order_x):
            z = z + coefficient_matrix[i3 + order_x*i2]*(x**(deg_x - i3))*(y**(deg_y - i2))
    return z

def using_P2(coefficient_matrix):
    x = sy.Symbol('x')
    y = sy.Symbol('y')
    z = P2(x, y, coefficient_matrix)
    # make the polynomial a logical sequence of monomials
    z = sy.expand(z)
    return z

order = 33
np.random.seed(2015)
coefficient_matrix = np.random.rand(order, order)
# coefficient_matrix = np.arange(1, order*order+1).reshape(order,order)

for func in (using_loop, using_P2, using_poly):
    t0 = tc()
    func(coefficient_matrix)
    t1 = tc()
    print('{:15s}: {:>5.2f}'.format(func.__name__, t1-t0))

这是一个小 coefficient_matrix:

输出的示例
In [277]: order = 3

In [278]: coefficient_matrix = np.arange(1, order*order+1).reshape(order,order)

In [279]: using_poly(coefficient_matrix)
Out[279]: Poly(9*x**2*y**2 + 8*x**2*y + 7*x**2 + 6*x*y**2 + 5*x*y + 4*x + 3*y**2 + 2*y + 1, x, y, domain='ZZ')

In [280]: using_P2(coefficient_matrix)
Out[280]: 9.0*x**2*y**2 + 8.0*x**2*y + 7.0*x**2 + 6.0*x*y**2 + 5.0*x*y + 4.0*x + 3.0*y**2 + 2.0*y + 1.0

In [281]: using_loop(coefficient_matrix)
Out[281]: 9*x**2*y**2 + 8*x**2*y + 7*x**2 + 6*x*y**2 + 5*x*y + 4*x + 3*y**2 + 2*y + 1