Lineare Optimierung

Das Folgende setzt ebenfalls Mengenlehre, Topologie und Nichtstandardanalysis voraus.

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\) zusammenstellen 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^m\), \(m \in {}^{\nu}\mathbb{N}^*\), \(ß := {}_27\) und \(A \in {}^{\nu}\mathbb{C}^{n \times n}\) beweist die GR die Verkürzung der Laufzeit \(T(n) = \mathcal{O}(n^{ß})\)¹ 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 = mn\) bestimmen (s. Matrixmultiplikation).\(\square\)

Satz: Zu allen regulären \(C\;(= A^HA)\) mit \(n\) wie zuvor lässt sich \(C^{-1} \in {}^{\nu}\mathbb{C}^{n\times n}\) über das Schur-Komplement² und Einheitsmatrizen \(I\) multipliziert mit \(\varepsilon \in \mathcal{O}(\tilde{\nu})\) in der Laufzeit \(T(n) = \mathcal{O}(p^2)\) bestimmen, da \[ T(n) = \hat{T}(\check{n}) + \mathcal{O}(p^2) = 4T(\tilde{4}n) +\mathcal{O}(p^2 + \tilde{2}p^2) = \dots = \mathcal{O}(p^2) \] aufgrund der GR gilt.\(\square\)

Korollar: Jedes lösbare lineare Gleichungssystem \(Ax = b\) lässt sich für reguläres oder durch \(\varepsilon\) wie oben gestörtes singuläres \(A^HA\) mit einer Genauigkeit aus \(]0, 1[\) in der Laufzeit \(\mathcal{O}(p^2)\) lösen.\(\square\)

Satz: Das Intexverfahren löst jedes lösbare LP mit Inter-/Extrapolationen in \(\mathcal{O}(p^2)\).

Beweis und Algorithmus: Sei \(z := m + n\) und \(d \in [0, 1]\) die Dichte von \(A\). Zuerst werden \[ b^{T}y – c^{T}x \le 0,\quad Ax \le b,\quad A^{T}y \ge c \] normiert und skaliert. Es sei \[ P_r := \{(x, y)^T \in {}^{\nu}\mathbb{R}_{\ge 0}^{z} : b^{T}y – c^{T}x \le r \in [0, \rho],\; Ax – b \le r{\upharpoonright}_m,\; c – A^{T}y \le r{\upharpoonright}_n\} \] mit Radius \[\rho := s|\min \{b_1, …, b_m, -c_1, …, -c_n\}|\] und Skalierungsfaktor \(s \in [1, 2]\). Dann gilt \(0{\upharpoonright}_z \in \partial P_{\rho}\). Nach dem starken Dualitätssatz³ löst das LP \[ \min \{ r \in [0, \rho] : (x, y)^T \in P_r\} \] ebenso die LPs \[ \max \{{c}^{T}x : c \in {}^{\nu}\mathbb{R}^{n}, x \in {P}_{\ge 0}\} \quad\text{und}\quad \min \{{b}^{T}y : y \in {}^{\nu}\mathbb{R}_{\ge 0}^{m}, A^{T}y \ge c\}. \] Seine Lösung ist der geometrische Schwerpunkt \(g\) des Polytops \(P_0\). Mit \[ p_k^* := \min \check{p}_k + \max \check{p}_k,\quad k = 1, …, \grave{z}, \] wird \(g\) durch \(p_0 := (x_0, y_0, r_0)^T\) approximiert, bis \(||\Delta p||_1\) hinreichend klein ist. Die Lösung \(t^o(x^o, y^o, r^o)^T\) des zweidimensionalen LPs \[ \min \{ r \in [0, \rho] : t \in {}^{\nu}\mathbb{R}_{> 0},\; t(x_0, y_0)^T \in P_r\} \] approximiert \(g\) besser. Es wird \(r \le \check{\rho}\) erreicht und mit \(t^o(x^o, y^o)^T\) wiederholt, bis ggf. \(g \in P_0\) in \(\mathcal{O}({}_2\rho^2 dmn)\) berechnet ist. Das Lösen aller zweidimensionalen LPs \(\min_k r_k\) durch Bisektion für \(r_k \in {}^{\nu}\mathbb{R}_{\ge 0}\), \(k = 1, …, z\), in jeweils \(\mathcal{O}(p^2)\) ermittelt \(q \in {}^{\nu}\mathbb{R}^k\) mit \( q_k := \Delta p_k \Delta r_k/r\) und \(r := \min_k \Delta r_k. \) Vereinfacht sei \(|\Delta p_1| = … = |\Delta p_{z}|\). Hierbei wäre \(\min r_{\grave{z}}\) für \(p^* := p + wq\) mit \(w \in {}^{\nu}\mathbb{R}_{\ge 0}\) ebenso zu lösen. Folgt \(\min_k \Delta r_k\, r = 0\), wird die Berechnung gestoppt, andernfalls wiederholt, bis \(\min r = 0\) oder \(\min r > 0\) feststeht.\(\square\)

