在多项式中替换矩阵时出现的问题

Question

例子：让

M = Matrix([[1,2],[3,4]]) # and 
p = Poly(x**3 + x + 1)    # then
p.subs(x,M).expand()

给出错误：

类型错误：无法添加和

这是非常合理的，因为前两项变成矩阵但最后一项（常数项）不是矩阵而是标量。为了解决这种情况，我将多项式更改为

p = Poly(x**3 + x + x**0)    # then

同样的错误仍然存在。我是否必须手动输入表达式，将 x 替换为 M？在这个例子中，多项式只有三个项，但实际上我遇到了（多元多项式）几十个项。

Answer 1

所以我觉得问题主要围绕Matrix polynomial:

的概念

（其中P是多项式，A是矩阵）

我认为这是说自由项是一个数字，它不能与矩阵的其余部分相加，实际上这两种类型之间的加法运算是未定义的。

TypeError: cannot add <class'sympy.matrices.immutable.ImmutableDenseMatrix'> and <class 'sympy.core.numbers.One'>

但是，这可以通过定义一个计算特定矩阵的矩阵多项式的函数来规避。这里的区别在于我们使用的是所需形状的 matrix exponentiation, so we correctly compute the free term of the matrix polynomial a_0 * I where I=A^0 is the identity matrix：

from sympy import *
x = symbols('x')
M = Matrix([[1,2],[3,4]])
p = Poly(x**3 + x + 1)

def eval_poly_matrix(P,A):
    res = zeros(*A.shape)
    for t in enumerate(P.all_coeffs()[::-1]):
        i, a_i = t
        res += a_i * (A**i)
    return res

eval_poly_matrix(p,M)

输出：

In this example the polynomial has only three terms but in reality I encounter (multivariate polynomials with) dozens of terms.

上面的函数 eval_poly_matrix 可以通过使用 .monoms() method 扩展为适用于多元多项式，如下所示：

from sympy import *
x,y = symbols('x y')
M = Matrix([[1,2],[3,4]])
p = poly( x**3 * y + x * y**2 + y )

def eval_poly_matrix(P,*M):
    res = zeros(*M[0].shape)
    for m in P.monoms():
        term = eye(*M[0].shape)
        for j in enumerate(m):
            i,e = j
            term *= M[i]**e
        res += term
    return res

eval_poly_matrix(p,M,eye(M.rows))

注意：一些健全性检查、边缘情况处理和优化是可能的：

多项式中存在的变量数量与作为参数传递的矩阵数量有关（前者永远不应大于后者，如果它低于需要存在的某些逻辑来处理它，我'我们只处理了两者相等的情况）
根据 matrix polynomial
关于 Horner's rule features in the comments of this question 的多元版本的讨论。这可能有助于减少矩阵乘法的次数。
处理矩阵多项式 x*y 与 y*x 不同的事实，因为矩阵乘法 is non-commutative . Apparently poly functions in sympy do not support 变量，但您可以使用 commutative=False 定义符号，似乎成为前进的道路

关于上面的第 4 点，SymPy 中支持 Matrix expressions，这对这里有帮助：

from sympy import *
from sympy.matrices import MatrixSymbol

A = Matrix([[1,2],[3,4]])
B = Matrix([[2,3],[3,4]])

X = MatrixSymbol('X',2,2)
Y = MatrixSymbol('Y',2,2)
I = eye(X.rows)

p = X**2 * Y + Y * X ** 2 + X ** 3 - I
display(p)

p = p.subs({X: A, Y: B}).doit()
display(p)

输出：

有关此功能的更多开发，请关注 #18555

在多项式中替换矩阵时出现的问题

Problems when substituting a matrix in a polynomial

sympy