Vorkonditionierung — warum GMRES hier n Iterationen braucht

Der LeSim-Solver arbeitet matrixfrei. Er kennt keine Systemmatrix, sondern nur einen Operator: BmAx(x) = b − A·x. Genau ein Aufruf davon lässt in LeSimStepCtrl.KPA_solve die komplette Blockphysik aller Blöcke durchlaufen:

const Ax = function(KnotenSpgVektor) {
    for (...) { Knoten[k].Spannung = KnotenSpgVektor[k]; Knoten[k].StromSum = 0; }
    for (let block of opt.Blocks) block.KPA_solve(opt);
    return opt.Knoten.map(kn => kn.StromSum);
};
const res = Solver.run(Ax, x0, SolverDebugInfo);

Die Zahl der Krylov-Iterationen ist damit nicht bloß eine Kennzahl des Lösers, sondern unmittelbar der Preis der Simulation. Bei einem Reihenstapel aus N Zellen kostet ein Aufruf O(N) — braucht der Löser zusätzlich O(N) Iterationen, landet man bei O(N²).

Diese Seite wertet dazu den vorhandenen Bestand an Debugmatrizen aus (./json/, dieselben Matrizen wie im Solver-Benchmark). Alle Zahlen werden beim Öffnen live im Browser gerechnet, mit dem echten sgmres.js dieser Seite — nichts ist fest eingetragen.

Matrizen

—

ausgewertet

Median it/n

—

Iterationen je Unbekannte

Korrelation

—

Iterationen gegen n

Singulär

—

Rangdefekt ≥ 1

Ein Krylov-Verfahren soll nach wenigen Iterationen abbrechen — das ist sein ganzer Sinn. Auf diesem Bestand passiert das nirgends. GMRES terminiert nach spätestens n Schritten exakt; es läuft hier also praktisch bis zum Ende durch und verhält sich wie ein direktes Verfahren, nur mit dem Aufwand eines iterativen.

n	#	BmAx	it/n

BmAx-Aufrufe über der Systemgröße. Die gestrichelte Linie ist „Aufrufe = n". Die Punkte liegen darauf, nicht darunter.

Matrix

Besetzungsmuster. Jeder Punkt ist ein Eintrag über 0,1 % der Zeilennorm.

Die Knotenmatrix eines Reihenstapels ist eine Kette — jeder Knoten hängt nur an seinen zwei Nachbarn, geschlossen über die äußere Beschaltung zu einem Ring. Daraus folgt eine harte Untergrenze, die nichts mit Rundung zu tun hat:

$$ A^k \text{ hat Bandbreite } k $$

Der Krylov-Raum span{b, Ab, A²b, …} transportiert Information also genau einen Knoten je Iteration. Um ein Kettenende mit dem anderen zu verbinden, braucht es mindestens so viele Matrix-Vektor-Produkte, wie der Graph Abstand hat. Das gilt unabhängig von Toleranz, Startwert und Skalierung.

—

Die Matrizen sind zudem fast alle singulär — Knotenpotentialverfahren ohne Bezugsknoten, jede schwimmende Komponente lässt sich um eine Konstante verschieben. Das ist aber nicht die Ursache der hohen Iterationszahl: einen Knoten je Komponente festzuhalten ändert sie praktisch nicht.

Bringt Umsortieren etwas? Die beiden Einträge rechts oben und links unten sind keine Störung, sondern die Kette selbst: die Knotennummerierung setzt die äußeren Klemmen an den Anfang (0, 1) und die versteckten Kettenknoten dahinter (4 … n−1). Physikalisch läuft die Reihe 0 → 4 → 5 → … → n−1 → 1, und in dieser Reihenfolge nummeriert wird daraus eine saubere Diagonale.

rechnet …

Für den Löser ändert das aber nichts. Eine symmetrische Umnummerierung wirkt auf den Krylov-Raum genauso wie auf alles andere:

$$ (P^{T}AP)^{k}\,(P^{T}b) \;=\; P^{T}\,(A^{k}b) $$

Der permutierte Krylov-Raum ist der permutierte alte, die Iterierten sind dieselben Vektoren in anderer Reihenfolge. Umnummerieren allein ist für ein Krylov-Verfahren eine Nulloperation. Nützlich wird es erst für den Vorkonditionierer: ein Band-LU muss wissen, wo das Band liegt — sonst müsste er statt einer Bandbreite 1 eine von n−2 faktorisieren.

Den Ausreißer einfach wegzulassen hilft ebenfalls nicht — er ist die physikalische Verbindung der unteren Primärklemme zum letzten Kettenknoten. Ohne ihn rechnet man einen anderen Schaltplan.

