在 Sympy 中因式分解多项式

Question

我正在进行一个非常简单的概率计算，从A-Z集合中获取X、Y、Z的子集（相应的概率为x、y、z）。

由于公式非常复杂，为了处理它们，我尝试使用sympy来简化（或收集或因式分解-我不知道确切的定义）这些多项式表达式。

所以...有了这个（一个非常简单的概率计算表达式，从A-Z集合中获取X、Y、Z的子集，并具有相应的概率x、y、z）。

import sympy as sp

x, y, z = sp.symbols('x y z')

expression = (
    x * (1 - x) * y * (1 - x - y) * z +
    x * (1 - x) * z * (1 - x - z) * y +

    y * (1 - y) * x * (1 - y - x) * z +
    y * (1 - y) * z * (1 - y - z) * x +

    z * (1 - z) * y * (1 - z - y) * x +
    z * (1 - z) * x * (1 - z - x) * y
)

我希望能得到类似于这样的东西

x * y * z * (6 * (1 - x - y - z) + (x + y) ** 2 + (y + z) ** 2 + (x + z) ** 2)

一个多项式，重写为尽可能少的操作（+、-、*、**等）。

---

我尝试使用 factor()、collect()、simplify()。但结果与我的预期不同。通常我得到的是：

2*x*y*z*(x**2 + x*y + x*z - 3*x + y**2 + y*z - 3*y + z**2 - 3*z + 3)

我知道sympy可以将多项式合并为简单形式：

我知道sympy可以合并多项式成为简单形式：

sp.factor(x**2 + 2*x*y + y**2)  # gives (x + y)**2

但是如何使sympy将上面的表达式中的多项式合并？

---

如果在sympy中无法完成此任务，也许还有其他选择？

- akaRem

Answer 1

将一些方法组合起来，这次偶然得到了一个不错的答案。有趣的是看看这种策略在你所生成的方程中是否经常奏效，或者像其名称所示，这只是这一次的幸运结果。

def iflfactor(eq):
    """Return the "I'm feeling lucky" factored form of eq."""
    e = Mul(*[horner(e) if e.is_Add else e for e in
        Mul.make_args(factor_terms(expand(eq)))])
    r, e = cse(e)
    s = [ri[0] for ri in r]
    e = Mul(*[collect(ei.expand(), s) if ei.is_Add else ei for ei in
        Mul.make_args(e[0])]).subs(r)
    return e

>>> iflfactor(eq)  # using your equation as eq
2*x*y*z*(x**2 + x*y + y**2 + (z - 3)*(x + y + z) + 3)
>>> _.count_ops()
15

顺便提一下，factor_terms和gcd_terms之间的区别在于factor_terms会更加努力地提取共同项，同时保留表达式的原始结构，很像你手动进行的操作（即在添加中寻找可以提取的公共项）。

>>> factor_terms(x/(z+z*y)+x/z)
x*(1 + 1/(y + 1))/z
>>> gcd_terms(x/(z+z*y)+x/z)
x*(y*z + 2*z)/(z*(y*z + z))

说句实话，

克里斯

- smichr

---

不错的组合，但很难理解。你能解释一下算法/思路吗？顺便说一句，我发现了一个明显的错误：x * (1 - x) * y * (1 - x - y) * z + ... -> x / (1 - x) * y / (1 - x - y) * z + ...，对于这样的方程，你的组合不起作用（我想这是因为显而易见的事情，但由于我不知道算法...）。 - akaRem

不确定这里的根本原因是什么，但当我有大量变量时，iflfactor最终无法执行所有必要的替换以返回原始变量。我必须在末尾添加几个e.subs(r)才能获取所有原始变量。 - Him

Answer 2

据我所知，目前还没有能够完全实现该功能的函数。我认为这是一个非常困难的问题。请参见Reduce the number of operations on a simple expression以获取更多讨论。

然而，在SymPy中有很多简化函数可供尝试。其中一个你没有提到的是gcd_terms，它可以因式分解出一个符号gcd而不进行扩展。它给出的结果如下：

>>> gcd_terms(expression)
x*y*z*((-x + 1)*(-x - y + 1) + (-x + 1)*(-x - z + 1) + (-y + 1)*(-x - y + 1) + (-y + 1)*(-y - z + 1) + (-z + 1)*(-x - z + 1) + (-z + 1)*(-y - z + 1))

另一个实用的函数是.count_ops，它可以计算表达式中操作的数量。例如：

>>> expression.count_ops()
47
>>> factor(expression).count_ops()
22
>>> e = x * y * z * (6 * (1 - x - y - z) + (x + y) ** 2 + (y + z) ** 2 + (x + z) ** 2)
>>> e.count_ops()
18

请注意，e.count_ops()与您自己计算的不同，因为SymPy会自动将6*(1-x-y-z)分配给6-6*x-6*y-6*z。

其他有用的函数：

• cse：对表达式执行公共子表达式消除。有时您可以简化各个部分，然后将其重新组合。这也有助于避免重复计算。
• horner：对多项式应用霍纳方案。如果多项式是单变量的，则最小化操作次数。
• factor_terms：类似于gcd_terms。实际上我并不完全清楚它们之间的区别。

请注意，默认情况下，simplify会尝试几种简化方法，并返回由count_ops最小化的结果。

- asmeurer

Answer 3

我曾经遇到过类似的问题，在我发现这个解决方案之前，我已经实现了自己的解决方案。我的解决方案似乎在减少操作数量方面做得更好。然而，我的解决方案也是通过对所有变量组合进行暴力样式的集合操作。因此，它的运行时间在变量数量上呈超指数增长。另一方面，我已经设法在不合理但远非实时的时间内运行它，可以处理7个变量的方程。

可能有一些方法来修剪这里的某些搜索分支，但我没有费心去做。欢迎进一步优化。

def collect_best(expr, measure=sympy.count_ops):
    # This method performs sympy.collect over all permutations of the free variables, and returns the best collection
    best = expr
    best_score = measure(expr)
    perms = itertools.permutations(expr.free_symbols)
    permlen = np.math.factorial(len(expr.free_symbols))
    print(permlen)
    for i, perm in enumerate(perms):
        if (permlen > 1000) and not (i%int(permlen/100)):
            print(i)
        collected = sympy.collect(expr, perm)
        if measure(collected) < best_score:
            best_score = measure(collected)
            best = collected
    return best

def product(args):
    arg = next(args)
    try:
        return arg*product(args)
    except:
        return arg

def rcollect_best(expr, measure=sympy.count_ops):
    # This method performs collect_best recursively on the collected terms
    best = collect_best(expr, measure)
    best_score = measure(best)
    if expr == best:
        return best
    if isinstance(best, sympy.Mul):
        return product(map(rcollect_best, best.args))
    if isinstance(best, sympy.Add):
        return sum(map(rcollect_best, best.args))

为了说明性能，this paper（需要付费，抱歉）列出了7个公式，这些公式是7个变量的5次多项式，含有最多29个项和158个操作。在应用rcollect_best和@smichr的iflfactor后，这7个公式中的操作数量如下：

[6, 15, 100, 68, 39, 13, 2]

和

[32, 37, 113, 73, 40, 15, 2]

分别为。其中一个公式中，iflfactor 的操作次数比 rcollect_best 多433%。此外，扩展后的公式中操作次数如下：

[39, 49, 158, 136, 79, 27, 2]

- Him