  Gaussian process (GP) regression with 1D inputs can often be performed in
linear time via a stochastic differential equation formulation. However, for
non-Gaussian likelihoods, this requires application of approximate inference
methods which can make the implementation difficult, e.g., expectation
propagation can be numerically unstable and variational inference can be
computationally inefficient. In this paper, we propose a new method that
removes such difficulties. Building upon an existing method called
conjugate-computation variational inference, our approach enables linear-time
inference via Kalman recursions while avoiding numerical instabilities and
convergence issues. We provide an efficient JAX implementation which exploits
just-in-time compilation and allows for fast automatic differentiation through
large for-loops. Overall, our approach leads to fast and stable variational
inference in state-space GP models that can be scaled to time series with
millions of data points.
