Probleme beim Ersetzen einer Matrix in einem Polynom

Jan 12 2021

Beispiel: let

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

gibt den Fehler:

TypeError: <class'sympy.matrices.immutable.ImmutableDenseMatrix '> und <class' sympy.core.numbers.One '> können nicht hinzugefügt werden

Dies ist sehr plausibel, da die beiden ersten Terme zu Matrizen werden, der letzte Term (der konstante Term) jedoch keine Matrix, sondern ein Skalar ist. Um diese Situation zu beheben, habe ich das Polynom in geändert

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

Der gleiche Fehler bleibt bestehen. Muss ich den Ausdruck von Hand eingeben und x durch M ersetzen? In diesem Beispiel hat das Polynom nur drei Terme, aber in Wirklichkeit begegne ich Dutzenden von Termen (multivariate Polynome mit).

Antworten

1 wsdookadr Jan 12 2021 at 12:40

Ich denke, die Frage dreht sich hauptsächlich um das Konzept des Matrixpolynoms :

(wobei P ein Polynom und A eine Matrix ist)

Ich denke, dies bedeutet, dass der freie Term eine Zahl ist und nicht mit dem Rest addiert werden kann, der eine Matrix ist. Tatsächlich ist die Additionsoperation zwischen diesen beiden Typen undefiniert.

TypeError: <class'sympy.matrices.immutable.ImmutableDenseMatrix '> und <class' sympy.core.numbers.One '> können nicht hinzugefügt werden

Dies kann jedoch umgangen werden, indem eine Funktion definiert wird, die das Matrixpolynom für eine bestimmte Matrix auswertet. Der Unterschied besteht darin, dass wir die Matrixexponentiation verwenden , sodass wir den freien Term des Matrixpolynoms korrekt berechnen, a_0 * Iwobei I=A^0sich die Identitätsmatrix der erforderlichen Form befindet:

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)

Ausgabe:

In diesem Beispiel hat das Polynom nur drei Terme, aber in Wirklichkeit begegne ich Dutzenden von Termen (multivariate Polynome mit).

Die eval_poly_matrixobige Funktion kann erweitert werden, um für multivariate Polynome zu arbeiten, indem die .monoms()Methode zum Extrahieren von Monomen mit Koeffizienten ungleich Null verwendet wird , wie folgt:

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

Hinweis: Einige Sanity Checks, Edge Cases Handling und Optimierungen sind möglich:

  1. Die Anzahl der im Polynom vorhandenen Variablen bezieht sich auf die Anzahl der als Parameter übergebenen Matrizen (die erstere sollte niemals größer als die letztere sein, und wenn sie niedriger ist als eine Logik, die vorhanden sein muss, um dies zu handhaben, habe ich nur den Fall behandelt wenn die beiden gleich sind)
  2. Alle Matrizen müssen gemäß der Definition des Matrixpolynoms quadratisch sein
  3. Eine Diskussion über eine multivariate Version der Horner-Regel findet sich in den Kommentaren dieser Frage . Dies kann nützlich sein, um die Anzahl der Matrixmultiplikationen zu minimieren.
  4. Behandeln Sie die Tatsache , dass in einer Matrix Polynom x*yunterscheidet sich von y*xda Matrixmultiplikation ist nichtkommutative . Anscheinend unterstützen Poly-Funktionen in Sympy keine nicht kommutativen Variablen, aber Sie können Symbole mit definieren, commutative=Falseund es scheint einen Weg nach vorne zu geben

Über den 4. Punkt oben gibt es Unterstützung für Matrix-Ausdrücke in SymPy, und das kann hier helfen:

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)

Ausgabe:

Weitere Entwicklungen zu dieser Funktion finden Sie unter # 18555