Lineare Algebra

Durchmessersatz für Polytope und Polyeder: Jeder Durchmesser eines durch \(m\) Restriktionen gegebenen \(n\)-dimensionalen Polytops bzw. Polyeders mit \(m,n\in{}^{\omega}\mathbb{N}_{\ge 2}\) ist maximal \(2(m + n – 3)\).

Beweis: Maximal \(\acute{m}\) Hyperebenen lassen sich zu einem unvollständigen Zyklus der Dimension \(2\) zusam-
menstellen und es gibt höchstens \(n-2\) Ausweichmöglichkeiten zur Seite in den restlichen Dimensionen. Das Überwinden jeder minimalen Strecke benötigt maximal zwei Kanten und liefert den Faktor \(2.\square\)

Satz zum Strassen-Algorithmus: Für hinreichend großes \(n := 2^{\ell}\), \(\ell \in {}^{\nu}\mathbb{N}^*\), \(\text{ß} := {}_27\) und \(A \in {}^{\nu}\mathbb{C}^{n \times n}\) beweist die GR die Verkürzung der Laufzeit \(T(n) = \mathcal{O}(n^{\text{ß}})\)¹ um ca. \(\tilde{3}\) zur Berechnung von \[ AA^H = \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix} \begin{pmatrix} A_{11}^H & A_{21}^H \\ A_{12}^H & A_{22}^H \end{pmatrix} = \begin{pmatrix} A_{11}A_{11}^H+A_{12}A_{12}^H & A_{11}A_{21}^H+A_{12}A_{22}^H \\ A_{21}A_{11}^H+A_{22}A_{12}^H & A_{21}A_{21}^H+A_{22}A_{22}^H \end{pmatrix}.\;\square \] Satz zur schnellen Matrixmultiplikation: Mit \(A\) wie \(B = (b_{ij}) \in {}^{\nu}\mathbb{C}^{n \times n}\) und \(n\) wie zuvor lässt sich \(AB\) für \[ s := \min\{2^k > \max_{i,j} (|a_{ij}|^2, |b_{ij}|^2)\hat{n} : k \in {}^{\nu}\mathbb{Z}\} \] aus \[ A_{11}(s B_{11} + B_{12}),\; A_{21}(s B_{11} + B_{12}),\; A_{12}(s B_{21} + B_{22}),\; A_{22}(s B_{21} + B_{22}) \] durch Rechnen modulo \(s\) in der Laufzeit \(\mathcal{O}(p^2)\) mit \(p = \ell n\) bestimmen (s. Matrixmultiplikation).\(\square\)

Revised INLIF und Farkas-Zertifikate

Das vollständige INLIF-Tableau macht alle aktuellen Koeffizienten explizit sichtbar. Dies ist für Herleitung, Nachvollziehbarkeit und Tableau-Spur vorteilhaft. Für große oder dünn besetzte Systeme ist es jedoch günstiger, nicht das gesamte transformierte Tableau zu speichern, sondern nur die aktuelle Basisstruktur sowie die für den nächsten Schritt benötigten Zeilen und Spalten zu berechnen. Diese Darstellung wird im Folgenden als revised INLIF bezeichnet.

Ausgangspunkt sei das System \(Ax\le b\). Durch Schlupfvariablen \(s\ge0\) entsteht \(Ax+s=b\). Mit \(H=(A\; I)\) und \(q=(x,s)^\top\) lautet das erweiterte Gleichungssystem \(Hq=b\). Die ursprünglichen Variablen \(x_j\) sind im freien INLIF zunächst freie Variablen; die Schlupfvariablen \(s_i\) sind Kegelvariablen. Zu jedem Zeitpunkt stehen \(m\) Variablen auf der Zeilenseite. Die zugehörigen Spalten von \(H\) bilden die aktuelle Basismatrix \(\mathcal B\). Die übrigen Spalten bilden die aktuelle Spaltenmatrix \(\mathcal N\). Damit gilt \(\mathcal B q_R+\mathcal N q_C=b\) und folglich \(q_R=\mathcal B^{-1}b-\mathcal B^{-1}\mathcal N q_C\). Setzt man \(\beta:=\mathcal B^{-1}b\) und \(G:=\mathcal B^{-1}\mathcal N\), so lautet das vollständige Tableau \(q_R=\beta-Gq_C\).

Das revised INLIF bildet die Matrix \(G\) nicht vollständig. Stattdessen wird nur \(\mathcal B\) bzw. eine Faktorisierung oder Produktformdarstellung von \(\mathcal B\) gespeichert. Die aktuelle rechte Seite wird aus \(\mathcal B\beta=b\) gewonnen. Eine aktive Prüfzeile ist verletzt, wenn die zugehörige Zeilenvariable eine Kegelvariable ist und \(\beta_k<0\) gilt. Zeilenvariablen, die aus ursprünglich freien Variablen hervorgegangen sind, sind dagegen reine Rückrechnungszeilen und werden nicht auf Nichtnegativität geprüft.

Für die Pivotwahl in einer verletzten Zeile \(k\) wird nicht die gesamte Tableauzeile berechnet. Stattdessen löst man das transponierte Basissystem \(\mathcal B^\top \lambda_k=e_k\). Dann gilt \(\lambda_k^\top=e_k^\top \mathcal B^{-1}\). Für eine aktuelle Spaltenvariable mit ursprünglicher Spalte \(h_j\) in \(H\) ergibt sich der zugehörige Tableaukoeffizient durch \(G_{kj}=e_k^\top\mathcal B^{-1}h_j=\lambda_k^\top h_j\). Dies ist die zentrale revised-INLIF-Formel:

Eine freie Spaltenvariable ist als Pivot zulässig, wenn \(G_{kj}\ne0\) gilt. Eine Kegelspaltenvariable ist nur dann als Pivot zulässig, wenn \(G_{kj}<0\) gilt. Damit bleibt die Pivotlogik des vollständigen INLIF-Tableaus unverändert; nur die Berechnung der Koeffizienten erfolgt implizit.

Ist eine Pivotspalte \(p\) gewählt, so wird die transformierte Spalte durch \(\mathcal B u=h_p\) berechnet. Dann gilt \(u=\mathcal B^{-1}h_p\) und insbesondere \(u_k=G_{kp}\). Die neue Basis entsteht, indem die Spaltenvariable \(q_{C_p}\) die bisherige Zeilenvariable \(q_{R_k}\) ersetzt. Schreibt man \(\mathcal B’=\mathcal B E\), wobei \(E\) die Einheitsmatrix ist, deren \(k\)-te Spalte durch \(u\) ersetzt wurde, so folgt \(\mathcal B’^{-1}=E^{-1}\mathcal B^{-1}\). Die rechte Seite kann daher ohne vollständiges Tableau aktualisiert werden: \(\beta’=E^{-1}\beta\). Für \(i\ne k\) folgt aus der Struktur von \(E\) explizit \(\beta’_k=\beta_k/u_k\) und \(\beta’_i=\beta_i-u_i\beta’_k\). Damit ersetzt ein Eta-Update die vollständige Tableau-Aktualisierung.

Der Vorteil liegt nicht darin, ein zweites duales System aufzubauen. Revised INLIF benötigt keine getrennte Führung von \(A\) und \(A^\top\) als zwei eigenständigen Ungleichungssystemen. Die transponierten Gleichungen \(\mathcal B^\top\lambda_k=e_k\) und, im LP-Fall, \(\mathcal B^\top y=\breve c_R\), beziehen sich auf dieselbe aktuelle Basis \(\mathcal B\). Die duale Information entsteht also aus der Kreuztableau-Struktur selbst und nicht aus einer separaten INLIF-Behandlung des dualen Problems.

Für ein LP

wird der erweiterte Kostenvektor \(\breve c=(c,0)^\top\) verwendet. Sind \(\breve c_R\) die Kosten der aktuellen Zeilenvariablen, so wird der Zielfunktionsmultiplikator \(y\) aus \(\mathcal B^\top y=\breve c_R\) bestimmt. Für eine Spaltenvariable \(q_j\) mit Spalte \(h_j\) lautet der reduzierte Koeffizient \(\overline c_j=\breve c_j-h_j^\top y\). Die dualen Ungleichungen werden also nicht als separates System gespeichert, sondern aus den reduzierten Koeffizienten der aktuellen Basis gelesen. Für die ursprünglichen \(x\)-Spalten ergibt sich \(\overline c_{x_j}=c_j-A_{\cdot j}^\top y\), und für die ursprünglichen Schlupfspalten \(\overline c_{s_i}=-y_i\). Die Bedingungen \(\overline c_{x_j}\le0\) und \(\overline c_{s_i}\le0\) entsprechen daher \(A^\top y\ge c\) und \(y\ge0\). Zusammen mit einem primal zulässigen \(x\) und \(b^\top y-c^\top x=0\) liegt ein primal-duales Optimalitätszertifikat vor.

Besonders natürlich ist im revised INLIF die Gewinnung eines Farkas-Zertifikats. Sei \(k\) eine verletzte aktive Prüfzeile mit \(\beta_k<0\). Für diese Zeile wurde zur Pivotprüfung bereits \(\mathcal B^\top\lambda_k=e_k\) gelöst. Damit gilt

Ist die Zeile blockiert, existiert also weder ein freier Pivot noch ein zulässiger Kegelpivot, so kann \(\lambda_k\) unmittelbar gegen das Originalsystem geprüft werden.

Für freie ursprüngliche Variablen muss in einer blockierten Zeile \(A^\top\lambda_k=0\) gelten. Denn jeder von null verschiedene Koeffizient einer freien Spaltenvariable würde einen freien Pivot ermöglichen. Für Kegelvariablen dürfen in der blockierten Zeile nur nichtnegative Koeffizienten verbleiben. Da die Schlupfspalten die Einheitsmatrix bilden, sind diese Koeffizienten gerade die Komponenten von \(\lambda_k\). Damit folgt \(\lambda_k\ge0\). Wird also verifiziert, dass

gilt, so ist

Tatsächlich würde aus einem hypothetischen \(x\) mit \(Ax\le b\) der Widerspruch folgen:

Bei Systemen mit zusätzlicher Nichtnegativitätsbedingung \(x\ge0\) lautet die entsprechende Farkas-Prüfung

Denn für jedes hypothetische \(x\ge0\) mit \(Ax\le b\) ergäbe sich dann der Widerspruch

Damit liefert revised INLIF die folgende zertifizierte Alternative:

Im LP-Fall kommt als dritter Abschluss ein primal-duales Optimalitätszertifikat hinzu:

Unbeschränktheit wird dagegen nicht allein aus einer offenen Verbesserungsspalte geschlossen. Sie wird erst akzeptiert, wenn ein Strahl \(d\) mit \(d\ge0\), \(Ad\le0\) und \(c^\top d>0\) im Originalsystem verifiziert ist.