Bemerkung: Die Laufzeit für \(m, n \in {}^{\nu}\mathbb{N}^*\) ist \(\mathcal{O}(dmn)\). Hierbei ist \(d\) die Dichte der Matrix \(A\). Absichtliche Rundungsfehler in \(b\) und \(c\) können eine Lösung eindeutig machen. Wird das Verfahren für verteiltes Rechnen in \({}^{\nu}\mathbb{R}^{\nu}\) optimiert, beträgt sie nur \(\mathcal{O}(1)\). Gilt \(\acute{v}_0 \le \lfloor v_j \rfloor – v_j\) für ganzzahlige Lösungen \(v_j\), bewährt es sich gut für (gemischt-)ganzzahlige Probleme und die (nicht-)konvexe (Pareto-)Optimierung.

Bemerkung: Das Verfahren nutzt das Zusammenfallen des berechneten geometrischen Schwerpunkts und des Inkugelmittelpunkts in \(P_0\). Das (paarweise) Drehen der Koordinatenachsen des ersteren wie von \(c\) durch eine einfache Rotationsmatrix sichert dies in vergleichbarer Laufzeit ab. Falls erforderlich schwächt ein gleicher kleiner Betrag die Restriktionen vorübergehend ab. Zahlen der Länge \(\mathcal{O}(\nu)\) lassen sich nur in \(\mathcal{O}(p)\) abarbeiten. Auch das Folgende lässt sich einfach ins Komplexe überführen:

Folgerung: Das LP \( \max \{{||x – {x}^{o}||}_1 : {c}^{T}x = {c}^{T}{x}^{o},\; Ax \le b,\; x – {x}^{o} \in {[-1, 1]}^{n},\; x \in {}^{\nu}\mathbb{R}_{\ge 0}^{n}\} \) kann zur ersten Lösung \({x}^{o}\) ggf. eine zweite in \(\mathcal{O}(p^2)\) bestimmen, wobei \({y}^{o}\) analog behandelt werden kann.\(\square\)

Folgerung: Das LP \( \max \{\nu|\lambda_j| + ||x_j||_1 : Ax_j = \lambda_j x_j,\; |\lambda_j| \in [0, r_j],\; x_j \in {[-1, 1]}^{n*}\} \) kann für jedes \(j = 1, …, n\) den Eigenwert \(\lambda_j \in {}^{\nu}\mathbb{R}\) und den Eigenvektor \(x_j \in {}^{\nu}\mathbb{R}^{n}\) der Matrix \(A \in {}^{\nu}\mathbb{R}^{n \times n}\) in \(\mathcal{O}(p^2)\) bestimmen.\(\square\)

Folgerung: Die LPs \(\max \{x_{j} : Ax = 0\}\) liefern alle Lösungen von \(Ax = b\) in \(\mathcal{O}(p^2)\). Die Matrix \(A\) ist genau dann regulär, wenn das LP \(\max \{||x||_{1} : Ax = 0\} = 0\) ist.\(\square\)

Folgerung: Sei \(\alpha_{j} := A_{.j}^{-1}\) für \(j = 1, …, n\) zur Matrix \(A^{-1} \in {}^{\nu}\mathbb{R}^{n \times n}\) und sei \(\delta_{ij}\) das Kronecker-Delta. Für reguläres \(A\) ist das Eigenprodukt \(\neq 0\) und jedes \( A \alpha_{j} = (\delta_{1j}, …, \delta_{nj})^{T} \) lässt sich in \(\mathcal{O}(p^2)\) lösen.\(\square\)

