Matrix-Balancing-Algorithmus

Ich habe eine Steuerungssystem-Toolbox von Grund auf neu und rein in Python3 geschrieben (schamloser Plug :) harold. Aus meiner bisherigen Forschung habe ich mich immer care.maus technischen / irrelevanten Gründen über den Riccati-Löser beschwert .

Daher habe ich meine eigenen Routinen geschrieben. Eine Sache, die ich nicht umgehen kann, ist, einen Hochleistungs-Ausgleichsalgorithmus zu erhalten, der mindestens so gut ist wie balance.m. Bevor Sie es erwähnen, wird die xGEBALFamilie in Scipy angezeigt und Sie können von Scipy aus wie folgt aufrufen: Angenommen, Sie haben ein Float-2D-Array A:

import scipy as sp
gebal = sp.linalg.get_lapack_funcs(('gebal'),(A,)) # this picks up DGEBAL
Ab, lo, hi, scaling , info = gebal(A, scale=1 , permute=1 , overwrite_a=0 )

Nun, wenn ich die folgende Testmatrix verwende

array([[ 6.      ,  0.      ,  0.      ,  0.      ,  0.000002],
       [ 0.      ,  8.      ,  0.      ,  0.      ,  0.      ],
       [ 2.      ,  2.      ,  6.      ,  0.      ,  0.      ],
       [ 2.      ,  2.      ,  0.      ,  8.      ,  0.      ],
       [ 0.      ,  0.      ,  0.000002,  0.      ,  2.      ]])

Ich bekomme

array([[ 8.      ,  0.      ,  0.      ,  2.      ,  2.      ],
       [ 0.      ,  2.      ,  0.000002,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  6.      ,  2.      ,  2.      ],
       [ 0.      ,  0.000002,  0.      ,  6.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      ,  0.      ,  8.      ]])

Wenn ich dies jedoch weitergebe balance.m, bekomme ich

>> balance(A)

ans =

    8.0000         0         0    0.0625    2.0000
         0    2.0000    0.0001         0         0
         0         0    6.0000    0.0002    0.0078
         0    0.0003         0    6.0000         0
         0         0         0         0    8.0000

Wenn Sie die Permutationsmuster überprüfen, sind sie gleich, die Skalierung ist jedoch deaktiviert. gebalgibt Einheitsskalierungen an, während matlab die folgenden Potenzen von 2 ergibt : [-5,0,8,0,2].

Anscheinend verwenden diese nicht die gleichen Maschinen. Ich habe verschiedene Optionen ausprobiert, wie Lemonnier, Van Dooren zweiseitige Skalierung, original Parlett-Reinsch und auch einige andere weniger bekannte Methoden in der Literatur, wie die dichte Version von SPBALANCE.

Ein Punkt, den ich vielleicht betonen möchte, ist, dass ich mir Benners Arbeit bewusst bin; insbesondere das Symplectic Balancing of Hamiltonian Matrices speziell für diesen Zweck. Beachten Sie jedoch, dass diese Art der Behandlung innerhalb gcare.m(generalisierter Riccati-Löser) erfolgt und der Ausgleich direkt über erfolgt balance.m. Daher würde ich mich freuen, wenn mich jemand auf die tatsächliche Implementierung hinweisen kann.

Offenlegung: Ich versuche wirklich nicht, Mathworks-Code zurückzuentwickeln: Ich möchte aus verschiedenen Gründen, einschließlich der Motivation dieser Frage, davon wegkommen, das heißt, ich weiß nicht, was es tut, was mich viel gekostet hat der Zeit zurück in den Tag. Meine Absicht ist es, einen zufriedenstellenden Ausgleichsalgorithmus zu erhalten, der es mir ermöglicht, CAREX-Beispiele so zu übergeben, dass ich Newton-Iterationsmethoden zusätzlich zum regulären Löser implementieren kann.

Ich habe eine ganze Weile gebraucht, um das herauszufinden, und wie immer wird es offensichtlich, nachdem Sie den Täter gefunden haben.

Nach Überprüfung der in David S. Watkins gemeldeten problematischen Fälle . Ein Fall, in dem das Balancieren schädlich ist. Elektron. Trans. Numer. Anal, 23: 1–4, 2006 und auch die Diskussion hier (beide zitiert in arXiv: 1401.5766v1 ) stellen sich heraus, dass matlab den Ausgleich verwendet, indem zuerst die diagonalen Elemente getrennt werden.

Mein erster Gedanke war, dass GEBAL dies gemäß der klassischen eingeschränkten Dokumentation zu LAPACK-Funktionen automatisch durchführte. Ich denke jedoch, was die Autoren unter Ignorieren der diagonalen Elemente verstehen, bedeutet nicht, sie aus den Zeilen- / Spaltensummen zu entfernen.

Wenn ich die Diagonale manuell aus dem Array entferne, stimmen beide Ergebnisse überein

import scipy as sp
gebal = sp.linalg.get_lapack_funcs(('gebal'),(A,)) # this picks up DGEBAL
Ab, lo, hi, scaling , info = gebal(A - np.diag(np.diag(A)), scale=1 , permute=1 , overwrite_a=0 )  

ergibt das gleiche Ergebnis wie balance.m(natürlich ohne die diagonalen Einträge).

Wenn ein Fortran-versierter Benutzer dies durch Überprüfen von dgebal.f bestätigen kann , wäre ich dankbar.

EDIT: Das obige Ergebnis bedeutet nicht, dass dies der einzige Unterschied ist. Ich habe auch verschiedene Matrizen konstruiert, bei denen GEBAL und balance.m auch nach dem Trennen der Diagonalen unterschiedliche Ergebnisse liefern.

Ich bin ziemlich neugierig, was der Unterschied sein könnte, aber es scheint, dass es keine Möglichkeit gibt, dies zu wissen, da es sich um einen in Matlab integrierten und daher geschlossenen Code handelt.

EDIT2 : Es stellte sich heraus, dass matlab eine ältere Version von LAPACK (wahrscheinlich vor 3.5.0) verwendete und bis 2016b auf die neuere Version aktualisiert zu sein scheint. Jetzt stimmen die Ergebnisse überein, soweit ich sie testen kann. Ich denke, das löst das Problem. Ich hätte es mit älteren LAPACK-Versionen testen sollen.