Rank-one Riemannian Subspace Descent for Nonlinear Matrix Equations

Problemstellung
Die Arbeit befasst sich mit der rechnerischen Herausforderung, symmetrisch positive definite (SPD) Lösungen zu nichtlinearen Matrixgleichungen (NMEs) zu finden, die in der Regelungstechnik, dynamischen Programmierung und stochastischen Filterung auftreten. Diese Gleichungen, zu denen die kontinuierliche und diskrete algebraische Riccati-Gleichung (CARE/DARE) sowie Lyapunov-Gleichungen gehören, sind oft als das Finden von $X \in \mathcal{P}_n$ mit $G(X, A, C, \dots) = 0$ formuliert. Dies wird als ein nichtlineares Least-Squares-Problem umformuliert: $\min_{X \in \mathcal{P}_n} \|G(X, A, C, \dots)\|_F^2$.

Das primäre Hindernis ist die Skalierbarkeit. Standard-Solver für dichte Matrixgleichungen verursachen typischerweise $O(n^3)$ Gleitkommaoperationen (flops) pro Iteration. Während $O(n^2)$-Methoden existieren, beruhen diese im Allgemeinen auf Annahmen über die Sparsität oder den Niedrigrang der Lösungen, die für allgemeine dichte Probleme jedoch nicht zutreffen. Infolgedessen werden Standard-Solver für dichte Matrizen bei großen Zustandsdimensionen (z. B. $n=10.000$) aufgrund von Zeit- und Speicherbeschränkungen rechnerisch unpraktikabel.

Methodik: Rank-one Riemannian Subspace Descent (R1RSD)
Die Autoren schlagen einen Rank-one Riemannian Subspace Descent (R1RSD) Algorithmus vor, der direkt auf der Mannigfaltigkeit der SPD-Matrizen ($\mathcal{P}_n$) operiert, die mit einer affinen-invarianten Riemannschen Metrik ausgestattet ist.

1. Subraumselektion mittels Spektralzerlegung: Im Gegensatz zu traditionellen Subraum-Descent-Methoden, die Gradienten auf feste Koordinatensubräume projizieren, zerlegt R1RSD den transformierten Riemannschen Gradienten in eine Summe von Rang-eins-Eigenkomponenten. Der Algorithmus wählt die dominante Eigenkomponente (entsprechend dem größten Eigenwertbetrag), um die Abstiegrichtung zu definieren.
2. Potenzverfahren zur Richtungsidentifikation: Um den $O(n^3)$-Kosten einer vollständigen Eigenzerlegung zu entgehen, wird die dominante Komponente mittels des Potenzverfahrens identifiziert. Die Autoren zeigen, dass eine geringe Anzahl von Potenziterationen (speziell $r = O(\log n)$) ausreicht, um die dominante Richtung effektiv zu approximieren.
3. Rang-eins-Updates und Komplexität: Die Update-Regel wird unter Verwendung der Matrixexponentialfunktion einer Rang-eins-Matrix abgeleitet, welche ein geschlossenes Rang-eins-Update zulässt. Durch die Aufrechterhaltung des Cholesky-Faktors $B_t$ von $X_t$ und dessen Inversen mittels Rang-eins-Updates (unter Verwendung der Woodbury-Identität und cholupdate) vermeidet der Algorithmus dichte Matrixmultiplikationen.
4. Gradient-Vektor-Produkte: Für Zielfunktionen, die als Kompositionen Standard-Matrixoperationen (Addition, Multiplikation, Invertierung, Log-Determinante) definiert sind, wird das Gradient-Vektor-Produkt $\nabla f(X)v$ in $O(n^2)$ Zeit berechnet, indem der Rechengraph durchlaufen und Matrix-Vektor-Produkte angewendet werden, wodurch die explizite Bildung der vollen Gradientenmatrix vermieden wird.
5. Schrittweitenwahl: Die Rang-eins-Struktur ermöglicht in vielen Fällen eine exakte Line Search bei minimalen zusätzlichen Kosten, obwohl auch feste oder Armijo-Schrittweiten praktikabel sind.

Wesentliche Beiträge