Folgerung: Das Polynom \(c_{\acute{n}}x^{\acute{n}} + \dots + c_1x + c_0\) lässt sich in \(\mathcal{O}(p^2)\) mit Messungen \((x_j, y_j) \in {}^{\nu}\mathbb{R}^2\) und Residuen \(r_j \in {}^{\nu}\mathbb{R}\) für \(j = 1, …, m\) sowie \((c, x)^T \in {}^{\nu}\mathbb{R}^{\grave{n}}\) aus \[ y = r + Xc,\quad X^Tr = 0,\quad X = (x_j^k)_{j,k},\ k = 0, …, \acute{n} \] bestimmen (auch unter zusätzlichen Nebenbedingungen).\(\square\)

Folgerung: Die diskrete Fouriertransformation (s. Wissenschaftliches Rechnen) kann die Padé-Approximation \(\frac{a(x)}{b(x)} = {\LARGE{\textbf{+}}}_{m=0}^n a_m x^m / {\LARGE{\textbf{+}}}_{m=0}^n b_m x^m\) von \(f(x)\) für \(x \in {}^{\nu}\mathbb{R}\) in \(\mathcal{O}(p^2)\) bestimmen.\(\square\)

Korollar: Intex- und zweidimensionale Bisektions- oder Newtonverfahren lösen die meisten lösbaren konvexen Programme \[ \min \{f_{1}(x) : x \in {}^{\nu}\mathbb{R}^{n},\ (f_{2}(x), …, f_{m}(x))^{T} \le 0\} \] mit konvexen Funktionen \(f_{k} \in {}^{\nu}\mathbb{R}\) für \(k = 1, …, m\) in polynomieller Laufzeit, wenn die Anzahl der Operanden \(x_{j}\) der \(f_{k}\) \(\le {\nu}^{\nu-3}\) ist und ein \(x\) existiert mit \(f_{k}(x) < 0\) für alle \(k > 1.\square\)

Korollar: Das LP \( \max \{x : y_{\acute{m}} – xy_m = b_{\acute{m}},\; y_n = y_0 = 0,\; x \le x_0 \in {}^{\nu}\mathbb{R}\} \) kann für \(m = 1, …, n\) nach dem Horner-Schema ein \(b_{\acute{n}}x^{\acute{n}} + \dots + b_1x + b_0 = 0\) lösendes \(x \in {}^{\nu}\mathbb{R}\) in \(\mathcal{O}(p^2)\) bestimmen.\(\square\)

Das Intexverfahren löst die meisten lösbaren LPs in \(\mathcal{O}(p^2)\), indem es Inter- und Extrapolationen mit Schwerpunktprojektionen verbindet. Es nutzt die Geometrie des zulässigen Polyeders und seiner Schnitte in den Projektionen und arbeitet unabhängig von Pivotregeln. In Situationen, in denen zusätzliche algebraische Struktur vorliegt (Toeplitz, Hankel, blockzirkulant oder allgemein faltungsähnlich), ist eine geometrische Suche zwar möglich, aber typischerweise weniger effizient als eine direkte strukturierte Faktorisierung. Für diese Fälle erweisen sich die FTD und die FSR-Verfahren als komplementär zum Intexverfahren.

Die FTD nutzt die Beobachtung, dass eine reguläre Matrix \(A\in{}^{\nu}\mathbb{R}^{n\times n}\), deren Einträge entlang der Diagonalen nur langsam variieren oder durch Faltungsoperatoren beschreibbar sind, sich gut als Produkt zweier Toeplitz-Matrizen \(T_{1},T_{2}\) approximieren lässt. Zu Beginn wird \(x_{0}=\max|A_{ij}|\) gewählt; anschließend werden die Folgen \(x_{k},y_{k}\) durch lokale arithmetische Mittelwerte benachbarter Einträge gebildet. Die diagonalen Summen \[ z_{1}=x_{0}y_{0}+x_{1}y_{1}+\ldots,\qquad z_{2}=x_{0}y_{-1}+x_{1}y_{0}+\ldots \] werden rekursiv berechnet, bis \(A\approx T_{1}T_{2}\) erreicht ist. Die Normierung der Spalten von \(T_{1}\) und \(T_{2}\) stabilisiert die Rekursion. Durch Rückwärtseinsetzen entlang der Diagonalen wird eine eindeutige Approximation erreicht. Die Anzahl der Schritte pro Diagonale beträgt \(\mathcal{O}(m)\). Werden anschließend \(T_{1}^{-1}\) und \(T_{2}^{-1}\) berechnet, folgt wegen der rein strukturierten Operationen \[ T(n)=\mathcal{O}(p^{2}),\qquad p=mn. \] Für symmetrische Systeme entsteht aus \(H=T_{1}^{T}T_{1}\) die Faktorisierung \(H\approx P^{T}DP\) mit diagonaler Matrix \(D\) und block-Toeplitz-Matrix \(P\), die im LUX-Verfahren (Linear Universal Exact) zur Lösung strukturierter LPs verwendet wird. Sie ersetzt klassische Cholesky- oder QR-Verfahren und erhält Toeplitz-, Hankel- und Zirkulantstrukturen exakt.

