Problemas ao substituir uma matriz em um polinômio

Question

Exemplo: deixe

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

dá o erro:

TypeError: não é possível adicionar <class'sympy.matrices.immutable.ImmutableDenseMatrix '> e <class' sympy.core.numbers.One '>

o que é muito plausível, uma vez que os dois primeiros termos se tornam matrizes, mas o último termo (o termo constante) não é uma matriz, mas um escalar. Para remediar esta situação, mudei o polinômio para

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

o mesmo erro persiste. Devo digitar a expressão manualmente, substituindo x por M? Neste exemplo, o polinômio tem apenas três termos, mas na realidade encontro (polinômios multivariados com) dezenas de termos.

Answer 1

Então, acho que a questão gira principalmente em torno do conceito de polinômio de matriz :

(onde P é um polinômio e A é uma matriz)

Acho que isso quer dizer que o termo livre é um número, e não pode ser adicionado com o resto que é uma matriz, efetivamente a operação de adição é indefinida entre esses dois tipos.

TypeError: não é possível adicionar <class'sympy.matrices.immutable.ImmutableDenseMatrix '> e <class' sympy.core.numbers.One '>

No entanto, isso pode ser contornado definindo uma função que avalia o polinômio da matriz para uma matriz específica. A diferença aqui é que estamos usando a exponenciação da matriz , então calculamos corretamente o termo livre do polinômio da matriz, a_0 * Ionde I=A^0é a matriz de identidade da forma necessária:

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)

Resultado:

Neste exemplo, o polinômio tem apenas três termos, mas na realidade encontro (polinômios multivariados com) dezenas de termos.

A função eval_poly_matrixacima pode ser estendida para funcionar para polinômios multivariados usando o .monoms()método para extrair monômios com coeficientes diferentes de zero , assim:

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))

Nota: Algumas verificações de integridade, manipulação de casos extremos e otimizações são possíveis:

1. O número de variáveis ​​presentes no polinômio está relacionado ao número de matrizes passadas como parâmetros (as primeiras nunca devem ser maiores que as últimas, e se for menor do que alguma lógica precisa estar presente para lidar com isso, eu apenas tratei do caso quando os dois são iguais)
2. Todas as matrizes precisam ser quadradas de acordo com a definição do polinômio da matriz
3. Uma discussão sobre uma versão multivariada das características da regra de Horner nos comentários desta questão . Isso pode ser útil para minimizar o número de multiplicações de matrizes.
4. Lide com o fato de que em um polinômio Matrix x*yé diferente y*xporque a multiplicação da matriz não é comutativa . Aparentemente, as funções poli em sympy não suportam variáveis não comutativas , mas você pode definir símbolos com commutative=Falsee parece haver um caminho a seguir aí

Sobre o 4º ponto acima, há suporte para expressões Matrix no SymPy, e isso pode ajudar aqui:

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)

Resultado:

Para mais desenvolvimentos neste recurso, siga # 18555