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

Technische Samenvatting: Schaalbare Gaussische Processen voor Geïntegreerde en Overlappende Metingen via Geaugmenteerde Toestandsruimtemodellen

Probleemstelling
Astronomische metingen worden vaak verkregen door signalen te integreren over eindige blootstellingstijden (bijvoorbeeld het accumuleren van foto-elektronen in een CCD). Wanneer de tijdschaal van de onderliggende astrofysische variabiliteit vergelijkbaar is met de blootstellingsduur, worden standaard Gaussische Proces (GP) modellen voor momentane metingen vertekende schatters. Het correcte rekening houden met deze integratie vereist dat de meting wordt gemodelleerd als een gemiddelde van het latente proces over het blootstellingsinterval.

Hoewel GP-regressie een standaardinstrument is voor het modelleren van stochastische fenomenen (bijvoorbeeld steroscillaties in radiale snelheidsdata), vormt het verwerken van geïntegreerde metingen aanzienlijke computationele uitdagingen. Voor niet-overlappende metingen met constante blootstellingstijden kunnen gespecialiseerde kernels (zoals die in celerite) een quasi-scheidbare structuur behouden, wat oplossingen van $O(N)$ mogelijk maakt. Echter, in het algemene geval met variabele blootstellingstijden en, cruciaal, overlappende blootstellingen (veelvoorkomend bij het combineren van data van meerdere instrumenten met verschillende cadans), verliest de covariantiematrix zijn quasi-scheidbare structuur. Dit dwingt tot het gebruik van dichte matrixoperaties, wat resulteert in geheugenkosten van $O(N^2)$ en rekentijden van $O(N^3)$. Bestaande methoden voor geïntegreerde kernels (bijvoorbeeld Luhn et al. 2026) vertrouwen op het construeren van de volledige dichte covariantiematrix, wat computationeel onuitvoerbaar wordt voor grote datasets ($N \sim 10^5$) die typerend zijn voor moderne multi-instrument surveys.

