Neural Inference of Gaussian Processes for Time Series Data of Quasars

The study of quasar light curves poses two problems: inference of the power spectrum and interpolation of an irregularly sampled time series. A baseline approach to these tasks is to interpolate a time series with a Damped Random Walk (DRW) model, in which the spectrum is inferred using Maximum Likelihood Estimation (MLE). However, the DRW model does not describe the smoothness of the time series, and MLE faces many problems in terms of optimization and numerical precision. In this work, we introduce a new stochastic model that we call $\textit{Convolved Damped Random Walk}$ (CDRW). This model introduces a concept of smoothness to a DRW, which enables it to describe quasar spectra completely. We also introduce a new method of inference of Gaussian process parameters, which we call $\textit{Neural Inference}$. This method uses the powers of state-of-the-art neural networks to improve the conventional MLE inference technique. In our experiments, the Neural Inference method results in significant improvement over the baseline MLE (RMSE: $0.318 \rightarrow 0.205$, $0.464 \rightarrow 0.444$). Moreover, the combination of both the CDRW model and Neural Inference significantly outperforms the baseline DRW and MLE in interpolating a typical quasar light curve ($\chi^2$: $0.333 \rightarrow 0.998$, $2.695 \rightarrow 0.981$). The code is published on GitHub.


Introduction
Quasars are astronomical sources powered by the accretion of matter falling onto the black hole in a galaxy's center. The accretion onto the black hole is a turbulent process, so the brightness varies intrinsically over time. This variability not only allows for studying the region around the black hole [1] but also enables one to trace the evolution of the Universe through strong gravitational lensing [2]. The former requires an interpretable model of radiation variability, whereas the latter requires interpolation of a time series with irregular time sampling.
Due to the accretion turbulence, the variability of quasar radiation is usually modelled by a stochastic Gaussian process called "Damped Random Walk" [DRW;3]. This model has two main problems: the first problem is that the DRW can model the quasar spectrum on time scales of years but not on the scales of months, where the quasar spectrum steepens.  and CDRW (orange) introduced in this paper. Error bars and χ 2 are derived from MCMC sampling [11; 12]. CDRW better describes the power spectrum at high frequencies. [4; 5; 6; 7]. The second problem is associated with the maximum likelihood estimation (MLE), which is used to infer the parameters of a DRW: this method has convergence issues, high computation complexity, and, sometimes, poor 2σ constraints on the inferred parameters.
In this work, we aim to solve these two problems of conventional DRW modelling of quasar light curves. In Section 2.1, we introduce a new physically-interpretable model Convolved Damped Random Walk (CDRW) that well describes the quasar spectrum on time scales, from days to decades. In Section 2.2, we propose the Neural Inference method that aims to solve problems of maximum likelihood estimation. Our approach is task-specific in contrast to the generalized approach in [8]. In Section 3, we present the data. In Section 4, we present the experiments. We conclude in Section 5.

Convolved Damped Random Walk model of quasar radiation
The conventional physical model of quasar light curve F (t) is given by the lamp-post radiation g LP and thin-disk reverberation P TD : where τ and ρ are the characteristic times of radiation from accretion onto the black hole and from its re-radiation by the outer disk, respectively [13]. Usually, F (t) is modelled as a DRW, but this model fails to describe the power spectrum at high frequencies (Fig. 2). This discrepancy can be mitigated if the driving lamp-post radiation g LP (t) is assumed to be DRW instead. In this case, convolution with thin-disk reverberation kernel P T D suppresses the power spectrum at high frequencies. The spectrum of this kernel is well fitted by a softened powerlaw. Such a derivation brings us to a new model of quasar radiation that we call Convolved Damped Random Walk (CDRW) with power spectrum: where ω is the circular frequency 2πf . Given τ ≥ ρ, τ is correlation time. In turn, the reverberation time ρ defines the smoothness of the curve. This model well describes the power spectra of real quasars on all frequencies (Fig. 2). The newly introduced smoothness of ρ results in better interpolation performance of the model (Fig. 1). Moreover, in border cases, it reduces to the main types of stochastic processes: white noise, Matern ν = 1.5 noise, DRW, and Random Walk (red noise).

Neural inference of Gaussian processes
Without prior knowledge of the quasar flux F observed at times t obs , the baseline approach to the inference of a Gaussian process is MLE, that is, maximization of likelihood p(F (t obs )|µ, σ 2 , τ, ρ). In the context of numerical methods, MLE requires matrix inversion that has computational complexity O(N 3 ). For light curves with a high number of observations N , such an inference requires much computational power, moreover, the numerical errors in matrix inversion can jeopardize the inference results. In the context of the theory of optimization, MLE faces convergence issues due to orders of magnitude differences between gradients on different parameters. Additionally, there are degeneracies between the parameters and biases of modes of the posterior distribution -e.g. Bessel correction for inverse Gamma posterior of sample variance [14].
To mitigate these negative aspects of Gaussian process inference, we introduce a novel method that we call Neural Inference. We use a neural network to predict the final correlation time τ and an initial guess of reverberation time ρ. In this case, maximization of the likelihood p(F (t obs ), τ |µ, σ 2 , ρ) leads to significantly more accurate parameters. Moreover, these estimations become more stable with respect to the stochastic realisation of the studied Gaussian process.