Die Schrittzahlschranke ergibt sich aus der Monotonie der Absorptionen. Jede freie Absorption entfernt eine noch freie Variable aus dem aktiven Restsystem. Kegelabsorptionen entstehen nur aus den bereits aktivierten Ursprungslagen. Unter der Extremalregel und der Nichtreaktivierung aktivierter Hyperebenen gilt daher in den stabilen Reduktionsklassen

Für LPs mit \(n\) primalen Entscheidungsvariablen und \(m\) Schlupf- bzw. Zeilenvariablen wird entsprechend die strukturelle Grenze \(N_{\mathrm{Absorption}}\le m+n\) verwendet, sofern jede Tableauvariable höchstens einmal in die berechnete Rolle absorbiert wird. Das Erreichen dieser Grenze ist jedoch für sich allein kein Unzulässigkeitsbeweis. Danach muss entweder ein primal-duales Optimalitätszertifikat oder ein verifiziertes Farkas- bzw. Unbeschränktheitszertifikat vorliegen. Andernfalls ist die Monotonie- oder Pivotinvariante verletzt, und der Fall bleibt algorithmisch offen.

Bezeichne \(s\) die Anzahl der Absorptionsschritte und \(F_t\) den Aufwand einer Basislösung im Schritt \(t\). Werden in Schritt \(t\) im Mittel \(q_t\) Kandidatenspalten geprüft und besitzt die Matrix eine wirksame Dichte \(\delta\), so ergibt sich für revised INLIF die Arbeitsabschätzung

Bei vollständiger Spaltenprüfung ist \(q_t=\mathcal O(n)\). Für \(m=\mathcal O(n)\) und \(s=\mathcal O(n)\) ergibt sich dann der Scananteil \(\mathcal O(L\delta n^3)\). Bei lokaler Pivotprüfung mit \(q_t=\mathcal O(1)\) sinkt dieser Anteil formal auf \(\mathcal O(L\delta n^2)\). Die Kosten der Basislösungen bleiben zusätzlich zu berücksichtigen. Im dichten Worst Case kann die bekannte kubische Ordnung bestehen bleiben. Der entscheidende Gewinn des revised INLIF liegt daher in geringerem Speicherbedarf, Sparse-Ausnutzung und darin, dass das dichte Tableau \(G=\mathcal B^{-1}\mathcal N\) nicht vollständig gespeichert werden muss. \(\square\)

Problemstellung und Bedeutung von \(\operatorname{pro}\).

Die Problemstellung ist wie beim INLIF-Verfahren. Der Parameter \(\operatorname{pro}\in\{1,-1\}\) legt die Interpretation der Eingabedaten fest.

Für \(\operatorname{pro}=1\) wird das eingegebene Problem als maximierendes primales Problem gelesen: \[ \max \; e+c^\top x \qquad\text{unter}\qquad Ax\le b,\qquad x\ge 0. \] Für \(\operatorname{pro}=-1\) wird die interne primal-duale Orientierung vertauscht. Das Verfahren arbeitet dann mit der transponiert-negierten Tableauform, so dass die interne Zielfunktion dieselbe geometrische Rolle spielt wie im maximierenden Fall. Der mathematische Zielwert darf dabei jedoch nicht mit falschem Vorzeichen ausgegeben werden. Deshalb wird der intern skalierte Zielfunktionswert bei der Rücktransformation mit \(\operatorname{pro}\) versehen: \[ f_{\mathrm{orig}} = e+\operatorname{pro}\,\kappa f_{\mathrm{int}}, \qquad \kappa>0. \] Der Faktor \(\operatorname{pro}\) gehört also zur Rücktransformation des internen Zielwertes, nicht zu den primalen und dualen Zulässigkeitsbedingungen selbst. Die Prüfgrößen \[ Ax\le b,\qquad A^\top y\ge c,\qquad b^\top y-c^\top x\ge 0 \] werden immer in der jeweils zurücktransformierten Originalorientierung ausgewertet.

Satz.
Das hybride Intex–Merit-Verfahren bestimmt für jedes lösbare LP der obigen Form eine primal-duale Lösung \((x^o,y^o)\), sofern die verwendete Rechengenauigkeit, die Relaxationsfolge und der abschließende CO die aktive Struktur des Problems auflösen. Der Intex-Kern arbeitet durch Inter-/Extrapolationen. Für eine Matrixdichte \(d\in[0,1]\) und eine geeignete binäre Genauigkeitsskala \({}_2(\tilde\alpha\rho)\) besitzt im unstrukturierten Fall die geometrische Intex-Phase die Ordnung \[ \mathcal O\!\left({}_2(\tilde\alpha\rho)^2\,dmn\right). \] Wird die lineare Algebra durch FTD/FSR strukturiert ausgeführt, tritt an die Stelle von \(dmn\) der jeweilige strukturierte Operatoraufwand. Sei \(z:=m+n.\)

Intex-Relaxation.
Statt unmittelbar \[ Ax\le b,\qquad A^\top y\ge c,\qquad b^\top y-c^\top x=0 \] zu erzwingen, betrachtet das Intexverfahren mit \(r\in[0,\rho]\) eine Familie relaxierter primal-dualer Polytope \[ P_r := \left\{ (x,y)^\top\in{}^{\nu}\mathbb R_{\ge 0}^{z} : b^\top y-c^\top x\le r,\; Ax-b\le r_{m\times 1},\; c-A^\top y\le r_{n\times 1} \right\}. \] Der Anfangsradius wird so gewählt, dass der Nullpunkt \(0\) sicher in \(P_\rho\) enthalten ist. Eine mögliche Wahl ist \[ \rho := s\left|\min\{b_1,\dots,b_m,-c_1,\dots,-c_n\}\right|, \qquad s\in]1,2], \] gegebenenfalls nach vorheriger Skalierung und Verschiebung der Daten.

Das ursprüngliche LP ist genau dann gelöst, wenn ein Paar \((x,y)\in P_0\) gefunden wird. Wegen starker Dualität² ist dies äquivalent zu \[ Ax\le b,\qquad x\ge0,\qquad A^\top y\ge c,\qquad y\ge0,\qquad b^\top y-c^\top x=0. \] Damit lösen dieselben Variablen zugleich \[ \max\{c^\top x:x\in{}^{\nu}\mathbb R_{\ge0}^n,\ Ax\le b\} \] und \[ \min\{b^\top y:y\in{}^{\nu}\mathbb R_{\ge0}^m,\ A^\top y\ge c\}. \]

Normierung und Skalierung.
Vor der geometrischen Iteration werden Zielfunktion, rechte Seiten und Restriktionsmatrix normiert. Ziel ist nicht eine Veränderung des Problems, sondern eine balancierte Repräsentation der drei Fehlerarten \[ Ax-b,\qquad c-A^\top y,\qquad b^\top y-c^\top x. \] Die Skalierung soll vermeiden, dass eine der drei Größen numerisch dominiert. Um zu garantieren, dass die Zielfunktionskopplung im Gradientenabstieg der anschließenden Merit-Politur weder verschwindet noch die Hinge-Fehler numerisch überdeckt, wird initial eine skalare Lineartransformation angewendet. Konkret werden der Vektor \(c\) sowie der konstante Offset \(e\) mit dem Faktor \(\lambda = \|c\|^{-1}\) (für \(c \neq 0\)) skaliert.

Unter dieser Transformation bleibt die primale Lösung \(x\) invariant, während die duale Lösung \(y\) und die absoluten Zielfunktionswerte streng linear mit \(\lambda\) skalieren. Am Ende des Gesamtverfahrens erfolgt eine exakte Rücktransformation in den ursprünglichen physikalischen Raum. Nach der Rechnung werden alle Variablen und Zielwerte in die ursprüngliche Orientierung zurücktransformiert. Insbesondere wird bei \(\operatorname{pro}=-1\) der intern entstandene Zielwert mit dem korrekten Vorzeichen ausgegeben.

Geometrischer Intex-Schritt.
Innerhalb eines festen Polytops \(P_r\) wird ein innerer Schwerpunkt approximiert. Für jede Koordinate wird ein zulässiges Intervall bestimmt. Für \(v \in P_r\) wird näherungsweise \[ v_k^* := \min\check v_k+\max\check v_k, \qquad k=1,\dots,z, \] gebildet. So entsteht ein geometrischer Mittelpunkt \(v=(x,y)^\top\), der innerhalb der aktuellen Relaxation bleibt und als stabiler Ausgangspunkt für die nächste Radiusverkleinerung dient.

Aus den letzten Schwerpunktbewegungen wird anschließend extrapoliert. Dazu wird mit \(w\ge 0\) eine Richtung \(\Delta v\) gebildet und ein maximal zulässiger Schritt entlang \(v(w)=v+w\Delta v\) bestimmt. Der neue Radius ergibt sich aus der kleinsten noch notwendigen Relaxation \[ r_{\mathrm{neu}} = \min\{r\ge0:v(w)\in P_r\}. \] Gelingt \(r_{\mathrm{neu}}=0,\) so ist ein primal-dual optimales Paar erreicht. Bleibt \(r_{\mathrm{neu}}>0\) stabil von null getrennt, deutet dies auf Unlösbarkeit, fehlende Genauigkeit oder eine nicht aufgelöste aktive Struktur hin.

Intex als zulässiger Startpunkt.
Der wesentliche Vorteil der Intex-Phase besteht darin, dass sie nicht blind im Außenraum sucht. Sie arbeitet mit relaxierten zulässigen Bereichen und liefert typischerweise einen Punkt, der primal-dual bereits gut orientiert ist: \[ Ax-b\lesssim r_{m \times 1},\qquad c-A^\top y\lesssim r_{n \times 1},\qquad b^\top y-c^\top x\lesssim r. \] Dieser Punkt ist noch nicht notwendig optimal, aber häufig nahe an einer aktiven Fläche des optimalen Polytops. Genau dort setzt die Merit-Politur an.

Primal-duale Merit-Bereinigung und algebraischer Crossover

Residuen.
Für \(x\ge0\) und \(y\ge0\) seien \[ r^p(x):=\max(Ax-b,0)\in{}^{\nu}\mathbb R^m, \] \[ r^d(y):=\max(c-A^\top y,0)\in{}^{\nu}\mathbb R^n \] die primalen und dualen Verletzungen. Der signierte Dualitätsspalt ist \[ g(x,y):=b^\top y-c^\top x. \] Für ein maximierendes primales Problem ist nur \[ g_+(x,y):=\max(g(x,y),0) \] als positiver Spalt zulässigkeitsrelevant. Im exakten Optimum gilt \[ r^p(x)=0,\qquad r^d(y)=0,\qquad g(x,y)=0. \] Die unskalierten Kontrollgrößen lauten \[ p_\infty=\|r^p(x)\|_\infty, \qquad d_\infty=\|r^d(y)\|_\infty, \qquad g_{\mathrm{abs}}=|b^\top y-c^\top x|. \] Die primalen und dualen Zielfunktionswerte sind \( p^*=c^\top x\) und \( d^*=b^\top y. \) Die konstante Verschiebung \(e\) wird erst bei der Ausgabe addiert.

