  Kalman filtering and smoothing are the foundational mechanisms for efficient
inference in Gauss-Markov models. However, their time and memory complexities
scale prohibitively with the size of the state space. This is particularly
problematic in spatiotemporal regression problems, where the state dimension
scales with the number of spatial observations. Existing approximate frameworks
leverage low-rank approximations of the covariance matrix. Since they do not
model the error introduced by the computational approximation, their predictive
uncertainty estimates can be overly optimistic. In this work, we propose a
probabilistic numerical method for inference in high-dimensional Gauss-Markov
models which mitigates these scaling issues. Our matrix-free iterative
algorithm leverages GPU acceleration and crucially enables a tunable trade-off
between computational cost and predictive uncertainty. Finally, we demonstrate
the scalability of our method on a large-scale climate dataset.
