Scalable Gaussian Processes for Integrated and Overlapping Measurements Via Augmented State Space Models

Technische Zusammenfassung: Skalierbare Gaußsche Prozesse für integrierte und überlappende Messungen mittels erweiterter Zustandsraummodelle

Problemstellung
Astronomische Messungen werden häufig durch die Integration von Signalen über endliche Belichtungszeiten gewonnen (z. B. die Akkumulation von Photoelektronen in einem CCD). Wenn die Zeitskala der zugrunde liegenden astrophysikalischen Variabilität mit der Belichtungsdauer vergleichbar ist, werden Standard-Gaußsche-Prozess-Modelle (GP) für instantane Werte zu verzerrten Schätzern. Die korrekte Berücksichtigung dieser Integration erfordert die Modellierung der Messung als Mittelwert des latenten Prozesses über das Belichtungsintervall.

Während die GP-Regression ein Standardwerkzeug zur Modellierung stochastischer Phänomene ist (z. B. Sternoszillationen in Radialgeschwindigkeitsdaten), stellt die Handhabung integrierter Messungen erhebliche rechnerische Herausforderungen dar. Für nicht überlappende, konstante Belichtungen können spezialisierte Kernel (wie die in celerite) eine Quasi-Trennbarkeit aufrechterhalten, was $O(N)$-Lösungen ermöglicht. Im allgemeinen Fall mit variablen Belichtungszeiten und, entscheidend, überlappenden Belichtungen (häufig bei der Kombination von Daten mehrerer Instrumente mit unterschiedlichen Abtastraten) verliert die Kovarianzmatrix jedoch ihre quasi-trennbare Struktur. Dies erzwingt die Verwendung dichter Matrixoperationen, was zu Speicherbedarf von $O(N^2)$ und Laufzeitkosten von $O(N^3)$ führt. Bestehende Methoden für integrierte Kernel (z. B. Luhn et al. 2026) beruhen auf der Konstruktion der vollen dichten Kovarianzmatrix, was für große Datensätze ($N \sim 10^5$), wie sie für moderne Multi-Instrumenten-Übersichtssurveys typisch sind, rechnerisch prohibitiv wird.

Methodik
Die Autoren schlagen vor, die Äquivalenz zwischen Gaußschen Prozessen und linearen Gaußschen Zustandsraummodellen (SSM) zu erweitern, um integrierte und überlappende Beobachtungen unter Beibehaltung der Skalierbarkeit zu behandeln.

1. Erweiterte Zustandsraumformulierung:
   Die Kerninnovation besteht in der Erweiterung des Standard-Zustandsvektors des SSM. Für einen latenten Prozess $x(t)$, der durch eine lineare stochastische Differentialgleichung (SDE) gesteuert wird, führen die Autoren einen Integralzustand $z(t)$ ein, sodass $dz/dt = x(t)$ gilt.

   • Der erweiterte Zustandsvektor wird zu $\tilde{x} = [x, z]^T$.
   • Die Systemdynamik wird aktualisiert, um die Integration des latenten Zustands einzubeziehen.
   • Reset-Mechanismus: Zu Beginn jeder Belichtung $t_s$ wird der Integralzustand $z$ auf Null zurückgesetzt. Am Ende der Belichtung $t_e$ wird die Beobachtung $y$ als normalisiertes Integral des Zustands über das Intervall modelliert: $y = z(t_e)/\delta + \epsilon$, wobei $\delta$ die Belichtungsdauer ist.

2. Handhabung von Überlappungen:
   Um überlappende Belichtungen mehrerer Instrumente zu berücksichtigen, weist das Framework jedem Instrument $i$ einen eigenen Integralzustand $z_i$ zu. Der erweiterte Zustandsvektor wird erweitert, um alle instrumentenspezifischen Integrale einzuschließen. Die Reset-Matrizen werden selektiv angewendet: $z_i$ wird nur zu Beginn einer Belichtung zurückgesetzt, die zu Instrument $i$ gehört, was die gleichzeitige Akkumulation von Integralen für überlappende Instrumente ermöglicht. Diese Konstruktion kodiert die für überlappende Daten erforderlichen Kreuzkorrelationen auf natürliche Weise, ohne dass dichte Matrixoperationen erforderlich sind.

