Problèmes lors de la substitution d'une matrice dans un polynôme

Question

Exemple: let

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

donne l'erreur:

TypeError: impossible d'ajouter <class'sympy.matrices.immutable.ImmutableDenseMatrix '> et <class' sympy.core.numbers.One '>

ce qui est très plausible puisque les deux premiers termes deviennent des matrices mais le dernier terme (le terme constant) n'est pas une matrice mais un scalaire. Pour remédier à cette situation, j'ai changé le polynôme en

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

la même erreur persiste. Suis-je obligé de taper l'expression à la main, en remplaçant x par M? Dans cet exemple, le polynôme n'a que trois termes mais en réalité je rencontre (polynômes multivariés avec) des dizaines de termes.

Answer 1

Je pense donc que la question tourne principalement autour du concept de polynôme Matrix :

(où P est un polynôme et A est une matrice)

Je pense que cela veut dire que le terme libre est un nombre, et il ne peut pas être ajouté avec le reste qui est une matrice, effectivement l'opération d'addition est indéfinie entre ces deux types.

TypeError: impossible d'ajouter <class'sympy.matrices.immutable.ImmutableDenseMatrix '> et <class' sympy.core.numbers.One '>

Cependant, cela peut être contourné en définissant une fonction qui évalue le polynôme de la matrice pour une matrice spécifique. La différence ici est que nous utilisons l' exponentiation matricielle , nous calculons donc correctement le terme libre du polynôme matriciel a_0 * Ioù I=A^0est la matrice d'identité de la forme requise:

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)

Production:

Dans cet exemple, le polynôme n'a que trois termes mais en réalité je rencontre (polynômes multivariés avec) des dizaines de termes.

La fonction eval_poly_matrixci-dessus peut être étendue pour fonctionner pour les polynômes multivariés en utilisant la .monoms()méthode pour extraire les monômes avec des coefficients différents de zéro , comme ceci:

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

Remarque: certaines vérifications de cohérence, la gestion des cas marginaux et des optimisations sont possibles:

1. Le nombre de variables présentes dans le polynôme est lié au nombre de matrices passées en paramètres (la première ne doit jamais être supérieure à la seconde, et si elle est inférieure à une certaine logique doit être présente pour gérer cela, je n'ai traité que le cas quand les deux sont égaux)
2. Toutes les matrices doivent être carrées selon la définition du polynôme matriciel
3. Une discussion sur une version multivariée des caractéristiques de la règle de Horner dans les commentaires de cette question . Cela peut être utile pour minimiser le nombre de multiplications matricielles.
4. Gérez le fait que dans un polynôme Matrix x*yest différent de y*xparce que la multiplication matricielle est non commutative . Apparemment, les fonctions poly dans sympy ne prennent pas en charge les variables non commutatives , mais vous pouvez définir des symboles avec commutative=Falseet il semble y avoir un moyen d'y aller

À propos du 4e point ci-dessus, il existe un support pour les expressions Matrix dans SymPy, et cela peut aider ici:

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)

Production:

Pour plus de développements sur cette fonctionnalité, suivez # 18555