Merit-Funktion.
Die Merit-Funktion fasst primale Verletzung, duale Verletzung und Dualitätsspalt zusammen: \[ \widehat{\Phi}(x,y) = \|r^p(x)\|_2^2 + \|r^d(y)\|_2^2 + \gamma^2\,g_+(x,y)^2, \qquad \gamma>0. \] In skalierten Variablen \( x=S_x\breve x\) und \(y=S_y\breve y \) und mit diagonalen Residualgewichten \( D_p=\operatorname{diag}(d_p), \) und \( D_d=\operatorname{diag}(d_d) \) lautet die skalierte Form

\( \widehat\Phi(\breve x,\breve y) = \|D_p\max(AS_x\breve x-b,0)\|_2^2\) \( + \|D_d\max(c-A^\top S_y\breve y,0)\|_2^2 + d_g^2 \max(b^\top S_y\breve y-c^\top S_x\breve x,0)^2. \)

Optional kann ein kleiner Regularisierungsterm \( \varepsilon\left(\|\breve x\|_2^2+\|\breve y\|_2^2\right) \) hinzugenommen werden.

Gradient.
Mit \[ w_i^p=d_{p,i}^2\,r_i^p(\breve x), \qquad w_j^d=d_{d,j}^2\,r_j^d(\breve y), \qquad w_g=d_g^2\,g_+(\breve x,\breve y) \] ergeben sich die Gradienten \[ \frac{{\downarrow}\Phi}{{\downarrow}\breve x_j} = s_{x,j} \left( {\LARGE{\textbf{+}}}_{i=1}^m A_{ij}w_i^p – c_j w_g \right), \qquad j=1,\dots,n, \] \[ \frac{{\downarrow}\Phi}{{\downarrow}\breve y_i} = s_{y,i} \left( b_i w_g – {\LARGE{\textbf{+}}}_{j=1}^n A_{ij}w_j^d \right), \qquad i=1,\dots,m. \] Der Term \(b_iw_g\) ist wesentlich: Der Gradient bezüglich \(y\) muss den Anteil von \(b^\top y\) im Dualitätsspalt enthalten. Ohne diesen Beitrag wird die duale Seite falsch gesteuert.

Zulässige Schrittführung.
Ausgehend vom Intex-Punkt wird ein Abstiegsschritt berechnet. In der einfachen projizierten Variante gilt \[ \breve x‘ = \Pi_{\mathbb R_+^n} \left( \breve x-\alpha\nabla_{\breve x}\widehat\Phi \right), \qquad \breve y‘ = \Pi_{\mathbb R_+^m} \left( \breve y-\alpha\nabla_{\breve y}\widehat\Phi \right). \] In der zulässigen hybriden Variante wird auf eine harte Orthantprojektion während der Merit-Phase verzichtet. Stattdessen wird ein Schritt verworfen oder verkürzt, wenn er die zulässige Relaxation oder die Nichtnegativität verlassen würde. Dadurch bleibt der Charakter der Intex-Relaxation erhalten.

Ein Schritt wird nur akzeptiert, wenn eine Armijo-artige Bedingung erfüllt ist: \[ \widehat\Phi(\breve x‘,\breve y‘) \le \widehat\Phi(\breve x,\breve y) – c_1\alpha \left( \|\breve x‘-\breve x\|_2^2+ \|\breve y‘-\breve y\|_2^2 \right). \] Andernfalls wird \(\alpha\) reduziert.

Barzilai-Borwein-Schrittweite.
Nach einem akzeptierten Schritt werden \[ s_k= \begin{pmatrix} \breve x^{\grave{k}}-\breve x^k\\[0.3em] \breve y^{\grave{k}}-\breve y^k \end{pmatrix}, \qquad q_k= \begin{pmatrix} \nabla_{\breve x}\widehat{\Phi}(\breve x^{\grave{k}},\breve y^{\grave{k}})-\nabla_{\breve x}\widehat{\Phi}(\breve x^k,\breve y^k)\\[0.3em] \nabla_{\breve y}\widehat{\Phi}(\breve x^{\grave{k}},\breve y^{\grave{k}})-\nabla_{\breve y}\widehat{\Phi}(\breve x^k,\breve y^k) \end{pmatrix} \] gebildet. Für \( \eta_k:=s_k^\top q_k>0 \) wird \( \alpha_{\grave{k}} = \tilde{\eta}_k\, s_k^\top s_k \) gesetzt und anschließend auf ein Intervall \( [\alpha_{\min},\alpha_{\max}] \) beschränkt. Ist \(\eta_k\le 0\) wird die bisherige brauchbare Schrittweite beibehalten oder eine konservative Startschrittweite verwendet.

Zweck.
Die Merit-Phase ist ein Verfahren erster Ordnung. Sie reduziert Residuen und Dualitätsspalt, liefert aber nicht immer exakt die aktive Basis. Deshalb folgt bei Bedarf ein algebraischer CO. Er versucht aus dem Näherungspunkt \((x,y)\) eine exakt zulässige primal-duale Lösung zu rekonstruieren.

Aktive Kandidaten.
Primal aktive Restriktionen haben solche Indizes \(i\), für die \(a_i^\top x\approx b_i\) gilt. Aktive Nichtnegativitätsbedingungen haben solche Indizes \(j\), für die \(x_j\approx0\) gilt. Aus diesen Kandidaten werden Gleichungssysteme der Größe \(n\) gebildet. Ein primaler Kandidat wird akzeptiert, wenn das rekonstruierte \(x\ge0\) auch \(Ax\le b\) erfüllt. Anschließend wird ein dualer Kandidat \(y\ge 0\) gesucht mit \(A^\top y\ge c\). Nur wenn zusätzlich \(|b^\top y-c^\top x|\) klein genug ist, wird das Paar als optimale Lösung akzeptiert.

Zertifizierter CO.
Ein isoliert guter primaler Kandidat reicht nicht aus. Ebenso reicht ein isoliert guter dualer Kandidat nicht aus. Übernommen wird nur ein primal-duales Paar, das gleichzeitig \[ p_\infty=0,\qquad d_\infty=0,\qquad g_{\mathrm{abs}}=0 \] bis zur geforderten Toleranz erfüllt. Dadurch wird vermieden, dass ein nur formal zulässiger dualer Kandidat mit falschem Zielwert die Lösung ersetzt.

Gestufter CO.
Der gestufte CO begrenzt die kombinatorischen Kosten maßgeblich.

In Stufe 1 ergeben sich die wahrscheinlich aktiven Restriktionen aus kleinen Slacks \(b_i-a_i^\top x\) und Variablen \(x_j.\) Diese Stufe ist schnell, kann aber bei schlecht getrennten aktiven Mengen scheitern.

In Stufe 2 wird der Kandidatensatz erweitert. Zusätzliche Restriktionen mit mittleren Slacks werden zugelassen. Dadurch steigt die Erfolgswahrscheinlichkeit, aber auch die Zahl der getesteten Basen.

In Stufe 3 wird ein vollständiger oder nahezu vollständiger Fallback zugelassen. Diese Stufe ist teuer, aber robust. Sie wird nur verwendet, wenn vorherige Stufen keine zertifizierte primal-duale Lösung liefern.

Der CO kann außerdem vollständig abgeschaltet werden. Dann gibt das Verfahren den besten bis dahin gefundenen Intex–Merit-Punkt aus, ohne algebraische Basisrekonstruktion. Flankiert wird diese algorithmische Begrenzung durch eine strikt speichereffiziente Implementierung auf Softwareebene. Die innere Schleife der kombinatorischen Kandidatensuche arbeitet vollständig allokationsfrei.

Sämtliche Vektoren und Index-Arrays werden vorab präallokiert und pro Iteration lediglich überschrieben. Dadurch wird kein dynamischer Speicherdruck erzeugt und die Garbage Collection der Laufzeitumgebung wird in der kritischen Phase nicht ausgelöst. Durch diese mechanische Effizienz verliert die theoretische Worst-Case-Komplexität der Basissuche in der Praxis ihren Schrecken, da das Durchlaufen der Permutationen lediglich minimale CPU-Zyklen beansprucht.

Timeout.
Wird ein Zeitlimit erreicht, bricht das Verfahren kontrolliert ab und liefert den besten bis dahin bekannten Zustand: \[ x_{\mathrm{best}},\qquad y_{\mathrm{best}},\qquad p_\infty,\qquad d_\infty,\qquad g_{\mathrm{abs}},\qquad p^*,\qquad d^*. \] Ist bis dahin keine zertifizierte Lösung gefunden worden, wird dies ausdrücklich angezeigt. Ein Timeout ist daher kein mathematischer Fehlschluss, sondern ein unvollständiger Rechenstand.

Ein primal unzulässiges LP hat kein \(x\ge0\) mit \(Ax\le b\): Weder Intex kann den Radius bis \(r=0\) senken, noch Merit \(p_\infty=0\) erzwingen. Typischerweise bleibt \( p_\infty>0\) oder die Radiusreduktion stagniert. Ist das duale LP unzulässig, dann existiert kein \(y\ge0\) mit \(A^\top y\ge c.\) Dann bleibt gewöhnlich \(d_\infty>0.\) Primal entspricht dies häufig einem unbeschränkten primalen Problem.

Ist das Problem zwar formal lösbar, aber numerisch schlecht skaliert oder degeneriert, können Intex und Merit stagnieren, obwohl eine Lösung existiert. Dann entscheidet der CO, ob die aktive Struktur trotzdem algebraisch rekonstruiert werden kann. Gelingt dies nicht, liefert das Verfahren keinen falschen Optimalitätsnachweis, sondern einen nichtzertifizierten Näherungspunkt.

Damit unterscheidet das Verfahren praktisch vier Fälle: \[ \begin{array}{ll} p_\infty=d_\infty=g_{\mathrm{abs}}=0 &\text{zertifiziert gelöst},\\[0.3em] p_\infty>0 &\text{primale Zulässigkeit nicht erreicht},\\[0.3em] d_\infty>0 &\text{duale Zulässigkeit nicht erreicht},\\[0.3em] p_\infty,d_\infty\approx0,\ g_{\mathrm{abs}}>0 &\text{zulässig, aber Optimalität nicht erreicht}. \end{array} \] Ein vollständiges Unlösbarkeitszertifikat im Sinne eines Farkas-Zeugen ist damit noch nicht automatisch enthalten. Es kann jedoch als zusätzliche Diagnosestufe ergänzt werden.

Strukturierte lineare Algebra durch FTD und FSR

