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

Question

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

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

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

...不过，我也不反对这种格式，只要它能最大限度地减少计算时间。我想要但不一定需要的只是一个看起来像这样的输出(在这种情况下，可以通过简单地说 z = 7*x + 2*y if x 和 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 更快地将项与系数匹配并将它们作为单项式序列返回？

Answer 1

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