Appendix A provides derivations supporting Section 3 in the main paper. In Appendix B, we explain our experimental setup, including dataset preparation and model implementation, in more detail. Finally, Appendix C provides additional results supporting our claims regarding the scalability of our method, together with additional results from the experiments presented in Section 4. In this section we provide detailed derivations of the ST-DGMRF joint distribution, for both firstorder transition models (Section A.1) and higher-order transition models (Section A.2). A.1 Joint distribution The LDS (see Section 2.2 and 3.1 in the main paper) defines a joint distribution over system states First, note that Eq. (1) can be written as a set of linear equations Moving all xk-terms to the left-hand side, we can rewrite this as a matrix-vector multiplication I F1 I F2 I ...... FKI | {z} Empty positions in F represent zero-blocks. Now, we can express x as an affine transformation of ϵ x = F 1c+F 1ϵ, (3) where F 1 exists because det(F) = 1. Since ϵ is distributed as ϵ N(0,Q 1) with Q = diag(Q0,Q1,...,QK), and c is deterministic, we can use the affine property of Gaussian distributions to obtain the joint distribution This reduces both computations and memory requirements. In contrast, the information vector η = Ωµcan be expressed compactly as η = FTQFF 1c = FTQc, (8) which can be computed efficiently using sparse and parallel matrix-vector multiplications on a GPU.

Appendix A provides derivations supporting Section 3 in the main paper. In this section we provide detailed derivations of the ST -DGMRF joint distribution, for both first-order transition models (Section A.1) and higher-order transition models (Section A.2). A.1 Joint distribution The LDS (see Section 2.2 and 3.1 in the main paper) defines a joint distribution over system states First, note that Eq. (1) can be written as a set of linear equations x We make use of this property in the DGMRF formulation and in the conjugate gradient method. Eq. 11 is converted into a discrete-time dynamical system by approximating ρ We consider two ST -DGMRF variants that capture different amounts of prior knowledge. DGMRF transition matrices can be parameterized accordingly. The air quality dataset is based on hourly PM2.5 measurements obtained from [ The raw PM2.5 measurements are log-transformed and standardized to zero mean and unit Ca. 50% of the nodes are masked out (purple nodes within We use a simple MLP with one hidden layer of width 16 with ReLU activations and no output non-linearity. The DGMRF parameters are not shared across time, allowing for dynamically changing spatial covariance patterns.