Motivation.
Intex, Merit und CO benötigen wiederholt Operationen der Form \[ Ax,\qquad A^\top y,\qquad b^\top y,\qquad c^\top x. \] Für eine allgemeine dünn besetzte Matrix kostet dies \( \mathcal O(\operatorname{nnz}(A)).\) Bei einer verwertbaren inneren Struktur von \(A\) als Toeplitz-, Hankel-, blockzirkulanter, faltungsartiger oder im Inneren nahezu translationsinvarianter Matrix lassen sich die Operatorauswertungen durch FTD und FSR ersetzen oder zumindest wesentlich komprimieren.

FTD.
Die FTD kann in der Praxis oftmals die Matrizen mit \(mn\) Einträgen durch wenige Generatoren, Diagonalfolgen, Verschiebungsdefekte oder Faltungsdaten beschreiben und damit operatoriell betrachten. Eine geeignete Modellform ist\[ A\approx K_s+UV^\top+E, \qquad K_s=T_1T_2\cdots T_s. \] Dabei sind \(T_1,\dots,T_s\) Toeplitz-, Hankel-, blockzirkulante oder verwandte Faltungsoperatoren, \[ U,V\in{}^{\nu}\mathbb R^{n \times r}, \qquad r\ll n, \] beschreiben eine niedrig-rangige Rand- oder Defektkorrektur, und \(E\) bezeichnet den verbleibenden Approximationsfehler. Der Fall \(s=1\) liefert einen einfachen Toeplitz- oder Faltungskern. Der Fall \(s>1\) erlaubt eine mehrfach iterierte Generierung des strukturierten Kerns.

Der Nutzen dieser Mehrfachiterierung besteht darin, dass der verbleibende numerische Rang der Korrektur sinken kann. Formal lässt sich dies durch \[ r_\varepsilon(s) = \min_{T_1,\dots,T_s} \operatorname{rank}_\varepsilon \left( A-T_1T_2\cdots T_s \right) \] wiedergeben. Hierbei bezeichnet \(r_\varepsilon(s)\) den numerischen Restrang nach Wahl von \(s\) strukturierten Faktoren. Ein zusätzlicher Faktor lohnt sich insbesondere dann, wenn die Verringerung von \(r_\varepsilon(s)\) den zusätzlichen Aufwand seiner Auswertung überwiegt.

Die Faktoren werden aus lokalen Mittelwerten, diagonalen Summen, Verschiebungsdefekten und rekursiven Ausgleichsschritten gewonnen. Statt der vollen Matrix werden nur Generator-, Diagonal-, Faltungs- und Korrekturdaten gespeichert. Matrix-Vektor-Produkte, transponierte Produkte und ausgewählte Projektionen können dadurch strukturiert ausgeführt werden.

Mit der effektiven Strukturbreite \(\ell\) und \( p=\ell n \) ergibt sich für die strukturierte Kernrechnung typischerweise ein Aufwand der Form \( \mathcal O(p^2) \) anstelle dichter unstrukturierter Kosten. Bei FFT-fähigen Faltungskernen können Matrix-Vektor-Produkte in der Form \( \mathcal O(s\,n{}_2n+nr) \) organisiert werden, wobei \(s\) die Zahl der strukturierten Faktoren und \(r\) den Rang der Korrektur bezeichnet.

FSR.
Die FSR ergänzt die FTD arithmetisch. Während die FTD eine strukturierte Darstellung des Operators liefert, organisiert FSR die tatsächliche Rechnung rekursiv und wortlängenschonend. Mit \(\ll v\) als Verschiebung um \(v\) Stellen lässt sich eine Matrix oder ein Operator in hoch- und niedriggewichtige Anteile zerlegen: \[ A=(A_{\mathrm{hi}}\ll v)+A_{\mathrm{lo}}, \] Produkte, Schur-Komplemente, inverse Teilsysteme und Rekonstruktionen werden rekursiv aus kleineren Anteilen aufgebaut. Somit bleiben Zwischenwerte kontrolliert und die BigInt- beziehungsweise BigFloat-Arithmetik wächst nicht unnötig an. Die Rekursion ist hier gut parallelisierbar. Daher eignet sich \ac{fsr} besonders für hohe Präzision, exakte Ganzzahlarithmetik und große strukturierte lineare Systeme.

Niedrigrangige Korrekturen.
Ist \( A\approx K_s+UV^\top \) und sind die Kernlösungen mit \(K_s\) strukturiert verfügbar, kann die niedrig-rangige Korrektur zusätzlich durch die Woodbury-Formel abgespalten werden: \[ (K_s+UV^\top)^{-1} = K_s^{-1} – K_s^{-1}U \left(I+V^\top K_s^{-1}U\right)^{-1} V^\top K_s^{-1}. \] Dies reduziert ein großes lineares Gleichungssystem auf strukturierte Kernlösungen und ein kleines Korrektursystem. Für den einfachen Toeplitz-Kern \( A=T+UV^\top \) ergibt sich darüber hinaus die Resolventenform \[ \det(\lambda I-A) = \det(\lambda I-T) \det\!\left(I_r-V^\top(\lambda I-T)^{-1}U\right). \] Das zweite Eigenprodukt besitzt nur die Größe \(r\times r\). Eigenwertanteile, die durch die Rand- oder Defektkorrektur verursacht werden, lassen sich daher als Nullstellen der kleinen rationalen Gleichung \[ \det\!\left(I_r-V^\top(\lambda I-T)^{-1}U\right)=0 \] untersuchen. So lassen sich Ausreißer-Eigenwerte von dem durch \(T\) bestimmten Kernspektrum trennen.

Zusammenspiel.
liefert die strukturierte Näherung \[ A\approx K_s+UV^\top+E, \qquad K_s=T_1T_2\cdots T_s, \] während die \ac{fsr} die daraus entstehenden Produkte, Inversionen, Schur-Komplement-Schritte und Rücktransformationen rekursiv ausführt.\\ \\ Die Parameter \( s, r, \|E\| \) beschreiben dabei drei verschiedene Aspekte der Struktur: \(s\) ist die Tiefe des struk- turierten Kerns, \(r\) der Rang der verbleibenden Rand- oder Defektkorrektur, und \(E\) der nicht erfasste Rest. Die algorithmisch günstige Zone liegt bei \(\qquad s\ll n,\qquad r\ll n,\qquad \|E\|\approx 0. \)

Einbindung in Intex–Merit.
Die FTD/FSR-Struktur beschleunigt dreifach:

Erstens die Intex-Phase, weil dort wiederholt Schranken, Schnitte und Schwerpunktbewegungen durch Operationen mit \(A\) und \(A^\top\) berechnet werden.

Zweitens bdie Merit-Bereinigung, weil Funktionswert und Gradient wesentlich bestehen aus \[ Ax-b,\qquad c-A^\top y,\qquad b^\top y-c^\top x. \]

Drittens ggf. den CO. Das Lösen eines dichten aktiven Systems der Größe \(n\) ohne Struktur kostet \(\mathcal O(n^3)\). Eine Toeplitz-, Hankel-, blockzirkulante oder faltungsartige Struktur ermöglicht FTD/FSR eine Reduktion auf strukturierte Kernoperationen, Schur-Komplement-Schritte und kleine Korrektursysteme.

Damit bilden Intex, Merit, CO, FTD und FSR keine getrennten Bausteine, sondern eine zusammenhängende Architektur: Intex = geometrische Reduktion,\(\quad\) Merit = numerische Bereinigung,\(\quad\) CO = algebraische Zertifizierung\(\quad\) FTD/FSR = strukturierte lineare Algebra.

Bezeichne \(\mathcal M_A\) den Aufwand einer vollständigen strukturerhaltenden Operatorauswertung mit \(A\) und \(A^\top\). Je nach Darstellung von \(A\) ergibt sich \[ \mathcal M_A= \begin{cases} \mathcal O(\operatorname{nnz}(A)), &\text{für dünn besetzte Matrizen},\\[0.3em] \mathcal O(mn), &\text{für dichte unstrukturierte Matrizen},\\[0.3em] \mathcal O(p^2), &\text{für} \mathrm{FTD}/\mathrm{FSR}\text{-strukturierte Kernoperationen},\\[0.3em] \mathcal O(s\,n {}_2n+nr), &\text{für} \mathrm{FFT}\text{-fähige Faltungsfaktoren mit Rang-}r\text{-Korrektur}. \end{cases} \] Hier bezeichnet \(p=\ell n\) die effektive strukturierte Problemgröße bei Strukturbreite \(\ell\), \(s\) die Zahl der strukturierten Faktoren im Kern \( K_s=T_1T_2\cdots T_s, \) und \(r\) den Rang der niedrig-rangigen Rand- oder Defektkorrektur. Die Intex-Phase besitzt damit die Ordnung \( \mathcal O\!\left({}_2(\tilde\alpha\rho)^2\,\mathcal M_A\right). \) Die Merit-Bereinigung benötigt bei \(K\) akzeptierten Iterationen und \(B\) zusätzlichen Backtracking-Auswertungen den Aufwand \( \mathcal O\!\left((K+B)\mathcal M_A\right). \)

Für den CO gilt ohne zusätzliche Struktur \( \mathcal C=\mathcal O(Tn^3), \) wenn \(T\) aktive Basiskandidaten tatsächlich getestet werden. Da für \(Q\) Kandidatengleichungen im ungünstigen Fall \( T\le \binom{Q}{n} \) gelten kann, ist ein ungezügelter CO potenziell kombinatorisch. Das Stufenmodell reduziert daher zunächst \(Q\) und aktiviert größere Kandidatenmengen erst, wenn kleinere aktive Systeme keine hinreichende Zertifizierung liefern.

Besitzt das aktive Teilsystem dieselbe Toeplitz-, Hankel-, blockzirkulante oder faltungsartige Struktur wie der Ausgangsoperator, kann der CO statt durch dichte Faktorisierungen durch strukturierte Kernoperationen, Schur-Komplement-Schritte und kleine Korrektursysteme ausgeführt werden. In diesem Fall ergibt sich näherungsweise \( \mathcal C=\mathcal O(Tp^2) \) beziehungsweise bei FFT-fähigen Produktkernen und niedrig-rangiger Korrektur eine entsprechende Abhängigkeit von \( s\,n {}_2n+nr. \)

Damit lautet die Gesamtordnung des hybriden Verfahrens \[ \mathcal O\!\left( {}_2(\tilde\alpha\rho)^2\mathcal M_A + (K+B)\mathcal M_A + \mathcal C \right). \] Die tatsächliche Laufzeit wird also durch drei Größen bestimmt: die Anzahl der geometrischen Intex-Schritte, die Anzahl der Merit-Auswertungen und den Aufwand der algebraischen Zertifizierung im CO.

Der Speicherbedarf beträgt neben der Matrix- oder Generatorrepräsentation im Intex–Merit-Kern \( \mathcal O(m+n). \) Bei dichter Speicherung ergibt sich insgesamt \( \mathcal O(mn+m+n)=\mathcal O(mn). \) Bei strukturierter Speicherung hängt der Speicherbedarf von der Generatorbreite, der Zahl der Faktoren und dem Rang der Korrektur ab. Für eine Darstellung \( A\approx K_s+UV^\top+E \) liegt er näher bei \( \mathcal O(p+nr+m+n), \) wobei \(p\) die strukturierte Kerndarstellung und \(nr\) die Speicherung der niedrig-rangigen Korrektur beschreibt.