Data
We carry out experiments on synthetic data that mimic real observations. The observations of quasars have 135 ± 25 days-long annual gaps due to occultation by the Sun and bad weather. Within the observation season in the n-th year, we model time sampling as a stochastic process: where P oissonP is a discrete Poisson process [15], that models short observation gaps, and is time in the night. The statistical form and parameters of the observation schedule correspond to the COSMOGRAIL dataset [16].
Finally, τ min ≈ 11 to a minimal relevant DRW break reported in the work [19]. The curves are normalized to (-1,1). A 2-year-long sample from a mock light curve is presented in Fig. 1.  Neural network. For training, we simulate a dataset of 300,000 labels log(τ ) and ρ according to eq. 4 and split them in train/validation/test in ratio 0.8/0.1/0.1. Both stochastic observation times t obs and radiation F (t obs ) are sampled in the runtime of training to avoid overfitting of random realisations. The architecture of the network is presented in Fig. 3. In the neural network, the data is interpolated with B-splines on the grid (10,365), then processed by three ResNet blocks [20] with 64,128,256 filters. It is followed by linear reduction of the filters to 64 and the attention layer, which is applied by the Learnable Dictionary Encoding (LDE) [21; 22] with 64 clusters. Finally, the encodings are linearly transformed to predictions of log(τ ) and ρ.

Results
The neural network is trained with MSE loss for four epochs until convergence. It reaches RMSE 0.15 for log(τ ) and 0.98 for ρ on the test set. On Tesla P100, an epoch takes half an hour, given that sampling of the data is done on the GPU. We carry out data sampling using Jax [23], the network itself is implemented in PyTorch [24].
Experiments. To compare our methods to baseline approaches, we simulate 1000 mock quasar light curves with different stochastic realisations of t obs and F (t obs ) (eq's. 3,4), but the same spectral parameters τ and ρ for all these curves. First, we pick τ = 540 as the median correlation time for the observed quasars [19]. Second, we pick ρ = 4.77, so it results in the break on a month scale, which is midway between DRW ρ = 0 and the maximal reported break of two months.
For the first experiment, we compare the performance of Neural Inference (Fig 3) and baseline MLE optimization [TNC; 25; 26] for the inference of log(τ ) and ρ. The resulting metrics are presented in Table 1. In this experiment, Neural Inference outperforms baseline MLE in every metric.
For the second experiment, we compare performance in the interpolation of the same curves on a grid with a regular cadence of one day. In this experiment, the baseline approach to quasars is interpolation with DRW, inferred with MLE. In turn, the introduced approach is interpolation with CDRW, inferred using the neural inference technique. The resulting metrics are presented in Table 2.
The difference between interpolation means µ(t) (MAE, RMSE) is marginal, but the χ 2 of the new method is significantly better and predicts much more reliable uncertainties than the baseline.

Conclusion
In this work, we introduced the new Convolved Damped Random Walk (CDRW) model of quasar radiation (Sec. 2.1) that well describes spectra of quasars on all the time scales. Additionally, we introduced a new method of inference of Gaussian processes, which we call Neural Inference (Sec. 2.2). This method uses neural networks to mitigate several problems of conventional MLE optimization. A combination of these two novelties results in significant improvement in inference and interpolation of a model of a typical quasars light curve (Tables 1, 2).

Impact statement
Impact on Astrophysics: The new CDRW model decouples the temporal scales of accretion onto a black hole from that of reverberation of the disk surrounding it. The introduced reverberation time ρ is the direct tracer of the characteristic radius of the disk R λ , that traces the mass of the black hole, Eddington luminosity, and accretion efficiency. Therefore, the CDRW model introduces an explicit constraint on those parameters, which is designed to enhance our understanding of galaxy physics.
On the other hand, ρ defines the short-scale characteristic time of quasar evolution. The distribution of this parameter across quasars can help optimize the observation schedule to maximize the obtained information while minimizing the load on the telescope.
In turn, the Neural Inference approach significantly elevates the quality of inference and interpolation of mock quasar light curves. The improvement in inference results in more accurate information about the physics of quasars. The improvement in interpolation quality reveals itself in the studies of gravitationally lensed quasars. The combination of CDRW and Neural Inference provides better uncertainties that result in the more accurate inference of lensing-induced time delays and, therefore more strict constraints on the Hubble rate.
Impact on other domains: Time series phenomena permeate our everyday lives -from finance and healthcare to climate science and astrophysics, since the mathematical formulation of the task is often similar. The Gaussian process, especially DRW, is a baseline approach in many fields. Moreover, people face the same problems with MLE-based inference in all these fields. Therefore the introduced Neural Inference technique can be applied out of the box to constrain the kernel of a Gaussian process for many different types of data. Additionally, the convolutional term in the CDRW model can be well-adjusted to any data in order to account for domain-specific short-term interaction.