Minimalpolynom einer quadratischen Matrix

[JETMAN]

Hallo,

ich bin auf der Suche nach einer Routine die mir das Minimalpolynom einer quadratischen Matrix errechnet. Im Packet Sympy fand ich allerdings nur die Bestimmung des Charakteristischen Polynoms (charpoly()).
Gibt es überhaupt eine solche Routine ähnlich der in Maxima "minimalpoly()".

Schöne Grüße
Roland

[JETMAN]

Für alle die es interessiert. Ich habe mal angefangen eine Funktion zu schreiben, die das Minimalpolynom einer Matrix bestimmt. Die Funktion ist noch nicht ganz fertig, funktioniert aber schon für alle Matrizen über Z:

Code: Alles auswählen

def minipoly(matrix, var = 'T'):
    [P, JCells] = matrix.jordan_cells()
    li_len = len(JCells)
    ews = {}

    for i in range(li_len):
        jc = JCells[i].eigenvals()
        gefunden = False

        for k in ews.keys():
            if k == list(jc.keys())[0]:
                gefunden = True
                z = ews.get(k) < jc.get(k)
                if z:
                    ews.update(jc)
                    print(z, k, ews.get(k), jc.get(k))

        if not gefunden:
            ews.update(jc)

    poly = ''
    poly_pur = ''
    tm = ''
    for k in ews.keys():
        if len(poly) > 0:
            poly = poly + ' * '
            poly_pur = poly_pur + ' * '
        if k > 0:
            tm = '-' + str(abs(k))
        elif k < 0:
            tm = '+' + str(abs(k))
        else:
            tm = ''
        poly = poly + '(' + str(var) + tm + ')**' + str(ews.get(k))
        poly_pur = poly_pur + '(' + str(var) + tm + '*Id)**' + str(ews.get(k))

    return eval(poly), poly_pur

Die Funktion liest die Jordanblöcke und wertet sie entsprechend für das Minimalpolynom aus.
Z.B.:

Code: Alles auswählen

	A = Matrix([
        	[3, 1, 1,-2],
        	[0, 2, 1, 0],
        	[0, 0, 2, 0],
        	[1, 1, 1, 0]])
        	
    	I_n = Matrix(identity(A.shape[0]))
    	[a, p] = minipoly(A)
    	pol = lambda T, Id: eval(p)
    	pol(A, I_n)

Das Ergebnis ist wie erwartet eine Nullmatrix!

Schöne Grüße
Roland

[__blackjack__]

@JETMAN: Da fehlen ein Importe damit der Code funktioniert. Ich vermute mal `Matrix` soll aus `sympy` kommen? Und `identity()` aus `numpy`?

In der Funktion wird bei `eval()` ein `T` verwendet was nirgends definiert ist. `eval()` ist sowieso eine blöde Idee. Insbesondere wenn man sowieso schon `sympy` hat und sich damit dynamisch einen Ausdruck aufbauen kann. Da `poly` sowieso nicht verwendet wird in Deinem Test, habe ich das jetzt erst einmal einfach rausgeworfen.

Um Gleichheitszeichen in Argumentlisten werden normalerweise keine Leerzeichen gesetzt.

Man kann auf der linken Seite einer Zuweisung tatsächlich eine Liste mit Namen schreiben, das macht aber niemand. Da steht normalerweise ein Tupel und meistens auch ohne irgendwelche Klammern, solange das nicht auf mehr als eine Zeile aufgeteilt werden muss.

Namen werden klein_mit_unterstrichen geschrieben, ausser Konstanten (KOMPLETT_GROSS) und Klassen (MixedCase). Kryptische Abkürzungen sind doof. Das soll man lesen und verstehen können ohne das man Rätselraten muss. `li_len`? Klingt nach einem asiatischen Namen.

Der Name ist aber sowieso unnötig, genau wie das `i`, denn man kann auch direkt über die Elemente von `JCells` iterieren, ohne den Umweg über einen Index.

Der Wert von `z` muss nicht an einen Namen gebunden werden, denn an der einzigen Stelle wo das ausser beim ``if`` noch verwendet wird, hat es *immer* den Wert `True`, sonst wäre der Programmablauf ja nicht in den ``if``-Zweig gegangen.

``list(jc.keys())[0]`` liefert einen beliebigen Schlüssel aus `jc` – ist das gewollt? Was ist mit den anderen Schlüsseln? Das sieht mir extrem nach einem Programmierfehler aus.

Anstelle der ganzen `get()`-Aufrufe würde man einen Index/Schlüsselzugriff machen, denn es ist ja bei jedem der Zugriffe garantiert dass es den Schlüssel gibt, womit `get()` keinen Sinn macht.

Zumindest im innersten ``if`` kann man die umschliessende Schleife abbrechen, denn eine weitere Abarbeitung kann am Zustand/Ergebnis nichts mehr ändern.