Korollar aus dem INLIF-Verfahren.
Jedes LP kann in \(\mathcal O(Ln^3)\) gelöst oder als unlösbar zertifiziert werden.\(\square\)

Das LINEAL-Verfahren (Linear Interpolation, Nearing Extrapolation, and Amending Lifting) verbindet Intex und INLIF zu einem einheitlichen Lösungsansatz für lineare Ungleichungssysteme und lineare Programme. Intex nähert sich der Lösung durch lineare Interpolation und Extrapolation innerhalb schrittweise verkleinerter Relaxationen; verbleibende oder stagnierende Verletzungen können anschließend durch das korrigierende Lifting von INLIF absorbiert und zurückgerechnet werden. Dabei bleibt Intex ein eigenständig erfolgreiches Primärverfahren, während INLIF bei Bedarf als algebraischer Abschluss dient.

Zweite Lösung auf derselben Zielfläche.
Ist \(x^o\) eine optimale Lösung, so kann das LP \[ \max \left\{ \|x-x^o\|_1: c^\top x=c^\top x^o,\; Ax\le b,\; x-x^o\in[-1,1]^n,\; x\in{}^{\nu}\mathbb R_{\ge0}^n \right\} \] eine zweite optimale Lösung bestimmen, sofern die optimale Fläche mehr als einen Punkt enthält. Der duale Vektor \(y^o\) kann analog behandelt werden.

Lineare Gleichungssysteme.
Die äquivalente Darstellung \(Ax\le b,\ -Ax\le -b\) macht ein lineares Gleichungssystem \(Ax=b\) zum INLIF-System, das sich über die Merit-Funktion \(\widehat{\Phi}(x)=\|Ax-b\|_2^2\) behandeln lässt. Im regulären Fall führt dies auf die Lösung des Systems, im singulären oder inkonsistenten Fall nur im Kleinste-Quadrate-Sinn. Der Gradient ist \(\nabla\Phi(x)=A^\top(Ax-b)\). Zusätzliche Nebenbedingungen wie \(x_j\ge0\) lassen sich durch Projektion oder durch LP-Formulierung ergänzen.

Damit wird Gauß-Elimination nicht ersetzt, sondern in eine zulässigkeitsorientierte Reduktionssicht eingebettet. Der Vorteil kann in der frühen Terminierung liegen: INLIF kann stoppen, sobald der transformierte Ursprung zulässig wird oder eine Restkegelzeile ein Unlösbarkeitszertifikat liefert. Für reguläre Gleichungssysteme bleibt Gauß das natürliche Direktverfahren; für überbestimmte, inkonsistente oder mit Ungleichungen gekoppelte Systeme ist INLIF dagegen als strukturelle Alternative plausibel.

Die Reduktionslogik von INLIF ist verletzungsgesteuert. Das aktuelle Tableau unterscheidet zwischen verletzten Aktivierungszeilen und erfüllten passiven Zeilen. Erfüllte Zeilen bleiben Bestandteil des Systems, erzwingen aber keinen Absorptionsschritt. Dadurch kann das Verfahren die aktuelle Reduktion auf die verletzten Restriktionen beschränken und nicht verletzte Ungleichungen bis zu ihrer möglichen späteren Aktivierung zurückstellen.

Regularität.
Die Matrix \(A\) ist genau dann regulär, wenn das homogene System \(Ax=0\) nur die Nulllösung besitzt. Dies kann ein LP der Form \[ \max\{\|x\|_1:Ax=0,\;x\in[-1,1]^n\} \] testen. Ist der optimale Wert null, so ist der Kern trivial.

Inverse Spalten.
Für reguläres \(A\in{}^{\nu}\mathbb R^{n\times n}\) kann jede Spalte \(\alpha_j\) von \(A^{-1}\) durch \(A\alpha_j=e_j\) bestimmt werden, wobei \(e_j\) der \(j\)-te Einheitsvektor ist. Ohne Struktur ist dies ein klassisches lineares Gleichungssystem; mit FTD/FSR kann die Struktur von \(A\) in der Lösung ausgenutzt werden.

Eigenwertprobleme.
Eigenwertgleichungen \(Ax=\lambda x\) lassen sich in Nebenbedingungen eines erweiterten Optimierungsproblems überführen. Praktisch interessiert dies dann, wenn sich zusätzliche Nebenbedingungen, Normierungen oder Strukturbedingungen an \(x\) und \(\lambda\) stellen. Für allgemeine dichte Matrizen bleiben klassische Spektralverfahren meist geeigneter; bei strukturierter oder nebenbedingter Formulierung kann die Merit-LP-Sicht jedoch vorteilhaft sein.

Polynom- und Padé-Approximation.
Auch Ausgleichsprobleme der Form \[ y=r+Xc,\qquad X^\top r=0 \] lassen sich als lineare oder quadratische Merit-Probleme formulieren. Zusätzliche Nebenbedingungen an die Koeffizienten \(c\) können direkt eingebaut werden. Für Padé-Approximationen kann die diskrete Fouriertransformation strukturierte Faltungsanteile liefern, die FTD/FSR anschließend verarbeiten.

Konvexe Programme.
Für konvexe Programme \(\quad \min \left\{ f_1(x): x\in{}^{\nu}\mathbb R^n,\; (f_2(x),\dots,f_m(x))^\top\le0 \right\}\quad\) kann dieselbe Idee auf nichtlineare Residuen übertragen werden. Dann werden die linearen Terme \(Ax-b\) durch die Verletzungen \(f_k(x)\) ersetzt. Die Grundstruktur bleibt: \[ \text{Zulässigkeit}+\text{Dualität oder Optimalitätsbedingung} \quad\longrightarrow\quad \text{Merit-Funktion}. \] Eine polynomielle Laufzeitaussage hängt hier jedoch zusätzlich von der Auswertbarkeit, Glattheit, Kondition und Struktur der Funktionen \(f_k\) ab.

Shift-polynomiale FSR-Multiplikation für ganzzahlige Matrizen

Die klassische Multiplikation zweier Matrizen über \(\mathbb Z\) ist bei großen Einträgen nicht nur durch die Anzahl der Skalarprodukte, sondern wesentlich auch durch die Kosten der Ganzzahlmultiplikation bestimmt. Mit Einträgen großer Bitlänge für \( A,B\in \mathbb Z^{N\times N} \) kann daher eine rekursive Blockstruktur vorteilhaft sein, wenn sie die Speicherlokalität verbessert, die Teilprodukte in eine günstige Akkumulationsordnung bringt und temporäre Produktmatrizen vermeidet. Im Folgenden werde \(N\) zunächst durch \(3\) teilbar vorausgesetzt. Schreibe \[ A= \begin{pmatrix} A_{11}&A_{12}&A_{13}\\ A_{21}&A_{22}&A_{23}\\ A_{31}&A_{32}&A_{33} \end{pmatrix}, \qquad B= \begin{pmatrix} B_{11}&B_{12}&B_{13}\\ B_{21}&B_{22}&B_{23}\\ B_{31}&B_{32}&B_{33} \end{pmatrix}, \] wobei alle Blöcke dieselbe Größe besitzen. Der shift-polynomiale Ansatz kodiert die drei Blockspalten von \(A\) und die drei Blockzeilen von \(B\) formal durch Potenzen eines zentralen Shiftsymbols \(s\). Für die drei inneren Kopplungen definiert man \[ P_1 = (A_{11}+sA_{21}+s^2A_{31}) (B_{11}+s^3B_{12}+s^6B_{13}), \] \[ P_2 = (A_{12}+sA_{22}+s^2A_{32}) (B_{21}+s^3B_{22}+s^6B_{23}), \] \[ P_3 = (A_{13}+sA_{23}+s^2A_{33}) (B_{31}+s^3B_{32}+s^6B_{33}). \] Dann ist \[ P=P_1+P_2+P_3 \] ein formales Matrixpolynom in \(s\), dessen Koeffizienten unmittelbar die Ergebnisblöcke von \(AB\) enthalten. Genauer gilt \[ (AB)_{pq} = [s^{\,p-1+3(q-1)}]P \qquad (p,q\in\{1,2,3\}). \] Die Ordnung der Koeffizienten ist also \[ s^0,s^1,s^2,s^3,s^4,s^5,s^6,s^7,s^8 \] und entspricht den Blockpositionen \[ 11,21,31,12,22,32,13,23,33. \]

Satz zur Korrektheit.
Seien \(A,B\in R^{N\times N}\) Matrizen über einem assoziativen, im Allgemeinen nicht kommutativen Ring \(R\). Das Shiftsymbol \(s\) sei zentral, d. h. \( sX=Xs \) für alle auftretenden Blöcke \(X\). Dann liefert die obige Koeffizientenextraktion exakt das Matrixprodukt \(AB\). Beweis. Für \(r\in\{1,2,3\}\) enthält der Ausdruck \[ (A_{1r}+sA_{2r}+s^2A_{3r}) (B_{r1}+s^3B_{r2}+s^6B_{r3}) \] den Term \( A_{pr}B_{rq} \) genau beim Exponenten \( (p-1)+3(q-1). \) Nach Summierung über \(r=1,2,3\) ergibt der Koeffizient von \( s^{p-1+3(q-1)} \) somit \( {\LARGE{\textbf{+}}}_{r=1}^3 A_{pr}B_{rq} = (AB)_{pq}. \) Die Reihenfolge der Faktoren bleibt dabei erhalten; insbesondere wird keine Kommutativität der Matrixblöcke vorausgesetzt.\(\square\)

Rekursive FSR-Fassung.
Wird die obige Konstruktion rekursiv auf die Blockprodukte angewendet, entsteht eine \(3\)-äre FSR-Multiplikation. Die Rekursion wird nicht notwendig bis zur Skalarebene durchgeführt. Für die Blattgröße \(L\) ist eine hybride Variante günstiger: \[ N \longrightarrow \frac{N}{3} \longrightarrow \frac{N}{3^2} \longrightarrow \cdots \longrightarrow L. \] Auf der Blattebene wird zur klassischen Ganzzahl-Matrixmultiplikation zurückgekehrt. Damit ist die Rekursionstiefe \( d={}_{3}(N/L), \) sofern \(N/L\) eine Potenz von \(3\) ist. Für \(N=729\) und \(L=27\) ergibt sich beispielsweise \( d=3. \) Die Methode ist jedoch nicht auf reine Dreierpotenzen beschränkt. Es genügt, dass \( N=3^dL \) mit einer praktisch geeigneten Blattgröße \(L\) gilt. So ist etwa \( 1026=3^3\cdot38, \) und \(1053=3^3\cdot39.\) Die Wahl der Blattgröße ist dabei wesentlich. Eine zu kleine Blattgröße führt zu tiefer Rekursion, vielen kleinen Blattaufrufen und erhöhtem Verwaltungs- und Speicheraufwand. Die praktische Operation lautet nicht \( \qquad T=AB,\qquad C\leftarrow C+T,\qquad \) sondern direkt \( C\leftarrow C+AB. \) Dies vermeidet temporäre Produktmatrizen. Für große \(\mathbb Z\)-Einträge ist dies wesentlich, da die Speicherverwaltung sonst einen erheblichen Teil der Laufzeit ausmachen kann.