Methodologie
De auteurs stellen voor de equivalentie tussen Gaussische Processen en Lineaire Gaussische Toestandsruimtemodellen (SSM's) uit te breiden om geïntegreerde en overlappende observaties te verwerken, terwijl de schaalbaarheid behouden blijft.

1. Geaugmenteerde Toestandsruimteformulering:
   De kerninnovatie is de uitbreiding van de standaard SSM-toestandsvector. Voor een latente proces $x(t)$ dat wordt beheerst door een lineaire stochastische differentiaalvergelijking (SDE), introduceren de auteurs een integraaltoestand $z(t)$ zodanig dat $dz/dt = x(t)$.

   • De geaugmenteerde toestandsvector wordt $\tilde{x} = [x, z]^T$.
   • De systeemdynamica worden bijgewerkt om de integratie van de latente toestand te omvatten.
   • Resetmechanisme: Aan het begin van elke blootstelling $t_s$ wordt de integraaltoestand $z$ gereset naar nul. Aan het einde van de blootstelling $t_e$ wordt de observatie $y$ gemodelleerd als de genormaliseerde integraal van de toestand over het interval: $y = z(t_e)/\delta + \epsilon$, waarbij $\delta$ de blootstellingsduur is.

2. Verwerking van Overlappingen:
   Om overlappende blootstellingen van meerdere instrumenten te accommoderen, kent het kader een distincte integraaltoestand $z_i$ toe aan elk instrument $i$. De geaugmenteerde toestandsvector wordt uitgebreid om alle instrument-specifieke integralen te bevatten. De resetmatrices worden selectief toegepast: $z_i$ wordt alleen gereset aan het begin van een blootstelling die behoort tot instrument $i$, wat gelijktijdige accumulatie van integralen voor overlappende instrumenten mogelijk maakt. Deze constructie codeert de voor overlappende data vereiste kruiscorrelaties op een natuurlijke manier zonder dat dichte matrixoperaties nodig zijn.

3. Schaalbare Inference:
   De auteurs leiden gemodificeerde Kalman-filter en Rauch–Tung–Striebel (RTS) Smoother-algoritmen af voor dit geaugmenteerde systeem.

   • Sequentiële Complexiteit: De filter- en smoothingstappen opereren op de geaugmenteerde toestandsdimensie ($d + N_{inst}$), waarbij een computationele complexiteit van $O(N)$ op een CPU wordt behouden.
   • Parallelle Complexiteit: Door gebruik te maken van de associativiteit van de filter- en smoothingverdelingen, wordt de methode herformuleerd als een all-prefix-sums-probleem. Dit maakt parallelle implementatie op GPU's mogelijk, wat de wandkloktijd reduceert tot $O(N/T + \log T)$, waarbij $T$ het aantal parallelle workers is.

4. Software-implementatie:
   De auteurs introduceren smolgp, een open-source Python/JAX-pakket. Het biedt drop-in-compatibiliteit met de bestaande tinygp-bibliotheek, met ondersteuning voor zowel standaard als blootstellingsbewust GP-modelleren. Het implementeert de hierboven beschreven sequentiële en parallelle oplossers.

Belangrijkste Resultaten

• Numerieke Equivalentie: De auteurs valideren de geaugmenteerde SSM-aanpak tegenover de volledige dichte GP-implementatie (Luhn et al. 2026) voor zowel single-instrument als multi-instrument (overlappende) scenario's. De resultaten voor geconditioneerde gemiddelden, varianties en log-likelihoods komen overeen tot op machineprecisie (verschillen $\sim 10^{-14}$ tot $10^{-15}$), wat de theoretische equivalentie bevestigt.
• Prestatiebenchmarks:
  • Geheugen: De SSM-aanpak is aanzienlijk geheugenefficiënter dan de dichte GP-methode. Waar dichte GPs $O(N^2)$ geheugen vereisen (meer dan 100 GB voor $N \sim 10^5$), schaalt de SSM-methode lineair, waardoor grote datasets haalbaar zijn op standaardhardware.
  • Rekentijd: Voor geïntegreerde data, waar geen quasi-scheidbaar (QSM) alternatief bestaat, biedt de SSM-methode een transformatieve snelheidswinst. Op een CPU is de sequentiële SSM-oplosser de snelste voor voorspelling en concurrerend voor likelihood-berekening. Op een GPU bereikt de parallelle SSM-oplosser de snelste conditioneringstijden, die QSM-methoden voor likelihoods evenaren.
  • Voorspelling: De methode ondersteunt efficiënte voorspelling op willekeurige testpunten met behulp van een specifiek algoritme (Algoritme 1) dat embarrassingly parallel en geheugenvriendelijk is in vergelijking met dichte matrixvermenigvuldigingen.

Betekenis en Claims
Het artikel claimt dat dit werk massale multi-instrument kruisvergelijkingen mogelijk maakt waarbij blootstellingsoverlapping alomtegenwoordig is, een scenario dat eerder beperkt werd door computationele kosten. Door een complexiteit van $O(N)$ (of $O(\log N)$ op GPU's) te behouden, ontsluit de methode het potentieel voor:

• Het analyseren van datasets met $N > 10^5$ observaties (bijvoorbeeld dagelijkse zonnebewaking van meerdere instrumenten zoals KPF, NEID, EXPRES en HARPS-N).
• Het toepassen van complexere modellen en hogere dimensionale datasets die eerder onbereikbaar waren vanwege de $O(N^3)$ bottleneck van geïntegreerde GPs.
• Het uitbreiden van schaalbare GP-modellering naar kernels die geen quasi-scheidbaarheid bezitten (bijvoorbeeld quasiperiodieke kernels) wanneer geïntegreerde metingen betrokken zijn.

De auteurs benadrukken dat het kader niet beperkt is tot tijdreeksen, maar van toepassing is op elke gebinde data met een sorteerbare onafhankelijke coördinaat, en dat het op natuurlijke wijze afgeleide observaties ondersteunt (bijvoorbeeld FF'-modellen) vanwege de opname van toestandsafgeleiden in de SDE-formulering. Het werk positioneert smolgp als een schaalbaar, flexibel instrument voor de tijdsdomeinanalyse van de astronomie van de volgende generatie.