Parallel hierzu beruhen die FSR-Verfahren auf rekursiven Verschiebungen und Modulo-Rekonstruktionen. Eine Matrix \(A\) wird in niedrig- und hochgewichtige Anteile zerlegt, \[ A=(A_{\mathrm{lo}}\ll s)+A_{\mathrm{hi}}, \] wobei \(\ll s\) eine generative Verschiebung bedeutet. Die Produkte der Einzelanteile werden rekursiv berechnet und anschließend modulo eines gut gewählten \(s\) zu einem exakten Ergebnis zusammengesetzt. Dadurch bleiben alle Zwischenwerte dynamisch klein, die Rekursion ist vollständig parallelisierbar, und BigInt- bzw. BigFloat-Arithmetik führen nicht zu Überlaufproblemen. Auf diese Weise gilt für Matrixmultiplikation, Inversion und Schur-Komplement-Schritte dieselbe Komplexität wie bei der FTD, \[ T_{\mathrm{FSR}}(n)=\mathcal{O}(p^{2}). \] Die drei Verfahren ergänzen einander in natürlicher Weise. Intex löst LPs durch Reduktion auf die Bestimmung eines Schwerpunkts im zulässigen Polyeder und eignet sich besonders für Probleme mit schwacher oder fehlender algebraischer Struktur. FTD löst strukturierte lineare Gleichungssysteme und symmetrische Faktorisierungen exakt und benötigt nur Faltungsoperatoren sowie diagonale Aktualisierungen. FSR führt die arithmetischen Operationen rekursiv und speichersparend aus und eignet sich für sehr große Ganzzahlen und hohe Präzision. Liegt in einem LP oder einem LGS eine Toeplitz-, Hankel- oder allgemeine Faltungsstruktur vor, so reduziert die Verkettung
\(A\approx T_{1}T_{2},\qquad T_{1}^{-1},T_{2}^{-1}\)    rekursiv durch FSR,     Projektion oder LP-Schnitt durch Intex,
den Gesamtaufwand auf \(\mathcal{O}(p^{2})\) für die strukturierten LGS und die LPs. Dadurch lassen sich auch große Probleme in den generativen Zahlräumen \({}^{\nu}\mathbb{R}^{n\times n}\) stabil und exakt lösen. Die Verfahren arbeiten unabhängig von der Wortlänge der Einträge und skalieren für Zahlen vom Umfang \(\mathcal{O}(\nu)\) nahezu linear.

Die Kombination von FTD, FSR und Intex bildet so einen universellen Werkzeugkasten für strukturierte lineare Algebra und Optimierung. Sie verbindet Struktur (Toeplitz, Hankel, Zirkulant), Rekursion (Verschiebung, Faltung) und Geometrie (Inter-/Extrapolation) zu einem Verfahren, das sowohl Regularität als auch exakte Lösungen in kurzer Laufzeit gewährleistet. Damit bietet sie einen durchgängigen Zugang zu Problemen, die bislang nur getrennt über numerische oder geometrische Ansätze gelöst wurden.

Die folgenden Größenordnungen fassen die Schnittpunkte zwischen klassischen und neuen Verfahren sowie die dadurch neu berechenbaren Bereiche zusammen. Die genannten Schwellen sollen heuristisch zu verstehen sein und hängen von Implementierung, Hardware und gewünschter Genauigkeit ab; sie markieren jedoch die Bereiche, in denen sich die neuen Verfahren systematisch gegenüber den klassischen Ansätzen durchsetzen könnten.

Programmtext des Simplexverfahrens

© 2008-2025 by Boris Haase

Seitenbeginn

1. s. Golub, Gene H.; van Loan, Charles F.: Matrix Computations; 3rd Ed.; 1996; Johns Hopkins University Press; Baltimore, S. 31 ff.
2. s. Knabner, Peter; Barth, Wolf: Lineare Algebra; 2., überarb. und erw. Aufl.; 2018; Springer; Berlin, S. 223
3. Vanderbei, Robert J.: Linear Programming; 3rd Ed.; 2008; Springer; New York., S. 60 – 65