1
RDDMPI: Residual Denoising Diffusion Model for Probabilistic Multivariate Time Series Imputation
arXiv:2609.11648v1 [cs.LG] 10 Sep 2026
Ramiro Valdes Jara, David Chapman, Adam Meyers
Early approaches relied primarily on statistical and traditional machine learning techniques, including interpolation methods, matrix and tensor factorization, and Gaussian process models. Although these methods provide interpretable solutions, their ability to capture complex nonlinear temporal dynamics and cross-variable dependencies present in modern multivariate datasets is limited. To overcome these limitations, deep learningbased approaches learn rich temporal and cross-variable representations directly from data, leading to substantial improvements in reconstruction accuracy and transforming the MTSI landscape. Existing deep learning approaches can be broadly categorized into deterministic and probabilistic methods. Deterministic methods estimate a single value for each missing entry and have achieved strong performance by capturing temporal dynamics and cross-variable interactions. However, they do not explicitly model uncertainty, which is problematic when multiple plausible imputations are consistent with the observed data. Probabilistic methods address this limitation by modeling a conditional distribution over missing values rather than a single point estimate. Among them, diffusion models have recently emerged as a powerful framework for multivariate time series imputation. By learning an iterative reverse denoising process conditioned on observed entries, diffusion models Index Terms—Multivariate time series imputation, residual can generate multiple plausible imputations, thereby enabling diffusion, diffusion model, uncertainty quantification. uncertainty estimation. Despite their success, existing diffusion-based methods typically reconstruct the full missing signal directly from I. I NTRODUCTION noise. Consequently, the denoising network must jointly learn ULTIVARIATE time series data are widely collected global temporal structure, cross-variable dependencies, local in real-world domains such as healthcare monitoring, dynamics, and stochastic variability, resulting in a highly traffic systems, finance, energy, and climate science [1], complex learning problem. This observation is particularly [2], [3], [4], [5]. In these settings, each variable evolves relevant because modern deterministic imputers are already over time while interacting with other variables. However, capable of recovering much of the underlying signal structure. missing observations are common due to sensor failures, In many cases, the remaining prediction error contains both irregular sampling, communication failures, or data acquisition systematic, potentially reducible, errors in the reconstruction limitations. Since missing values can degrade downstream tasks and irreducible conditional variability arising when multiple such as forecasting, anomaly detection, decision support, and missing-value realizations are consistent with the observed system monitoring, multivariate time series imputation (MTSI) data. has become a fundamental problem in time series analysis. This raises a central question: can applying diffusion to residThe goal of MTSI is to recover missing entries by exploit- ual correction, rather than full-signal reconstruction, simplify ing both temporal continuity and cross-variable correlations. the generative learning problem and improve reconstruction accuracy and uncertainty estimation? We argue that residual Ramiro Valdes Jara, Adam Meyers are with the Department of Industrial modeling leads to a simpler and more targeted generative and Systems Engineering, University of Miami, United States of America (e-mail: {rjv71, axm8336}@miami.edu). problem. By reformulating probabilistic imputation in residual David Chapman is with the Department of Computer Science, University space, diffusion can focus on correcting systematic baseline of Miami, United States of America (e-mail: [email protected]). errors and modeling the remaining conditional variability rather Corresponding author: Adam Meyers. The code can be accessed at https://github.com/ramirovaldesjara/RDDMPI. than reconstructing the full signal from scratch. However,
Abstract—Multivariate time series imputation (MTSI) aims to recover missing values in temporal data composed of multiple interdependent variables. This problem is central to real-world applications such as healthcare monitoring, traffic networks, and energy systems. Recent diffusion-based approaches have shown strong potential for probabilistic imputation by learning to generate missing values through iterative denoising. However, most existing approaches perform diffusion directly in the original data space, requiring the denoising network to simultaneously capture global structure, temporal dynamics, and stochastic variability. This makes the generative task unnecessarily complex, especially when modern deterministic imputers can already provide accurate initial reconstructions. To address this limitation, we propose RDDMPI, a conditional residual diffusion framework that operates directly in residual space. Instead of modeling the full missing signal directly, we reformulate probabilistic imputation as a baseline-residual decomposition, where a pretrained model captures the dominant signal and a diffusion process models the residual uncertainty. To better exploit deterministic guidance, RDDMPI conditions the reverse denoising process on both the baseline-completed signal and its latent representation, while a reliability-aware conditioning mechanism adaptively controls the influence of baseline information during residual generation. This formulation simplifies the diffusion learning objective, enabling it to focus on structured correction terms rather than reconstructing the full signal. Experiments on multiple benchmark datasets demonstrate that RDDMPI consistently improves both reconstruction accuracy and uncertainty quantification.
M
2
leveraging deterministic predictions as conditioning information be categorized into deterministic (Section II-A) and probabilisintroduces an additional challenge. The quality of baseline tic approaches (Section II-B). In this section, we review these imputations can vary across datasets, missingness patterns, and two research directions, with particular emphasis on recent regions of the time series. While accurate baseline predictions diffusion-based methods that motivate the proposed residual provide valuable guidance, unreliable estimates may propagate diffusion framework. errors through the diffusion process. Therefore, conditioning information should be incorporated adaptively according to its A. Deterministic Time Series Imputation estimated reliability. Deterministic imputation methods aim to recover missing Motivated by these observations, we propose RDDMPI, a values by exploiting temporal dependencies and cross-variable Residual Denoising Diffusion Model for Probabilistic Mul- correlations present in the observed data. These methods tivariate Time Series Imputation. RDDMPI decomposes the produce a single imputed estimate for each missing entry and imputation task into two stages. First, a pretrained deterministic have demonstrated strong performance across a wide range of imputation model produces a baseline-completed signal and a applications. structured latent representation. Second, a conditional diffusion Traditional approaches include interpolation and smoothing model is trained in residual space, where it learns to generate methods [6], as well as low-rank matrix and tensor completion correction terms only for the missing regions. The diffusion methods [7], [2]. These methods are efficient and interpretable, process is conditioned on both the baseline reconstruction exploiting local smoothness or shared low-dimensional structure and its latent representation, allowing the model to focus on across variables and time. However, their smoothness and probabilistic residual refinement rather than reconstructing the low-rank assumptions limit their ability to capture complex full signal from scratch. In addition, RDDMPI incorporates a nonlinear dependencies. reliability-aware conditioning mechanism to adaptively control Early deep learning approaches, such as GRU-D [8] and the influence of baseline information and reduce the propagation BRITS [9], primarily relied on recurrent neural networks of unreliable deterministic estimates. (RNNs) to model temporal dynamics while explicitly handling The main contributions of this work are summarized as missing observations. However, sequential processing can follows: accumulate reconstruction errors and limit parallelization across timesteps during training and inference. More recently, • We reformulate probabilistic multivariate time series imputation as a baseline-residual decomposition, separat- attention-based architectures have become increasingly popular ing deterministic signal reconstruction from probabilistic due to their ability to capture long-range temporal dependencies and complex interactions among variables. SAITS [10] employs residual uncertainty modeling. diagonally masked self attention and combines two recon• We propose RDDMPI, a conditional residual diffusion framework that performs diffusion directly in residual struction stages, whereas NRTSI [11] treats observations as space and leverages both baseline-completed signals and permutation-equivariant sets and progressively imputes missing values without recurrent processing. Other methods jointly latent baseline representations to guide denoising. model temporal and feature-level dependencies through multi• We introduce a reliability-aware conditioning mechanism that adaptively modulates the contribution of deterministic dimensional, global-local, or structured attention mechanisms baseline information during the reverse diffusion process. [12], [13]. Recent methods can be divided into two broad categories. • We provide extensive empirical evidence across five One adapts general-purpose time series architectures, and the benchmark datasets showing that RDDMPI achieves other develops specialized architectures for time series imputastate-of-the-art performance in most evaluated settings, tion. General-purpose architectures based on temporal convoluconsistently outperforming deterministic and probabilistic tions [14], two-dimensional temporal variation modeling [15], baselines in both reconstruction accuracy and uncertainty inverted attention [16], and multiscale mixing [17] have been quantification adapted to imputation through masked reconstruction objectives. The remainder of this paper is organized as follows. These approaches benefit from advances in general time series Section II reviews existing deterministic and probabilistic representation learning but are not designed exclusively for approaches for multivariate time series imputation. Section III missing-value reconstruction. formulates the imputation problem. Section IV presents RDIn contrast, specialized imputation architectures introduce DMPI, including the residual diffusion formulation, the conmechanisms tailored to partially observed data. ImputeFormer ditional denoising network, the reliability-aware conditioning [18] employs low-rankness-induced attention to exploit the mechanism, and the theoretical motivation. Section V reports latent low-dimensional structure of spatiotemporal data and experimental results and ablation studies. Finally, Section VI improve generalization across sensors and missingness patterns. concludes the paper. T1 [19] introduces one-to-one channel-head binding, assigning each variable to a dedicated attention head to strengthen the II. R ELATED W ORK modeling of variable-specific dynamics while retaining crossResearch on multivariate time series imputation (MTSI) has variable interactions. Beyond architectural design, another evolved considerably, progressing from deterministic recon- line of research incorporates explicit inductive biases about struction methods toward probabilistic generative approaches signal structure or distributional similarity. Decompositioncapable of modeling uncertainty. Existing methods can broadly based methods separate trend, periodic, and local components
3
to facilitate the reconstruction of long-range and recurring FGTI [31] introduces spectral information to better reconstruct patterns [20]. Optimal-transport approaches instead recover periodic and high-frequency components. SSSD [32] replaces missing values by aligning distributions of time series patches. conventional denoising backbones with structured state-space In particular, PSW-I [21] incorporates spectral information layers designed to capture long-range temporal dependencies. into a Wasserstein-based discrepancy, enabling nonparametric Beyond diffusion, flow matching has recently emerged as an imputation under temporal and distributional shifts without alternative generative formulation for probabilistic imputation. training a separate predictive model. GiFlow [33] constructs a graph-informed source distribution For spatiotemporal data, relationships among variables can from the observed signal and learns a continuous transport proalso be encoded by a graph. Graph-recurrent architectures cess toward the target distribution, enabling efficient generation propagate information across both time and connected variables, while explicitly modeling spatial and temporal dependencies. while sparse spatiotemporal attention improves scalability to However, it is specifically designed for settings in which larger sensor networks [22], [23]. These methods are effective meaningful graph structure is available. when a meaningful spatial or relational graph is available. Despite their differences, existing diffusion-based imputation However, general multivariate time series may not provide methods share a common characteristic: diffusion is performed such a graph, requiring relationships among variables to be directly in the original data space. Consequently, the denoising learned directly from the observed data. network must learn both the dominant signal structure and the Despite their strong reconstruction performance, determinis- residual uncertainty simultaneously while reconstructing the tic methods inherently produce only a single imputed trajectory. entire missing signal from noise. This motivates the exploration Consequently, they do not explicitly represent the conditional of alternative formulations that simplify the generative task variability of missing values, which can be problematic when while retaining the uncertainty modeling benefits of diffusionmultiple plausible completions exist for a partially observed based imputation. sequence. This limitation has motivated the development of probabilistic approaches that model distributions over missing C. Residual Modeling and Position of This Work values rather than single-point estimates. Residual modeling has been widely adopted in machine learning as a strategy for simplifying complex prediction tasks. B. Probabilistic Time Series Imputation For example, gradient boosting incrementally improves an Probabilistic imputation methods seek to estimate the con- initial predictor by fitting successive models to its remaining ditional distribution of missing values given the observed errors [34]. More recently, residual formulations have been indata. Unlike deterministic approaches, these methods provide corporated into diffusion models. RDDM [35] separates image uncertainty estimates alongside imputations, allowing them restoration into a directional residual diffusion process and a to represent multiple plausible reconstructions and better stochastic noise diffusion process, distinguishing restoration characterize ambiguity in the missing regions. from sample diversity. This idea has recently gained attention Early probabilistic approaches primarily relied on latent- in time series forecasting, where RDIT [36] combines a point variable models and Bayesian formulations. Methods such as forecaster with conditional diffusion over its residuals, allowing V-RIN [24] combine recurrent architectures with variational the deterministic model to capture the central forecast while inference to learn uncertainty-aware latent representations, diffusion models the remaining predictive distribution. while Gaussian process-based approaches, including Multi-task RDDMPI proposes a residual diffusion formulation for GP [25] and GP-VAE [26], provide principled uncertainty quan- probabilistic multivariate time series imputation. Unlike RDDM, tification through probabilistic priors. However, these methods which addresses image restoration, and RDIT, which predicts struggle to capture highly complex temporal dependencies. future values, RDDMPI models residuals only at missing More recently, diffusion models have emerged as a pow- positions. A pretrained deterministic imputer first reconstructs erful probabilistic framework for time series imputation. By the dominant signal structure, after which diffusion models learning to reverse a gradual noising process, diffusion models the conditional distribution of its remaining errors. The can approximate complex conditional distributions without denoising process is conditioned on the baseline-completed strong distributional assumptions. CSDI [27] established the signal, its latent representation, and the observation mask, foundation for diffusion-based time series imputation by while reliability-aware conditioning adaptively controls the formulating the task as a conditional score-based diffusion influence of potentially inaccurate baseline information. Thus, process, demonstrating that diffusion models can generate RDDMPI applies residual diffusion specifically to missingmultiple plausible imputations while providing uncertainty value reconstruction without requiring the diffusion model to estimates. generate the complete missing signal from noise. Subsequent work has improved diffusion-based imputation through more informative structural conditioning, conIII. P ROBLEM D EFINITION sistency constraints, and specialized denoising architectures. Spatiotemporal diffusion models incorporate graph and geoIn this section, we present the mathematical setup and graphic relationships into the conditioning process [28], while problem definition. Let X0 ∈ RV ×L denote a multivariate consistency-aware formulations encourage coherent estimates time series with V variables and L timesteps, and let M ∈ across variables, timesteps, or complementary views [29], [30]. {0, 1}V ×L be the corresponding observation mask, where
4
( 1, Mv,ℓ = 0,
if X0,v,ℓ is observed, otherwise.
Using the mask, the observed and missing components of the time series can be written as X0ob = M ⊙ X0 ,
(1)
X0mi = (1 − M ) ⊙ X0 ,
(2)
where ⊙ denotes element-wise multiplication. The objective of multivariate time series imputation is to recover the missing component X0mi from the available observations X0ob . Formally, an imputation model seeks to learn a mapping X̂0 = f (X0ob , M ),
(3)
where X̂0 denotes the reconstructed complete time series. From a probabilistic perspective, the goal is to estimate the conditional distribution
deterministic baseline exclusively within the unobserved region. These corrections capture both systematic baseline errors and the remaining conditional variability. The proposed residual diffusion model is conditioned on the baseline-completed signal, the baseline latent representation, and the observation mask. We collectively denote this conditioning information as c = X̃ base , H base , M . (7) After sampling a residual correction R̂0mi , the final imputed series is obtained by preserving observed entries and correcting the baseline only on missing positions as X̂0 = M ⊙ X0 + (1 − M ) ⊙ X̃ base + R̂0mi . (8) Figure 1 illustrates the overall RDDMPI training process. B. Conditional Residual Diffusion
To model uncertainty over missing values, RDDMPI employs a conditional Denoising Diffusion Probabilistic Model (DDPM) [37] operating in residual space. DDPMs are latent-variable which characterizes all plausible values of the missing entries generative models that learn a data distribution p(X0 ) through given the observed context. Deterministic methods approximate a sequence of latent variables {Xt }Tt=1 . These latent variables this distribution through a single point estimate, whereas correspond to progressively noisier versions of the original probabilistic approaches seek to model the full conditional data, where X0 denotes a clean sample and XT approaches distribution and provide uncertainty-aware imputations. RDan isotropic Gaussian distribution after repeated perturbations. DMPI belongs to the latter category and models this conditional The diffusion framework consists of two Markov processes: distribution through a residual diffusion process. a forward diffusion process that gradually corrupts data by adding Gaussian noise and a reverse diffusion process that IV. M ETHODOLOGY learns to recover clean samples from noisy observations. In this section, we present RDDMPI, a residual diffusion In the proposed framework, we adapt this formulation to framework for probabilistic multivariate time series imputa- the residual imputation setting by performing diffusion on the tion. We begin in Section IV-A by introducing a baseline- missing-region residual Rmi . Let {Rmi }T denote the noisy t 0 t=1 residual decomposition that reformulates imputation as residual residual states generated by the forward diffusion process. modeling around a deterministic reconstruction. Section IV-B Given the conditioning information c defined in Eq. (7), the then presents a conditional diffusion process that learns the objective is to learn the conditional residual distribution p(Rmi | 0 distribution of residual corrections over missing entries. Next, c), which characterizes plausible residual corrections around Section IV-C describes the proposed conditional denoising the deterministic baseline reconstruction. network, including the reliability-aware conditioning mecha1) Forward Diffusion Process: The forward diffusion pronism used to incorporate baseline information during denoising. cess progressively perturbs the clean residual target with Finally, Section IV-D provides a theoretical motivation for Gaussian noise over T diffusion steps. Formally, the forward performing diffusion in residual space rather than directly in Markov chain is defined as the original data space. T Y mi mi q(R1:T | R0mi ) = q(Rtmi | Rt−1 ), (9) A. Baseline-Residual Decomposition t=1 √ mi mi q(Rtmi | Rt−1 ) = N Rtmi ; αt Rt−1 , (1 − αt )I , (10) Given the observed component X0ob and mask M , a pretrained deterministic imputation model first produces a baselineAt each diffusion step, a small amount of Gaussian noise is completed signal X̃ base , obtained by filling the missing entries added while preserving part of the original residual signal. As with the baseline model’s deterministic predictions, and its the process progresses, the residual state becomes increasingly base latent representation H as corrupted and gradually loses its structure. The distribution at X̃ base , H base = fbase (X0ob , M ). (5) an arbitrary diffusion step t can be written in closed form as √ Instead of modeling the missing signal directly, we define q(Rtmi | R0mi ) = N Rtmi ; ᾱt R0mi , (1 − ᾱt )I , (11) the missing-region residual target as t Y αt := 1 − βt , ᾱt = αs . (12) mi base R0 = (1 − M ) ⊙ X0 − X̃ . (6) s=1 p(X0mi | X0ob , M ),
(4)
The residual target is defined only at missing positions, ensuring that the diffusion model learns corrections to the
where βt denotes the variance schedule and ᾱt represents the cumulative signal retention coefficient up to diffusion step
5
Noise
Masked Noise
Minimize
Residual Noisy Residual
Input Time Series Pretrained Deterministic Baseline
Inject Missingness
Baseline Completed Signal
Partially Observed Series
Diffusion Step t
Conditional Denoising Network
Noise Prediction
Conditioning Information c Baseline Latent Representation Mask M
Fig. 1: Overview of the residual diffusion training process. Missingness is injected into the input time series to obtain a partially observed series, which is processed by a pretrained baseline to produce a baseline-completed signal, a latent representation and a residual target. Noise is added only to the missing-region residual, and the diffusion model is conditioned on the baseline-completed signal, the baseline latent representation, the noise timestep, and the mask to predict the masked noise. Training minimizes the discrepancy between the true masked noise and the predicted noise.
t. Using the reparameterization property of DDPMs, a noisy residual at any diffusion step can be sampled directly as Rtmi =
√ √ ᾱt R0mi + 1 − ᾱt ϵ,
ϵ ∼ N (0, I).
(13)
As t increases, the residual state converges toward a Gaussian noise distribution. A complete derivation of the residual forward process is provided in Appendix B. 2) Reverse Diffusion Process: The objective of the reverse process is to progressively remove noise from Rtmi while conditioning on c. Starting from RTmi ∼ N (0, I), the model learns a conditional reverse Markov chain mi pθ (R0:T | c) = p(RTmi )
T Y
mi pθ Rt−1 | Rtmi , t, c ,
(14)
t=1 mi pθ Rt−1 | Rtmi , t, c
=N
mi Rt−1 ; µθ (Rtmi , t, c), σt2 I
.
where µθ (Rtmi , t, c) denotes the learned conditional mean of mi Rt−1 under the reverse transition, and σt2 is the corresponding reverse-process variance, fixed according to the diffusion noise ᾱt−1 schedule. We set σt2 = β̃t , where β̃t = 1− 1−ᾱt βt . Following the standard DDPM parameterization, the reverse mean is expressed through the noise prediction network ϵθ , whose output ϵ̂θ = ϵθ (Rtmi , t, c) denotes the predicted noise, as 1 1 − αt mi mi µθ (Rt , t, c) = √ Rt − √ ϵ̂θ . (15) αt 1 − ᾱt 3) Training Objective: The reverse process introduced above depends on the noise prediction network ϵθ , which is trained to estimate the injected Gaussian noise ϵ used to construct Rtmi at an arbitrary diffusion step t in Eq. (13) by minimizing h i 2 Lθ = ERtmi , ϵ, t ϵ − ϵ̂θ Rtmi , t, c ⊙ (1 − M ) 2 , (16)
Importantly, the sampled noise ϵ is not provided directly as an input to the denoising network. Instead, it serves as the supervision target in the training objective, while ϵθ receives the noisy residual Rtmi , the diffusion step t, and the conditioning information c. The mask restricts optimization to missing entries, preventing the model from being penalized on observed values already fixed by the conditioning signal. The complete self-supervised training procedure is provided in Appendix C. 4) Inference: After training, the denoising network ϵθ parameterizes the reverse diffusion process. The overall inference procedure of RDDMPI is summarized in Algorithm 1. The baseline model first provides a deterministic baseline reconstruction and latent representation. Then, the residual diffusion model performs N independent reverse sampling trajectories from t = T to t = 1, producing a set of residual samples {R0mi,n }N n=1 . These samples approximate the conditional residual distribution around the deterministic baseline. Each sample produces a plausible imputation by adding the sampled residual to the baseline at missing positions. The complete sample set is retained for probabilistic evaluation and uncertainty quantification, whereas its element-wise median is used as the point estimate for reconstruction-accuracy evaluation. The architecture used to parameterize ϵθ (Rtmi , t, c) is described next in Section IV-C. C. Conditional Denoising Network The diffusion process is parameterized by a conditional denoising network ϵθ (Rtmi , t, c). As illustrated in Figure 2, the network consists of three main components: a reliabilityaware input fusion mechanism, a FiLM-based conditioning module using the baseline latent representation, and a stack of
6
Diffusion-Step Embedding
Diffusion Step t Mask M Baseline-Completed Signal
Temporal Convolutional Module
Sigmoid
Reliability Map A
Initial Hidden Representation
FiLMConditioned Representation
Residual Denoising Block (1)
Residual Denoising Block (L)
Output Projection
Noise Prediction
Input Projection
Noisy Residual
Input Projection Residual Denoising Block (Internal Structure)
TemporalPosition Embedding
Variable Embedding
Mask M
Baseline Latent Representation
Reliability Map A Skip Connections
Diffusion-Step Embedding Block Input
Side Info Temporal Transformer
Variable Transformer
Gated Activation Unit
Block Output
Conditional Denoising Network
Fig. 2: Architecture of the conditional denoising network ϵθ , which receives the noisy residual Rtmi , diffusion step t, and e base , H base , M ), and outputs the noise prediction b conditioning information c = (X ϵθ . The first residual denoising block receives (0) (b−1) the FiLM-conditioned representation, Ut = Z̃t , and each subsequent block receives the output of the preceding block, Ut . A complete mathematical description of the residual denoising blocks is provided in Appendix E.
Algorithm 1: Imputation (Sampling) with RDDMPI. Input : Partially observed sample X0 ; mask M ; pretrained deterministic baseline fbase ; trained denoising network ϵθ ; number of stochastic samples N . Output : Imputed sequence X̂0 . Compute observed context: X0ob ← M ⊙ X0 ; Obtain deterministic baseline-completed signal and latent representation: X̃ base , H base ← fbase (X0ob , M ); Define conditioning information: c ← (X̃ base , H base , M ); for n = 1 to N do Initialize missing-region residual with Gaussian noise: mi,(n) RT ∼ N (0, I); for t = T to 1 do Predict noise using the conditional residual diffusion model: ϵ̂θ ← ϵθ (Rtmi,n , t, c); Compute DDPM reverse mean: t µθ ← √1αt Rtmi,n − √1−α ϵ̂ ; 1−ᾱt θ Sample reverse step: mi,n Rt−1 ← µθ + σt z , z ∼ N (0, I); end end Estimate final residual correction using the element-wise n oN mi,n mi median: R̂0 ← median R0 ; n=1
Recover the missing signal by adding the residual correction to the baseline: X̂0mi ← (1 − M ) ⊙ X̃ base + R̂0mi ; Return the full imputed sequence: X̂0 ← X0ob + X̂0mi ;
residual denoising blocks. Together, these components allow the denoiser to adaptively incorporate information from the deterministic baseline while estimating the noise contained in the current residual state. 1) Reliability-Aware Conditioning: The initial input to the denoising network is the noisy missing-region residual Rtmi ∈ RV ×L , where V is the number of variables and L
e base is the sequence length. The baseline-completed signal X is incorporated separately through a reliability-aware fusion mechanism. Since the quality of deterministic baseline predictions can vary across timesteps and variables, conditioning information should not be incorporated uniformly. To address this issue, we introduce a learned reliability gate that adaptively controls the influence of the baseline reconstruction on the denoising representation. The reliability map is computed as h i . (17) A = σ tconv X̃ base , M Here, tconv (·) is a lightweight temporal convolutional module, σ(·) denotes the sigmoid function, and A ∈ [0, 1]V ×L acts as a learned reliability gate over variables and timesteps. The reliability module applies short- and long-range depthwise temporal convolutions independently to each variable, followed by pointwise projections and a sigmoid activation. The resulting map acts as a learned gate rather than an explicitly supervised estimate of baseline error. The noisy residual and baseline-completed signal are projected independently to dh hidden channels. The initial hidden representation is constructed as e base ), Zt = ϕres (Rtmi ) + A ⊙ ϕctx (X
(18)
where ϕres and ϕctx are the independent input projections. The reliability map is broadcast across the hidden channels dh . Thus, the noisy residual defines the primary denoising state, while the deterministic baseline is incorporated through reliabilityweighted fusion. The resulting representation Zt ∈ RV ×L×dh serves as the initial hidden representation, which is subsequently modulated by latent conditioning (Section IV-C2) before being processed by the residual denoising blocks. The same reliability map A is passed to later layers as additional conditioning information.
7
2) FiLM-Based Latent Conditioning: The latent representation extracted from the deterministic baseline is additionally used to condition the denoising process. Because its original shape depends on the selected baseline architecture, it is first aligned with the variable and temporal dimensions of the denoising representation. A learned projection then produces multiplicative and additive FiLM [38] tensors as
D. Theoretical Motivation for Residual Diffusion
(19)
A natural question is whether performing diffusion in residual space provides a principled advantage over directly modeling the missing signal. Let X0mi denote the missing target and let f (c) be a deterministic baseline prediction computed from conditioning information c. We define the residual variable on the missing entries as R0mi = (1 − M ) ⊙ X0mi − f (c) . (21)
where galign (·) denotes the baseline-specific reshaping and temporal alignment operation, and gFiLM (·) is a learned linear projection. After dimension permutation, γH , δH ∈ RV ×L×dh , FiLM conditioning is applied as
This transformation corresponds to a conditional reparameterization rather than a fundamentally different inference problem. In particular, the conditional distributions satisfy pR (r | c) = pX r + f (c) | c , (22)
(γH , δH ) = gFiLM galign H base
,
which implies that residual diffusion is equivalent to standard conditional diffusion under a change of variables. The practical advantage of this formulation arises when the where Zt is the initial hidden representation. This operation deterministic baseline removes a substantial fraction of the injects structured information from the baseline latent repre- conditional mean of the missing signal. Let m(c) = E[X0mi | c] sentation into the denoising process by adaptively scaling and denote the conditional mean. By adding and subtracting m(c), shifting the hidden features independently across channels, the residual can be decomposed as variables, and timesteps. The result is the FiLM-conditioned R0mi = X0mi − m(c) + m(c) − f (c) . (23) hidden representation Z̃t ∈ RV ×L×dh , which is provided as input to the residual denoising blocks. Through this Taking squared norms and expectations yields conditioning, the denoiser can exploit the temporal and crossh i h i h i 2 2 2 E R0mi = E X0mi − m(c) + E ∥m(c) − f (c)∥ variable dependencies encoded by the deterministic baseline while retaining the noisy residual as its primary denoising state. + 2E X0mi − m(c), m(c) − f (c) , (24) | {z } 3) Residual Denoising Blocks: As illustrated in Figure 2, =0 the FiLM-conditioned hidden representation Zet is processed by a stack of residual denoising blocks inspired by DiffWave [39]. where the cross term vanishes due to conditional centering, mi We denote the representation entering the b-th block by since E[X0 −m(c) | c] = 0. A complete derivation is provided (b−1) (0) et . Each block produces an updated in Appendix A-A. This decomposition shows that the residual Ut , with Ut = Z representation that is passed to the next block, together with a consists of two components: the intrinsic conditional irreducible uncertainty and the squared bias of the deterministic baseline skip representation used by the final output path. Within each block, a projected diffusion-step embedding model. Consequently, when f (c) closely approximates m(c), is first added to the incoming representation, allowing the the diffusion model only needs to learn the remaining residual denoiser to adapt to the current noise level. Temporal and variation rather than the conditional mean structure underlying variable transformer layers then process the representation the missing signal. This interpretation is particularly relevant in score-based sequentially. The temporal transformer models dependencies diffusion. Let sX (xt , c, t) = ∇xt log p(xt | c) denote the across timesteps separately for each variable, while the variable conditional score in data space. When diffusion is performed transformer models interactions among variables separately in residual space, the practical score induced by the baseline, at each timestep. The resulting features are combined with denoted s , is given by f projected side information consisting of temporal-position √ embeddings, variable embeddings, the observation mask, and sf (rt , c, t) = sX rt + ᾱt f (c), c, t , (25) the learned reliability map A. The conditioned representation is passed through a gated whereas the ideal score centered at the conditional mean, activation and projected into a residual update and a skip denoted sm , is √ representation. The residual update is combined with the sm (rt , c, t) = sX rt + ᾱt m(c), c, t . (26) block input and forwarded to the next block, while the skip representations from all blocks are aggregated by the Under a standard regularity assumption that sX (·, c, t) is Lt output path. Two final pointwise projections, with a ReLU Lipschitz in its first argument, we obtain √ activation between them, produce b ϵθ = ϵθ (Rtmi , t, c) ∈ RV ×L , ∥sf (rt , c, t) − sm (rt , c, t)∥ ≤ Lt ᾱt ∥f (c) − m(c)∥. (27) mi which has the same variable and temporal dimensions as Rt and represents the Gaussian noise predicted for the reverse Therefore, if the deterministic baseline accurately approximates diffusion process. A complete mathematical description of the the conditional mean, the induced residual score remains close block operations, residual and skip connections, and output to the ideal centered score. This yields a lower-energy and aggregation is provided in Appendix E. more correction-oriented target distribution, resulting in a Z̃t = (1 + γH ) ⊙ Zt + δH ,
(20)
8
more favorable learning problem for finite-capacity denoising networks in practice. Complete derivations of the score-based analysis are provided in Appendix A-B. V. E VALUATION In this section, we evaluate RDDMPI on multivariate time series imputation under different missingness scenarios. We first describe the datasets, evaluation protocol, baseline methods, metrics, and implementation details. We then report quantitative results for reconstruction accuracy and uncertainty quantification, followed by ablation studies. A. Evaluation Setup
are adapted here using masked reconstruction, whereas others were designed explicitly for imputation. We briefly summarize the role of each baseline below. The general-purpose time series models include: DLinear [44]: A linear baseline that decomposes temporal dynamics through linear projections, providing a strong low-complexity benchmark. • ModernTCN [14]: A modern convolutional architecture based on large-kernel temporal convolutions, designed to capture multi-scale temporal patterns efficiently. • iTransformer [16]: A Transformer variant that models multivariate time series through variable-wise tokenization, allowing for capturing cross-variable interactions. • TimesNet [15]: A temporal convolutional model that captures multi-period temporal variation by transforming time series into structured 2D representations.
•
1) Datasets: We evaluate RDDMPI on five multivariate time series benchmark datasets: ETTh1 [40], ETTh2 [40], Weather [41], Exchange [42], and Illness [43]. These datasets cover The specialized deterministic imputation methods include: diverse domains, including energy, climate, finance, and public • SAITS [10]: A self-attention-based imputation model health. For all datasets, we use an input window length of 96. specifically designed for multivariate time series, with ETTh1 and ETTh2 contain hourly electricity transformer masked training objectives for missing-value reconstrucmeasurements, each with 7 correlated variables related to tion. oil temperature and load covariates. Weather contains 21 • ImputeFormer [18]: A transformer-based imputation meteorological variables collected at 10-minute intervals from architecture tailored for generalizable spatiotemporal the Max Planck Institute weather station. Exchange contains imputation through low-rank structured attention. 8 daily international exchange rate series, representing a non• T1 [19]: A CNN-transformer hybrid imputation model stationary financial time series setting. Illness contains 7 weekly that combines temporal convolutional feature extraction influenza-related variables from the U.S. Centers for Disease with selective cross-variable information transfer through Control and Prevention (CDC) surveillance system. During channel-head binding. testing, non-overlapping windows are used to avoid repeated evaluation over the same temporal regions. The probabilistic imputation methods include: 2) Experimental Design: All experiments are conducted • GP-VAE [26]: A latent-variable model that combines a using five random seeds: 2, 102, 202, 302, and 402. Training variational autoencoder with a Gaussian-process prior to is performed on an NVIDIA RTX 6000 GPU. To evaluate gencapture temporal correlations. eralization under different observability regimes, we consider • CSDI [27]: A conditional score-based diffusion model both point-wise and structured block missingness. for probabilistic time series imputation, representing a In the point missing setting, we test the model under strong diffusion-based baseline with uncertainty-aware missing ratios of 0.2, 0.4, 0.6, and 0.8, corresponding to generation. 20%, 40%, 60%, and 80% missing entries, with missing • FGTI [31]: A frequency-aware diffusion model that positions sampled independently and uniformly at random. extracts spectral information from the observed values In the block missing setting, we simulate realistic sensor failure to guide the denoising process. patterns by combining two types of corruption: (i) a point All baseline implementations are based on established missing component with rate 0.05, where 5% of entries are frameworks including Time-Series Library1 , PyPOTS [45], and randomly removed, and (ii) a sequence missing component with Awesome-Imputation [46] repositories to ensure reproducibility rate 0.0015, meaning that approximately 0.15% of positions and fair comparison. initiate a missing block. For every initiated block, its length 4) Evaluation Metrics: Let M denote the set of artificially is sampled uniformly from the integers between 24 and 96, and the corresponding consecutive entries of that variable are masked positions used for evaluation, consisting of elements masked. Blocks extending beyond the end of a window are (v, ℓ), where v indexes the variable and ℓ indexes the timestep. truncated, and overlapping blocks are combined. This mixed Its cardinality is denoted by |M|. At position (v, ℓ), the protocol produces localized measurement dropouts together ground-truth and imputed values are denoted by yv,ℓ and x̂v,ℓ , with extended contiguous gaps, yielding a more challenging respectively. For reconstruction accuracy, we use mean absolute error (MAE) and mean squared error (MSE), computed over and practically relevant evaluation scenario. 3) Baselines: We compare the proposed method against the positions in M as ten representative baselines spanning three model categories: X 1 MAE = |x̂v,ℓ − yv,ℓ | , (28) general-purpose time series models, specialized deterministic |M| (v,ℓ)∈M imputation methods, and specialized probabilistic imputation methods. This distinction is important because some baselines 1 https://github.com/thuml/Time-Series-Library were originally developed for general time series analysis and
9
TABLE I: Imputation performance on five benchmark datasets under point and block missing scenarios. Results are averaged across four point missing ratios (0.2, 0.4, 0.6, 0.8). Dataset abbreviations: ETTh1/2 = Eh1/2, Exchange = Exch, Illness = Illn, and Weather = Wthr. Best results are marked in bold and second-best results in underlined. SAITS MSE MAE
ImputeFormer MSE MAE
TimesNet MSE MAE
T1 MSE
MAE
GP-VAE MSE MAE
CSDI MSE MAE
FGTI MSE MAE
Eh1
iTransformer MSE MAE
Point 0.0548 Block 0.0168
0.1365 0.0851
0.2184 0.2906 0.1145 0.2111 0.1646 0.2569 0.1351 0.2113 0.2856 0.3044 0.1634 0.2567 0.0807 0.1679 0.6246 0.5860 0.0641 0.1470 0.0718 0.1520 0.1728 0.2844 0.0592 0.1735 0.1016 0.2097 0.0257 0.1075 0.0594 0.1543 0.0920 0.2117 0.0258 0.1077 0.2842 0.4231 0.0179 0.0896 0.0170 0.0877
Eh2
ModernTCN MSE MAE
Point 0.0442 Block 0.0203
0.1123 0.0744
0.0801 0.1851 0.0576 0.1524 0.0712 0.1745 0.4370 0.4100 0.6667 0.4457 0.0749 0.1799 0.0442 0.1264 1.4139 0.8712 0.0735 0.1444 0.0982 0.1696 0.0770 0.1890 0.0477 0.1414 0.0546 0.1552 0.1471 0.2711 0.2654 0.2672 0.0533 0.1592 0.0292 0.1014 0.6001 0.5663 0.0561 0.1070 0.0948 0.1321
Exch
DLinear MSE MAE
Point 0.0022 Block 0.0043
0.0203 0.0181
0.0053 0.0453 0.0095 0.0664 0.0033 0.0339 0.2311 0.3815 0.0693 0.1098 0.0034 0.0335 0.0019 0.0218 0.6415 0.6944 0.0485 0.1135 0.0056 0.0409 0.0063 0.0557 0.0054 0.0498 0.0034 0.0336 0.1864 0.3330 0.1217 0.1233 0.0036 0.0362 0.0115 0.0312 0.5100 0.6290 0.2237 0.1458 0.0177 0.0457
Illn
RDDMPI (Ours) MSE MAE
Point 0.0236 Block 0.0390
0.0714 0.1282
0.2175 0.3046 0.0835 0.1736 0.1211 0.2201 0.2873 0.3044 0.3394 0.3423 0.0969 0.2054 0.0252 0.0849 0.6355 0.5177 0.0806 0.1479 0.0860 0.1464 0.2603 0.3691 0.1521 0.2647 0.3345 0.3598 0.1475 0.2373 0.2697 0.3100 0.1368 0.2535 0.0874 0.1792 0.3270 0.4122 0.0977 0.1895 0.0516 0.1283
Wthr
Models Metric
Point 0.0310 Block 0.0233
0.0315 0.0254
0.0472 0.0882 0.0425 0.0799 0.0920 0.1452 0.0453 0.0647 0.0475 0.0594 0.0471 0.0903 0.0363 0.0582 0.1768 0.2578 0.0321 0.0316 0.0335 0.0320 0.0495 0.1053 0.0371 0.0827 0.1003 0.1480 0.0251 0.0328 0.0401 0.0461 0.0390 0.0855 0.0248 0.0400 0.0627 0.1349 0.0234 0.0254 0.0236 0.0256
MSE =
1 |M|
X
2
(x̂v,ℓ − yv,ℓ ) .
(29)
(v,ℓ)∈M
MAE measures the average reconstruction error, while MSE penalizes large deviations more strongly. For probabilistic imputation, we use the continuous ranked probability score (CRPS), which jointly evaluates the accuracy and distributional quality of the predictive distribution. Given a predictive cumulative distribution function Fv,ℓ , CRPS is defined as X Z ∞ 1 2 CRPS = (Fv,ℓ (z) − I(z ≥ yv,ℓ )) dz. |M| −∞ (v,ℓ)∈M
(30) In this paper, the predictive distribution is approximated using (n) N = 100 generated samples {x̃v,ℓ }N n=1 at each masked position, with empirical cumulative distribution N
F̂v,ℓ (z) =
1 X (n) I x̃v,ℓ ≤ z . N n=1
(31)
Following the common evaluation protocol for diffusion-based time series imputation, we compute a normalized quantilebased approximation. Let Q = {0.05, 0.10, . . . , 0.95}, and let (q) x̂v,ℓ be the empirical q-quantile of the generated samples. Then, P (q) 1 X 2 (v,ℓ)∈M ρq yv,ℓ − x̂v,ℓ P CRPS = , (32) |Q| (v,ℓ)∈M |yv,ℓ | q∈Q
where the quantile loss is ρq (u) = u(q − I(u < 0)). 5) Implementation Details: All methods are evaluated using identical train, validation, and test splits under the experimental protocol described above. For general-purpose time series models DLinear [44], ModernTCN [14], iTransformer [16], and TimesNet [15], we employ a unified training setup using 0.4 point-wise masking during training. Optimization is performed using Adam with learning rate 10−3 , batch size 16, and a maximum of 300 epochs with early stopping. For the specialized imputation models SAITS [10], ImputeFormer [18], T1 [19], GP-VAE [26], CSDI [27], and FGTI [31] we preserve the original training procedures and hyperparameters reported in their official implementations whenever available. An analysis of parameter count, training cost, and inference efficiency is provided in Appendix D.
RDDMPI adopts T1 [19] as the pretrained deterministic backbone, providing both the baseline-completed signal and the latent representation used to condition the residual diffusion model. The deterministic backbone is pretrained and subsequently frozen during diffusion training. For the diffusion component, we employ a DDPM formulation with T = 50 diffusion steps. The variance schedule increases quadratically between the initial and final noise levels, allocating smaller noise increments to the earlier diffusion steps. Its architectural hyperparameters are selected separately for each dataset. These include the learning rate, number of residual denoising blocks, and number of hidden channels, among other hyperparameters. The complete dataset-specific configurations are reported in Appendix F. During training, a self-supervised masking strategy is adopted, where a masking ratio is randomly sampled and applied to the observed entries to construct reconstruction targets. During inference, N = 100 residual samples are generated through the reverse diffusion process. The median prediction is used for deterministic evaluation, while the complete sample set is used to compute probabilistic metrics such as CRPS and to construct predictive intervals. B. Experimental Results In this section, we evaluate the proposed method from two complementary perspectives: deterministic reconstruction accuracy and probabilistic uncertainty quantification. Since the goal of probabilistic imputation is not only to recover missing values accurately but also to provide well-calibrated predictive distributions, we report both types of metrics separately for clarity. The tables in the main paper report results averaged across the four point missing ratios for conciseness. Full results for each individual missing ratio, together with the corresponding standard deviations and block missing results, are provided in Appendix G. 1) Reconstruction Accuracy: We first assess deterministic imputation quality using MAE and MSE. The corresponding quantitative results across all benchmark datasets are reported in Table I. RDDMPI achieves the best result in 18 of the 20 comparisons, indicating that residual diffusion provides a robust improvement over both deterministic and probabilistic baselines in reconstruction accuracy. Relative to its deterministic backbone, T1, our proposed model improves in 19 of the 20 comparisons. The improvement is particularly consistent
10
RDDMPI CRPS
GP-VAE CRPS
CSDI CRPS
FGTI CRPS
Eh1
Point Block
0.1313 0.0811
0.7376 0.5298
0.1406 0.0867
0.1454 0.0832
Eh2
Point Block
0.0641 0.0412
0.6395 0.4076
0.0826 0.0596
0.0979 0.0751
Exch
TABLE II: Probabilistic imputation performance in terms of CRPS on five benchmark datasets under point and block missing scenarios. Results for the point missing setting are averaged across four missing ratios (0.2, 0.4, 0.6, 0.8). Dataset abbreviations: ETTh1/2 = Eh1/2, Exchange = Exch, Illness = Illn, and Weather = Wthr. Best results are marked in bold, and second-best results are marked in underlined.
Point Block
0.0179 0.0153
0.7848 0.6903
0.0990 0.1238
0.0354 0.0391
Models Metric
Figure 3 provides a qualitative comparison with CSDI. Both methods recover the local signal morphology, but the median imputations produced by RDDMPI more closely follow the ground-truth targets, particularly along the descending segment in ETTh2 and the short-term fluctuations in Exchange. RDDMPI also produces narrower predictive intervals that remain centered around the reconstructed trajectory and contain most of the target values. In these examples, its generated imputations are therefore less dispersed and maintain close agreement with the ground truth. Additional qualitative results under different missing ratios are provided in Appendix H. C. Ablation Studies
Wthr
Illn
As a framework, RDDMPI is designed to extend a strong deterministic imputation baseline with probabilistic residual Point 0.0728 0.6668 0.1558 0.1445 Block 0.1238 0.4793 0.1873 0.1173 refinement. We first examine whether its benefits depend on the selected deterministic backbone. We then evaluate the Point 0.0420 0.4501 0.0421 0.0465 Block 0.0338 0.2383 0.0336 0.0348 contributions of baseline conditioning and reliability-aware fusion. under block missingness, where RDDMPI outperforms T1 in 1) Effect of the Deterministic Backbone: To determine both MSE and MAE on all five datasets. These results suggest whether the proposed framework depends on the use of T1, that the benefits of residual diffusion are especially consistent we instantiate RDDMPI with either T1 or ImputeFormer as when the deterministic baseline is less accurate. the pretrained deterministic backbone. Table IV compares each 2) Uncertainty Quantification: We next evaluate probabilis- complete model with its corresponding deterministic baseline. tic imputation quality using the continuous ranked probability RDDMPI reduces both MSE and MAE for both backbones score (CRPS), which measures the compatibility between the across every evaluated point and block missingness setting. For predicted distribution and the observed target values. Since T1, the gains generally increase with the point missing ratio, deterministic methods do not model predictive distributions, indicating that residual correction becomes increasingly useful CRPS comparisons are reported against GP-VAE, CSDI, and as the baseline reconstruction problem becomes more difficult. FGTI baselines in Table II. The improvements obtained with ImputeFormer are larger RDDMPI obtains the lowest CRPS in eight of the ten dataset because its baseline errors leave more substantial correction and missingness settings. In particular, it obtains the best result terms for the residual diffusion model to recover. These results on all five datasets under point missingness and on ETTh1, demonstrate that the proposed formulation can generalize to ETTh2, and Exchange under block missingness. FGTI performs different deterministic imputation architectures, although the best on Illness under block missingness, while CSDI achieves magnitude of the improvement depends on the quality and a marginally lower CRPS on Weather under block missingness. error structure of the selected deterministic backbone. 2) Component Ablation: To isolate the contributions of the Taken together, these results show that the advantages of RDDMPI are not limited to the median imputation used to proposed components, we compare three residual diffusion compute MSE and MAE. These results are consistent with the variants on the Exchange dataset that progressively incorporate hypothesis that modeling corrections around a deterministic the proposed design choices. reconstruction provides a more focused generative objective RDDMPI without baseline conditioning. In this variant, than reconstructing the complete missing signal directly from the model performs diffusion in residual space but does noise. not receive any conditioning derived from the deterministic
Fig. 3: Select examples of probabilistic time series imputation on the ETTh2 and Exchange datasets (zoomed-in views). The red crosses denote observed values, and the blue circles denote the ground-truth imputation targets. For each method, the median imputation is shown as a solid line, while the shaded region represents the 5% and 95% predictive quantiles.
11
TABLE III: Component ablation of RDDMPI on Exchange dataset under point missing and block missing settings. Lower is better. Point Missing Method
Block Missing
MSE
0.2 MAE
CRPS
MSE
0.4 MAE
CRPS
MSE
0.6 MAE
CRPS
MSE
0.8 MAE
CRPS
MSE
MAE
CRPS
RDDMPI w/o baseline conditioning
0.00184
0.01436
0.01322
0.00129
0.01536
0.01511
0.00193
0.02045
0.01824
0.00309
0.02995
0.02622
0.00948
0.02665
0.02151
RDDMPI w/o reliabilityaware fusion
0.00184
0.01408
0.01283
0.00129
0.01646
0.01462
0.00183
0.02006
0.01765
0.00309
0.02920
0.02585
0.01007
0.03077
0.02495
RDDMPI (full)
0.00210
0.01465
0.01312
0.00148
0.01705
0.01511
0.00193
0.02055
0.01814
0.00300
0.02883
0.02521
0.00426
0.01805
0.01530
TABLE IV: Comparison between deterministic baselines and RDDMPI on the ETTh1 dataset under point missing and block missing settings. For each deterministic model, RDDMPI uses that model as its pretrained baseline. IF denotes ImputeFormer. ∆ denotes the percentage improvement of RDDMPI over its corresponding baseline, where higher ∆ indicates greater improvement. Lower is better for MSE and MAE. Best results are marked in bold. Method
T1
MSE
MSE
Point Missing 0.4 0.6 MAE MSE MAE
Block Missing MSE
0.8 MAE
MSE
MAE
Base 0.0296 0.1112 0.0397 0.1276 0.0657 0.1632 0.1878 0.2694 0.0258 0.1077 RDDMPI 0.0228 0.0955 0.0336 0.1122 0.0480 0.1362 0.1144 0.2020 0.0168 0.0850 ∆ (%)
IF
0.2 MAE
22.97
14.12
15.37
12.07
26.94
16.54
39.08
25.02
34.88
21.08
Base 0.0797 0.1697 0.1546 0.2232 0.2931 0.3137 0.6148 0.5108 0.0594 0.1542 RDDMPI 0.0294 0.1091 0.0402 0.1263 0.0625 0.1566 0.1441 0.2292 0.0208 0.0973 ∆ (%)
63.11
35.71
74.00
43.41
78.68
50.08
76.56
55.13
64.98
36.90
that the deterministic reconstruction provides useful temporal and cross-variable guidance at moderate missingness, while its errors are not yet sufficiently severe to require adaptive gating. However, this behavior changes in the more difficult regimes. The complete RDDMPI obtains the best MSE, MAE, and CRPS at point-0.8 and under block missingness. At point-0.8, the limited number of observed values makes the deterministic reconstruction less dependable. Under block missingness, contiguous gaps additionally remove local temporal context. Incorporating baseline information uniformly can therefore propagate inaccurate estimates into the denoiser. The reliability-aware mechanism mitigates this effect by adapting the contribution of baseline-derived information across variables and timesteps. Overall, the ablation results do not indicate that additional conditioning is uniformly beneficial. Instead, they show that residual diffusion alone can be sufficient in simpler regimes, baseline conditioning becomes useful as reconstruction difficulty increases, and reliability-aware fusion provides its clearest benefit when baseline errors have a greater effect on the final imputation.
backbone. In particular, neither the baseline-completed signal nor the associated latent features are provided to the denoiser. As a result, conditioning mechanisms such as FiLM modulation and reliability-aware fusion are entirely removed. The model therefore relies solely on the noisy residual and diffusion step. This setting isolates the effect of residual diffusion alone. RDDMPI without reliability-aware fusion. Here, baseline VI. C ONCLUSIONS conditioning is available, including both the completed signal and latent features from the deterministic backbone. However, the proposed reliability-aware fusion mechanism is removed, In this paper, we presented RDDMPI, a conditional residual meaning that all conditioning information is treated uniformly. diffusion framework for probabilistic multivariate time series This variant corresponds to a conditional residual diffusion for- imputation. By decomposing the imputation problem into a demulation and allows us to evaluate whether naive conditioning terministic initial reconstruction and a residual refinement stage, is sufficient. the proposed approach allows the diffusion model to focus on Full RDDMPI. The full model combines residual diffusion modeling structured correction terms rather than the full missing with baseline conditioning and the proposed reliability-aware signal directly. In addition, the framework leverages both latent fusion mechanism. In this setting, the learned reliability gate representations extracted from the pretrained baseline and adaptively modulates the contribution of baseline-derived reliability-aware conditioning mechanisms to guide denoising information, rather than incorporating the deterministic context under diverse missingness patterns. Across five benchmark datasets, RDDMPI achieves the best reconstruction result in uniformly during denoising. Results analysis. Table III shows that the contributions of 18 of the 20 aggregated MSE and MAE comparisons and the baseline conditioning and reliability-aware fusion depend on lowest CRPS in eight of the ten probabilistic comparisons. the missingness regime. At point missing ratios 0.2 and 0.4, The main limitation of RDDMPI is its computational cost. the variants without baseline conditioning or without reliability- The framework requires a separately pretrained deterministic aware fusion obtain the best results. When observations backbone, and generating N stochastic imputation samples remain sufficiently distributed throughout the sequence, residual through T reverse steps entails N T denoiser evaluations. Future diffusion can recover the correction term with limited reliance work will therefore investigate accelerated or reduced-step on baseline-derived guidance. sampling, and joint training of the deterministic and diffusion At point-0.6, baseline conditioning without reliability-aware components. Extending the framework to other downstream fusion achieves the best MSE, MAE, and CRPS. This suggests tasks also represents an important direction for future research.
12
R EFERENCES [1] X. Yi, Y. Zheng, J. Zhang, and T. Li, “St-mvl: filling missing values in geosensory time series data,” in Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence, 2016, p. 2704–2710. [2] H. Tan, G. Feng, J. Feng, W. Wang, Y.-J. Zhang, and F. Li, “A tensor-based method for missing traffic data completion,” Transportation Research Part C: Emerging Technologies, vol. 28, pp. 15–27, 2013. [3] F. V. Nelwamondo, S. Mohamed, and T. Marwala, “Missing data: A comparison of neural network and expectation maximization techniques,” Current Science, pp. 1514–1521, 2007. [4] A. T. Hudak, N. L. Crookston, J. S. Evans, D. E. Hall, and M. J. Falkowski, “Nearest neighbor imputation of species-level, plot-scale forest structure attributes from lidar data,” Remote Sensing of Environment, vol. 112, no. 5, pp. 2232–2245, 2008. [5] S. Van Buuren and K. Groothuis-Oudshoorn, “mice: Multivariate imputation by chained equations in r,” Journal of Statistical Software, vol. 45, no. 3, p. 1–67, 2011. [6] X. Yi, Y. Zheng, J. Zhang, and T. Li, “St-mvl: filling missing values in geo-sensory time series data,” in Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence, ser. IJCAI’16, 2016, p. 2704–2710. [7] H.-F. Yu, N. Rao, and I. S. Dhillon, “Temporal regularized matrix factorization for high-dimensional time series prediction,” in Proceedings of the 30th International Conference on Neural Information Processing Systems. Red Hook, NY, USA: Curran Associates Inc., 2016, p. 847–855. [8] Z. Che, S. Purushotham, K. Cho, D. Sontag, and Y. Liu, “Recurrent neural networks for multivariate time series with missing values,” 2016. [9] W. Cao, D. Wang, J. Li, H. Zhou, Y. Li, and L. Li, “Brits: bidirectional recurrent imputation for time series,” in Proceedings of the 32nd International Conference on Neural Information Processing Systems, 2018, p. 6776–6786. [10] W. Du, D. Côté, and Y. Liu, “Saits: Self-attention-based imputation for time series,” Expert Systems with Applications, vol. 219, p. 119619, 2023. [11] S. Shan, Y. Li, and J. B. Oliva, “Nrtsi: Non-recurrent time series imputation,” in ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2023, pp. 1–5. [12] P. Bansal, P. Deshpande, and S. Sarawagi, “Missing value imputation on multidimensional time series,” Proc. VLDB Endow., vol. 14, no. 11, p. 2533–2545, 2021. [13] Q. Suo, W. Zhong, G. Xun, J. Sun, C. Chen, and A. Zhang, “Glima: Global and local time series imputation with multi-directional attention learning,” in 2020 IEEE International Conference on Big Data (Big Data), 2020, pp. 798–807. [14] L. Donghao and W. Xue, “ModernTCN: A modern pure convolution structure for general time series analysis,” in The Twelfth International Conference on Learning Representations, 2024. [15] H. Wu, T. Hu, Y. Liu, H. Zhou, J. Wang, and M. Long, “Timesnet: Temporal 2d-variation modeling for general time series analysis,” in International Conference on Learning Representations, 2023. [16] Y. Liu, T. Hu, H. Zhang, H. Wu, S. Wang, L. Ma, and M. Long, “itransformer: Inverted transformers are effective for time series forecasting,” in The Twelfth International Conference on Learning Representations, 2024. [17] S. Wang, H. Wu, X. Shi, T. Hu, H. Luo, L. Ma, J. Y. Zhang, and J. ZHOU, “Timemixer: Decomposable multiscale mixing for time series forecasting,” in International Conference on Learning Representations (ICLR), 2024. [18] T. Nie, G. Qin, W. Ma, Y. Mei, and J. Sun, “Imputeformer: Low ranknessinduced transformers for generalizable spatiotemporal imputation,” in Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, ser. KDD ’24. New York, NY, USA: Association for Computing Machinery, 2024, p. 2260–2271. [19] D. Park, H. Ryu, S. Bae, K. Park, and H.-S. Kim, “T1: One-to-one channel-head binding for multivariate time-series imputation,” in The Fourteenth International Conference on Learning Representations, 2026. [20] S. LIU, X. Li, G. Cong, Y. Chen, and Y. JIANG, “Multivariate time-series imputation with disentangled temporal representations,” in The Eleventh International Conference on Learning Representations, 2023. [21] H. Wang, zhengnan li, H. Li, X. Chen, M. Gong, BinChen, and Z. Chen, “Optimal transport for time series imputation,” in The Thirteenth International Conference on Learning Representations, 2025. [22] A. Cini, I. Marisca, and C. Alippi, “Filling the g ap s: Multivariate time series imputation by graph neural networks,” in International Conference on Learning Representations, 2022.
[23] I. Marisca, A. Cini, and C. Alippi, “Learning to reconstruct missing data from spatiotemporal graphs with sparse observations,” arXiv preprint arXiv:2205.13479, 2022. [24] A. W. Mulyadi, E. Jun, and H.-I. Suk, “Uncertainty-aware variationalrecurrent imputation network for clinical time series,” IEEE Transactions on Cybernetics, vol. 52, no. 9, pp. 9684–9694, 2022. [25] E. V. Bonilla, K. Chai, and C. Williams, “Multi-task gaussian process prediction,” in Advances in Neural Information Processing Systems, vol. 20. Curran Associates, Inc., 2007. [26] V. Fortuin, D. Baranchuk, G. Raetsch, and S. Mandt, “Gp-vae: Deep probabilistic time series imputation,” in Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, ser. Proceedings of Machine Learning Research, vol. 108. PMLR, 2020, pp. 1651–1661. [27] Y. Tashiro, J. Song, Y. Song, and S. Ermon, “Csdi: conditional score-based diffusion models for probabilistic time series imputation,” in Proceedings of the 35th International Conference on Neural Information Processing Systems, 2021. [28] M. Liu, H. Huang, H. Feng, L. Sun, B. Du, and Y. Fu, “Pristi: A conditional diffusion framework for spatiotemporal imputation,” in 2023 IEEE 39th International Conference on Data Engineering (ICDE), 2023, pp. 1927–1939. [29] X. Wang, H. Zhang, P. Wang, Y. Zhang, B. Wang, Z. Zhou, and Y. Wang, “An observed value consistent diffusion model for imputing missing values in multivariate time series,” in Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, ser. KDD ’23. New York, NY, USA: Association for Computing Machinery, 2023, p. 2409–2418. [30] J. Zhou, J. Li, G. Zheng, X. Wang, and C. Zhou, “Mtsci: A conditional diffusion model for multivariate time series consistent imputation,” in Proceedings of the 33rd ACM International Conference on Information and Knowledge Management, ser. CIKM ’24. New York, NY, USA: Association for Computing Machinery, 2024, p. 3474–3483. [31] X. Yang, Y. Sun, X. Yuan, and X. Chen, “Frequency-aware generative models for multivariate time series imputation,” in The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024. [32] J. L. Alcaraz and N. Strodthoff, “Diffusion-based time series imputation and forecasting with structured state space models,” Transactions on Machine Learning Research, 2023. [33] Z. Zhang, A. Einizade, J. H. Giraldo, and O. Fink, “Spatiotemporal imputation with graph-informed flow matching,” 2026. [34] J. H. Friedman, “Greedy function approximation: A gradient boosting machine.” The Annals of Statistics, vol. 29, no. 5, pp. 1189 – 1232, 2001. [35] J. Liu, Q. Wang, H. Fan, Y. Wang, Y. Tang, and L. Qu, “Residual denoising diffusion models,” 2024. [36] C.-Y. Lai, Y.-C. Ning, and D. S. Boning, “Rdit: Residual-based diffusion implicit models for probabilistic time series forecasting,” 2025. [37] J. Ho, A. Jain, and P. Abbeel, “Denoising diffusion probabilistic models,” in Advances in Neural Information Processing Systems, vol. 33. Curran Associates, Inc., 2020, pp. 6840–6851. [38] E. Perez, F. Strub, H. de Vries, V. Dumoulin, and A. Courville, “Film: visual reasoning with a general conditioning layer,” in Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence and Thirtieth Innovative Applications of Artificial Intelligence Conference and Eighth AAAI Symposium on Educational Advances in Artificial Intelligence, 2018. [39] Z. Kong, W. Ping, J. Huang, K. Zhao, and B. Catanzaro, “Diffwave: A versatile diffusion model for audio synthesis,” in International Conference on Learning Representations, 2021. [40] H. Zhou, S. Zhang, J. Peng, S. Zhang, J. Li, H. Xiong, and W. Zhang, “Informer: Beyond Efficient Transformer for Long Sequence TimeSeries Forecasting,” Proceedings of the AAAI Conference on Artificial Intelligence, vol. 35, pp. 11 106–11 115, 2021. [41] Wetterstation, “Weather.” [Online]. Available: https://www.bgc-jena.mpg. de/wetter/ [42] G. Lai, W.-C. Chang, Y. Yang, and H. Liu, “Modeling long- and short-term temporal patterns with deep neural networks,” in The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval. Association for Computing Machinery, 2018, p. 95–104. [43] CDC, “Illness.” [Online]. Available: https://gis.cdc.gov/grasp/fluview/ fluportaldashboard.html [44] A. Zeng, M. Chen, L. Zhang, and Q. Xu, “Are transformers effective for time series forecasting?” 2023. [45] W. Du, “PyPOTS: A Python Toolkit for Data Mining on PartiallyObserved Time Series,” KDD 2023 MiLeTS, 2023.
13
[46] W. Du, J. Wang, L. Qian, Y. Yang, F. Liu, Z. Wang, Z. Ibrahim, H. Liu, Z. Zhao, Y. Zhou, W. Wang, K. Ding, Y. Liang, B. A. Prakash, and Q. Wen, “Tsi-bench: Benchmarking time series imputation,” arXiv preprint arXiv:2406.12747, 2024. [47] J. Song, C. Meng, and S. Ermon, “Denoising diffusion implicit models,” 2022.
14
A PPENDIX A T HEORETICAL A NALYSIS OF R ESIDUAL D IFFUSION A. Residual Energy Decomposition Let m(c) = E[X0mi | c] denote the conditional mean. We introduce and subtract m(c) to obtain R0mi = X0mi − m(c) + m(c) − f (c) . We now analyze the second moment of the residual. Expanding the squared norm yields ∥R0mi ∥2 = ∥X0mi − m(c)∥2 + ∥m(c) − f (c)∥2 + 2⟨X0mi − m(c), m(c) − f (c)⟩. Taking the conditional expectation with respect to c, we obtain E ∥R0mi ∥2 | c = E ∥X0mi − m(c)∥2 | c + ∥m(c) − f (c)∥2 + 2 E ⟨X0mi − m(c), m(c) − f (c)⟩ | c . Vanishing of the cross term. We now show that the cross term vanishes. Since m(c) = E[X0mi | c], we have E[X0mi − m(c) | c] = 0. Moreover, m(c) − f (c) is deterministic given c. Therefore, using linearity of expectation, E ⟨X0mi − m(c), m(c) − f (c)⟩ | c = E[X0mi − m(c) | c], m(c) − f (c) = ⟨0, m(c) − f (c)⟩ = 0. Thus, E ∥R0mi ∥2 | c = E ∥X0mi − m(c)∥2 | c + ∥m(c) − f (c)∥2 . Taking the expectation over c and applying the law of total expectation, we obtain E ∥R0mi ∥2 = E ∥X0mi − m(c)∥2 + E ∥m(c) − f (c)∥2 . B. Score-Based Diffusion Perspective Let xt denote the noisy variable at diffusion step t, and let pt (xt | c) denote its conditional marginal distribution induced by the forward diffusion process. The corresponding conditional score function is defined as sX (xt , c, t) = ∇xt log pt (xt | c). When performing diffusion in residual space, the noisy residual is given by √ √ rt = ᾱt R0mi + 1 − ᾱt ϵ, ϵ ∼ N (0, I). This induces a corresponding noisy variable in data space through the transformation √ xt = rt + ᾱt f (c). The score function associated with this shifted representation is therefore √ sf (rt , c, t) = sX rt + ᾱt f (c), c, t . In contrast, the ideal centered residual representation would shift by the conditional mean m(c) = E[X0mi | c], leading to the score √ sm (rt , c, t) = sX rt + ᾱt m(c), c, t . Assume that for each diffusion step t, the score function sX (·, c, t) is Lt -Lipschitz in its first argument, i.e., ∥sX (u, c, t) − sX (v, c, t)∥ ≤ Lt ∥u − v∥, Applying the Lipschitz property with u = rt + we obtain ∥sf (rt , c, t) − sm (rt , c, t)∥ = sX (rt +
√
√
ᾱt f (c),
v = rt +
√
ᾱt f (c), c, t) − sX (rt +
∀u, v.
ᾱt m(c),
√
√ ᾱt m(c), c, t) ≤ Lt ᾱt ∥f (c) − m(c)∥.
This inequality shows that, under regularity conditions, the discrepancy between the practical residual score sf and the ideal centered score sm is controlled by the baseline approximation error ∥f (c) − m(c)∥. In particular, if f (c) provides an accurate approximation of the conditional mean m(c), then the induced residual score remains close to the optimal centered score. √ Moreover, the discrepancy is modulated by the diffusion coefficient ᾱt , implying that the influence of baseline error diminishes as the diffusion process progresses toward higher noise levels. Consequently, residual diffusion yields a score field that is both stable and closer to the ideal centered representation, leading to a simpler and more favorable learning problem for a finite capacity denoising network.
15
A PPENDIX B M ATHEMATICAL D ETAILS OF F ORWARD P ROCESS D ERIVATION OF R ESIDUAL D IFFUSION This appendix derives the closed-form expression for the noisy residual Rtmi given in Eq. (13), starting from the one-step forward diffusion transition in Eq. (10). The one-step forward transition is √ √ mi Rtmi = αt Rt−1 + 1 − αt ϵt , ϵt ∼ N (0, I). Similarly, mi Rt−1 =
√
mi αt−1 Rt−2 +
p
1 − αt−1 ϵt−1 ,
ϵt−1 ∼ N (0, I).
(33)
Substituting Eq. (33) into the one-step forward transition, we obtain √ p √ √ mi Rtmi = αt αt−1 Rt−2 + 1 − αt−1 ϵt−1 + 1 − αt ϵt p √ √ mi = αt αt−1 Rt−2 + αt (1 − αt−1 )ϵt−1 + 1 − αt ϵt . Continuing this expansion recursively, we obtain √ Rtmi = αt αt−1 · · · α1 R0mi p + αt αt−1 · · · α2 (1 − α1 )ϵ1 p + αt αt−1 · · · α3 (1 − α2 )ϵ2 + ··· p √ + αt (1 − αt−1 )ϵt−1 + 1 − αt ϵt . Qt √ Now define ᾱt = i=1 αi . Thus, the first term becomes ᾱt R0mi . Now consider the noise terms. Since ϵi ∼ N (0, I) are i.i.d., each scaled term is p αt · · · αi+1 (1 − αi )ϵi ∼ N (0, αt · · · αi+1 (1 − αi )I) . Since linear combinations of independent Gaussian variables remain Gaussian, the sum of these terms is also Gaussian with variance given by the sum of the individual variances; that is, t X
αt · · · αi+1 (1 − αi ) = 1 − ᾱt .
i=1
Therefore, the total noise term becomes
√
Finally, we obtain Rtmi =
√
1 − ᾱt ϵ,
ᾱt R0mi +
√
ϵ ∼ N (0, I). 1 − ᾱt ϵ,
ϵ ∼ N (0, I).
This yields the closed-form forward noising expression used in Eq. (13). A PPENDIX C RDDMPI T RAINING P ROCEDURE Following the conditional residual diffusion formulation described in Section IV, let X0 ∈ RV ×L denote a training window sampled from the data distribution, with V variables and sequence length L. We use v ∈ {1, . . . , V } to index variables and ℓ ∈ {1, . . . , L} to index timesteps. Let M ∈ {0, 1}V ×L be the binary observation mask, where Mv,ℓ = 1 indicates an observed value, whereas Mv,ℓ = 0 indicates a missing value. During training, a subset of originally observed entries is randomly masked and treated as missing targets. The residual diffusion process is then applied only to the missing-region residual R0mi , while the observed mask, the baseline-completed signal, and the latent baseline representation remain fixed and serve as conditioning information. At diffusion step t, Gaussian noise ϵ ∼ N (0, I) is injected into the missing-region residual √ √ Rtmi = ᾱt R0mi + 1 − ᾱt ϵ, Qt where {αt }Tt=1 is a predefined variance schedule and ᾱt = s=1 αs . The model is trained to predict the injected noise using a conditional noise predictor ϵθ Rtmi , t, c . Since diffusion is applied only to the missing-region residual, the denoising loss is computed exclusively on missing positions: " # V X L X 1 2 L(θ) = E P (1 − Mv,ℓ ) (ϵv,ℓ − ϵθ (·)v,ℓ ) . v,ℓ (1 − Mv,ℓ ) v=1 ℓ=1
16
The self-supervised training procedure of the proposed residual diffusion model is summarized in Algorithm 2. Algorithm 2: Training of RDDMPI 1: Input: training data distribution q(X0 ), pretrained deterministic baseline fbase , number of optimization iterations Niter , Qt noise schedule {αt }Tt=1 with ᾱt = s=1 αs 2: Output: trained conditional noise predictor ϵθ 3: for i = 1 to Niter do 4: Sample a training window X0 ∼ q(X0 ) and a diffusion timestep t ∼ Uniform({1, . . . , T }) 5: Sample a mask by randomly hiding a subset of currently observed entries, producing M 6: Compute observed and missing components: X0ob ← M ⊙ X0 , X0mi ← (1 − M ) ⊙ X0 ob 7: Obtain baseline-completed signal and latent representation: X̃ base, H base ← f base (X0 , M ) mi base 8: Construct the missing-region residual target: R0 ← (1 − M ) ⊙ X0 − X̃ 9: Sample Gaussian noise ϵ ∼ N (0, I) with the same shape R0mi √ √ as mi mi 10: Apply forward noising to the residual target: Rt ← ᾱt R0 + 1 − ᾱt ϵ 11: Predict noise using the residual diffusion model: ϵ̂ ← ϵθ Rtmi , t, c PV PL 2 1 (1 − M ) (ϵ − ϵ̂ ) 12: Take a gradient step minimizing ∇θ P (1−M v,ℓ v,ℓ v,ℓ v=1 ℓ=1 v,ℓ ) v,ℓ 13: end for A PPENDIX D E FFICIENCY A NALYSIS Table V presents computational efficiency and performance metrics across RDDMPI, DLinear [44], ModernTCN [14], iTransformer [16], TimesNet [15], SAITS [10], ImputeFormer [18], T1 [19], GP-VAE [26], CSDI [27] and FGTI [31]. TABLE V: Computational efficiency and performance comparison on the ETTh1 and Illness datasets. Params (M): parameters in millions; Train Speed: ms per iteration; Inference Speed: ms per sample; MAE: Mean Absolute Error averaged over point missing ratios 0.2, 0.4, 0.6, and 0.8 (lower is better). Dataset Model
Parameters (M) Train Speed (ms/iter) Inference Speed (ms/sample)
MAE
RDDMPI (Ours) DLinear ModernTCN iTransformer TimesNet ETTh1 SAITS ImputeFormer T1 GP-VAE CSDI FGTI
1.23 0.02 1.72 0.22 0.59 5.27 1.37 0.54 0.11 1.19 1.96
213.96 4.02 7.42 7.27 22.06 51.68 21.83 12.80 23.04 188.98 45.72
2220.24 0.16 0.50 0.08 1.28 0.28 0.47 0.33 19.60 2315.40 3696.13
0.1365 0.2906 0.2111 0.2569 0.2567 0.2113 0.3044 0.1679 0.5860 0.1470 0.1520
RDDMPI (Ours) DLinear ModernTCN iTransformer TimesNet SAITS ImputeFormer T1 GP-VAE CSDI FGTI
0.44 0.02 0.76 0.22 4.69 5.27 1.37 0.54 0.11 0.41 0.52
313.97 5.48 22.10 18.56 52.15 83.61 39.84 15.24 23.19 268.52 87.67
12614.79 1.07 3.46 1.70 15.81 6.32 5.79 5.47 19.24 13997.47 18958.83
0.0714 0.3046 0.1736 0.2201 0.2054 0.3044 0.3423 0.0849 0.5177 0.1479 0.1464
Illness
Table V shows that RDDMPI has a substantially higher inference cost than deterministic baselines, although it achieves better imputation accuracy on both datasets. This difference is primarily caused by iterative probabilistic sampling rather than the model size. Deterministic methods produce a point estimate through a single forward pass, whereas generating one stochastic imputation sample with a diffusion model requires a complete reverse trajectory of T sequential denoising steps. At each step, the current noisy state is updated according to 1 1 − αt xt−1 = √ xt − √ ϵθ (xt , t, cond) + σt z, z ∼ N (0, I). αt 1 − ᾱt
17
Also, a single reverse trajectory produces one imputed sample, i.e., one possible estimate of the missing values. To approximate the predictive distribution, RDDMPI generates N independent samples, each requiring its own T -step reverse trajectory. Inference therefore requires N T denoising-network evaluations in total. The resulting samples are aggregated to obtain the reported point imputation and are also used to quantify predictive uncertainty. Although the N trajectories may be evaluated in parallel through batching, the T denoising steps within each trajectory remain sequential. Consequently, inference cost depends on both the number of diffusion steps T and the number of generated samples N , explaining the higher inference times of RDDMPI, CSDI, and FGTI relative to single-pass deterministic methods. a) Accelerated Diffusion Sampling.: To mitigate this limitation, we explore faster sampling strategies such as Denoising Diffusion Implicit Models (DDIM) [47], which reduce the number of required steps while maintaining performance. DDIM introduces a non-Markovian deterministic sampling process: p √ xt−1 = ᾱt−1 x̂0 + 1 − ᾱt−1 ϵθ (xt , t), where x̂0 is the predicted clean signal. By selecting a subset of timesteps {t1 , . . . , tK } with K ≪ T , inference can be significantly accelerated. We evaluate DDIM-based sampling with reduced step counts and observe that RDDMPI maintains competitive performance even with substantially fewer inference steps, while achieving notable speedups. This highlights the potential of accelerated diffusion methods to bridge the efficiency gap with deterministic models.
RDDMPI (T = 50) MSE MAE CRPS
RDDMPI (T = 10) MSE MAE CRPS
MSE
∆ (in %) MAE
CRPS
ETTh1
0.2 0.4 0.6 0.8
0.0229 0.0337 0.0480 0.1145
0.0956 0.1123 0.1362 0.2020
0.0924 0.1080 0.1310 0.1940
0.0236 0.0345 0.0495 0.1216
0.0970 0.1139 0.1390 0.2102
0.0957 0.1117 0.1358 0.2031
+3.06 +2.37 +3.13 +6.20
+1.46 +1.42 +2.06 +4.06
+3.57 +3.43 +3.66 +4.69
ETTh2
0.2 0.4 0.6 0.8
0.0228 0.0305 0.0437 0.0797
0.0771 0.0924 0.1150 0.1649
0.0438 0.0528 0.0657 0.0941
0.0223 0.0295 0.0422 0.0760
0.0786 0.0935 0.1157 0.1646
0.0455 0.0540 0.0666 0.0941
-2.19 -3.28 -3.43 -4.64
+1.95 +1.19 +0.61 -0.18
+3.88 +2.27 +1.37 0.00
Exchange
0.2 0.4 0.6 0.8
0.0022 0.0015 0.0019 0.0030
0.0147 0.0171 0.0206 0.0288
0.0131 0.0151 0.0181 0.0252
0.0022 0.0042 0.0035 0.0030
0.0150 0.0185 0.0219 0.0290
0.0138 0.0169 0.0198 0.0257
0.00 +180.00 +84.21 0.00
+2.04 +8.19 +6.31 +0.69
+5.34 +11.92 +9.39 +1.98
Illness
0.2 0.4 0.6 0.8
0.0057 0.0081 0.0132 0.0673
0.0449 0.0518 0.0664 0.1225
0.0460 0.0526 0.0680 0.1248
0.0059 0.0092 0.0151 0.0720
0.0468 0.0547 0.0702 0.1355
0.0494 0.0563 0.0727 0.1370
+3.51 +13.58 +14.39 +6.98
+4.23 +5.60 +5.72 +10.61
+7.39 +7.03 +6.91 +9.78
Weather
TABLE VI: Results of RDDMPI under point missing ratios (0.2, 0.4, 0.6, 0.8). This table compares the default setting (T = 50) with a reduced number of steps (T = 10), along with the relative degradation ∆ in %.
0.2 0.4 0.6 0.8
0.0246 0.0271 0.0313 0.0409
0.0246 0.0271 0.0317 0.0425
0.0327 0.0362 0.0423 0.0573
0.0233 0.0270 0.0321 0.0412
0.0258 0.0289 0.0336 0.0441
0.0361 0.0401 0.0462 0.0599
-5.28 -0.37 +2.56 +0.73
+4.88 +6.64 +5.99 +3.76
+10.40 +10.77 +9.22 +4.54
Models Metric
Future work should focus on improving inference efficiency through advanced sampling techniques, such as adaptive timestep selection or one-step diffusion models. These directions could further reduce computational cost while preserving the strong performance benefits of diffusion-based imputation. A PPENDIX E R ESIDUAL D ENOISING B LOCK A RCHITECTURE This appendix provides a detailed description of the residual denoising blocks used in the conditional noise prediction network ϵθ (Rtmi , t, c). Let V denote the number of variables, L the sequence length, dh the number of hidden channels, and Nblk the number of residual denoising blocks. The reliability-aware fusion and FiLM-based conditioning mechanisms described in Section IV-C produce the FiLM-conditioned (0) hidden representation Z̃t , which is the input to the residual denoising stack and initialized as Ut = Z̃t ∈ RV ×L×dh . This representation is passed sequentially through Nblk residual denoising blocks. For block b ∈ {1, . . . , Nblk }, we define (b,0)
Ut (b−1)
(b−1)
= Ut
,
where Ut is the output of the preceding block. Each block then applies diffusion-step conditioning, temporal attention, variable attention, side-information conditioning, a gated activation, and residual and skip projections. The complete sequence of operations is summarized in Algorithm 3 and described below.
18
A. Diffusion-Step Conditioning The discrete diffusion step t is represented by a fixed sinusoidal embedding, followed by two learned linear projections with SiLU activations as (0) et = SiLU We,2 SiLU We,1 et + be,1 + be,2 , (0)
where et is the fixed sinusoidal representation, ddiff denotes the dimension of the diffusion-step embedding, and et ∈ Rddiff is the resulting diffusion-step embedding. Within block b, this embedding is projected to the hidden dimension as (b)
dt
(b)
(b)
= Wdiff et + bdiff ∈ Rdh .
The projected embedding is broadcast across variables and timesteps and added to the incoming representation: (b,1) (b,0) (b) Ut = Ut + Broadcast dt . (b,1)
Thus, Ut
∈ RV ×L×dh contains both the hidden residual features and information about the current diffusion noise level.
B. Temporal and Cross-Variable Modeling The representation is first processed by a temporal transformer independently for each variable as (b,1) (b,2) (b) Ut [v, :, :] = Ttime Ut [v, :, :] , v = 1, . . . , V. For each variable, the corresponding L × dh slice is treated as a temporal sequence. Therefore, the temporal transformer models dependencies across timesteps without mixing the variable dimension. The resulting representation is subsequently processed by a variable transformer independently at each timestep as (b,3) (b,2) (b) Ut [:, ℓ, :] = Tvar Ut [:, ℓ, :] , ℓ = 1, . . . , L. At each timestep, the corresponding V × dh slice is treated as a sequence over variables. The variable transformer therefore models cross-variable interactions while preserving the temporal positions. Both transformer modules use multi-head self-attention, residual connections, layer normalization, and a position-wise feed-forward network with a GELU activation. For an input sequence X, an individual attention head is computed as ! (XWjQ )(XWjK )⊤ √ XWjV , headj (X) = softmax dk and the multi-head output is MHA(X) = Concat (head1 (X), . . . , headH (X)) W O , where H is the number of attention heads and dk is the dimension of each query and key head. C. Side-Information Conditioning Let S ∈ RV ×L×dside denote the side-information tensor provided to every residual denoising block. The side information consists of a temporal-position embedding, a learned variable embedding, the observation mask, and the learned reliability map. For each temporal coordinate τℓ , we construct a dtime -dimensional sinusoidal embedding. Its components are given by τℓ Etime (ℓ, 2j) = sin , 2j/dtime 10000τ ℓ Etime (ℓ, 2j + 1) = cos . 100002j/dtime The same temporal embedding is replicated across all variables, resulting in Etime ∈ RV ×L×dtime , where dtime denotes the dimension of the temporal-position embedding. Each variable v is assigned a learned embedding vector ev ∈ Rdvar , where dvar denotes the dimension of the variable embedding, through an embedding lookup table: ev = Embvar (v). This embedding is replicated across all timesteps to obtain Evar ∈ RV ×L×dvar . The observation mask M ∈ {0, 1}V ×L indicates which entries are available as conditioning information. The learned reliability map A ∈ [0, 1]V ×L is computed by the reliability-aware mechanism described in Section IV-C. After treating M and A as one-dimensional feature channels, the complete side-information tensor is S = Concat (Etime , Evar , M, A) ∈ RV ×L×dside , where dside denotes the total dimension of the side-information features and is given by dside = dtime + dvar + 2, with the additional two dimensions corresponding to the observation mask M and reliability map A. The reliability map therefore affects
19
the network twice: first during the reliability-weighted fusion of the baseline-completed signal and again as side information within every residual denoising block. The output of the variable transformer and the side-information tensor are projected independently to 2dh hidden dimensions. (b) (b) (b) Let ψmid and ψcond denote the learned pointwise projections in block b, implemented as 1 × 1 convolutions. Specifically, ψmid (b) maps dh to 2dh dimensions, whereas ψcond maps dside to 2dh dimensions. Their outputs are combined as (b,4) (b) (b,3) (b) = ψmid Ut + ψcond (S) , Ut (b,4)
where Ut ∈ RV ×L×2dh . Because both operators are pointwise, the same learned transformation is applied independently at every variable-timestep position, without mixing information across variables or timesteps. D. Gated Activation The conditioned representation is split equally along its hidden dimension as (b,4) (b,4) (b,4) Ut,gate , Ut,filter = Split Ut , where both components belong to RV ×L×dh . A gated activation is then applied as (b,5) (b,4) (b,4) Ut = σ Ut,gate ⊙ tanh Ut,filter . The sigmoid component controls how much information is propagated, while the hyperbolic-tangent component provides the (b,5) nonlinear candidate features. Their element-wise product gives Ut ∈ RV ×L×dh . E. Residual and Skip Outputs (b)
The gated representation is projected pointwise from dh to 2dh hidden dimensions. Let ψout denote the learned block-output (b) projection, implemented as a 1 × 1 convolution: ψout : RV ×L×dh → RV ×L×2dh . The projected representation is therefore (b,6) (b) (b,5) Ut = ψout Ut ∈ RV ×L×2dh . Because the projection is pointwise, it transforms the hidden features independently at every variable-timestep position, without mixing information across variables or timesteps. The resulting representation is divided equally along the hidden dimension into a residual update and a skip connection as (b) (b) (b,6) Ut,res , Ut,skip = Split Ut , √ (b) (b) where Ut,res , Ut,skip ∈ RV ×L×dh . The residual output is added to the block input and scaled by 1/ 2 (b,0)
(b)
Ut
=
Ut
(b)
+U √ t,res . 2
√ (b) The resulting tensor Ut becomes the input to the next residual denoising block. The factor 1/ 2 helps control the magnitude (b) of the hidden activations across the sequence of blocks. The skip output Ut,skip is retained for the final noise prediction path. F. Skip Aggregation and Noise Prediction After all residual denoising blocks have been evaluated, their skip outputs are summed and normalized as Ut,skip = √
N blk X 1 (b) U . Nblk b=1 t,skip
The aggregated skip representation is processed by two learned pointwise projections. The first projection, ψout,1 , is a 1 × 1 convolution that preserves the hidden dimension: ψout,1 : RV ×L×dh → RV ×L×dh . It is followed by a ReLU activation as Utout = ReLU (ψout,1 (Ut,skip )) . The second projection, ψout,2 , is a 1 × 1 convolution that maps the dh hidden features to one predicted noise value at every variable-timestep position: ψout,2 : RV ×L×dh → RV ×L . The final noise prediction is b ϵθ = ψout,2 Utout ∈ RV ×L . Both output projections operate independently at each variable-timestep position. The temporal and cross-variable interactions have already been incorporated by the transformer layers within the residual denoising blocks. Therefore, b ϵθ = ϵθ Rtmi , t, c has the same variable and temporal dimensions as the noisy residual and represents the Gaussian noise predicted for the reverse diffusion process.
20
Algorithm 3: Sequential Operations of the Residual Denoising Network (0) Input: FiLM-conditioned representation Ut = Zet ; diffusion step t; side-information tensor S; number of residual blocks Nblk Output: Predicted noise b ϵθ = ϵθ (Rtmi , t, c) Compute the diffusion-step embedding et (Subsection E-A) ; Initialize an empty collection of skip representations; for b = 1 to Nblk do (b,0) (b−1) Set the block input: Ut ← Ut ; (b,1) (b,0) (b) (b) Project and inject the diffusion-step embedding: Ut ← Ut + Broadcast(Wdiff et + bdiff ); (b,2) (b) (b,1) Apply temporal self-attention independently to each variable: Ut ← Ttime (Ut ); (b,3) (b) (b,2) Apply cross-variable self-attention independently at each timestep: Ut ← Tvar (Ut ); (b,4) (b) (b,3) (b) Apply the pointwise hidden and side-information projections: Ut ← ψmid (Ut ) + ψcond (S); (b,4) (b,4) (b,4) Split Ut into Ut,gate and Ut,filter ; (b,5) (b,4) (b,4) Apply the gated activation: Ut ← σ(Ut,gate ) ⊙ tanh(Ut,filter ); (b,6) (b) (b,5) Apply the pointwise output projection: Ut ← ψout (Ut ); (b,6) (b) (b) Split Ut into the residual update Ut,res and skip representation Ut,skip ; √ (b) (b,0) (b) Update the residual path: Ut ← (Ut + Ut,res )/ 2; (b) Store Ut,skip ; PNblk (b) √ Aggregate and normalize the skip representations: Ut,skip ← b=1 Ut,skip / Nblk ; Apply the first pointwise output projection: Utout ← ReLU(ψout,1 (Ut,skip )); Apply the final pointwise output projection: b ϵθ ← ψout,2 (Utout ); return b ϵθ ;
A PPENDIX F DATASET-S PECIFIC RDDMPI H YPERPARAMETERS Table VII reports the selected dataset-specific optimization and architectural hyperparameters used to train RDDMPI. These include the batch size, learning rate, number of residual denoising blocks, number of hidden channels, and number of attention heads. The deterministic T1 backbone is pretrained separately for each dataset and remains frozen during residual diffusion training. TABLE VII: Dataset-specific hyperparameters for RDDMPI. Nblk denotes the number of residual denoising blocks, and dh denotes the number of hidden channels. Dataset ETTh1 ETTh2 Exchange Illness Weather
Batch Size
Learning rate
Nblk
dh
Number of heads
16 16 16 64 16
10−3
4 4 4 4 4
128 128 128 64 64
8 8 8 8 8
10−3 10−3 10−3 10−4
A PPENDIX G F ULL R ESULTS This section reports the complete quantitative results underlying the aggregated comparisons presented in the main paper. For point missingness, results are shown separately for missing ratios of 0.2, 0.4, 0.6, and 0.8, rather than averaging across these ratios as in the main tables. We also report the results under block missingness. All experiments are conducted using five random seeds: 2, 102, 202, 302, and 402. For each seed, the models are trained independently and evaluated on the corresponding artificially masked test data. MAE, MSE, and CRPS are first computed for each seed, after which we report their mean and standard deviation across the five runs. The mean tables summarize average performance, while the standard deviation tables quantify sensitivity to model initialization, training stochasticity, and the generated missingness patterns. We provide both the disaggregated results and their variability for transparency, reproducibility, and a more complete assessment of model robustness.
21
TABLE VIII: Full results under point missing ratios (0.2, 0.4, 0.6, 0.8) across datasets. Best results are marked in bold, and second-best results are marked in underlined.
Weather
Illness
Exchange
ETTh2
ETTh1
Models RDDMPI (Ours) Metric MSE MAE
DLinear MSE MAE
ModernTCN MSE MAE
iTransformer MSE MAE
SAITS MSE MAE
ImputeFormer MSE MAE
TimesNet MSE MAE
T1 MSE
MAE
GP-VAE MSE MAE
CSDI MSE MAE
FGTI MSE MAE
0.0229 0.0337 0.0480 0.1145
0.0956 0.1123 0.1362 0.2020
0.1158 0.2320 0.0533 0.1594 0.0900 0.2033 0.0335 0.1180 0.0797 0.1697 0.0916 0.2075 0.0296 0.1113 0.3531 0.4602 0.0266 0.1016 0.0253 0.1009 0.1166 0.2232 0.0584 0.1632 0.1054 0.2172 0.0568 0.1482 0.1547 0.2232 0.0998 0.2123 0.0398 0.1277 0.5075 0.5365 0.0408 0.1225 0.0391 0.1220 0.2179 0.2927 0.0943 0.2020 0.1526 0.2533 0.1162 0.2097 0.2932 0.3138 0.1401 0.2434 0.0658 0.1632 0.6978 0.6235 0.0610 0.1512 0.0634 0.1537 0.4234 0.4143 0.2520 0.3198 0.3105 0.3537 0.3337 0.3691 0.6149 0.5109 0.3220 0.3637 0.1878 0.2694 0.9399 0.7237 0.1281 0.2128 0.1594 0.2315
Avg 0.0548
0.1365
0.2184 0.2906 0.1145 0.2111 0.1646 0.2569 0.1351 0.2113 0.2856 0.3044 0.1634 0.2567 0.0807 0.1679 0.6246 0.5860 0.0641 0.1470 0.0718 0.1520
0.0228 0.0305 0.0437 0.0797
0.0771 0.0924 0.1150 0.1649
0.0562 0.1578 0.0390 0.1269 0.0495 0.1463 0.1442 0.2685 0.1490 0.2279 0.0524 0.1542 0.0255 0.0954 0.6106 0.5741 0.0415 0.1008 0.0402 0.1149 0.0557 0.1560 0.0403 0.1286 0.0560 0.1553 0.1788 0.2950 0.2451 0.2845 0.0559 0.1577 0.0310 0.1057 0.8039 0.6763 0.0532 0.1206 0.0502 0.1313 0.0799 0.1871 0.0540 0.1501 0.0693 0.1744 0.3431 0.3897 0.6173 0.4458 0.0729 0.1787 0.0437 0.1283 1.4455 0.9409 0.0712 0.1475 0.0726 0.1594 0.1287 0.2396 0.0973 0.2043 0.1099 0.2219 1.0818 0.6869 1.6554 0.8247 0.1182 0.2291 0.0764 0.1764 2.7956 1.2937 0.1280 0.2089 0.2299 0.2728
Avg 0.0442
0.1123
0.0801 0.1851 0.0576 0.1524 0.0712 0.1745 0.4370 0.4100 0.6667 0.4457 0.0749 0.1799 0.0442 0.1264 1.4139 0.8712 0.0735 0.1444 0.0982 0.1696
0.2 0.4 0.6 0.8
0.0147 0.0171 0.0206 0.0288
0.0029 0.0364 0.0053 0.0512 0.0016 0.0262 0.1106 0.2863 0.0115 0.0401 0.0017 0.0258 0.0007 0.0151 0.5095 0.6226 0.0315 0.0889 0.0029 0.0250 0.0025 0.0317 0.0075 0.0602 0.0021 0.0278 0.1480 0.3225 0.0134 0.0493 0.0021 0.0266 0.0012 0.0174 0.5551 0.6523 0.0382 0.1008 0.0021 0.0279 0.0053 0.0453 0.0109 0.0717 0.0032 0.0328 0.2328 0.3882 0.0451 0.0951 0.0033 0.0321 0.0020 0.0217 0.6376 0.6981 0.0472 0.1147 0.0034 0.0355 0.0105 0.0680 0.0144 0.0823 0.0064 0.0488 0.4330 0.5289 0.2070 0.2549 0.0067 0.0496 0.0038 0.0332 0.8640 0.8046 0.0772 0.1496 0.0142 0.0753
Avg 0.0022
0.0203
0.0053 0.0453 0.0095 0.0664 0.0033 0.0339 0.2311 0.3815 0.0693 0.1098 0.0034 0.0335 0.0019 0.0218 0.6415 0.6944 0.0485 0.1135 0.0056 0.0409
0.0057 0.0081 0.0132 0.0673
0.0449 0.0518 0.0664 0.1225
0.1392 0.2691 0.0236 0.1078 0.0722 0.1724 0.1344 0.2193 0.2115 0.2625 0.0434 0.1468 0.0079 0.0561 0.2975 0.3760 0.0394 0.1055 0.0288 0.0918 0.0790 0.1842 0.0261 0.1079 0.0860 0.1985 0.1772 0.2448 0.2690 0.3014 0.0394 0.1390 0.0095 0.0597 0.5426 0.4767 0.0506 0.1160 0.0362 0.1045 0.1781 0.2829 0.0451 0.1447 0.0965 0.2093 0.2600 0.3040 0.3376 0.3512 0.0670 0.1818 0.0156 0.0759 0.7072 0.5588 0.0772 0.1525 0.0660 0.1452 0.4736 0.4823 0.2394 0.3339 0.2298 0.3001 0.5777 0.4493 0.5394 0.4542 0.2377 0.3541 0.0676 0.1479 0.9947 0.6593 0.1550 0.2178 0.2128 0.2441
Avg 0.0236
0.0714
0.2175 0.3046 0.0835 0.1736 0.1211 0.2201 0.2873 0.3044 0.3394 0.3423 0.0969 0.2054 0.0252 0.0849 0.6355 0.5177 0.0806 0.1479 0.0860 0.1464
0.0246 0.0271 0.0313 0.0409
0.0246 0.0271 0.0317 0.0425
0.0354 0.0745 0.0305 0.0632 0.0911 0.1436 0.0252 0.0308 0.0324 0.0355 0.0351 0.0729 0.0249 0.0375 0.0806 0.1582 0.0269 0.0253 0.0255 0.0252 0.0316 0.0585 0.0302 0.0563 0.0873 0.1431 0.0272 0.0359 0.0329 0.0392 0.0348 0.0685 0.0268 0.0407 0.1201 0.2048 0.0289 0.0277 0.0285 0.0278 0.0449 0.0863 0.0401 0.0773 0.0900 0.1443 0.0387 0.0557 0.0409 0.0514 0.0444 0.0880 0.0330 0.0529 0.1859 0.2753 0.0319 0.0317 0.0346 0.0323 0.0769 0.1334 0.0692 0.1228 0.0998 0.1497 0.0902 0.1363 0.0839 0.1114 0.0743 0.1317 0.0606 0.1016 0.3206 0.3927 0.0408 0.0415 0.0457 0.0428
Avg 0.0310
0.0315
0.0472 0.0882 0.0425 0.0799 0.0920 0.1452 0.0453 0.0647 0.0475 0.0594 0.0471 0.0903 0.0363 0.0582 0.1768 0.2578 0.0321 0.0316 0.0335 0.0320
0.2 0.4 0.6 0.8
0.2 0.4 0.6 0.8
0.2 0.4 0.6 0.8
0.2 0.4 0.6 0.8
0.0022 0.0015 0.0019 0.0030
TABLE IX: The standard deviation of Table VIII.
Weather
Illness
Exchange
ETTh2
ETTh1
Models RDDMPI (Ours) Metric MSE MAE 0.2 0.4 0.6 0.8
DLinear MSE MAE
ModernTCN MSE MAE
iTransformer MSE MAE
SAITS MSE MAE
ImputeFormer MSE MAE
TimesNet MSE MAE
T1 MSE
MAE
GP-VAE MSE MAE
CSDI MSE MAE
FGTI MSE MAE
0.0025 0.0028 0.0023 0.0081
0.0008 0.0015 0.0015 0.0051
0.0071 0.0047 0.0026 0.0026 0.0038 0.0023 0.0046 0.0068 0.0087 0.0080 0.0064 0.0081 0.0014 0.0017 0.0086 0.0057 0.0022 0.0031 0.0017 0.0017 0.0068 0.0039 0.0035 0.0028 0.0033 0.0032 0.0137 0.0125 0.0184 0.0116 0.0028 0.0036 0.0024 0.0017 0.0130 0.0065 0.0058 0.0044 0.0045 0.0035 0.0167 0.0065 0.0070 0.0045 0.0112 0.0046 0.0318 0.0241 0.0421 0.0281 0.0133 0.0051 0.0036 0.0023 0.0119 0.0057 0.0034 0.0040 0.0058 0.0048 0.0141 0.0071 0.0156 0.0089 0.0195 0.0064 0.0501 0.0367 0.0613 0.0457 0.0206 0.0120 0.0074 0.0066 0.0136 0.0093 0.0218 0.0106 0.0171 0.0096
Avg 0.0039
0.0022
0.0111 0.0056 0.0072 0.0047 0.0095 0.0041 0.0251 0.0200 0.0326 0.0233 0.0108 0.0072 0.0037 0.0031 0.0118 0.0068 0.0083 0.0055 0.0073 0.0049
0.2 0.4 0.6 0.8
0.0012 0.0022 0.0044 0.0081
0.0034 0.0030 0.0049 0.0079
0.0030 0.0024 0.0022 0.0027 0.0027 0.0023 0.0071 0.0078 0.0126 0.0056 0.0029 0.0026 0.0012 0.0015 0.0730 0.0339 0.0052 0.0063 0.0037 0.0054 0.0014 0.0013 0.0003 0.0006 0.0003 0.0005 0.0163 0.0152 0.0429 0.0221 0.0011 0.0011 0.0012 0.0011 0.0784 0.0389 0.0085 0.0084 0.0016 0.0029 0.0010 0.0016 0.0008 0.0018 0.0016 0.0019 0.0937 0.0546 0.1761 0.0650 0.0014 0.0028 0.0015 0.0018 0.1417 0.0374 0.0112 0.0093 0.0035 0.0027 0.0035 0.0024 0.0041 0.0034 0.0035 0.0028 0.2958 0.1241 0.3109 0.0934 0.0021 0.0020 0.0041 0.0037 0.0974 0.0177 0.0285 0.0192 0.0792 0.0425
Avg 0.0040
0.0048
0.0022 0.0019 0.0019 0.0021 0.0020 0.0019 0.1032 0.0504 0.1356 0.0465 0.0019 0.0021 0.0020 0.0020 0.0976 0.0320 0.0133 0.0108 0.0220 0.0134
0.2 0.4 0.6 0.8
0.0019 0.0008 0.0005 0.0001
0.0003 0.0005 0.0003 0.0002
0.0002 0.0005 0.0004 0.0010 0.0002 0.0002 0.0154 0.0221 0.0084 0.0035 0.0002 0.0009 0.0002 0.0004 0.0400 0.0259 0.0383 0.0789 0.0019 0.0009 0.0008 0.0003 0.0008 0.0005 0.0009 0.0003 0.0056 0.0082 0.0054 0.0117 0.0007 0.0007 0.0008 0.0003 0.0300 0.0197 0.0548 0.0970 0.0011 0.0030 0.0006 0.0009 0.0006 0.0005 0.0006 0.0006 0.0106 0.0130 0.0417 0.0532 0.0004 0.0007 0.0006 0.0004 0.0373 0.0208 0.0698 0.1167 0.0006 0.0052 0.0004 0.0007 0.0003 0.0005 0.0004 0.0011 0.0229 0.0210 0.2059 0.1700 0.0005 0.0010 0.0006 0.0012 0.0497 0.0251 0.1148 0.1518 0.0059 0.0203
Avg 0.0008
0.0003
0.0005 0.0006 0.0005 0.0007 0.0005 0.0006 0.0136 0.0161 0.0653 0.0596 0.0005 0.0008 0.0005 0.0006 0.0392 0.0229 0.0694 0.1111 0.0024 0.0073
0.2 0.4 0.6 0.8
0.0007 0.0014 0.0020 0.0389
0.0029 0.0028 0.0036 0.0138
0.0159 0.0107 0.0043 0.0120 0.0200 0.0249 0.0327 0.0314 0.0447 0.0276 0.0052 0.0132 0.0018 0.0061 0.0642 0.0333 0.0075 0.0088 0.0045 0.0104 0.0098 0.0091 0.0031 0.0047 0.0077 0.0049 0.0385 0.0121 0.0768 0.0367 0.0050 0.0073 0.0013 0.0016 0.0805 0.0268 0.0090 0.0083 0.0088 0.0117 0.0514 0.0329 0.0089 0.0134 0.0228 0.0181 0.0445 0.0182 0.0384 0.0161 0.0135 0.0200 0.0022 0.0028 0.0734 0.0267 0.0153 0.0127 0.0145 0.0169 0.0858 0.0268 0.0644 0.0334 0.0728 0.0316 0.1379 0.0425 0.1117 0.0419 0.0426 0.0230 0.0387 0.0271 0.0599 0.0184 0.0290 0.0131 0.0891 0.0312
Avg 0.0108
0.0058
0.0407 0.0199 0.0202 0.0159 0.0308 0.0199 0.0634 0.0260 0.0679 0.0306 0.0166 0.0159 0.0110 0.0094 0.0695 0.0263 0.0152 0.0107 0.0292 0.0175
0.2 0.4 0.6 0.8
0.0018 0.0020 0.0020 0.0014
0.0003 0.0006 0.0014 0.0022
0.0026 0.0013 0.0010 0.0055 0.0040 0.0013 0.0022 0.0012 0.0039 0.0034 0.0014 0.0037 0.0017 0.0025 0.0079 0.0096 0.0023 0.0006 0.0014 0.0022 0.0007 0.0004 0.0006 0.0012 0.0018 0.0005 0.0016 0.0018 0.0039 0.0036 0.0011 0.0026 0.0022 0.0032 0.0081 0.0087 0.0022 0.0010 0.0027 0.0028 0.0015 0.0003 0.0025 0.0040 0.0019 0.0002 0.0022 0.0033 0.0004 0.0047 0.0005 0.0028 0.0027 0.0058 0.0101 0.0109 0.0024 0.0012 0.0027 0.0032 0.0004 0.0005 0.0020 0.0046 0.0042 0.0013 0.0079 0.0116 0.0161 0.0219 0.0011 0.0015 0.0043 0.0106 0.0258 0.0244 0.0007 0.0014 0.0056 0.0048
Avg 0.0018
0.0011
0.0013 0.0006 0.0015 0.0038 0.0030 0.0008 0.0035 0.0045 0.0061 0.0084 0.0010 0.0026 0.0027 0.0055 0.0130 0.0134 0.0019 0.0010 0.0031 0.0033
TABLE X: Full results under the block missing scenario across datasets. Best results are marked in bold, and second-best results are marked in underlined. Dataset
RDDMPI (Ours) MSE MAE
ETTh1 0.0168 ETTh2 0.0203 Exchange 0.0043 Illness 0.0390 Weather 0.0233
0.0851 0.0744 0.0181 0.1282 0.0254
DLinear MSE MAE
ModernTCN MSE MAE
iTransformer MSE MAE
SAITS MSE MAE
ImputeFormer MSE MAE
TimesNet MSE MAE
T1 MSE
MAE
GP-VAE MSE MAE
CSDI MSE MAE
FGTI MSE MAE
0.1728 0.2844 0.0592 0.1735 0.1016 0.2097 0.0257 0.1075 0.0594 0.1543 0.0920 0.2117 0.0258 0.1077 0.2842 0.4231 0.0179 0.0896 0.0170 0.0877 0.0770 0.1890 0.0477 0.1414 0.0546 0.1552 0.1471 0.2711 0.2654 0.2672 0.0533 0.1592 0.0292 0.1014 0.6001 0.5663 0.0561 0.1070 0.0948 0.1321 0.0063 0.0557 0.0054 0.0498 0.0034 0.0336 0.1864 0.3330 0.1217 0.1233 0.0036 0.0362 0.0115 0.0312 0.5100 0.6290 0.2237 0.1458 0.0177 0.0457 0.2603 0.3691 0.1521 0.2647 0.3345 0.3598 0.1475 0.2373 0.2697 0.3100 0.1368 0.2535 0.0874 0.1792 0.3270 0.4122 0.0977 0.1895 0.0516 0.1283 0.0495 0.1053 0.0371 0.0827 0.1003 0.1480 0.0251 0.0328 0.0401 0.0461 0.0390 0.0855 0.0248 0.0400 0.0627 0.1349 0.0234 0.0254 0.0236 0.0256
22
TABLE XI: Standard deviations corresponding to the block missing results reported in Table X. Dataset
RDDMPI (Ours) MSE MAE
ETTh1 0.0020 ETTh2 0.0024 Exchange 0.0074 Illness 0.0384 Weather 0.0041
0.0025 0.0043 0.0091 0.0675 0.0024
DLinear MSE MAE
ModernTCN MSE MAE
iTransformer MSE MAE
SAITS MSE MAE
ImputeFormer MSE MAE
TimesNet MSE MAE
T1 MSE
MAE
GP-VAE MSE MAE
CSDI MSE MAE
FGTI MSE MAE
0.0050 0.0024 0.0066 0.0060 0.0345 0.0129 0.0061 0.0100 0.0155 0.0130 0.0085 0.0077 0.0051 0.0049 0.0114 0.0067 0.0022 0.0046 0.0020 0.0041 0.0099 0.0143 0.0117 0.0163 0.0102 0.0132 0.0168 0.0148 0.1646 0.0641 0.0070 0.0109 0.0078 0.0129 0.0693 0.0356 0.0324 0.0202 0.1091 0.0341 0.0012 0.0027 0.0020 0.0060 0.0018 0.0049 0.0852 0.0602 0.1218 0.0673 0.0014 0.0042 0.0192 0.0204 0.0608 0.0364 0.3921 0.1583 0.0322 0.0363 0.1112 0.0515 0.0985 0.0867 0.2518 0.1288 0.0663 0.0893 0.1174 0.0934 0.1194 0.1058 0.0885 0.1076 0.1543 0.1103 0.0984 0.1226 0.0719 0.1089 0.0062 0.0041 0.0050 0.0040 0.0140 0.0059 0.0061 0.0040 0.0023 0.0025 0.0047 0.0020 0.0052 0.0031 0.0084 0.0102 0.0051 0.0026 0.0032 0.0039
TABLE XII: Full CRPS comparison of RDDMPI, GP-VAE, CSDI, and FGTI under point missing and block missing scenarios across datasets. For point missingness, results are reported at missing ratios 0.2, 0.4, 0.6, and 0.8. Best results are marked in bold, and second-best results are marked in underlined. Lower is better. (a) Point missing (b) Block missing RDDMPI CRPS
GP-VAE CRPS
CSDI CRPS
0.2 0.4 0.6 0.8
0.0924 0.1080 0.1310 0.1940
0.5829 0.6739 0.7856 0.9080
0.0970 0.0966 0.1165 0.1161 0.1446 0.1469 0.2044 0.2222
Avg
0.1313
0.7376
0.1406 0.1454
0.2 0.4 0.6 0.8
0.0438 0.0528 0.0657 0.0941
0.4202 0.4967 0.6921 0.9488
0.0570 0.0682 0.0848 0.1203
Avg
0.0641
0.6395
0.0826 0.0979
0.2 0.4 0.6 0.8
0.0131 0.0151 0.0181 0.0252
0.7039 0.7360 0.7887 0.9106
0.0766 0.0869 0.0998 0.1328
0.0216 0.0239 0.0304 0.0658
Avg
0.0179
0.7848
0.0990
0.0354
0.2 0.4 0.6 0.8
0.0460 0.0526 0.0680 0.1248
0.4830 0.6118 0.7244 0.8479
0.1111 0.0917 0.1209 0.1026 0.1617 0.1412 0.2296 0.2427
Avg
0.0728
0.6668
0.1558
0.2 0.4 0.6 0.8
0.0327 0.0362 0.0423 0.0573
0.2764 0.3575 0.4809 0.6855
0.0335 0.0365 0.0368 0.0400 0.0421 0.0468 0.0555 0.0626
Avg
0.0420
0.4501
0.0421 0.0465
Weather
Illness
Exchange
ETTh2
ETTh1
Models Metric
FGTI CRPS
0.0657 0.0756 0.0928 0.1576
0.1445
Dataset ETTh1 ETTh2 Exchange Illness Weather
RDDMPI GP-VAE CSDI 0.0811 0.0412 0.0153 0.1238 0.0338
0.5298 0.4076 0.6903 0.4793 0.2383
0.0867 0.0596 0.1238 0.1873 0.0336
FGTI 0.0832 0.0751 0.0391 0.1173 0.0348
23
TABLE XIII: Standard deviations of the CRPS results reported in Table XII. (a) Point missing (b) Block missing RDDMPI CRPS
GP-VAE CRPS
CSDI CRPS
FGTI CRPS
0.2 0.4 0.6 0.8
0.0012 0.0017 0.0017 0.0064
0.0073 0.0050 0.0056 0.0098
0.0023 0.0046 0.0045 0.0103
0.0021 0.0031 0.0044 0.0091
Avg
0.0027
0.0069
0.0054
0.0047
0.2 0.4 0.6 0.8
0.0022 0.0019 0.0027 0.0050
0.0261 0.0269 0.0264 0.0134
0.0037 0.0044 0.0061 0.0116
0.0030 0.0019 0.0010 0.0229
Avg
0.0030
0.0232
0.0065
0.0072
0.2 0.4 0.6 0.8
0.0004 0.0005 0.0003 0.0002
0.0253 0.0220 0.0220 0.0298
0.0652 0.0798 0.0944 0.1247
0.0011 0.0025 0.0043 0.0182
Avg
0.0004
0.0248
0.0910
0.0065
0.2 0.4 0.6 0.8
0.0038 0.0028 0.0030 0.0149
0.0058 0.0250 0.0206 0.0195
0.0099 0.0077 0.0094 0.0155
0.0079 0.0128 0.0171 0.0380
Avg
0.0061
0.0177
0.0106
0.0190
0.2 0.4 0.6 0.8
0.0005 0.0008 0.0020 0.0034
0.0165 0.0152 0.0186 0.0421
0.0008 0.0012 0.0017 0.0020
0.0032 0.0039 0.0043 0.0066
Avg
0.0017
0.0231
0.0014
0.0045
Weather
Illness
Exchange
ETTh2
ETTh1
Models Metric
Dataset ETTh1 ETTh2 Exchange Illness Weather
RDDMPI GP-VAE CSDI 0.0015 0.0021 0.0070 0.0546 0.0032
0.0080 0.0326 0.0188 0.0783 0.0194
0.0038 0.0107 0.1331 0.1111 0.0029
FGTI 0.0027 0.0185 0.0320 0.0969 0.0058
24
A PPENDIX H A DDITIONAL Q UALITATIVE R ESULTS We provide additional qualitative imputation visualizations under varying point wise missing ratios (20%, 40%, 60%, and 80%) and block missingness for ETTh2. These figures illustrate median predictions and 90% predictive intervals, highlighting reconstruction accuracy and uncertainty calibration.
Fig. 4: Visualization of probabilistic imputation results on ETTh2 (20% missingness). The results are for a time series sample with all 7 features. The median and 90% predictive intervals are shown in green, observed points in red and targets in blue.
25
Fig. 5: Visualization of probabilistic imputation results on ETTh2 (40% missingness). The results are for a time series sample with all 7 features. The median and 90% predictive intervals are shown in green, observed points in red and targets in blue.
26
Fig. 6: Visualization of probabilistic imputation results on ETTh2 (60% missingness). The results are for a time series sample with all 7 features. The median and 90% predictive intervals are shown in green, observed points in red and targets in blue.
27
Fig. 7: Visualization of probabilistic imputation results on ETTh2 (80% missingness). The results are for a time series sample with all 7 features. The median and 90% predictive intervals are shown in green, observed points in red and targets in blue.
28
Fig. 8: Visualization of probabilistic imputation results on ETTh2 (block missingness). The results are for a time series sample with all 7 features. The median and 90% predictive intervals are shown in green, observed points in red and targets in blue.