Vor der zweiten Schleife in der Funktion macht es keinen Sinn `tm` an einen Wert zu binden, der dann nirgends verwendet wird.

Die Schleife kann auch gleich über `ews.items()` gehen, denn man braucht im Schleifenkörper grundsätzlich Schlüssel *und* Wert.

Das wiederholte ”erweitern” einer Zeichenkette mit ``+`` ist ineffizient. Man sammelt deshalb die Teilzeichenketten in einer Liste und setzt sie dann am Ende mit der `str.join()`-Methode zusammen. Hier kann man als Trennzeichen dann die Multiplikation wählen. Das macht den Code auch einfacher, weil man am Anfang nicht mehr den Sonderfall hat bei dem Kein ' * ' hinzugefügt werden muss.

`var` ist bereits eine Zeichenkette. Da macht es keinen Sinn noch einmal `str()` mit aufzurufen. Überhaupt ist das Zusammenstückeln von Zeichenketten und Werten mit ``+`` und `str()` nicht wirklich Python sondern eher BASIC. In Python gibt es dafür Zeichenkettenformatierung.

Code: Alles auswählen

def minipoly(matrix, var='T'):
    _, jordan_cells = matrix.jordan_cells()
    ews = dict()
    for cell in jordan_cells:
        eigenvals = cell.eigenvals()
        gefunden = False
        m = list(eigenvals.keys())[0]
        for k in ews:
            if k == m:
                gefunden = True
                if ews[k] < eigenvals[k]:
                    ews.update(eigenvals)
                    print(True, k, ews[k], eigenvals[k])
                    break

        if not gefunden:
            ews.update(eigenvals)

    poly_pur = list()
    for k, value in ews.items():
        if k == 0:
            tm = ''
        else:
            tm = '{}{}'.format('-' if k > 0 else '+', abs(k))
        poly_pur.append('({}{}*Id)**{}'.format(var, tm, value))

    return ' * '.join(poly_pur)

Aber wie gesagt möchte man eher nicht evil `eval()` verwenden, insbesondere wenn man in der letzten Schleife stattdessen eine Sympy-Objekt erstellen kann (`add` und `sub` aus dem `operator`-Modul):

Code: Alles auswählen

    T, Id = symbols('T Id')
    terms = list()
    for k, value in ews.items():
        terms.append((add if k < 0 else sub)(T, k * Id)**value)
    return Mul(*terms)

[JETMAN]

@__blackjack__
Mein Code waren die ersten Überlegungen zu dem Thema "das Minimalpolynom einer Matrix" aus den Jordanblöcken in Python zu realisieren.
Ich möchte mich herzlichst für die excellenten Ausführungen bedanken. Da ich "Anfänger" in Sachen "Python" bin werde ich mir deine Hinweise genau anschauen um hieraus zu lernen.

Ein Frage bleibt noch, wenn ich das erzeugte Polynom auf eine Matrix anwenden möchte, brauche ich doch die Funktion eval()?!?

Zum Beispiel (nehme ich Deine Funktion "minipoly"):

Code: Alles auswählen

A = Matrix([
        [3, 1, 1,-2],
        [0, 2, 1, 0],
        [0, 0, 2, 0],
        [1, 1, 1, 0]])
        
idn = Matrix([
        [1,0,0,0],
        [0,1,0,0],
        [0,0,1,0],
        [0,0,0,1]])

mp = minipoly(A)

Ausgabe.
(T-1*Idn)**1 * (T-2*Idn)**2

Code: Alles auswählen

mpol = lambda T, Idn: eval(mp)
mpol(A, idn)

Ausgabe:
⎡0 0 0 0⎤
⎢0 0 0 0⎥
⎢0 0 0 0⎥
⎣0 0 0 0⎦

Wenn ich nun eval() weglasse, bekomme ich als Ausgabe das Polynom und nicht das Ergebnis.

Schöne Grüße

[__blackjack__]

@JETMAN: Jain, wenn Du da tatsächlich noch eine Zeichenkette zusammenbastelst, brauchst Du `eval()`. Aber ich würde, wenn ich sowieso schon `sympy` an Bord habe, wie gesagt, gar keine Zeichenkette mit Quelltext zusammenbasteln, sondern ein Sympy-Objekt das den Ausdruck repräsentiert. Da kommt dann auch nicht '(T-1*Id)**1 * (T-2*Id)**2' heraus, sondern '(-2*Id + T)**2*(-Id + T)', wenn man das für die Ausgabe in eine Zeichenkette umwandeln lässt. Ist natürlich das selbe, aber eben schon etwas vereinfacht. Ausrechnen kann man das dann beispielsweise über die `subs()`-Methode in dem man konkrete Werte für die Platzhalter angibt.