Vorzeichen.
Da die rekursive Fassung direkt in die Zielmatrix akkumuliert und keine Koeffizientenfenster durch ganzzahlige Divisionen zurückgewinnt, ist die Methode nicht auf nichtnegative Einträge beschränkt. Für gemischt vorzeichenbehaftete BigInt-Einträge gilt derselbe Korrektheitssatz. Die Vorzeichen werden vollständig durch die gewöhnliche Multiplikation und Addition in \(\mathbb Z\) behandelt. Dies unterscheidet diese Fassung von reinen Packungsverfahren, bei denen mehrere Koeffizienten in einer großen ganzen Zahl kodiert und anschließend durch Shifts oder Divisionen wieder getrennt werden. Bei gemischten Vorzeichen erfordert eine solche Rückgewinnung zusätzliche Balancierungs- oder Divisionsschritte. Die hier verwendete FSR-Fassung vermeidet diese Rekombination.

Empirischer Befund.
Die folgenden Messungen sind empirische Laufzeitbefunde für dichte BigInt- Matrizen. Sie hängen von Prozessorarchitektur, Threadzahl, Julia-Version, GMP-Anbindung, Blattgröße und Speicherverwaltung ab. Sie zeigen jedoch deutlich, dass der rekursive Verwaltungsaufwand bei hinreichend großen Matrizen überkompensiert werden kann. Für dichte Matrizen der Größe \( 2187\times2187 \) mit nichtnegativen BigInt-Einträgen von \(128\) Bit und hybrider Rekursion bis zur Blattgröße \( L=81 \) ergab sich auf einer \(24\)-Thread-Workstation ein serieller Speedup von etwa \( 3{,}54 \) gegenüber klassischer GMP-basierter Matrixmultiplikation. In der parallelen Variante wurde gegenüber der klassisch parallelisierten Multiplikation ein Speedup von etwa \( 1{,}38 \) gemessen. Für dichte Matrizen der Größe \( 729\times729 \) mit nichtnegativen BigInt-Einträgen von \(256\) Bit und hybrider Rekursion bis zur Blattgröße \( L=27 \) ergab sich auf derselben Workstation ein serieller Speedup von etwa \( 3{,}36 \) gegenüber klassischer GMP-basierter Matrixmultiplikation. In der parallelen Variante wurde ein Speedup von etwa \( 1{,}33 \) gegenüber der klassisch parallelisierten Multiplikation gemessen. Auch für gemischt vorzeichenbehaftete Einträge bleibt der Effekt erhalten. Für dichte Matrizen der Größe \( 729\times729 \) mit gemischten \(256\)-Bit-BigInt-Einträgen und Blattgröße \( L=27 \) wurde auf dem \(24\)-Thread-System ein serieller Speedup von etwa \( 3{,}80 \) und ein paralleler Speedup von etwa \( 1{,}41 \) gegenüber der entsprechenden klassischen GMP-basierten Multiplikation gemessen. Für die Größe \( 1026\times1026 \) mit gemischten \(256\)-Bit-BigInt-Einträgen gilt \( 1026=3^3\cdot38. \) Die Blattgröße \( L=38 \) ergab einen seriellen bzw. parallelen Speedup von etwa \(4{,}22\) bzw. \(1{,}32\). Somit beschränkt sich das Verfahren nicht auf reine Dreierpotenzen, sondern bleibt für Größen der Form \(3^dL\) sinnvoll einsetzbar. Die optimale Blattgröße kann auch kleiner sein: Kleinere Blätter erhöhen zwar den Rekursions- und Speicheraufwand, können durch günstigere Blattmultiplikationen dennoch schneller sein.

Vergleich mit binärem FastShift.
Ein benachbarter praktischer Vergleich mit einem binären FastShift-Verfahren ergab für gemischte \(256\)-Bit-BigInt-Einträge der Größe \( 1024\times1024 \) auf \(24\) Threads eine Laufzeit von etwa \( 24{,}84\ \mathrm{s}. \) Die \(3\)-äre shift-polynomiale FSR-Fassung berechnete dagegen das nahe benachbarte Produkt der Größe \( 1026\times1026 \) in etwa \( 5{,}30\ \mathrm{s}. \) Obwohl die Matrixgrößen nicht identisch sind, deutet dieser Vergleich in der betrachteten Einstellung auf einen praktischen Vorteil der \(3\)-ären Fassung hin. Ein wesentlicher Grund liegt darin, dass das binäre Packungsverfahren mit breiteren Zwischenzahlen und anschließender Rekombination arbeitet, während die \(3\)-äre FSR-Fassung direkt in die Zielblöcke akkumuliert.

Einordnung.
Die shift-polynomiale FSR-Multiplikation ist nicht als Ersatz für hochoptimierte Gleitkomma-BLAS-Multiplikation zu verstehen. Ihr natürlicher Anwendungsbereich liegt vielmehr bei exakter Ganzzahlarithmetik, rationaler Arithmetik und strukturierten Rekursionen, bei denen die Kosten der einzelnen arithmetischen Operationen hoch sind. Der Vorteil entsteht nicht aus einer bloßen Umformung der Matrixmultiplikation, sondern aus der Kombination von \[ \text{großer Matrixgröße} + \text{großer Ganzzahlen} + \text{hybrider }3\text{-ärer Rekursion} + \text{direkter In-Place-Akkumulation}. \] Besonders wichtig ist die Wahl der Blattgröße. Eine zu tiefe Rekursion, etwa mit sehr kleiner Blattgröße, erhöht den Verwaltungs- und Speicheraufwand. Eine zu flache Rekursion nutzt die Struktur dagegen nicht ausreichend aus. Die optimale Blattgröße ist daher ein algorithmischer Parameter und hängt von Matrixgröße, Bitlänge, Prozessorarchitektur und Parallelisierung ab.

Dreizweigige Projektorstruktur für lineare Gleichungssysteme

Grundblock.
Für \(a,c,d\in\mathbb K\), \(\mathbb K\in\{\mathbb R,\mathbb C\}\), sei \[ M(a,c,d) = \begin{pmatrix} \check{(a+d)}&0&\check{(a-d)}\\ 0&c&0\\ \check{(a-d)}&0&\check{(a+d)} \end{pmatrix}. \] Zu einem Vektor \(x=(x_1,x_2,x_3)^\top\) werden die lokalen Zweigkoordinaten \[ s=\check{(x_1+x_3)},\qquad m=x_2,\qquad q=\check{(x_1-x_3)} \] eingeführt. Die Rücktransformation lautet \[ x_1=s+q,\qquad x_2=m,\qquad x_3=s-q. \] Dann wirkt \(M(a,c,d)\) auf den drei Zweigen diagonal: \[ s\mapsto as,\qquad m\mapsto cm,\qquad q\mapsto dq. \] Insbesondere gilt \[ M(a,c,d)M(e,f,g)=M(ae,cf,dg), \] und \[ M(1,1,1)=I_3. \] Die Klasse ist daher additiv und multiplikativ abgeschlossen und bildet eine kommutative Algebra, die zu \(\mathbb K^3\) isomorph ist. Für \(a,c,d\in\mathbb K^*\) gilt \[ M(a,c,d)^{-1}=M(\tilde a,\tilde c,\tilde d). \]

Blocksysteme der Größe \(3m\).
Seien \[ A^{(s)},A^{(m)},A^{(q)}\in\mathbb K^{m \times m}. \] Aus diesen drei Zweigmatrizen wird eine Matrix \[ A_0\in\mathbb K^{3m\times3m} \] gebildet, deren \(3\times3\)-Block an Position \((i,j)\) durch \[ (A_0)_{ij} = M\!\left(a^{(s)}_{ij},a^{(m)}_{ij},a^{(q)}_{ij}\right) \] gegeben ist. Damit besitzt \(A_0\) eine dreizweigige Projektorstruktur. Für ein LGS \[ A_0x=b \] zerlegt sich \(x\) und \(b\) blockweise in die Zweige \[ s_j=\check{(x_{3j-2}+x_{3j})},\qquad m_j=x_{3j-1},\qquad q_j=\check{(x_{3j-2}-x_{3j})}, \] sowie entsprechend \[ b^{(s)}_j=\check{(b_{3j-2}+b_{3j})},\qquad b^{(m)}_j=b_{3j-1},\qquad b^{(q)}_j=\check{(b_{3j-2}-b_{3j})}. \] Dann zerfällt das ursprüngliche \(3m\times3m\)-System exakt in drei gewöhnliche \(m\times m\)-Systeme: \[ A^{(s)}s=b^{(s)},\qquad A^{(m)}m=b^{(m)},\qquad A^{(q)}q=b^{(q)}. \] Nach deren Lösung wird \(x\) durch \[ x_{3j-2}=s_j+q_j,\qquad x_{3j-1}=m_j,\qquad x_{3j}=s_j-q_j \] rekonstruiert. Damit gilt \[ A_0\text{ regulär} \quad\Longleftrightarrow\quad A^{(s)},A^{(m)},A^{(q)}\text{ regulär}. \]

Schur-Erweiterung.
Für Systemgrößen, die nicht vollständig durch die dreizweigige Struktur beschrieben werden, sei \[ A= \begin{pmatrix} A_0&R\\ R^\top&C \end{pmatrix}, \] wobei \[ A_0\in\mathbb K^{3m\times3m} \] die dreizweigige Projektorstruktur besitzt, \[ R\in\mathbb K^{3m\times r}, \qquad C\in\mathbb K^{r\times r} \] die Rest- und Randkopplungen beschreiben. Das LGS \[ \begin{pmatrix} A_0&R\\ R^\top&C \end{pmatrix} \begin{pmatrix} x\\y \end{pmatrix} = \begin{pmatrix} b_0\\b_1 \end{pmatrix} \] wird mit dem Schur-Komplement gelöst. Aus \[ A_0x+Ry=b_0 \] folgt \[ x=A_0^{-1}(b_0-Ry). \] Einsetzen in \[ R^\top x+Cy=b_1 \] ergibt \[ \left(C-R^\top A_0^{-1}R\right)y = b_1-R^\top A_0^{-1}b_0. \] Das Schur-Komplement lautet also \[ W=C-R^\top A_0^{-1}R. \] Nach Lösung des kleinen Systems \[ Wy=b_1-R^\top A_0^{-1}b_0 \] wird der strukturierte Anteil durch \[ x=A_0^{-1}(b_0-Ry) \] rekonstruiert. Dabei wird \(A_0^{-1}\) nicht explizit gebildet. Die Produkte \[ A_0^{-1}R,\qquad A_0^{-1}b_0 \] werden durch die drei Zweigsysteme \[ A^{(s)},\qquad A^{(m)},\qquad A^{(q)} \] berechnet.