• Algorithmische Innovation: Die Arbeit führt eine Methode ein, die Iterationen entlang einer einzelnen, sorgfältig ausgewählten Rang-eins-Richtung aktualisiert, die aus dem Riemannschen Gradienten abgeleitet ist, und sich damit von Methoden unterscheidet, die auf feste Subräume projizieren.
• Komplexitätsreduktion: Der vorgeschlagene Algorithmus erreicht einen Rechenaufwand von $O(n^2 \log n)$ pro Iteration (dominiert durch $r$ Potenziterationen der Ordnung $O(n^2)$) und benötigt unter Standardannahmen lediglich $O(n)$ Iterationen. Dies steht im Gegensatz zu State-of-the-Art Subraum-Descent-Methoden (z. B. [43]), die $O(n^2)$ pro Iteration, aber $O(n^2)$ Gesamtiterationen erfordern, sowie zu Standard-Solvern, die $O(n^3)$ pro Iteration benötigen.
• Theoretische Garantien: Die Autoren liefern eine Konvergenzanalyse sowohl für geodätisch glatte (nicht-konvexe) als auch für geodätisch stark konvexe Fälle. Sie etablieren eine $O(n)$-Iterationsbindung für den allgemeinen Fall und eine verbesserte Bindung für den stark konvexen Fall, mit Gesamtoprationszahlen von $O(n^3 \log n)$ bzw. $O(n^3 \kappa \log n)$.
• Praktische Implementierung: Die Arbeit enthält eine öffentlich verfügbare MATLAB-Implementierung, die Low-Rank-Updates nutzt, um die Effizienz zu gewährleisten, wodurch die $O(n^2)$-Komplexität pro Iteration in der Praxis realisiert wird.

Numerische Ergebnisse
Die Autoren vergleichen R1RSD mit State-of-the-Art-Solvern, einschließlich MATLABs icare, Structure-Preserving Doubling Algorithms (SDA) und der Burer-Monteiro-Faktorisierungs-basierten Koordinaten-Subraum-Descent-Methode (BMFC).

• Skalierbarkeit: Experimente an CARE, DARE und anderen NMEs mit Dimensionen bis zu $n=10.000$ zeigen, dass Standard-Solver (SDA, icare) und BMFC bei großen $n$ entweder keine Lösungen zurückgeben oder den Speicher überfordern. R1RSD löst diese Instanzen erfolgreich.
• Effizienz: Für $n=100$ konvergiert R1RSD in signifikant weniger Iterationen (1–2 Größenordnungen weniger) als BMFC.
• Robustheit: Der Algorithmus bleibt effektiv, selbst wenn die optimale Schrittweite stark variiert, und die gierige Subraumselektion (unter Verwendung von $\approx 10$ Potenziterationen) bietet ein starkes Gleichgewicht zwischen Konvergenzgeschwindigkeit und Aufwand pro Iteration.

Bedeutung und Behauptungen
Die Arbeit behauptet, dass R1RSD einen praktischen Ansatz für hochdimensionale und dichte SPD-beschränkte Matrixgleichungen bietet, bei denen traditionelle dichte Solver prohibitiv sind. Durch die Reduktion der Komplexität pro Iteration auf $O(n^2)$ bei gleichzeitiger Beibehaltung einer vernünftigen Iterationszahl schließt die Methode die Lücke zwischen theoretischer Skalierbarkeit und praktischer Anwendbarkeit für großskalige Regelungsprobleme. Die Autoren merken an, dass dieser Ansatz besonders relevant für zeitvariable und adaptive Regelszenarien ist (wie etwa State-Dependent Riccati Equation Frameworks), in denen das Lösen von Riccati-Gleichungen bei jedem Zeitschritt andernfalls unmöglich wäre. Die Arbeit betont, dass dieser Ansatz Rang-eins-Mannigfaltigkeits-Updates effektiv nutzen kann, um Probleme zu handhaben, die Standard-Solver nicht bewältigen können, ohne dabei auf Annahmen über die Sparsität oder den Niedrigrang der Lösungen oder Koeffizienten angewiesen zu sein.