Ausgangspunkt ist der im Beitrag Parameteridentifikation untersuchte Integraloperator. Dort wird eine unbekannte Parameterfunktion aus einer Integralgleichung rekonstruiert. Nach der Diskretisierung entsteht eine dichte Matrix. Mit wachsender Zahl der Stützstellen steigen nicht nur Speicherbedarf und Rechenzeit. Zugleich wird die inverse Aufgabe zunehmend schlecht konditioniert. Ein feineres Gitter verbessert die Rekonstruktion daher nicht automatisch.
In diesem Beitrag dient derselbe Operator als überschaubares Beispiel für eine Niedrigrangapproximation ohne Aufbau der vollständigen Matrix und für einen schnelleren Algorithmus zur regularisierten Lösung. Sein Kern \(\exp(xs)\) besitzt eine besondere Struktur: Er lässt sich sehr genau durch eine kurze Summe separierbarer Terme approximieren.
Man könnte versuchen, die Berechnung mit der vollständig aufgebauten Matrix durch zusätzliche Rechenkerne zu beschleunigen. Der wesentlich größere Hebel liegt hier jedoch in der Wahl eines geeigneteren numerischen Verfahrens. Kurz gesagt: Algorithmus schlägt Hardware. Wir nutzen die Struktur des Kerns und berechnen die Wirkung des Operators, ohne zuvor sämtliche Matrixeinträge zu speichern.
Vom Integraloperator zur dichten Matrix
Gesucht wird im Modellproblem die Funktion \(\lambda(s)=\exp(s)\) aus der Integralgleichung
\[ (A\lambda)(x)=\int_0^1 \exp(xs)\lambda(s)\,\mathrm ds. \]
Wir verwenden die Stützstellen \(x_i=i/n\) und die Gewichte \(w_j\) der summierten Trapezregel. Diese Gewichte werden in der diagonalen Gewichtsmatrix \(W=\operatorname{diag}(w_0,\ldots,w_n)\) zusammengefasst. Damit entsteht die diskrete Operatormatrix
\[ K_{ij}=\exp(x_i x_j),\qquad W=\operatorname{diag}(w_0,\ldots,w_n),\qquad A=KW. \]
Die Multiplikation mit \(W\) von rechts gewichtet die \(j\)-te Spalte der Kernmatrix \(K\) mit dem zugehörigen Quadraturgewicht \(w_j\). Somit gilt für die Einträge der Operatormatrix \(A_{ij}=w_j\exp(x_i x_j)\).
Die Matrix besitzt mit \(N=n+1\) Stützstellen insgesamt \(N^2\) Einträge. Ihre Speicherung und auch ein Matrix-Vektor-Produkt benötigen daher quadratischen Aufwand.
Was bedeutet separierbar?
Ein einzelner separierbarer Term hat die Form \(p(x)q(s)\): Der erste Faktor hängt nur von \(x\), der zweite nur von \(s\) ab. Eine endliche Summe solcher Produkte heißt separierbare Darstellung:
\[ k_r(x,s)=\sum_{k=0}^{r-1}p_k(x)q_k(s). \]
Für den Exponentialkern liefert die Taylor-Reihe diese Darstellung unmittelbar:
\[ \exp(xs) =\sum_{k=0}^{\infty}\frac{(xs)^k}{k!} =\sum_{k=0}^{\infty} \underbrace{\frac{x^k}{\sqrt{k!}}}_{p_k(x)} \underbrace{\frac{s^k}{\sqrt{k!}}}_{q_k(s)}. \]
Schon die ersten Terme \(1+xs+x^2s^2/2+\ldots\) machen die Trennung sichtbar: In jedem Summanden lassen sich die \(x\)- und \(s\)-Anteile voneinander trennen. Nach dem Abschneiden bei \(r\) Summanden besteht die Näherung aus \(r\) solchen Produkten.
Wertet man den \(k\)-ten Term an allen Stützstellen aus, erhält man die Vektoren \(\boldsymbol p^{(k)}=(p_k(x_i))_{i=0}^n\) und \(\boldsymbol q^{(k)}=(q_k(x_j))_{j=0}^n\). Der zugehörige Anteil der diskreten Matrix lautet
\[ \bigl(M^{(k)}\bigr)_{ij}=p_k(x_i)q_k(x_j), \qquad M^{(k)}=\boldsymbol p^{(k)} \bigl(\boldsymbol q^{(k)}\bigr)^\mathsf T. \]
Jede Spalte von \(M^{(k)}\) ist damit ein skalares Vielfaches desselben Vektors \(\boldsymbol p^{(k)}\). Der Spaltenraum ist eindimensional und \(M^{(k)}\) hat – sofern der Term nicht verschwindet – Rang eins. Die Summe von \(r\) solchen Matrizen besitzt höchstens Rang \(r\). Genau deshalb ist das Abschneiden der Taylor-Reihe zugleich eine Niedrigrangapproximation.
Auf dem Gebiet \(0\leq x,s\leq1\) gilt für den Taylor-Rest die einfache Abschätzung
\[ \max_{x,s\in[0,1]} \left|\exp(xs)-k_r(x,s)\right| \leq \frac{\mathrm e}{r!}. \]
Der Fehler fällt hier also sehr schnell mit dem Rang. Für \(r=16\) beträgt die theoretische Schranke etwa \(1{,}3\cdot10^{-13}\). In den Testrechnungen lag der tatsächlich beobachtete maximale Kernfehler bei ungefähr \(5\cdot10^{-14}\).
Rechnen ohne vollständige Matrix
Wir bilden die voll besetzte Faktormatrix \(P\in\mathbb R^{N\times r}\) mit \(r\ll N\). Sie besitzt \(N\) Zeilen, aber nur \(r\) Spalten:
\[ P_{ik}=\frac{x_i^k}{\sqrt{k!}}, \qquad K\approx PP^\mathsf T, \qquad A\approx PP^\mathsf T W. \]
Für einen Vektor \(v\) wird der Operator dann in drei kleinen Schritten angewendet:
\[ Av\approx P\bigl(P^\mathsf T(Wv)\bigr). \]
Die volle Matrix \(A\) wird dabei nicht aufgebaut. Zuerst wird mit den Quadraturgewichten multipliziert, anschließend auf \(r\) Koeffizienten verdichtet und zuletzt wieder auf die \(N\) Stützstellen abgebildet.
| Darstellung | Speicher | Operatoranwendung |
|---|---|---|
| dichte Matrix | \(\mathcal O(N^2)\) | \(\mathcal O(N^2)\) |
| separierbare Faktoren | \(\mathcal O(Nr)\) | \(\mathcal O(Nr)\) |
Von der symmetrisierten Darstellung zum reduzierten Tikhonov-System
Ausgangspunkt ist das konsistent diskretisierte Tikhonov-Gleichungssystem aus dem Beitrag Parameteridentifikation, Abschnitt „Zur Regularisierung“. Dort werden die adjungierte Abbildung und die korrekte Berücksichtigung der Quadraturgewichte hergeleitet. In der hier verwendeten Schreibweise lautet das System
\[ (\beta I_N+A^2)\alpha_\beta=Ay, \qquad \alpha_\beta\in\mathbb R^N. \]
Dieses lineare Gleichungssystem bestimmt die \(N\) gesuchten Werte der regularisierten Parameterfunktion an den Stützstellen. Dabei bezeichnet \(A^2=AA\) das gewöhnliche Matrixprodukt und keine elementweise Potenz. \(I_N\) ist die \(N\times N\)-Einheitsmatrix. Der Parameter \(\beta>0\) stabilisiert die schlecht konditionierte inverse Aufgabe.
Symmetrisierung mit den Quadraturgewichten
Wegen der Quadraturgewichte ist \(A=KW\) im gewöhnlichen euklidischen Skalarprodukt nicht symmetrisch. Mit
\[ B=W^{1/2}AW^{-1/2}=W^{1/2}KW^{1/2}. \]
wird \(A\) durch eine Ähnlichkeitstransformation in die symmetrische Matrix \(B\) überführt. Beide Matrizen beschreiben denselben diskreten Operator in unterschiedlich skalierten Koordinaten. Mit \(z=W^{1/2}\alpha_\beta\) und \(b=W^{1/2}y\) lautet das transformierte Tikhonov-System
\[ (\beta I_N+B^2)z=Bb. \]
Niedrigrangreduktion mit der Singulärwertzerlegung
Nun wird die Niedrigrangapproximation \(K\approx PP^\mathsf T\) einbezogen. Mit der \(N\times r\)-Faktormatrix \(Q=W^{1/2}P\) folgt
\[ B\approx W^{1/2}PP^\mathsf T W^{1/2} =QQ^\mathsf T. \]
Die Singulärwertzerlegung \(Q=U\Sigma V^\mathsf T\) liefert damit
\[ B\approx U\Sigma^2U^\mathsf T, \qquad B^2\approx U\Sigma^4U^\mathsf T. \]
Die Spalten von \(U\in\mathbb R^{N\times r}\) bilden eine orthonormale Basis des verwendeten Unterraums. Deshalb wird die Lösung in der Form \(z\approx Uc\) mit \(c\in\mathbb R^r\) angesetzt. Setzt man diese Darstellung und die Approximation für \(B\) in das transformierte Tikhonov-System ein, erhält man wegen \(U^\mathsf TU=I_r\)
\[ U(\beta I_r+\Sigma^4)c =U\Sigma^2U^\mathsf T b. \]
Die Multiplikation von links mit \(U^\mathsf T\) entfernt den gemeinsamen Faktor \(U\). Damit folgt das reduzierte Tikhonov-System
\[ (\beta I_r+\Sigma^4)c=\Sigma^2U^\mathsf T b. \]
Dieses System bestimmt nur noch die \(r\) Koeffizienten von \(c\). Da \(\Sigma\) diagonal ist, kann jeder Koeffizient mit einem skalaren Filterfaktor berechnet werden:
\[ c_k=\frac{\sigma_k^2}{\beta+\sigma_k^4}(U^\mathsf T b)_k, \qquad \alpha_\beta\approx W^{-1/2}Uc. \]
Filterwirkung der Regularisierung
Die Filterwirkung wird im Vergleich mit der unregularisierten Inversion sichtbar. Für \(\sigma_k>0\) lässt sich der Faktor schreiben als
\[ \frac{\sigma_k^2}{\beta+\sigma_k^4} = \frac{1}{\sigma_k^2} \underbrace{\frac{\sigma_k^4}{\beta+\sigma_k^4}} _{d_\beta(\sigma_k)}. \]
Der erste Faktor \(1/\sigma_k^2\) gehört zur unregularisierten Inversion. Der zweite Faktor \(d_\beta(\sigma_k)\) liegt zwischen null und eins und dämpft diese Inversion. Für \(\sigma_k^4\ll\beta\) liegt er nahe null. Für \(\sigma_k^4\gg\beta\) liegt er nahe eins. Ist \(\sigma_k=0\), ergibt die ursprüngliche Formel \(0/\beta=0\). Ein Ausdruck \(0/0\) tritt wegen \(\beta>0\) nicht auf.
Rangabschneiden und Tikhonov-Regularisierung sind deshalb zwei verschiedene Schritte. Der Rang \(r\) legt fest, welche durch \(U\) beschriebenen Richtungen überhaupt im Ansatz vorkommen. Innerhalb dieses Unterraums dämpft \(\beta\) vor allem die Beiträge zu kleinen Singulärwerten.
Wie genau ist die Approximation?
Für \(n=128\) und \(\beta=10^{-6}\) zeigt sich, dass die Niedrigranglösung die dichte regularisierte Lösung ab etwa \(r=12\) praktisch reproduziert:
| Rang \(r\) | \(\max|K-K_r|\) | rel. Abweichung zur dichten Tikhonov-Lösung | max. Fehler zu \(\exp(x)\) |
|---|---|---|---|
| 8 | \(2{,}79\cdot10^{-5}\) | \(1{,}29\cdot10^{-4}\) | \(2{,}31\cdot10^{-3}\) |
| 12 | \(2{,}26\cdot10^{-9}\) | \(6{,}76\cdot10^{-9}\) | \(2{,}95\cdot10^{-3}\) |
| 16 | \(5{,}02\cdot10^{-14}\) | \(5{,}40\cdot10^{-10}\) | \(2{,}95\cdot10^{-3}\) |
Dass \(r=8\) hier zufällig einen etwas kleineren Fehler zur exakten Funktion liefert, bedeutet nicht, dass der Kern genauer approximiert wird. Das stärkere Rangabschneiden wirkt zusätzlich regularisierend. Ab \(r=12\) begrenzt nicht mehr die Kernapproximation die Rekonstruktion, sondern vor allem die Wahl von \(\beta\) und die Diskretisierung.
Im rauschfreien Referenzfall verbessert eine Verringerung von \(\beta=10^{-6}\) auf \(10^{-7}\) den maximalen Rekonstruktionsfehler. Beispielsweise sinkt er bei \(n=1024\) von \(2{,}06\cdot10^{-3}\) auf \(7{,}25\cdot10^{-4}\), bei \(n=8192\) von \(2{,}04\cdot10^{-3}\) auf \(7{,}09\cdot10^{-4}\). Bei einer absoluten Toleranz von \(10^{-3}\) liegen mit \(\beta=10^{-7}\) alle untersuchten Fälle ab \(n=256\) innerhalb der Toleranz.
Dieser Vergleich ist keine allgemeine Empfehlung für die Wahl von \(\beta\). Er zeigt nur, dass \(\beta=10^{-7}\) im analytisch erzeugten Referenzfall einen kleineren Fehler liefert. Die allgemeine Wahl des Regularisierungsparameters ist eine eigene Fragestellung und wird hier nicht weiter behandelt.
Was bringt das praktisch?
Die folgenden Messungen vergleichen den regularisierten Arbeitsablauf mit vollständig aufgebauter Matrix und den Niedrigrangalgorithmus bei \(r=16\) und \(\beta=10^{-6}\). Jeder Wert ist der Median aus fünf unabhängigen Einthreadläufen auf demselben Rechner.
| Implementierung | \(n\) | dicht | Niedrigrang | Faktor |
|---|---|---|---|---|
| MATLAB | 1024 | 52,930 ms | 0,256 ms | 207 |
| Python/NumPy | 1024 | 55,298 ms | 0,623 ms | 89 |
| MATLAB | 8192 | 18,852 s | 2,244 ms | 8402 |
| Python/NumPy | 8192 | 17,763 s | 3,459 ms | 5136 |
Bei \(n=8192\), also \(N=8193\), benötigt bereits eine einzelne dichte Matrix in doppelter Genauigkeit etwa 512 MiB. Die beiden voll besetzten Faktoren \(P\) und \(Q\) besitzen dagegen jeweils nur 16 Spalten und benötigen zusammen ungefähr 2 MiB.
Noch deutlicher wird der Verfahrensvorteil in einem bewusst ungleichen Vergleich. Der Arbeitsablauf mit vollständig aufgebauter Matrix darf alle 40 physischen Kerne des Rechners verwenden. Der Niedrigrangalgorithmus läuft dagegen nur auf einem Kern.
| Implementierung | vollständige Matrix 40 Threads |
Niedrigrang 1 Thread |
Faktor |
|---|---|---|---|
| MATLAB | 2,439 s | 2,244 ms | 1087 |
| Python/NumPy | 4,127 s | 3,459 ms | 1193 |
Selbst unter diesen Bedingungen ist der Niedrigrangalgorithmus in beiden Implementierungen mehr als tausendmal schneller. Zusätzliche Threads helfen ihm in diesem Fall kaum, weil die Operationen mit Faktormatrizen von nur 16 Spalten zu wenig Arbeit bieten, um viele Rechenkerne effizient auszulasten. Der Wechsel des Verfahrens bewirkt daher wesentlich mehr als die Parallelisierung der Rechnung mit vollständig aufgebauter Matrix.
Fazit
Die entscheidende Einsicht ist die Niedrigrangstruktur
\[ K\approx PP^\mathsf T. \]
Sie ermöglicht es, auf den Aufbau der vollständigen Matrix zu verzichten. Bei \(n=8192\) und \(r=16\) benötigen die Faktormatrizen ungefähr \(2\) MiB statt \(512\) MiB für eine einzelne vollständige Matrix. Der Niedrigrangalgorithmus ist mit einem Rechenkern mehr als tausendmal schneller als die Berechnung mit vollständiger Matrix und 40 Rechenkernen.
Die Ausnutzung der mathematischen Struktur bewirkt damit wesentlich mehr als zusätzliche Parallelisierung. Algorithmus schlägt Hardware.
Hinweis: Der Beitrag wurde mit Unterstützung KI-gestützter Werkzeuge ausgearbeitet. Mathematische Aussagen, Implementierungen und Messergebnisse wurden fachlich geprüft und durch reproduzierbare Rechnungen validiert.