SPD-Fall.
Ist \[ A^{(s)},A^{(m)},A^{(q)} \] symmetrisch positiv definit und ist auch das Schur-Komplement \[ W=C-R^\top A_0^{-1}R \] symmetrisch positiv definit, so kann die gesamte Rechnung durch Cholesky-Faktorisierungen erfolgen: \[ A^{(s)}=L_sL_s^\top,\qquad A^{(m)}=L_mL_m^\top,\qquad A^{(q)}=L_qL_q^\top, \] sowie \[ W=L_WL_W^\top. \] Dies ist numerisch günstiger als eine allgemeine dichte Faktorisierung.

Komplexität.
Eine direkte dichte Lösung eines Systems der Größe \( 3m+r \) besitzt kubische Ordnung \[ \mathcal O\!\left((3m+r)^3\right). \] Die strukturierte Schur-Variante benötigt im SPD-Fall drei Faktorisierungen der Größe \(m\times m\), die Bildung des Schur-Komplements und die Lösung des kleinen \(r\times r\)-Systems. Damit ergibt sich näherungsweise \[ \mathcal O\!\left( 3m^3 + 3m^2r + mr^2 + r^3 \right). \] Für festes kleines \(r\) dominiert \( \mathcal O(3m^3) \) statt \( \mathcal O(27m^3). \) Der ideale kubische Strukturvorteil beträgt daher näherungsweise \[ \frac{(3m)^3}{3m^3}=9. \] Die Kopplungsterme \[ A_0^{-1}R,\qquad R^\top A_0^{-1}R \] verringern diesen Idealwert praktisch, ohne den asymptotischen Hauptvorteil für \(r\ll m\) aufzuheben.

Numerisches Beispiel.
In einem Julia-Test wurde ein symmetrisch positiv definites System der Größe \( 3m+r \) erzeugt. Der strukturierte Kern \(A_0\) bestand aus \(3\times3\)-Blöcken der Form \(M(a,c,d)\), während \(R\) eine kleine Randkopplung und \(C\) einen kleinen positiv definiten Restblock beschrieb. Für \( m=100,\qquad r=5 \) ergab sich \( 3m+r=305. \) Die direkte Cholesky-Lösung benötigte ungefähr \( 0{,}001232\,\mathrm{s}, \) die strukturierte Schur-Cholesky-Lösung ungefähr \( 0{,}000285\,\mathrm{s}. \) Die Fehler lagen im Bereich \[ \|x_{\mathrm{Schur}}-x_{\mathrm{true}}\|_2 \approx 7{,}909\cdot10^{-15}, \] \[ \|Ax_{\mathrm{Schur}}-b\|_2 \approx 5{,}414\cdot10^{-14}. \] Für \( m=300,\qquad r=5 \) ergab sich \( 3m+r=905. \) Die direkte Cholesky-Lösung benötigte ungefähr \( 0{,}006438\,\mathrm{s}, \) die strukturierte Schur-Cholesky-Lösung ungefähr \( 0{,}004157\,\mathrm{s}. \) Für \( m=500,\qquad r=5 \) ergab sich \( 3m+r=1505. \) Die direkte Cholesky-Lösung benötigte ungefähr \( 0{,}052780\,\mathrm{s}, \) die strukturierte Schur-Cholesky-Lösung ungefähr \( 0{,}008407\,\mathrm{s}. \) Damit wurde im größten Test ein praktischer Beschleunigungsfaktor von ungefähr \[ \frac{0{,}052780}{0{,}008407}\approx 6{,}28 \] erzielt. Der theoretische kubische Vergleich liefert für diesen Fall \[ \frac{1505^3}{3\cdot500^3+5^3}\approx9{,}09. \] Die Abweichung vom Idealwert erklärt sich durch Kopplungskosten, Speicherzugriffe und Implementierungskonstanten.

Einordnung.
Die Methode löst keine beliebigen dichten LGS schneller. Ihr Vorteil entsteht genau dann, wenn der führende Kern eine dreizweigige Projektorstruktur besitzt oder nach geeigneter Umordnung beziehungsweise Transformation näherungsweise eine solche Struktur annimmt. In diesem Fall wird ein großes System nicht als monolithisches dichtes LGS behandelt, sondern in mehrere kleinere Zweigsysteme und ein kleines Schur-Komplement zerlegt: \[ \boxed{ 3m+r \quad\longrightarrow\quad 3\cdot m+r. } \] Bei rekursiver Schachtelung der Zweigsysteme kann diese Zerlegung weitergeführt werden. Für \(k\) dreizweigige Ebenen entsteht formal eine Zerlegung in \( 3^k \) unabhängige Zweige. Damit passt die Methode natürlich zu rekursiven strukturierten Verfahren wie FSR.

Es sei

eine polynomiale Abbildung vom Grad höchstens \(2\). Für alle \((x,y)\in\mathbb{C}^2\) gelte

Der Satz behandelt den quadratischen Spezialfall der Jacobi-Vermutung. Es wird gezeigt, dass \(H\) als polynomialer Automorphismus kein echtes Doppelpunktepaar besitzt.

Lemma zur Trennung infinitesimaler Ordnungen

Sei \(\delta:=\tilde{\omega}\ne0\) infinitesimal und seien \(a_r,\ldots,a_m\) endliche, von \(\delta\) unabhängige Größen. Gilt

so gilt

Beweis.

Angenommen, nicht alle Koeffizienten seien null. Sei \(j\) der kleinste Index mit \(a_j\ne0\). Dann folgt aus Gleichung (1)

Wegen \(\delta\ne0\) ergibt sich

Die rechte Seite ist infinitesimal. Der von \(\delta\) unabhängige Koeffizient \(a_j\) ist dagegen entweder null oder nichtinfinitesimal. Wegen \(a_j\ne0\) entsteht ein Widerspruch. Daher sind alle Koeffizienten null.\(\square\)

Bedeutung des Ordnungslemmas

Das Lemma besagt insbesondere, dass Terme verschiedener infinitesimaler Ordnung einander nicht ausgleichen können. Ein Term der Ordnung \(\delta^3\) kann keinen von null verschiedenen Term der Ordnung \(\delta^2\) aufheben. Allgemein gilt

nur dann, wenn zunächst \(a_r=0\) gilt. Danach ist dieselbe Argumentation auf die nächsthöhere Ordnung anzuwenden. Dies ist die infinitesimale Form des Vergleichs homogener Polynomgrade.

Grad null

Hat \(H\) den Grad null, so ist \(H\) konstant. Dann gilt \( JH=0 \) und damit \( \operatorname{eig} JH=0. \) Dies widerspricht \(e\ne0\). Eine Keller-Abbildung vom Grad null existiert daher nicht.\(\square\)

Grad eins

Hat \(H\) den Grad eins, so besitzt \(H\) die Form \( H(z)=c+Az \) mit \( z=(x, y)^\top, c\in\mathbb{C}^2 \) und \(A\in\mathbb{C}^{2\times2}\).

Da \( JH=A \) und \( \operatorname{eig} A=e\ne0 \) gelten, ist \(A\) mit der Umkehrabbildung \(H^{-1}(u)=A^{-1}(u-c)\) invertierbar. Somit ist \(H\) ein affiner Automorphismus und besitzt kein echtes Doppelpunktepaar.\(\square\)

Normierung des quadratischen Falls

Sei nun \(H\) vom Grad genau \(2\). Dann besitzt \(H\) eine eindeutige Zerlegung

wobei \(Q\) ein homogenes quadratisches Vektorpolynom ist.

Am Ursprung gilt \( JH(0)=A. \) Daher folgt die Invertierbarkeit von \(A\) aus der Keller-Bedingung \( \operatorname{eig} A=e\ne0. \) Wird

Aus \( JK(z)=A^{-1}JH(z) \) folgt

Trennung des linearen und quadratischen Eigenproduktanteils

Da \(R\) homogen quadratisch ist, ist \( M(z):=JR(z) \) eine homogene lineare Matrixfunktion. Es gilt \( JK(z)=I+M(z). \) Für eine \(2\times2\)-Matrix \(M\) gilt

Die Keller-Bedingung liefert daher

Nun wird \(z\) durch \(\delta z\) ersetzt. Die Homogenität liefert \( M(\delta z)=\delta M(z) \) und somit folgt

Der erste Summand besitzt die Ordnung \(\delta\), der zweite die Ordnung \(\delta^2\). Nach dem Ordnungslemma können sie einander nicht ausgleichen. Die Division durch \(\delta\ne0\) ergibt nacheinander \( \operatorname{tr}M(z)=0 \) und \( \operatorname{eig} M(z)=0. \) Da \(z\) beliebig war, gelten die Polynomidentitäten \( \operatorname{tr}JR=0 \) und \( \operatorname{eig} JR=0. \) Dies berücksichtigt ausdrücklich, dass der quadratische infinitesimale Term den linearen nicht ausgleichen kann.

Struktur des quadratischen Anteils

Schreibe \( R=(P,S), \) wobei \(P\) und \(S\) homogene quadratische Polynome sind.

Aus

folgt die Existenz eines homogenen kubischen Polynoms \(h\) mit

Damit gilt

Mit \( \operatorname{eig} JR=0 \) folgt daher

Schreibe

Dann gilt

Nach Einsetzung in Gleichung (12) und Koeffizientenvergleich erhält man

Ist \(a\ne0\), so wird \( r:=\tilde a b \) gesetzt. Dann folgt \( c=ar^2 \) und \( d=ar^3. \) Somit gilt

Ist \(a=0\), so folgt aus \(ac=b^2\) zunächst \(b=0\) und aus \(bd=c^2\) anschließend \(c=0\). Daher gilt

In beiden Fällen existieren eine lineare Form \( \ell(x,y)=\alpha x+\beta y \) und ein \(\kappa\in\mathbb{C}\) mit \( h=\kappa\ell^3. \) Aus Gleichung (11) folgt

Mit \( w:=3\kappa(\beta,-\alpha) \) gilt \( R(z)=w\ell(z)^2. \) Außerdem ist

Damit besitzt jede normierte quadratische Keller-Abbildung die Form

Direkter Doppelpunktausschluss

Angenommen, es existierten \(z_1\ne z_2\) mit

Wird die lineare Form \(\ell\) auf \(K(z)\) angewendet, so folgt wegen \(\ell(w)=0\)

Aus Gleichung (19) folgt daher

Somit stimmen die quadratischen Korrekturterme überein:

Daher ergibt sich \( z_1=z_2 \) aus