Links-Vorkonditionierung braucht keine Änderung am Löser. Weil M⁻¹(b − A·x) = M⁻¹b − (M⁻¹A)·x genau die erwartete Form hat, ist sie eine Hülle um BmAx:

const pre   = makeBandPrecond(A);          // bandPrecond.js
const BmAx2 = wrapBmAx(BmAx, pre);
const res   = sgmres(BmAx2, x0, eps);      // unveraendert

Jacobi — nur die Diagonale. Billig, aber die schlechte Kondition sitzt hier nicht auf der Diagonale.
RCM + Band-LU — Knoten mit Reverse Cuthill-McKee so umsortieren, dass A möglichst schmalbandig wird, dann das Band exakt lösen. Für eine Kette wird die Bandbreite 1, der Vorkonditionierer damit exakt.

rechnet …

Matrix	n	Grad	bw	Modus	roh	Jacobi	Band	Faktor	res roh	res Band

Gezählt werden BmAx-Aufrufe, nicht Iterationen — das ist die Größe, die in der Simulation Rechenzeit kostet. Das Residuum wird immer mit dem unvorkonditionierten Operator nachgerechnet: bei Links-Vorkonditionierung sieht der Löser nur noch ein verzerrtes Residuum, sein eigenes „ok" sagt nichts über die Lösungsgüte.

Ein Vorkonditionierer lohnt nur, wenn er über viele Löserläufe wiederverwendet wird: ihn aufzubauen kostet die Matrix, und die kostet über OR.getABD n+1 BmAx-Aufrufe. Die Topologie ändert sich zwischen Zeitschritten nicht — die Zahlenwerte aber schon. Der Diagonalanteil enthält c/dt, ein Faktor 2 entspricht also einer Halbierung der Schrittweite.

Matrix	n	roh	dt gleich	×1,1	×1,5	×2	×4	×8

BmAx-Aufrufe, wenn M aus A gebaut und auf A mit skalierter Diagonale angewendet wird.

rechnet …

Dazu die Amortisation: bei N=60 kostet ein Neuaufbau 64 Aufrufe, während der Vorkonditionierer je Lauf etwa 54 spart. Jeden Schritt neu aufbauen ist damit schlechter als gar nichts zu tun; alle 10 Schritte ergibt Faktor ≈ 3,9, alle 50 Schritte ≈ 5,7 — aber dann ist M so veraltet, dass es kaum noch wirkt. Beides zusammen landet realistisch bei Faktor 2 bis 3.

Was daraus folgt

Der teure Teil ist nicht der Löser, sondern das Produkt aus Iterationszahl und Aufrufkosten. Es gibt drei Wege, es zu senken:

Vorkonditionieren. Allgemein und für jede Kettentopologie wirksam — Reihenstapel, MMC, Leitungsmodelle. Der Haken ist die Beschaffung von A: über getABD zu teuer für häufigen Neuaufbau. Sinnvoll erst, wenn die Blöcke ihren Leitwertbeitrag selbst stempeln können (eine KPA_stamp-Schnittstelle neben KPA_solve), denn dann kostet der Aufbau O(N) statt n+1 vollen Physikdurchläufen.
Die inneren Knoten gar nicht erst anlegen. Ein Block mit Reihenstruktur kann seine internen Knoten selbst eliminieren und dem globalen Löser nur seine Klemmen zeigen (Schur-Komplement). Bei eingefrorenem Schaltzustand ist jede Zelle exakt affin in ihren Klemmenspannungen — nachgewiesen auf Maschinengenauigkeit —, das Zweitor folgt also aus wenigen Sondierungen und die Iterationszahl wird von n unabhängig.
Eine geschlossene Lösung im Block. Was MultiCellDAB tut. Faktisch ein exakter, hartverdrahteter Vorkonditionierer für genau diese Struktur: er tauscht n Krylov-Iterationen gegen eine Formel. Schnell, aber eine zweite Herleitung derselben Physik — mit dem bekannten Risiko, dass beide auseinanderdriften.

Für den Löser selbst wäre der interessanteste Ausbau eine Strukturerkennung: maximalen Knotengrad bei relativer Schwelle bestimmen und bei ≤ 2 zyklisch-tridiagonal vorkonditionieren. Das hilft jedem Reihenstapel, ohne dass ein Block dafür eigene Physik braucht.

Rechenkern dieser Seite: ./bandPrecond.js (Adjazenz mit relativer Schwelle, Reverse Cuthill-McKee, Band-LU mit Diagonalregularisierung) und das unveränderte ./sgmres.js. Zurück zum Solver-Benchmark.