3. Skalierbare Inferenz:
   Die Autoren leiten modifizierte Kalman-Filter- und Rauch–Tung–Striebel (RTS)-Glättungs-Algorithmen für dieses erweiterte System ab.

   • Sequenzielle Komplexität: Die Filter- und Glättungsschritte operieren auf der Dimension des erweiterten Zustands ($d + N_{inst}$) und halten eine rechnerische Komplexität von $O(N)$ auf einer CPU aufrecht.
   • Parallele Komplexität: Unter Ausnutzung der Assoziativität der Filter- und Glättungsverteilungen wird die Methode als Problem aller Präfixsummen neu formuliert. Dies ermöglicht eine parallele Implementierung auf GPUs und reduziert die Wandzeit auf $O(N/T + \log T)$, wobei $T$ die Anzahl der parallelen Worker ist.

4. Software-Implementierung:
   Die Autoren stellen smolgp vor, ein Open-Source-Python/JAX-Paket. Es bietet Drop-in-Kompatibilität mit der bestehenden tinygp-Bibliothek und unterstützt sowohl Standard- als auch belichtungsbewusste GP-Modellierung. Es implementiert die oben beschriebenen sequenziellen und parallelen Löser.

Hauptergebnisse

• Numerische Äquivalenz: Die Autoren validieren den erweiterten SSM-Ansatz gegen die vollständige Implementierung dichter GPs (Luhn et al. 2026) für Szenarien mit einem Instrument sowie mit mehreren Instrumenten (überlappend). Die Ergebnisse für bedingte Mittelwerte, Varianzen und Log-Likelihoods stimmen bis auf Maschinengenauigkeit überein (Unterschiede $\sim 10^{-14}$ bis $10^{-15}$), was die theoretische Äquivalenz bestätigt.
• Leistungsbenchmarks:
  • Speicher: Der SSM-Ansatz ist deutlich speichereffizienter als die Methode dichter GPs. Während dichte GPs $O(N^2)$ Speicher benötigen (für $N \sim 10^5$ über 100 GB), skaliert die SSM-Methode linear, was große Datensätze auf Standardhardware machbar macht.
  • Laufzeit: Für integrierte Daten, für die keine quasi-trennbare (QSM) Alternative existiert, bietet die SSM-Methode eine transformative Beschleunigung. Auf einer CPU ist der sequenzielle SSM-Löser für Vorhersagen am schnellsten und für die Likelihood-Berechnung konkurrenzfähig. Auf einer GPU erreicht der parallele SSM-Löser die schnellsten Konditionszeiten und konkurriert bei Likelihoods mit QSM-Methoden.
  • Vorhersage: Die Methode unterstützt effiziente Vorhersagen an beliebigen Testpunkten mittels eines spezifischen Algorithmus (Algorithmus 1), der im Vergleich zu dichten Matrixmultiplikationen embarrassingly parallel und speicherschonend ist.

Bedeutung und Behauptungen
Die Arbeit behauptet, dass diese Forschung massive Multi-Instrumenten-Kreuzvergleiche ermöglicht, bei denen Belichtungsüberlappungen allgegenwärtig sind, ein Szenario, das zuvor durch Rechenkosten begrenzt war. Durch die Aufrechterhaltung einer Komplexität von $O(N)$ (oder $O(\log N)$ auf GPUs) erschließt die Methode das Potenzial für:

• Die Analyse von Datensätzen mit $N > 10^5$ Beobachtungen (z. B. tägliche Sonnenüberwachung durch mehrere Instrumente wie KPF, NEID, EXPRES und HARPS-N).
• Den Einsatz komplexerer Modelle und höherdimensionaler Datensätze, die zuvor aufgrund der $O(N^3)$-Flaschenhals-Problematik integrierter GPs unlösbar waren.
• Die Erweiterung skalierbarer GP-Modellierung auf Kernel, die keine Quasi-Trennbarkeit aufweisen (z. B. quasiperiodische Kernel), wenn integrierte Messungen involviert sind.

Die Autoren betonen, dass das Framework nicht auf Zeitreihen beschränkt ist, sondern auf alle gebinnten Daten mit einer sortierbaren unabhängigen Koordinate anwendbar ist, und dass es aufgrund der Einbeziehung von Zustandsableitungen in die SSM-Formulierung natürliche Unterstützung für Ableitungsbeobachtungen (z. B. FF'-Modelle) bietet. Die Arbeit positioniert smolgp als skalierbares, flexibles Werkzeug für die Analyse astronomischer Zeitreihen der nächsten Generation.