Dies widerspricht der Annahme \(z_1\ne z_2\). Folglich besitzt \(K\) kein echtes Doppelpunktepaar.\(\square\)

Explizite Umkehrabbildung

Die normierte Abbildung \( K(z)=z+w\ell(z)^2 \) besitzt die polynomiale Umkehrabbildung

Denn aus \(\ell(w)=0\) folgt

Daher gilt

Für die ursprüngliche Abbildung \( H(z)=c+A K(z) \) wird \( u:=A^{-1}(y-c) \) gesetzt. Dann ist die Umkehrabbildung \(H^{-1}\) polynomial wegen

Kanonische Normalform

Ist \(w=0\), so ist \(K\) die Identität. Sei daher \(w\ne0\). Wähle einen Vektor \(r\in\mathbb{C}^2\) mit \( \ell(r)=1. \) Wegen \(\ell(w)=0\) bilden \(w\) und \(r\) eine Basis von \(\mathbb{C}^2\). Mit \( z:=uw+vr \) gilt \( \ell(z)=v \) und daher

In den Koordinaten \((u,v)\) lautet die Abbildung

Bis auf affine Koordinatenänderungen besitzt jede quadratische Keller-Abbildung somit eine der beiden Formen

Die Keller-Klasse \(F=p(x+f(x,y))\) und \(G=q(y+f(x,y))\)

Es seien \( f\in\mathbb C[x,y] \) und \( p,q\in\mathbb C[t]. \) Betrachtet wird die polynomiale Abbildung

mit

und

Es werde vorausgesetzt, dass die Keller-Bedingung

mit einer Konstanten \( e\in\mathbb C^* \) für alle \((x,y)\in\mathbb C^2\) gilt.

Es wird gezeigt, dass die Keller-Bedingung die äußeren Polynome \(p\) und \(q\) zwingend auf affine Polynome reduziert und dass \(f\) vollständig bestimmt werden kann. Anschließend wird eine explizite polynomiale Umkehrabbildung angegeben.

Einführung der inneren Koordinaten

Setze

und

Die partiellen Ableitungen von \(u\) und \(v\) lauten

sowie

Berechnung der Jacobi-Matrix

Mit der Kettenregel folgt

und

Somit besitzt die Jacobi-Matrix die Form

Ihr Eigenprodukt ist

Damit lautet die grundlegende Eigenproduktidentität

Die Keller-Bedingung ist daher äquivalent zu

Einheitenargument im Polynomring

Der Polynomring \( \mathbb C[x,y] \) ist ein Integritätsbereich. Seine Einheiten sind genau die von null verschiedenen konstanten Polynome. Die rechte Seite von Gleichung (37) ist eine Einheit in \(\mathbb C[x,y]\). Folglich muss jeder der drei Faktoren auf der linken Seite ebenfalls eine Einheit sein. Es existieren daher Konstanten

mit

und

Dabei gilt notwendig

Nichtkonstanz der inneren Polynome

Das Polynom \(u=x+f(x,y)\) kann nicht konstant sein.

Angenommen, es gelte \( u=\gamma \) mit einer Konstanten \(\gamma\in\mathbb C\). Dann wäre \( F=p(\gamma) \) konstant. Folglich verschwände die erste Zeile der Jacobi-Matrix von \(H\), und es wäre \( \operatorname{eig} JH=0. \) Dies widerspricht der Voraussetzung \(e\ne0\). Somit ist \( u=x+f(x,y) \) nicht konstant. Entsprechend kann auch \( v=y+f(x,y) \) nicht konstant sein, da anderenfalls \(G=q(v)\) konstant wäre und wiederum \( \operatorname{eig} JH=0 \) folgte.

Die äußeren Polynome müssen affin sein

Aus \( {}^1p(u)=a \) und der Nichtkonstanz von \(u\) folgt, dass \({}^1p\) selbst konstant sein muss.

Tatsächlich ist der Einsetzungshomomorphismus

wegen der Nichtkonstanz von \(u\) injektiv. Aus

folgt daher die Polynomidentität

Mit einer Konstanten \(\alpha\in\mathbb C\) besitzt \(p\) somit die Form

Analog folgt mit einer Konstanten \(\beta\in\mathbb C\) aus \( {}^1q(v)=b \) die Darstellung

Da \(a,b\ne0\) gelten, sind \(p\) und \(q\) invertierbare affine Polynome.

Insbesondere schließt die Keller-Bedingung in dieser Klasse alle echt nichtlinearen äußeren Polynome aus:

Differentialgleichung für \(f\)

Aus Gleichung (41) folgt

Zur Lösung dieser polynomialen Differentialgleichung werden die linearen Koordinaten

und

eingeführt.

Dann gilt

Für die Differentialoperatoren erhält man

und

Daher ist

Gleichung (46) wird somit zu

Durch Integration nach \(s\) folgt

wobei \( h\in\mathbb C[t] \) ein beliebiges Polynom ist.

In den ursprünglichen Koordinaten lautet die allgemeine Lösung

Wegen

kann sie auch in der Form

geschrieben werden.

Vollständige Form der Keller-Abbildungen

Aus den Gleichungen (43), (44) und (53) folgt

und

Damit besitzt jede Keller-Abbildung der betrachteten Klasse die Form

und

wobei

und

gelten.

Umgekehrt erfüllt jede Abbildung der Form (58)–(60) die Keller-Bedingung

Beweis.

Aus Gleichung (53) folgt

und

Daher gilt

Außerdem ist \( {}^1p=a \) und \( {}^1q=b. \) Mit Gleichung (36) folgt

Damit ist auch die Rückrichtung bewiesen.\(\square\)

Darstellung in Summen- und Differenzkoordinaten

Setze erneut

Für die inneren Koordinaten

und

gilt zunächst

\\ Ferner folgt aus Gleichung (53)

Die innere Abbildung

besitzt somit in den Koordinaten \((r,s)\) die besonders einfache Dreiecksform

Die Differenzkoordinate \(r=x-y\) bleibt unverändert, während die Summenkoordinate \(s=x+y\) mit \(c\) skaliert und um ein Polynom in \(r\) verschoben wird.

Da \(c\ne0\), ist diese Dreiecksabbildung polynomial invertierbar.

Explizite Umkehrabbildung

Seien

die Bildkoordinaten.

Da

und

gelten, werden zunächst die affinen äußeren Abbildungen rückgängig gemacht:

und

Dann gilt

und

Durch Subtraktion folgt

Aus Gleichung (62) folgt

Daher ist

Mit

erhält man

Schließlich folgen aus Gleichung (49)

und

Somit lautet die explizite Umkehrabbildung

wobei \( u=\tilde a\,(U-\alpha) \) und \( v=\tilde b\,(V-\beta) \) gelten.

Da \(h\) polynomial ist und \(a,b,c\ne0\) gelten, ist auch \(H^{-1}\) polynomial.

Direkter Doppelpunktausschluss

Angenommen, es existierten zwei Punkte

mit

Da \(p\) und \(q\) affine Polynome mit von null verschiedenen linearen Koeffizienten sind, folgt aus Gleichung (71)

und

Durch Subtraktion der Gleichungen (72) und (73) erhält man

Somit stimmen die Differenzkoordinaten überein: \( r_1=r_2. \) Daher gilt auch

Durch Addition der Gleichungen (72) und (73) folgt

Mit Gleichung (75) ergibt sich

Wegen \(c\ne0\) folgt

Aus der Gleichheit der Summen- und Differenzkoordinaten folgt schließlich \( x_1=x_2 \) und \( y_1=y_2. \) Dies widerspricht der Annahme eines echten Doppelpunktepaares. Folglich ist \(H\) injektiv.\(\square\)

Der normierte Fall

Im normierten Keller-Fall \( e=1 \) gilt

Damit ist

Die allgemeine Form von \(f\) lautet

Sind zusätzlich die äußeren Polynome normiert durch \( p(t)=t \) und \( q(t)=t, \) so gilt \( a=b=1 \) und daher \( c=1. \) In diesem Fall reduziert sich \(f\) auf

Die Abbildung lautet dann

Ihre Umkehrabbildung ist

Tatsächlich bleibt \( F-G=x-y \) unverändert.

Scherungsdarstellung

Die innere Abbildung

besitzt in den Koordinaten

die Form

Sie ist somit eine Dreiecksabbildung, bestehend aus einer Skalierung der zweiten Koordinate und einer Scherung:

Die vollständige Abbildung \(H\) entsteht anschließend durch die affinen Transformationen

in den beiden Bildkomponenten.

Somit ist jede Keller-Abbildung dieser Klasse eine Verkettung aus

und

Zusammenfassung

Seien \( f\in\mathbb C[x,y] \) und \( p,q\in\mathbb C[t]. \) Definiere

und

Dann gilt

genau dann, wenn Konstanten

und ein Polynom \( h\in\mathbb C[t] \) existieren, sodass

und

gelten.

In diesem Fall besitzt \(H\) die explizite polynomiale Umkehrabbildung (70). Insbesondere ist \(H\) ein polynomialer Automorphismus und besitzt kein echtes Doppelpunktepaar.

Die Keller-Bedingung erzwingt innerhalb dieser umfassenden Klasse sowohl die Affinität der äußeren Polynome \(p\) und \(q\) als auch die vollständige Summen-Differenzen-Struktur des inneren Polynoms \(f\).\(\square\)

Bemerkung: Der Fall \(F(x,y) = x + f(x,y)\) und \(G(x,y) = y +p(F(x,y))\) ist einfacher zu beweisen.

Gesamtsatz

Jede zweikomponentige polynomiale Abbildung \( H:\mathbb{C}^2\longrightarrow\mathbb{C}^2 \) vom Grad höchstens \(2\) mit \( \operatorname{eig} JH=e\ne0 \) ist ein polynomialer Automorphismus. Genauer gilt:

Insbesondere existiert bis zum Grad \(2\) kein echtes Doppelpunktepaar.

Die infinitesimale Ordnungsrechnung liefert dabei im quadratischen Fall

Da ein Term der Ordnung \(\delta^2\) keinen von null verschiedenen Term der Ordnung \(\delta\) ausgleichen kann, folgt nacheinander \( \operatorname{tr}JR=0 \) und \( \operatorname{eig} JR=0. \) Die resultierende quadratische Abbildung ist eine Scherung und besitzt die explizite polynomiale Umkehrabbildung (25).\(\square\)

Reichweite des Beweises

Der vorstehende Beweis gilt vollständig für \( \deg H\le2. \) Für eine Abbildung höheren Grades trennt das Ordnungslemma zwar die verschiedenen homogenen Beiträge der Identität \( \operatorname{eig} JH=e, \) es beweist jedoch nicht allein, dass sämtliche höheren homogenen Bestandteile von \(H\) verschwinden. Eine Erweiterung auf beliebige Grade benötigt daher zusätzlich entweder eine zulässige Reduktion auf den quadratischen Fall oder einen gesonderten Ausschluss der höheren homogenen Bestandteile.\(\square\)