MissScore: High-Order Score Estimation in the Presence of Missing Data
Authors: Wenqin Liu
Haoze Hou
Erdun Gao
Biwei Huang
Qiuhong Ke
Howard Bondell
Mingming Gong
Venue: ICML 2025 poster (2025)
Keywords: Generative models; Score matching; Tabular data
Primary area: deep_learning->generative_models_and_autoencoders

Abstract: Score-based generative models are essential in various machine learning applications, with strong capabilities in generation quality. In particular, high-order derivatives (scores) of data density offer deep insights into data distributions, building on the proven effectiveness of first-order scores for modeling and generating synthetic data, unlocking new possibilities for applications. However, learning them typically requires complete data, which is often unavailable in domains such as healthcare and finance due to data corruption, acquisition constraints, or incomplete records. To tackle this challenge, we introduce MissScore, a novel framework for estimating high-order scores in the presence of missing data. We derive objective functions for estimating high-order scores under different missing data mechanisms and propose a new algorithm specifically designed to handle missing data effectively. Our empirical results demonstrate that MissScore accurately and efficiently learns the high-order scores from incomplete data and generates high-quality samples, resulting in strong performance across a range of downstream tasks.

MissScore: High-Order Score Estimation in the Presence of Missing Data
Wengin Liu'
Haoze Hou?
Erdun Gao
Biwei Huang"
Qiuhong Ke Howard Bondell' Mingming Gong 1
Abstract
medical imaging (Song et al., 2021; Chung & Ye, 2022). De- noising Score Matching (DSM) (Vincent, 201 1), an efficient method for estimating the score of the data density from core-based
Score-based generative models are essential in various machine learning applications, with strong capabilities in generation quality, In partic- ular, high-order derivatives (scores) of data den- sity offer deep insights into data distributions, building on the proven effectiveness of first-order scores for modeling and generating synthetic data, unlocking new possibilities for applications. How- ever, learning them typically requires complete data, which is often unavailable in domains such as healthcare and finance due to data corruption, acquisition constraints, or incomplete records. To tackle this challenge, we introduce MissScore, a novel framework for estimating high-order scores in the presence of missing data. We derive ob- jective functions for estimating high-order scores under different missing data mechanisms and pro- pose a new algorithm specifically designed to han- dle missing data effectively. Our empirical re- sults demonstrate that MissScore accurately and efficiently learns the high-order scores from in- complete data and generates high-quality samples, resulting in strong performance across range of downstream tasks.
samples, has become widely used in training generative models (Ho et al., 2020; Song & Ermon, 2020). Beyond the first-order score, high-order derivatives of the data density, which we refer to as high-order scores, offer more refined local approximations of the data distribution, such as its curvature, and enable new model capabilities. For instance, they can improve the mixing speed of sam- pling methods and provide insights into quantifying the uncertainty in denoising problem (Dalalyan & Karagulyan, 2019; Sabanis & Zhang, 2019; Meng et al., 2021). Addition- ally, Lu et al. (2022) empirically demonstrated that incorpo- rating high-order score matching improves the likelihood of score-based diffusion ordinary differential equations on both synthetic and real data, while maintaining high-quality generation. Furthermore, high-order scores have been uti- lized in recovering causal structures (Rolland et al., 2022; Sanchez et al., 2022; Liu et al., 2024).
Despite their promise, learning high-order scores usually requires training the model on complete data (Meng et al., 2021; Lu et al., 2022). However, in many real-world see- narios, such as healthcare, finance, and social networks, data often contain missing values due to privacy constraints or high sampling costs (Rubin, 1976; Shpitser, 2016). A straightforward approach to handling missing data is to im- pute the missing values and train the model on the imputed dataset. However, imputation methods often compromise data quality, and potentially leading biased results and sig- nificantly degrading the performance of downstream tasks (Ouyang et al., 2023). Furthermore, these methods fail to capture the inherent uncertainty associated with missing data, producing a distribution of imputed values that poorly represents the true underlying data distribution (Van Buuren & Groothuis -otbbuii--udsoooom.. 2011; Gondara & Wang, 2018; Wang et al., 2021). Some alternative approaches, such as using generative adversarial networks (GANs) vari- ational auto-encoders (VAEs) to directly approximate the data generation model from incomplete data (Li et al.. 2019; Gain & Shpitser, 2018), 8), require training additional networks, which can be computationally expensive and may also re- sult in model inconsistency. To address this issue, some works propose explicitly constrain the learning objective
1. Introduction
The first-order derivative of the log data density, also known as (Stein) score (Liu et al., 2016), plays an important role in various machine learning applications, including data synthesis (Song & Ermon, 2019; 2020; Kim et al., 2022), super-resolution (Li et al., 2022), and inverse problems in
'School of Mathematics and Statistics, The University of Mel- bourne, Australia IInstitute of Statistics and Big Data, Renmin Uni- versity of China, China AAustralian Institute for Machine Leamingg The University of Adelaide, Australia 'Halcroglu Data Science Institute, UC San Diego United States 'Faculty of Information Technology, Monash University. Australia 'Department of Ma- chine Learning, Mohamed bin Zayed University of Artificial In- telligence, United Arab Emirates. Correspondence to: Mingming Gong <mingming, gong @ uni.elb.edu..u.>.
Proceedings of the nnnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 20255 by the author(s).

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
within the model, which can enhance performance (Städler & Bühlmann, 2012; Gao al., 2022). These studies em- phasize the need for approaches that can directly and more efficiently handle missing data.
to capture data distributions more effectively. However, these methods are limited by their inability to model the inherent uncertainty in missing data and often struggle to outperform traditional methods in survey data (Wang et al., 2021).
Back to score estimation from incomplete data, Ouyang et al. (2023) adopts a similar way by introducing a diffusion- based framework that learns the first-order score directly from incomplete data. In principle, high-order scores could be estimated from learned first-order score model trained on incomplete data using automatic differentiation. How- ever, this approach becomes computationally impractical for high-dimensional data and large model sizes, particularly when using deep neural networks based models (Meng et al., 2021). Furthermore, automatic differentiation introduces additional estimation errors, as small errors does not always lead to small estimation error for high-order scores. More- over, methods based on GANs (Goodfellow et al., 2020) or VAES (Kingma, 2013) do not inherently capture score information, regardless of whether data is missing or com- plete. In contrast, score-based models naturally integrate this information (Li et al., 2019; Ho et al., 2020; Gain & Sh- pitser, 2018). These limitations highlight the need of using score-based models when estimating high-order scores in the presence of missing data.
Recent advances in generative modeling have explored learn- ing directly from incomplete data, moving beyond explicit imputations. Variational autoencoder-basee approaches han- dle complex tabular data through tailored likelihoods, hier- archical structures, or importance weighting (Nazabal et al., 2020; Mattei & Frellsen, 2019; Ma et al., 2020; Peis et al., 2022). Diffusion-based models offer an alternative by learn- ing score functions over partially observed data, supporting both data synthesis and imputation (Ouyang et al., 2023; Zhang et al., 2024; Chen et al., 2024; Zheng & Charoen- phakdee, 2022; Kotelnikov et al., 2023). However, many of these methods rely on auxiliary networks, make strong assumptions about the missingness mechanism, or are re- stricted to first-order scores. In contrast, our approach di- rectly estimates high-order score functions from incomplete data, enabling richer local information and broader applica- bility.
Score Matching. Score matching estimates the gradient of the log-density, known as the score function, of a data distribution (Hyvärinen & Dayan, 2005). It is particularly effective for complex, high-dimensional density models with intractable partition functions. Among its variants, DSM (Vincent, 2011) has emerged as a foundational method, learning the score of perturbed distribution by minimizing regression loss. While first-order scores are widely used across applications (Song & Ermon, 2019; Li et al., 2022), higher-order scores offer richer insights by capturing the curvature of the data distribution (Dalalyan & Karagulyan, 2019). These high-order derivatives
Contributions. In this work, we introduce MissScore, a novel score-based framework for learning high-order scores in the presence of missing data. We derive objective func- tions for estimating these scores under different missing data mechanisms, using DSM to recover the true score function. While our framework is general and applicable scores of any order, we focus second-order scores (the Hessian of the log density) in our experiments. Our results show that MissScore efficiently and accurately approximates high- order scores with missing data. In addition, we demonstrate that our model improves both sampling speed and data qual- ity in data generation tasks, with the quality of the generated samples validated across several downstream tasks.
2019; Sabanis & Zhang,2019), improve sampling efficiency and uncertainty quantification (Meng et al., 2021). However, learning high-order scores is computationally expensive and prone to inaccuracies, particularly in high-dimensional settings. Recent work by Meng et al. (2021) addresses this challenge by leveraging Tweedie's 's formula to estimate high-order derivatives elfi- ciently. Lu et al. (2022) further demonstrated that high-order score matching improves the likelihood of diffusion- models while maintaining high sample quality. Missecore builds on these insights, extending DSM to high-order score estimation in the presence of missing data, ensuring both computational efficiency and robustness.
2. Related Work
Missing Data Problem. Learning from incomplete observa- tions is a common challenge in real-world datasets, arising from factors such as data corruption or incomplete records (Rubin, 1976; Little& & Rubin, 2019). This issue has been extensively studied, with imputation being primary solu- tion (Poulos& & Valle, 2018; 8; Jäger et al., 2021; Shadbahr et al., 2023; Paterakis et al., 2024). Traditional imputa- tion methods, such as mean or median substitution, tend to compromise data diversity and introduce biases in down- stream tasks. More advanced techniques, including machine learning models and deep generative approaches like GAIN, MICE, and MIDA (Muzellec al., 2020; Van Buuren & Groothuis-Oudshoomm 2011; Gondara & Wang, 2018), aim
3. Estimating High-Order Scores with Missing Data
In this section, we provide an overview of high-order DSM and introduce our approach for handling missing data.
2

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
3.1. Background High-Order DSM
Theorem 3.2. (Meng al., 2021) Given score functions
Si (X), (8), ..0)),,.,.)1((88))))) SX), k-th order score model sx({;; 0) can be obtained by optimizing the following objective:
Consider a data distribution Paan(X) and a model distribu- tion p(x; 0) over R The first order score is the gradi- ent of log Pdama (x) with respect to x, denoted s (x) = Vx log (x). Correspondingly, the score function of p(x; 0) is denoted as 8) (x; 0) = Vx p(;; 0). In DSM, instead of directly estimating the score function from the original data, the method works by introducing noise from a predefined noise distribution qo (x(x) into the data. The objective is then estimate the score of the perturbed data distribution qo (x) ʃ go (x xx))Paan(x()
0' arg min Bpan (x) Bog (8/x)
|| 18*x x (8) (x;;; 0)]112]
where Sk is a polynomial of (8)., sx(()) such that:
To achieve so, DSM minimizes the following objective function,
fu (3, 81 (8)...... sx(x))
8 + (21(81(8)
ik = 1,
(OsM(()= - 2 E (x) Eg (81x))
1
||s| (x: 0)
0² 2 fk + 0² fk 1 (s1(x)+ +
if 2 2.
2
(3)
Vx Vx 2o (x|x)
(1)
We have sn (x; 0¹) = (x) for almost all X.
It has been shown that minimizing Eq. (1) is equivalent to minimizing the score matching loss between $ (x; 0) and 81 (x) under certain regularity conditions (Vincent, 20. When the noise distribution go (x(xx) is Gaussian, i.., (N (²)x) I)),, the objective simplifies to
3.2. High-Order DSM with Missing Data
Let x = (I1, (21, ,2....,,za E Rd be random vector sampled from an unknown data distribution Pdata (x), and m = (m1, m,...., ma) E (0, 1 be a binary mask where m, = indicates that x; is missing, and m, = 0 indi- cates that H, is observed. The observed data can be ex- press Xats = x O (1 m) na @ m, where O denotes element-wise multiplication, and na indicates the missing value. Perturbing x with Gaussian noise results in X x Nix, and the corresponding conditional density function is qa (x(x) := (2ma²) a exp( (K-x) where ø is a pre- specified constant.
LDSM (0) = 2 (x) By (81))))
0)+
(81x)
o
The learned score function implicitly learns how to "de- noise" the perturbed data x. guiding it back toward the true data distribution through the optimization of Eq. (2). By focusing on estimating the score of the noise-perturhed distribution qo (x) instead of the original data distribution (x).. DSM offers more efficient approach to score esti- mation compared to other techniques (Hyvärinen & Dayan, 2005; Song et al., 2020). Meng et al. (2021) provide derivation of DSM using Tweedie's formula (Efron, 201 1), and they generalize this approach to incorporate high-order moments of x based on x, allowing them to develop an objective function for learning high-order scores.
2
Rubin (1976) categorized missing data mechanisms into three types based on the dependency between the missing indicator (m) and the complete data (x): (1) Missing Com- pletely at Random (MCAR), where m is independent of the complete data x; (2) Missing at Random (MAR), where m depends only on the observed values Xots: and (3) Missing Not at Random (MNAR), where m depends on both the observed values Xabs and the missing values. Many previous efforts have concentrated on addressing the complexities of MNAR cases; however, M(C)AR cases have received less attention, as they allow for recovering the true distribution without the need for additional assumptions (Wang et al., 2019; Yang et al., 2019; Little& & Rubin, 2019).
Theorem 3.1. (Meng et al., 2021) X[
[x 81,....8.) where 0x E DDD denotes n-fold ten- sor multiplications,
is a polynomial .,,,,,,.s,)), (x), and su (x) represents the k-th order of x 8i (x), ...,,
(x),
score of qo (x) = Pdata (x) (x(x)dxx
The MCAR mechanism simplifies modeling by climinat- ing bias introduced by missingness but is rarely realistic, MAR, the other hand, is more applicable in real-world settings, as it allows missingness to depend on observed data. This dependency, however, can introduce bias in esti- mates if the missing data mechanism is not accounted for, making the problem considerably more challenging than in the MCAR framework. Inverse Probability Weighting
Theorem 3.1 shows that there exists an equality between high-order moments of the posterior distribution of x given x and high-order scores with respect to x. Leveraging The- orem 3.1 and the least squares estimation of [ the objectives for approximating the k-th order scores sy (xx) can be constructed as follows.
3

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
(IPW) has proven useful in correcting bias by assigning
If there exists unique O exists such that S2(X) <x 0²), then 0* =anssnn
weights to observed data points inversely proportional to
their probability of being observed (Wooldridge, 2007; Sea-
0
man & White, 2013). This ensures that observations more likely to be missing receive higher weights, enabling un- biased estimates that reflect the full data distribution and mitigating the bias introduced by the MAR mechanism.
We can conclude that the global optimum of the second- order DSM with a missing mask corresponds to the oracle second-order score. A detailed proof is provided in Ap- pendix A.2. We now extend to any desired order under M(C)AR mechanisms. Theorem 3.5 indicates that k-th or- der DSM with missing mask can learn the oracle k-th order score and the proof is provided in Appendix A.3.
The following theorem presents our first theoretical result, demonstrating that DSM with a missing indicator cover the oracle score under M(C)AR mechanisms, defined as the gradient of log Pdata (x) with respect to x.
Theorem 3.5. Given score functions 8y (i))).... : _k-- (5). and the missing probability of every element lies between O and 1, which is p(m; 1) E [0, 1) for all i E {1,2.... d]. If we correctly model the k-th order derivative sk (X), there exists O' such that 8& (x, 0*) = 8% (8), then
Theorem 3.3. Let the missing mechanism of be either MCAR or MAR, with the missing probability of every ele- ment lying between O and I, i.e., p(m; 1) E [0, 1) for all d EE (1.2..... d}. Define the objective function:
0r arg min Ex
||| 8* x - fx 51 S2 84-1).(...,,. (x),
(0) = Ex/x.. [[()((x)0)
0
Sx (xxx0)) 3 8 //wl// wwl?
+
2z
(x x)} 2 will?
where = 1 m if the missing mechanism of x is MCAR, and w P|mx=00x=x] if the missing mechanism of x is MAR.
1-m
where the weight w is defined as:
m.
MCAR,
W1=
m m
if MAR.
4. Training Score Models by High-Order DSM with Missing data
If there exists a unique 0 such that th (xx) = =s1(x; 0"), then
In this section, describe the training process for high- order score models in the presence of missing data and evaluate their empirical performance. While our analysis specifically focuses on the first- and second-order scores, the approach can be applied to any order of scores.
0" = arg minJuss
It follows that the global optimum of DSM with missing mask aligns with the oracle score under M(C)AR missing mechanisms. The detailed proof is provided in Appendix A. Extending this to the second-order score model, and building on Theorem 3.2. the second-order DSM with a missing mask can successfully recover the oracle second- order score, which corresponds the Hessian of log Pauta (x) with respect to
As derived in Theorem 3.3, the first-order score model si (X:; 0) learned by minimizing the following objective function using the observed data,
LDSM ( = Ex Nens. m
[| (Xabai 0)
Theorem 3.4. Suppose the first-order score sı (x) is given, and the missing mechanism of x is either MCAR or MAR, with p(m; 1) E [0, for all f1,2.....d}; Define the objective function:
1 + o² (Xots Xobs)) } 3 W1
Similarly, the second-order score model s2 0) is learned using the following objective, as derived in Theorem 3.4,
(O) = Exx E |x
[x2(((; 0) + (X)s] (x)
I -- zz
w2 12
(D)SM(0) = E Ex Nots. [1 (Xobs; 0) + (Xooss )s} (Xoba)
x2
T I Zabr Lons
where z x-x and W2 is defined as:
W2
(5)
(1 (1-m) m) -m)
if MCAR,
However, training the second-order score model, (xx 0). requires knowledge of the first-order score $s (x). Therefore,
W2=
(1-m)11m)
if MAR.
/f/mm =0
4

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
adopt a multi-task objective to train both 8 (x: 0) and 82 (x(x;; 0) simultaneously,
setting, 82 (Xobs: m) (1 m) m) while for MAR,
82 (Nobs: m) (1 (1-m)[1mm)
/Plmm [=mmm==0x=x =0|x=x] The formal analysis of
Gjoin (0) = (0) + w WLDgSM
(6)
variance reduction is provided in sections A.4 and A.5 of
the Appendix.
where w E IR+ is a tunable coefficient. P/m = 0|x = Xobs in MCAR and Pimm" = 0|x = Xots] in MAR are esti- mated using logistic regression models, where the response variable indicates whether the data missing or observed, and the predictors are the observed variables. In the experi- ments, missing values are handled by replacing them with O for continuous variables and creating new category for discrete variables. One-hot encoding is then applied to dis- crete variables. Element-wist multiplication with the mask naturally mitigates the impact of replacing missing values with zeros when computing the objective. The algorithm is provided in Appendix B.
We perform an empirical analysis to assess the impact of VR on training score models with DSM and DySM using incomplete data. The full data is generated from a 2-d Gaussian distribution and we simulate the incomplete data under MCAR, training $ (Katei 0) and 8 (Katta 0) using a joint learning objective Eq. (6). Using sample size of 1000 and missing ratio of 0.3, In Figure 1, we compare the estimated score and Hessian for the first dimension against the ground truth at noise levels = {0.1, 0.001], J. both with and without VR. The results indicate that VR is essential for accurate estimation in both DSM and D₂SM when o is close to zero, while its importance diminishes at higher values of σ. As a increases, both methods still achieve reasonable score estimates even without VR. Additionally, when using complete data and varying the missing ratio a {0.1, 0.3, 0.5) with DSM and VR at a 0.001, we observe that the first- and second-order score estimates remain close to the ground truth. While performance degrades as the proportion of missing data increases, the estimates remain generally accurate.
4.1. Improving stability with Variance Reduction
It important note that, in order match the score of the true distribution Paan(()) needs to be close to zero for both DSM and DySM, so that qa (x) closely approximates Pama(x), However, training score models using denoising methods can suffer from high variance when a approaches zero. This challenge motivates the use of variance reduction techniques. Building on existing variance reduction meth- ods for DSM (Song& & Kingma, 2021; Meng et al., 2021), we propose tailored variance reduction techniques specifically for training DSM with missing data, as follows.
4.2. Scalability and Numerical Stability
We show that the proposed method efficiently and accurately estimates second-order scores across different missing ra- tios, as summarized in Table 1. To achieve this, we generate 10 synthetic datasets with known ground truth, consisting of 100-dimensional correlated multivariate normal distri- butions that include varying levels of missing data under MCAR mechanism. The covariance matrix for these dis- tributions is constructed using eigenvalues t E {1, to vary degrees of correlation. We evaluate the performance of the estimated 8 and the diagonal of 82 by calculating the mean squared error (MSE) between the estimated scores and the ground truth scores derived from the complete data, across various values of a. Our results indicate that the jointly optimized 81 (Xals; 0) and (x,s; 0) achieve empir- ical performance close to the ground truth. As previously mentioned, when a approaches zero, the estimates with- out VR become unreliable for both the score and Hessian, likely due convergence issues. Although performance de- clines with an increase in missing data, the estimates remain reasonable.
(2 Tz
LDSM-VR (0) = (0) - BEx Eu..Nio..
Iz O gi (Xoba, m)
① gi (Xiins m) +
0²
where gi (Xobs: m) = 1 m under MCAR, and gi (Xobs, m) P/m=0
1-m under MAR.
For the second-order model with missing data, we imple- ment a variance reduction (VR) technique using antithetic sampling (James, 1985; Meng et al., 2021), which involves utilizing two negatively correlated sample vectors centered around x. The objective function is then formulated as
CD,SM-VR(()) Ex... mL
{
(3
)2 + 2b (x )2
I zz + o
3 2 g2 (Xbss m)
5. Sampling with Missing Data via Second-Order Score Models
(8)
where the antithetic samples are defined as x xts = Xohs +
In this section, we illustrate how our second-orde score model 82 s>(;; 0), trained on incomplete data, improves sam-
oz and x Sobg = Xobs OZ. Here, 2b = $2 + S1sT. and
T (xxx)) (xx 2xb((xtt))) (x Under the MCAR
5

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
Ground Truth
Ground Truth
Ground buth
Ground Iuuh
DSM (.11
DSM Cor 0.1)
DSM-VR ta = 0.01
DSMVV 0.0)
DSM (o=0.0003
DSM (or 0.001] (o
DaMAR m0.)
DSM-VR (a 0.1)
VR (a- 0:11
DSM 0
DSM VR
DS 0.001)
DSM-VA Co 001)
0.51
DSM-VR 0.5)
(a) Estimated $
(b) Estimated $2
(c) Estimated $1
(d) Estimated $2
Figure I. Comparison of estimated 8 and s3 under different conditions. (a) and (b) show estimates with DSM and D2SM varying the noise level with fixed missing ratio 0.3. (c) and (d) show estimates with DSM and DySM varying the missing ratio a.
Table 1. Mean squared error (MSE) between the estimated first-order and second-orde order scores and the ground truth is evaluated across 5,000 test samples. We vary the noise scales a and missing ratios a, with each configuration tested using 10 random seeds.
a 0.0
a = 0.1
= 0.3
a 0.5
Methods
0.1
== 0.01
0= 0.1
0= 0.01
0 0.1
0 0.01
0' 0.1
== 0.01
81
0.28 + 0.01
0.42
0.29 0.01
0.42 0.02
0.32 0.01
0.01
0.37 0.01
0.44 0.02
s, (VR)
0.07 +
0.07 0.00
0.09 0.00
0.13 2.34
0.09
0.15 +1000.01
0.23 0.00
0.27 = 0.01
82
0.16 -: 0.02
15.42 0.47
0.16 : 0.02
29.74
0.17 + 0.03
26.08 + 2.03
(2R)
0.04 1: 0.00
0.05 4 0.00
0.04 0.00
0.04 0.00
0.04 0.00
0.05 0.00
0.05 + 0.00
0.06 0.01
ple quality and enables faster convergence in terms of the number of sampling steps required. We demonstrate the ef- fectiveness of the proposed model through simulations and real-world datasets containing missing values, comparing its performance against various baseline methods.
generates comparable data to Langevin dynamics with fewer iterations, resulting in data that is more concentrated around the original distribution, while Langevin dynamics yields noisier and more dispersed results.
Following Kim al. (2022); Ouyang et al. (2023), we con- duct experiments on simulated Bayesian Network dataset and real Census dataset (Kohavi, 1996) to illustrate the efficiency and effectiveness of the data generated by our proposed model trained on missing data.
Langevin dynamics. Langevin dynamies samples data from (x) by utilizing the first-order score function 81 (x) (Bussi & Parrinello, 2007; Song & Ermon, 2019). Starting with a prior distribution xx). a fixed step size e > 0, and an initial value &o ~ (xx,, Langevin dynamics iteratively updates the samples as follows:
Baselines. We evaluate the proposed method using both Langevin and Ozaki sampling against several baseline tech- niques for synthetic data generation on datasets with miss- ing values. Specifically, we implement (1) a vanilla DSM model that removes rows with missing values (termed DSM- delete), and (2) STaSy (Kim et al., 2022), state-of-the-art score-based basre model, which significantly outperforms other approaches for tabular data (termed STaSy-mean). Since STaSy requires complete datasets for training, we apply mean imputation to handle any missing values in the train- ing data.
1 x, 2 cs (x,) + Vez.
where 2, ~ N(,, I) represents Gaussian noise.
Ozaki Sampling. Following Meng et al. (2021), Ozaki discretization improves data synthesis by integrating second- order information from (x() to precondition the sampling process. The updates for Ozaki sampling are performed as follows:
X = -x,-1 + (1-151(((-1) + 11/2 Ze,
(10)
Metrics. Following Kim et al. (2022); Ouyang et al. (2023), we employ two criteria. fuiteria,uueelity and machine learning (ML) efficiency. to assess the quality of the generated synthetic tabular data. For evaluating fidelity, we utilize the model- agnostic library SDMetries. The result ranges from 0 to 100%. A higher score indicates better overall quality of the synthetic data. To measure ML efficiency, we adopt the same pipeline as Kim et al. (2022); Kotelnikov al. (2023), training various models-including Decision Tree, AdaBoost, Logistic Regression, MLP Classifier, Random
where 2, ~ NO, I), M,- (022(R1 1) I, and Σ,
(62cx ())
Illustration. We use the Roll dataset to demonstrate the effectiveness of Ozaki Sampling, focusing on its conver- gence rate and quality of data generation through second- order information. Both methods employ a step size of e = 0.005 and a missing ratio of 0.5 under the MCAR missing mechanism, As shown Figure 2, Ozaki Sampling
6

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
Forest, and XGBoost on the synthetic data and testing them with real data. Our primary metric is classification accuracy, and we also report AUROC and Weighted-F1 scores in the Appendix C.4. All experimental results are based on three repetitions.
using CAM pruning to derive the final directed acyclic graph (DAG) (Bühlmann et al., 2014).
Baselines. We utilize the MissForest imputation method to address missing data, followed by the implementation of DiffAN (Sanchez. et al., 2022) (termed MissDiffAN) and DAGMA (Bello et al., 2022) (termed MissForest) for struc- ture learning. DiffAN serves as a diffusion diffusiin-bbsed adaptation of the approach proposed by (Rolland et al., 2022), ensur- ing a fair comparison. Furthermore, core with MíssDAG (Gao et al., 2022) and MissOTM (Vo et al., 2024), both of which are prominent methods that have shown superior performance in causal discovery with missing data relative to various other baselines.
Results. Figure 5 and Table 2 demonstrate the effectiveness of the proposed method on both the simulated Bayesian Network dataset and the Census data, showing superior performance in terms of fidelity and utility compared to other baselines. Specifically, these results confirm that the impute-then-ge seatte-thhe--genrrat approach introduces bias, whereas directly learning from missing data significantly improves the performance of the generative model. Furthermore, the advantages of the proposed model become more pronounced as the missing ratios increase. Additional details and results of the experiments can be found in Appendix C.
compare MissS-
Metries, All quantitative results are averaged over 10 ran- dom initiazations. For comparing the estimated DAG with the ground-truth one, we report commonly used metries: Order Divergence and Structural Hamming Distance (SHD). Order Divergence measures the number of errors in the or- dering, while SHD indicates the minimum number of edge additions, deletions, and reversals needed to convert the recovered DAG into the true one. Lower values for both metries are preferred. Order Divergence is calculated only for MissScore and MissDiffAN, as these are the only order- based methods.
Table Fidelity and ML efficiency evaluation of MissScore using
Langevin and Ozaki samplings, along with other baselines, on the Census dataset with a missing ratio of 0.3 under MCAR.
Langevin
Ozaki
DSM-delete
STaSy-mean
Fidelity 86%
88%
73%
82%
MI. efficiency
80%
81%
70%
77%
Simulations. We simulate synthetic datasets generating a ground-truth DAG from the graph model Erdōs-Rényi (ER). Each function f/ is constructed from a multi-layer perceptron (MLP) and a multiple index model (MIM) with random coefficients. We consider general scenario of non- euqal variances, sampling 1000 observations according to all missing mechansims: MCAR, MAR, MNAR at 10% and 30% missing rates and complete data. In the main text, we report the SHD and runtime in MCAR cases with missing ratio 0.1 varies with number of dimensions, using Gaussian noise. Specifically, the number of edges is set equal to the dimensionality.
6. Causal Discovery with Missing Data via Second-Order Scores
Background. Causal discovery aims to identify causal rela- tionship from purely observational data. However, the task is ill-posed without additional assumptions. Assuming an additive noise model (ANM) allows for the identification of causal structures. In this context, consider the ANM defined as x; = f((PA))) + 24. where ft is nonlinear function and 24 is a Gaussian noise. Rolland et al. (2022) proposed an order-based algorithm that uses the second order score of an ANM with a probability distribution Paat(x) to identify leaf nodes, and iteratively determine the topological order of the variables. However, the computation of the Hessian requires complete data, which poses challenges in real-world scenar- ios such as clinical trials, and biology, where missing data is common.
Results. In Figure 4, our approach demonstrates comparable performance to the state-of-the-art methods MissDAG and MissOTM, although it shows slightly lower performance in the high-dimensional scenario with d = 50. This may be attributed to the challenges of training a denoising score matching model with a limited number of samples. How- ever, MissScore significantly reduces computation time com- pared to MissOTM and MissDAG. offering an efficient alter- native while maintaining similar results. With MissForest, our performance is similar, possibly due to the additional constraints enforced during training by DAGMA. However, when comparing our method to those that rely on imputa- tion followed by causal discovery approach, we find that our approach outperforms MissDiffAN, This suggests that im- putation can introduce bias in downstream tasks. Additional results for various missing ratios, mechanisms, and order
A straightforward approach to address missing data problem is to first impute the incomplete entries using ohf-tf imputation methods and then apply existing causal discov- ery methods. However, this two-step approach can be sub- optimal, the imputation process may introduce bias for modeling the underlying data distribution. Our method mit- igates this issue by directly training second-order model with incomplete data, thereby reducing potential bias. Since the Hessian only provides information about variable order, we adopt a strategy similar to Rolland et al. (2022); Sanchez et al. (2022), first computing the topological order and then
7

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
(a) Data
(b) Langevin (100)
(c) Langevin (200)
(d) Langevin (300)
(c) Ozaki (50)
(f) Ozaki (100)
(g) Ozaki (200)
(h) Ozaki (300)
Figure Sampling Swiss Roll dataset with a step size of 0.005 using Langevin dynamics and Ozaki sampling. Ozaki sampling achieves higher-quality samples with fewer iterations compared to Langevin dynamics under MCAR with missing ratio of 0.5. Figure (a) displays the dataset; Figures (b)-(d) show the results from Langevin dynamics; and Figures (e)-(h) present the results from Ozaki sampling. The numbers in parentheses indicate the number of iterations taken.
11:0
1182
50
2180
Langevin
DSM-mean
Ozaki
4.
118
St
mean
Dimets
Missing Ratio
Missing Rabio
Pacicas
Mashorat
Mudcors
Natova
MCAR
MAR
(a) SHD
(b) Runtime(scc)
Figure 3. Fidelity evaluation of MissScore using Langevin and Ozaki samplings. along with other baselines, the Bayesian Network dataset with varying missing ratios a (0.1, 0.5, 0.7, 0.9) under different missing mechanisms.
Figure 4. Data is generated under MCAR with missing ratio of 0.1, varying dimensions The sample size is 1000, and fa corresponds to an MLP. Left: SHD; Right: Runtime. The shaded area indicates confidence.
{10, 20, 50) using ER graph model.
divergences are provided, along with further experimental details, in Appendix D.4.
competitive results against state-of-the-art approaches.
Despite these advancements, certain limitations remain. The effectiveness of high-order score estimation diminishes in low-noise environments without variance reduction, partic- ularly in scenarios with high levels of missingness. These issues also extend to downstream tasks, such as causal dis- covery and sampling, where performance tends to decline slightly with increasing dimensionality. Additionally, di- rect testing on real-world data poses a challenge due to the unavailability of ground truth.
7. Conclusion and Future Work
In this work, we present novel framework for directly estimating high-order data density scores in the presence of missing data using DSM under M(C)AR mechanisms. Unlike methods that rely on imputation or additional model training, our approach handles missing data directly and ef- ficiently. Distinct from prior work like MissDiff, MissScore focuses high-order score learning and introduces an ora- cle estimator under the MAR mechanism. Empirical results demonstrate the efficiency and accuracy of our framework in estimating second-orde orderr scores, with applications showcas- ing improved sampling quality through Ozaki discretization in sernarios with missing data. Furthermore, our causal discovery method scales with dimensionality and achieves
Future work could explore extending high-order score es- timation to domains such image and time-series data, while investigating the application of DSM directly han- dle MNAR data. Additionally, incorporating constraints into score-based models to enhance causal discovery with missing data holds significant potential for advancing score- based learning complex and missing data scenarios.
8

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
Acknowledgements
accelerated mri Medical image analysis, 80:102479, 2022.
This research was supported by the University Mel- bourne's Research Computing Services and the Petascale Campus Initiative. WL was supported by Melbourne Re- search Scholarship. EG was supported by the Centre for Augmented Reasoning (CAR). QK was supported by ARC grant DE250100030. MG was supported by ARC grants DE210101624 and DP240102088, as well as WIS-MBZUAI grant 142571.
Dalalyan, A. S. and Karagulyan, A. User-friendly guaran- tees for the langevin monte carlo with inaccurate gradi- ent. Stochastic Prucesses and their Applications, 129()2); 5278-5311, 2019.
Synthetic Data Metrics. DataCebo, Inc., 2023.
Dempster, A. P., Laird, N. M., and Rubin, D. B. Maximum likelihood from incomplete data via the em algorithm. Journal of the royal statistical society: series B (method- ological), 39(1):1-22, 1977.
Impact Statement
This paper aims to advance the field of machine learning by presenting novel framework for high-order score estima- tion in the presence of missing data. Our work has the po- tential to enhance applications in data synthesis, generative modeling and causal discovery, which are crucial in domains such as healtheare, finance, and social sciences. By directly addressing the challenges posed by missing data, our frame- work offers an alternative traditional inputation-basek methods, enabling more accurate and unbiased modeling of incomplete datasets. Our contribution aligns with ongoing efforts to improve the ability of machine learning systems to operate under realistic and challenging data conditions, promoting more reliable, efficient, and equitable solutions across diverse applications. Additionally, the focus on high- order score estimation introduces new opportunities for im- proving sampling quality, uncertainty quantification, and the recovery of underlying causal structures, making this work applicable to both foundational research and practical implementations.
Efron, B. Tweedie's 's formula and selection bias. Journal of the American Statistical Association, 166)4 1614, 2011.
Esteban, C., Hyland, S. L., and Rätsch, G. Real-valued (medical) time series generation with recurrent condi- tional gans. arXiv preprint 1706.0:7000000000 2017.
Gain, A. and Shpitser, I. Structure learning under miss- ing data. In International conference on probabilistic graphical models, pp. 121-132. Pmlr, 2018.
Gao, E., Ng, I. Gong, M., Shen, L., Huang, W., Liu, T., Zhang. K., and Bondell. H. Missdag: Causal discovery in the presence of missing data with continuous additive noise models. Advances in Neural Information Process- ing Systems, 35:5024-50033 2022.
Gondara, L. and Wang, K. Mida: Multiple imputation using denoising autoencoders. In Advances in Knowledge Dis- covery and Data Mining: 22nd Pacific-Asia Conference, PAKDD 2018, Melbourne, VIC, Australia, June 3-6, 2018, Proceedings, Part III 22, pp. 260-272. Springer, 2018.
References
Bello, K., Aragam, B., and Ravikumar, P. Dagma: Learn- ing dags via m-matrices and a log-determinan acyclicity characterization. Advances in Neural Information Pro- cessing Systems, 35:8226-888 8239, 2022.
Goodfellow, I., Pouget-Abadie. J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y. Generative adversarial networks. Communications of the ACM, 63(1 13:3111114444 2020.
Bühlmann, P., Peters, J., and Ernest, J. Cam: Causal additive models, high-dimensional order search and penalized regression. The Annals of Statistics, 42(6)) 2014.
Ho, J., Jain, A., and Abbeel, Denoising diffusion proba- bilistic models. Advances in neural information process- ing systems,
2020.
Bussi, G. and Parrinello, M. Accurate sampling using langevin dynamics. Physical Review E E-Statistical, Non- linear, and Sofi Matter Physics, 7555)0056707, 2007.
Hyvärinen, A. and Dayan, P. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
Chen, Z., Li, H., Wang, F., Zhang, O., Xu, H., Jiang, X., Song, Z.,, and Wang, H. Rethinking the diffusion models for missing data imputation: A gradient flow perspective. Advances in Neural Information Processing Systems, 37: 112500112103, 2024.
Jäger, S., Allhorn, A., and Bießmann, F. A benchmark for data imputation methods. Frontiers in big Data, 4: 693674, 2021.
James, B. Variance reduction techniques. Journal of the Operational Research Society, 36(6):525-530, 1985.
Chung, H. and Ye, C. Score-based diffusion models for

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
Volicoeur-Martinee A., Fatras, K., and Kachman, T. Gen- erating and imputing tabular data via diffusion and flow- based -bossted trees. In International Conference on Artificial Intelligence and Statistics, pp. 1288-1296, PMLR, 2024.
Meng, C., Song, Y., Li, W., and Ermon, S. Estimating high order gradients of the data distribution by denoising. Advances in Neural Information Processing Systems, 34: 25369, 2021.
Morales-Alvarez, P., Gong, W., Lamb, A., Woodhead, S., Peyton Jones, S., Pawlowski, N., Allamanis, M., and Zhang, C. Simultaneous missing value imputation and structure learning with groups. Advances in Neural Infor- mation Processing Systems, 35:2001 1-20024, 2022.
Kim, J., Lee, C., and Park, N. Stasy: score-based tabular data synthesis. arXiv preprint arXiv:2210.04011 2022.
Kingma, D. P. Auto-encoding variational bayes. arXiv 2013.
preprint arXiv:
Muzellec, B., Josse, J., Boyer, C., and Cuturi, M. Miss- ing data imputation using optimal transport. In Interna- tional Conference on Machine Learning, pp. 7130-7140. PMLR, 2020.
Kohavi, R. Census Income. UCI Machine Learning Reposi- tory, 1996. DOI: htpps://o..ogg//024333//CCC
Kotelnikov, A., Baranchuk, D., Rubachev, I., and Babenko, A. Tabddpm: Modelling tabular data with diffusion mod- els. In International Conference on Machine Learning, pp. 17564-17579. PMLR, 2023.
Nazabal, A., Olmos, P. M., Ghahramani, Z., and Valera, I. Handling incomplete heterogeneous data using vaes. Pattern Recognition, 107: 107:107501, 2020.
Li, H., Yang. Y., Chang, M., Chen, S., Feng, H., Xu, Z., Li, Q. and Chen, Y. Srdiff: Single image super-resolution with diffusion probabilistic models. Neurocomputing, 479:47-59, 2022.
Ouyang, Y., Xie, L.,, Li, C., and Cheng, G. Missdiff: Train- ing diffusion models on tabular data with missing values. arXiv preprint arXiv:2337.0044 00460000044,,. 2023.
Paterakis, G., Fafalios, S., Charonyktakis, P., Christophides, V., and Tsamardinos, I. Do we really need imputation in automl predictive modeling? ACM Transactions on Knowledge Discovery from Data, 18(6):1 64, 2024.
Li, S. C.-X., Jiang, B., and Marlin, B. Misgan: Learning from incomplete data with generative adversarial net- works. arXiv preprint arXiv: (o0p:9902.9999 2019.
Peis, I., Ma, C., and Hernández-Lobato, M. Missing data imputation and acquisition with deep hierarchical models and hamiltonian monte carlo. Advances in Neural Information Processing Systems, 35:35839-35851, 2022.
Little, R. J. and Rubin, D. B. Statistical analysis with missing data, volume 793. John Wiley Sons, 2019.
Liu, Q., Lee, J., and Jordan, M. A kernelized stein discrep- ancy for goodness-of-fit tests. In International conference on machine learning. pp. 276-284. PMLR, 2016.
Poulos, J. and Valle, R. Missing data imputation for su- pervised learning. Applied Artificial Intelligence, 32(2): 186-196, 2018
Liu, W., Huang, B., Gao, E., Kc, Q., Bondell, H., and Gong, M. Causal discovery with mixed linear and nonlinear additive noise models: scalable approach. In Causal Learning and Reasoning, pp. 1237-1263 PMLR, 2024.
Rolland, P., Cevher, V., Kleindessner, M., Russell, C., Janz- ing. D., Schölkopf, B., and Locatello, F. Score matching enables causal discovery of nonlinear additive noise mod- els. In International Conference on Machine Learning, pp. 18741-18753. PMLR, 2022.
Lu, C., Zheng, K., Bao, F., Chen, J., Li, C.,, and Zhu, Maximum likelihood training for score-based diffusion odes by high order denoising score matching. In Inter- national Conference on Machine Learning, pp. 14429- 14460. PMLR, 2022.
Rubin, D. B. Inference and missing data. Biometrika, 63(3): 581 581-592, 1976.
Sabanis, S. and Zhang, Y. Higher order langevin monte carlo algorithm. 2019.
Ma, C., Tschiatschek, S., Turner, R., Hernández-Lobato, J. M., and Zhang, C. Vaem: deep generative model for heterogeneous mixed type data. Advances in Neural Information Processing Systems, 33:1 13:17237-1247, 2020.
Sanchez, P., Liu, X., O'Neil, A. Q., and Tsaftaris, S. A. Diffusion models for causal discovery via topological ordering. arXiv preprint arXiv:2210.06201 2022.
Mattei, P.-. and Frellsen, J. Miwae: Deep generative modelling and imputation of incomplete data sets. In International conference on machine learning. pp. 4413- 4423. PMLR, 2019.
Seaman, S. R. and White, R. Review of inverse probabil- ity weighting for dealing with missing data. Statistical methods in medical research, 22(3):278-255 2013.
10

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
Shadbahr, T., Roberts, M., Stanczuk, J., Gilbey, J., Teare, P., Dittmer, S., Thorpe, M., Tomé, R. V., Sala, E.,, Lió, P., et al. The impact of imputation quality machine learning classifiers for datasets with missing values. Com- munications Medicine, 3(1):139, 2023.
Wang, Z., Akande, O., Poulos, J., and Li, Are deep learning models superior for missing data imputation in large surveys? evidence from an empirical comparison. arXiv preprint arXiv:2103.09316, 2021.
Wooldridge, M. Inverse probability weighted estimation for general missing data problems. Journal of economet- rics, 141(2): (28))))221113000 2007.
Shpitser, I. Consistent estimation of functions of data miss- ing non-monotonically and not at random. Advances in Neural Information Processing Systems, 29, 2016.
Yang, S., Wang. I.., and Ding, P. Causal inference with confounders missing not random. Biometrika, 106(4): 875-888, 2019.
Song. Y. and Ermon, S. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019.
Zhang, H., Fang. L.,, and Yu, S. Unleashing the potential of diffusion models for incomplete data imputation. arXiv preprint arXiv:2405.20604 2024.
Song. Y. and Ermon, Improved techniques for train- ing score- bosor- generative models. Advances in neural information processing systems, 33:12438-12444 2020.
Zheng, S. and Charoenphakdee, N. Diffusion models for
Song, Y. and Kingma, D. P. How train your energy-based models. arXiv preprint arXiv:2101.03288 2021.
missing value imputation in tabular data. arXiv preprint arXiv:2210.17133 2022.
Song. Y., Garg, S., Shi, J., and Ermon, S. Sliced score matching: A scalable approach to density and score es- timation. In Uncertainty in Artificial Intelligence, pp. 574-584. PMLR, 2020.
Zheng, X., Aragam, B., Ravikumar, P. K., and Xing, E. P. Dags with no tears: Continuous optimization for structure learning. Advances in neural information processing systems, 31, 2018.
Song, Y., Shen, L.., Xing, I.., and Ermon, S. Solving inverse problems in medical imaging with score-based generative models. arXiv preprint
2021.
Städler, N. and Bühlmann, Missing values: sparse in- verse covariance estimation and an extension sparse regression. Statistics and Computing. 22:219-235, 2012.
Stekhoven, D. and Bühlmann, P. Missforest-non non- parametric missing value imputation for mixed-type data. Bioinformatics, 121111111118,.111111111111111111 2012.
Tu, R., Zhang, C., Ackermann, P., Mohan, K., Kjellström, H., and Zhang, K. Causal discovery in the presence of missing data. The 22nd International Conference on Artificial Intelligence and Statistics, pp. 1762-1770. Pmlr, 2019.
Van Buuren, and Orusthuis-Oudshoonn,, mice: Multi- variate imputation by chained equations in Journal of statistical software, 45:1-67, 2011.
Vincent, A connection between score matching and de- noising autoencoders. Neural computation, 23(7):1661- 1674, 2011.
Vo, V., Zhao, H., Le, T.. Bonilla, E. V., and Phung, D. Optimal transport for structure learning under missing data. arXiv preprint arXiv:2402. 15255, 2024.
Wang, X., Zhang, R., Sun, Y., and Qi, J. Doubly robust joint learning for recommendation on data missing not at ran- dom. In International Conference on Machine Learning, pp. 6638 6638-6647. PMLR, 2019.

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
Appendix
Table of Contents
A Proofs
13
A.1 Proof of Theorem 3.3
13
A.2 Proof of Theorem
A.3 Proof of Theorem 3.5
16
A.4 Proof of Eq.(7)
A.5 Proof of Eq.(8)
17
B Data Generation under Different Missing Mechanisms
19
Additional Information on Sampling
20
C.I Dataset Description and Processing
20
C,2 Evaluation Methods
21
C33 Model Architecture
22
C.4 Experimental Results
22
C.S Broader Empirical Validation
22
Additional Information on Causal Discovery
D.1 Related Work
24
D.2 Evaluation Metrics...
25
D.3 Model Architecture
25
Experimental Results
25
12

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
A. Proofs
A.1. Proof of Theorem 3.3
Proof. Under the MCAR missing mechanism, we know that
Ex m Ex/x.
s(x: 0) +
xx - x)
1 (1 m
02
s(x: 0) + - -x) 0 (1 m
(((x;)) [ + J (x xx)
= xx Fm|
11 m
xxxx
s(x; 0)
+ 0 (x x) 1 Emix m
Exx
8 s(x;
0z (X
x)
2
Wi
0²
Denote lo (8|x) NNxx; )²), only need to show there exists some constant C independent of A such that
Exxx [[(s((::))
Vxqo (x[x)) @ wille = Exx [[[s((x: 0) s(x) ww +
(11)
For simplicity, we consider the case d For the right hand side of (q.(11),
R.H.S W1 (Bz,z (i):0)] F.. (s(;00)s())}} C.
For the left hand side of Eq.(11),
L.H.S = w1 EEzz [i: - 2B. (s(ē; Vz log ( + C.
Hence, we only need show
Ea (s(z;0)s(T)) Ez,x [s(2; V2 log do (f[z)]
(12)
We have,
Ex (s(z00)s({)) - / aapa
((;; (ple))a (F(x))d dē
/ / (tt;0)p(<) (g. (T(=z) dedā
/ ((;;0)p(z) / (iz:0)pa. V V (log go (z) dxdž = Eg,w [s(;;))V log go (ā|z))
drds
Hence, we get our desired result under the MCAR missing mechanism. By similar argument the proof for the MCAR missing mechanism, under the MAR missing mechanism, we have:
13

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
s(x: jz (x-x)
1 m
Ex fx/xm
P m = 0|x
= =Ex,x. x.m
0) 0) + 1
s(x;
02 1x -x] 3
m
m
0|x
Ex x Bm|&x
s((x 0) + -े (8 - x)
3
m
v PP m 0|x
s(x: + 0² 02(x²- (x-x)) x1
mm m)
Ex
P m 0|x x
{ s(x: 0) + //2 (X x)
1
Exxxx
where the last equality comes from Em x m) Em (1 m) P[mm 0 =x]. Then following the steps in the proof for the MCAR missing mechanism, taking w as vector of 1, we get our desired result.
A.2. Proof of Theorem 3.4
Proof. Under the MCAR missing mechanism, we know that
27
Ex.x.m
s2(x;;)) + (Xt)s (x)+ +
I zz
2 {11 -m)mm - m) []
σ²
Exxx
(x 0) + (x)s) (x)
zz
[m]xx (1 m)(1 m)
x2
27
F
82 (x 0) + (x)s
W2
r²
It is sufficient to show
{ (xx 0) + sx(ss (xx) + I- zz T
27
Exx
W22
(13)
= ExEx [1152((;0) - -s(())} o /W2122 [ +C.
with some constant C' independent of 0. For the right hand side of (13))
E Ex
[1(22 (x: 0) s2())) 2 vwall
d d
W2,65 (Ea.. [s (ā;; 0)] 2E,, [92.55 (ā: 0) 0)8.1; + Er.. [> (f)]
i=1 j=1
Er., [02.0 (z; 0)] - 2E, [52,1 (ā; (z; + Ea.. [3 + C
W2,65
where C1 some constant independent of 0.
14

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
For the left hand side of Eq.(13), note that I-ax
V
(x
Vzlo (x x)
VVx (x/x)) T then we know that
2
I zz
E x
S2 (x;
+ (X))) (x)
W2
x²
W2.35
Ez [02,,5 (z; 0)]
+ 2E
2g, (T; (T)s.())))
2E
|22,5 (z;
0)
V
lo
2E
$2.15 (ā: (x 90
Cz.
Comparing left and right hand side of Eq.(13), it sufficient to show
B (#; 82. [i) = E [22,, (ā: 0)V a; Va, log (a/2z)]
(14)
E [s2 (o
E (8))
We have
Ez,. [#2,1)(;; ((:;0))..;(()))
$2,15 (z: 02 (3) (T)ps (z)dī
V V2,p(ã)
}) ps
$2,1( (2; 0)
(x(()) (i)d
p(z)
/ (2,) 0)
V,, V(z)(T)
Va,p(T) Vayp(T)
pa (3)
p(ū)
p²(t)
/ (e(1:))V: Va,p(ūjdk
(15)
f (2(0:0)
V..((T) Va,p(a) dā.
(16)
For Eq.(1
/ 62(5; Va,p=(E)dā
(f } dē
82 02
//
(5: 0) (o
(xr)dddr
//
82 (2 01] v3. Va log 1o (1))))
qo (6 (x) (ddd
+ //
V2, 9o
2, 90 (z)x)
(2 0)
Pa (2r dāde
yo xr)
//
(a:
0)
V V.,
log aa (7)zz
Por (z, d)dädr
(olz)
(
+
(5 0)
Pa (x, FddTdo
(Z
Lo (2)xr)
=
E
(s2(;:)) 0) V 2. Va, log ( [(
+ E
(5)z) (2)
For E16(1
(2((;;)) 0) Ve,p(T) (())(
dā /
V Va,P(T) V2,p(ã)
82 s2(;;; 0))
p(o) p())dz
p(ā)
p())
E
[s2(;; ())
15

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
Combine all the results, get the desired result under the MCAR mechanism. By similar argument in the proof under the MCAR mechanism, we have
zz T a²
{(1 m) m) }
Exx.m
+ss))
+
P Pmmm = =0 x
E Em/xx( (1 - m) ml
- zz
= =xx,x
82 0)
(8)s)
(x)
0²
mm 0| x
{
I zz
-
82 & 0) +
(X)s}
x
r²
where the last equality comes from Em x.x [1
m1 - m) |
Empx((1 -m)(1 - m) ] F(mm = 0|x = x).
([11-)]]]
Then following the steps in section A.2, taking W2 in A.2 as 1, we get our desired result for the MAR mechanism.
A.3. Proof of Theorem 3.5
We first consider the MCAR case. Recall that w = 1 m
Exx.m [1
8 ww]]]
= Exx,m [I] f (
80 m|]
Ex,xE Em [I] [ x fx(( (3)
8* m)mm |2]
Exxx
[]] f (x, $1 ..., (x)
Em @* (1 (1-m)} |?
Noting that Em x (1 m) is constant. We can show that the solution for the weighted least square equation for any constant matrix c is
arg E - [[***xx())} o o<]]]]
f0 x [x]
Such equation holds since
E 1||(0'x - h(x)) oc|2] =E ||| [8"x x 3 c
+ E [x] Mx)) c
2 E
By Theorem 3.2, we know that E For the MAR case, recall that w "Im*=0Xx=x
x x
Sk (x)) (x)
sı
Hence, the desired result follows.
1-m
E Exx.m
[1 x
fa (xx ss
0 wl]]
Exxm
f8"x 18*x-1A0 - (1....s(().....s_113)...
m
P I 0
|11100x - (8) (x) Su X 0)) + 8*
(1- m)
ExxEmxx
P = 0|x =
Exx [Il (8*x- fx( si (8) (x: 0))} [13]
Thus, by Theorem 3.2, the solution to the criterion function is
arg min Exxxm #
[l] { xx fx(
Sk- (x; ((x))
80@"w/ [wwl3 -
16

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
A.4. Proof of Eq.(7)
Denote the oracle criterion function
|| [s((())0) + 0 (x - x) o Og m)
LDSM 0) == 1x m B RxIx m
where (x, m) =1 under MAR and gı (x)
0fm==x== x
we want to match the score of true data distribution pxx), should approximately zero for both DSM and D2SM so that 2o(X) (x) is close to p(x). According to Taylor expansion we have,
(x: 0) 0² (x ( o
(O) Ex, x.m
81 + oz: 0) [ o m)
= Ex,mB
81 +o 0/2
o (x. m) + +0(1)
Ex B.mmE.uu
(1(x;;
@ (2.
+ (111
-x.mB N(o.))
()
2 gi (x.
m)ll
+ (25((x))))) 2 (x m)
+ 0(1)
where z X-x where (1) is bounded as a approaches zero. However, when evaluating the expectation above from samples, the variances of [2Sg(x.m)) mill and a(x:0) both increase without bound as a nears zero, due to the terms involving a²
and o² in the denominator. This leads to a significant increase in the variance of the DSM loss, complicating optimization process. As consequence, DSM may become unstable and fail to converge when is small, highlighting the need for methods to reduce variance.
We have,
Oz
m) m)l 2
2
||g (x, m)ll
B (
(x;
2²
where d is the dimension of the data distribution p( p(x). Therefore, we can construct variable that is, for sufficiently small positively correlated with LOSM while having an expected value of zero:
co(x;2)= (2)1 + Oggi (x, m)
|z O gi m)]l
||g: m) |
d2
Subtracting it from LOSM will yield an estimator with reduced variance for DSM training with missing data:
LOSM vR (0) = (DSM(() Fx B E2-.())) ((=((x;0))
2z 3 gi (x, m) 1
+ @ggi m) +
σ²
Here omit the part || [x.)(x since it is independent of 9.
A.5. Proof of (8)
Similar to proof A.4, consider the oracle criterion function
| ss + (1)ss
I- zz
(D,SM(()) Eg.x.m
(x)
g2 (x, m)
o²
17

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
where g2 (Xabs, m) = (1 - m1 m) under MCAR, and g2 (Xobs, m)
(1 m)
under MAR.
B mm
x= X=xii
Denote (xx = 82 (Kahsi 0) + Si (Naba))) (Xrhs), we have
I zz
EDSSM( Exxxm
52 0) (x)()) (x)
2 gx m
o²
[[
{ 1/(xx -z: 0) +
I-zz zz
= =Rx,m 1
2 g2(2x, (x, m)
r²
{{ 1b((xx + 02: +
I nz
T
2
2z o²
= =x Ba
+ 2xb(x + 02: 0)
3 0z2(x, m)
2
Denote this (z;; 0) as the ijth term of (b(x;)) Oij = Lij 2:2g and Dij as the ijth term of ge g2 m) and according to Taylor expansion, we have,
E2 E.-Noo)
Mis( + 02:00² +
Φi + 22b2i 02;0) Pij
(t)
02
= E. (~((0)))
(z) (x;
Sij 2V
2 (x;
Pui
+c] @gyle,
+ 2000gyle
0(1)
o2
where: 2 Z-z with C²= = dig and (x + o2; representing the derivative of thiy(a + a2;0) with respect to æ.
can be treated constant that does not depend on 0, and ((1) remains bounded when - 0. However, when calculating the expectation from samples, the variances of * and Vehig (z; 0) * increase without bound as a + 0, due to the presence of a and o² in the denominator. This causes significant rise in variance for D2SM, making the optimization process more difficult. Consequently, D2SM can become unstable and fail to converge as a approaches zero, necessitating the use of variance reduction techniques.
012
In this case, for estimating Viht (a:; 0) n2. we utilize the antithetic sampling technique same as Meng et al. (2021) reduce variance.
employ the same variance reduction method outlined in the proof of A.4. However, to bypass the need
Denote 7. + and ã_ x + as the antithetic samples, according to Taylor expansion, the ijth term of the D2SM(() then becomes,
K (Mu(Er::))+
2 + (Wo(E_ 220 )"} @ (0( 0z)
(D.SM(()s E. E.B.NN(O..)
th (Z4;)0))) +
Wi (ā 0) + U (=_ +2 uh Ue (ā (a= 0) c
2 ())
0²
(Z):)0))) 4. (8_ 0) +2 24 0²
+ c () 0))
(1)( Ao
2 @
(1;
2y o²
4 0) 2V 2v + +c (5)
ps ()) + c o (3s
2 (vrr
where C' - (y
constant with respect to the optimization.
Therefore, we have the variance reduction for CD.SM-VR(( (0) which is equivalent to optimizing Eq. (3.4) up to control variate. Moreover, when approaches zero, optimizing Eq. (8) more stable.
LD,SM (0) = Ex F N(..)
{f
(1(xt)2 + 2x +
I zz
3 v()) @ @g((x,
(17)
18

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
Algorithm 1 MissScore
Input: Observed data Xobs, score models $1 (=; (2(;)), noise level o, coefficient w Infer the missingness mask m -aa repeat
Sample noise z ~~ I)
Compute perturbed data: Xobs Xobs if missing mechanism is MAR then
Estimate P|m = 0xx = Xoia] and IP P'mm"
= 0|x = Xabs. using fitted logistic models
end if
Update parameters using gradient descent on Ve (Eq.(4) (4) + w Eq.(5)) until convergence
where the antithetic samples are defined as x+ x + oz and x = x - oz. Here, sb = + s1s1 and L (ab(xt) + ++(()) - (x)(
B. Data Generation under Different Missing Mechanisms
Missing data mechanisms can vary significantly, but they are typically categorized into three main types defined by Rubin (1976): missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR). In our experiments, we simulate missing data based on these mechanisms as follows:
MCAR (Missing Completely at Random): Missing values are generated uniformly, with each data point having an equal probability being missing, determined by a predefined missing rate, a. Specifically, missing values are generated using a Bernoulli distribution, Ber(a:), where each entry missing independently with probability a.
MAR (Missing at Random): In this scenario, missing values are generated using logistic model. A random subset of the variables is selected to remain fully observed, while the remaining variables have missing values depending on the fully observed ones. The missingness is determined by a logistic model with random coefficients, scaled to achieve the target proportion of missing data for the variables influenced by the fully observed subset.
MNAR (Missing Not at Random): The MNAR mechanism modeled using logistic masking model. implements two mechanisms and in either case, weights are random and the intercept selected to attain the desired proportion of missing values.
Missing probabilities for each variable are determined by a logistic model that takes all the variables (including those with missing data) inputs;
Variables are split into two sets: a set of input variables for the logistic model and a set of variables whose missingness is determined by the logistic model. The input variables are masked using MCAR process, meaning the missingness in the second set depends on the missingness in the input set.
In all experiments, for MAR missing mechanism, use logistic regression to estimate the likelihood of each data point being observed. For MNAR, we utilize the same training objective as MCAR, while recognizing that this method may introduce some bias. The algorithm for training the first- and second-order models under different missing mechanisms is outlined in Algorithm 1.
Scalability and numerical stability for MAR We have included additional synthetic results for the MAR mechanism in Table following the same setup as described in Section 4.2. The results exhibit similar pattern, though slightly worse than MCAR, which may be attributed the inverse probability, potentially increasing variance and leading to less stable estimations. However, the estimates remain reasonable and demonstrate empirical performance closely aligned with the ground truth, along with variance reduction. Furthermore, to illustrate the impact of potential model misspecification by the logistic regression model in MAR, we conduct an experiment comparing ground truth p with the estimated ý from the logistic regression model on the same synthetic dataset, with varying missing ratios. When comparing the performance shown in Tables and 4, although potential model misspecification in the logistic regression model influences the estimation of missing probabilities, this impact kept within reasonable range in the variance-reduced version.
19

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
Algorithm 2 MissScore for Sampling with Missing Data
Input: Score models S (-; 0*), S2(: )" step size number of iterations ' Initialize: Xy ~ ((x) = 1 to T
Sample noise 2z (0,1)
if Ozaki sampling then
Compute M (0 1) (x,()111)
Compute -t (e (c20k--1 1 1) (x
Update X, X,-1 1 (1--) 54-122 2t
else
Update X, = 81-1 + [ (x1 + Vezi
end
Return: Xyr
Table 3. Mean squared error (MSE) between the estimated first-order and second-order scores and the ground truth is evaluated across 5,000 test samples. We vary the noise scales and missing ratios a, with each configuration tested using 10 random seeds. MAR with estimated missing probability á by logistic model.
a 0.0 (Complete data)
a 0.1
O 0.3
at 0.5
Methods
o 0.1
σ 0.01
0.1
0.01
0.1
o 0.01
0.1
0 0.01
$1
0.28 + 0.01
0.42 = 0.02
0.32 +: 0.00
0.51 + 0.02
0.36
0.02
0.52
0.00
0.45 = 0.01
0.56 + 0.03
8y (VR)
0.07 0.00
0.07 + 0.00
0.11 0.00
0.13 + 0.01
0.15
1.5+0.00 0.01
0.16 0.02
0.22 + 0.02
0.24 0.04
$2
0.16 + 0.02
15.42 = 0.47
0.28 0.02
31.22 1.08
0.30 = 0.04
34.24 4 5.14
0.36 = 0.04
35.41 : 1.03
(RR)
0.04 0.00
0.05 + 0.00
0.04 0.00
0.05 + 0.00
0.06 0.00
0.06 0.00
0.08 + 0.01
0.07 0.00
Table 4. Mean squared error (MSF) between the estimated first-order and second-order scores and the ground truth is evaluated across 5,000 test samples. We vary the noise scales and missing ratios o, with each configuration tested using 10 random seeds. MAR with ground truth missing probability p.
a 0.0 (Complete data)
Q 0.1
a 0.3
a 0.5
Methods
a = 0.1
σ 0.01
0 0= 0.1
= 0.01
0 0=1
a= = 0.01
0 = 0.1
a = 0.01
81
0.28 0.01
0.42 + 0.02
0.24 0.01
0.42 + 0.02
0.27 0.01
0.41 0.01
0.33 + 0.01
0.43 0.04
si (VR)
0.07 +: 0.00
0.07 = 0.00
0.06 0.00
0.06 + 0.00
0.09 -.090.00
0.10 0.00
0.19 = 0.00
0.22 + 0.01
82
0.16 + 0.02
15.42 + 0.47
0.24 0.01
16.66 = 4.34
0.27 - 0.01
41.47 +1777..5 7.28
0.33 + 0.01
39.89 4.03
R2(VR)
0.04 0.00
0.05 + 0.00
0.02 0.00
0.03 0.00
0.03 0.00
0.03 0.00
0.05 + 0.00
0.05 0.00
C. Additional Information on Sampling
In the sampling experiments with the Swiss-Roll dataset under MCAR, we use a small perturbation o = 0.01 and jointly optimize Eq. (7) and (8), where s (x) R 81 (x) and s2(x) R 82 (x). The sample size is set to 5000. Both s (x:; 0) and s2 (x; 0) are modeled using 3-layer MLP with a latent size of 128 and a Softplus activation function. We use learning rate of 0.001, batch size of 64 and train for 100 epochs, which takes approximately 4 minutes Gold 6448H CPU. The experiments in Section 4 also utilize the same model configuration for training. In the Ozaki sampling experiments, we only use the diagonal of s2 (x; 0) to avoid the computational costs associated with the inversion, exponentiation, and decomposition of 52(:; 0). The algorithm is presented in Algorithm 2.
Intel(R) Xeon(R)
The following sections provide experimental details on data generation with the simulated Bayesian Network and real Census data.
C.1. Dataset Description and Processing
Bayesian Network Details regarding the data generated from a Bayesian Network can be found in Section B. in (Ouyang et al., 2023).
20

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
Census Census dataset is binary classification dataset that predict whether income exceeds 50K/yr based census data (Kohavi, 1996). Also known Adult dataset.
The statistical information of datasets used in our experiments is in Table 5. #train, #test, #continuous, and #categorical mean the number of training data, testing data, continuous columns, and categorical columns, respectively.
Table S. Synthetic and Real-World Datasets Used in Experiments.
Dataset
#Train
#Test
#Categorical
#Continuous
Bayesian Network
2000
20000
3
2
Census
16000
4000
9
For data processing, we follow standard pre- and post-processing procedures for mixed-typp tabular data. Specifically, we apply min-max normalization to continuous variables and reverse this scaling during generation. For discrete variables, we use one-hot encoding and apply a rounding function after the softmax function during generation.
C.2. Evaluation Methods
We adopt the "train on synthetic, test on real (TSTR)" framework (Esteban et al., 2017). a widely used method for assessing the quality of sampling data from generative model (Kim et al., 2022; Ouyang et al., 2023; Li et al., 2019). The experimental results for sampling in this paper are calculated as follows:
1. We first download dataset and use its existing train-test split.
2. Then we generate synthetic records equal in number to the original training set using various synthetic data generation methods.
3. Using the synthetic training records from Step 2, we train base classifiers to make predictions. We conduct hyperparameter search for each classifier, considering Decision Tree, AdaBoost, Logistic Regression, MLP Classifiers, Random Forest, and XGBoost for the classification tasks. The hyperparameters and their candidate settings follow those described in (Kim et al., 2022; Ouyang et al., 2023), and are summarized in Table 26 of (Kim et al., 2022).
4. Finally, we evaluate the classifiers using the testing dataset, applying range of evaluation metrics for comprehensive assessment.
Steps 2 to repeated three times for each dataset, and the average scores for each method across all evaluation metries calculated. The detailed metrics used in experiment include:
1. Accuracy:
This is calculated using the accuracy score function from the sk wee.leenetreiic.... module.
2. Weighted-F1:
Weighted-F1 = Σ wis 1=
where N is the total number of classes. The weight for the i-th class, w, = |ー with p; representing the proportion of the i-th class's size relative to the total dataset. Here, is the Fl score for the i-th class, calculated using the One-vs-Rest strategy. This weighting approach is designed to prioritize the evaluation of synthesized tables by giving more importance to smaller classes, which are often prone to being overlooked by the model, thus addressing mode collapse.
3. AUROC: This is calculated using the roc_auc_score function from the sklearm rn. eelearn.merriaae as module.
4. SDMetrics: This metric evaluates synthetic data by comparing it against the real data, as described (Dat, 2023).
Among metrics, higher score indicates better overall quality of the synthetic data.
21

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
C.3. Model Architecture
We a perturbation of a 0.1 and jointly optimize Eq. (7) and Eq. (8). the Bayesian Network experiment, we follow the same configuration described earlier, with each rum taking approximately 20 minutes. For the census dataset, employ a simple MLP consisting of 5 Linear layers, LeakyReLU activation, Layer Normalization, and Dropout with probability of 0.2 in the first layer. The learning rate is set to 0.001. The first two layers use latent size of 128, while the last three layers use a latent size of 1024. We train with batch size of 256 for 250 epochs, with each experiment taking approximately 4 hours. All experiments are performed on an Intel(R) (R) Gold 6448H CPU. For the downstream classifier, we use the same base hyperparameters as listed Table 26 of (Kim et al., 2022).
C.4. Experimental Results
In the following experimental results, we use a missing data ratio of a = 0.3 and apply XGBoost for the downstream tasks, without delving into specific implementation details. Table 6 presents the utility evaluation MissScore using both Langevin and Ozaki samplings, compared other baseline methods, on the Census dataset with a missing ratio of 0.3. Additionally, Table 7 summarizes the Accuracy, AUROC, and Weighted-FI metrics as the missing ratio varies. Figure illustrates the fidelity evaluation of MissScore, again using Langevin and Ozaki samplings alongside other baselines, on the Bayesian dataset.
Table Utility evaluation of MissScore using Langevin and Ozaki samplings, along with other baselines, on Census dataset with missing ratio 0.3.
Criterion
Mechanism
Langevin
Ozaki
DSM-delete
DSM-mean
STaSy-mean
MCAR
0.80
0.81
0.70
0.75
0.77
Accuracy
MAR
0.82
0.82
0.69
0.77
0.74
MNAR
0.81
0.80
0.59
0.80
0.75
MCAR
0.84
0.84
0.57
0.67
0.62
AUROC
MAR
0.85
0.86
0.46
0.75
0.61
MNAR
0.86
0.86
0.52
0.76
0.63
MCAR
0.52
0.52
0.24
0.32
0..1
Weighted-Fl
MAR
0.61
0.60
0.32
0.38
0.38
MNAR
0.69
0.68
0.41
0.52
0.42
Table Evaluation of MissScore using Ozaki samplings the Census dataset with varying missing ratios (0.1, 0.5, 0.7.
Methods
0.1
0.3
0.5
0.7
0.9
MCAR
0.81
0.81
0.80
0.79
0.77
Accuracy
MAR
0.71
0.82
0.82
0.82
0.79
MNAR
0.80
0.80
0.83
0.74
0.72
MCAR
0.85
0.84
0.84
0.86
0.61
AUROC
MAR
0.85
0.86
0.87
0.85
0.83
MNAR
0.85
0.87
0.86
0.85
0.80
MCAR
0.54
0.52
0.41
0.66
0.22
Weighted-Fi
MAR
0.64
0.60
0.67
0.65
0.63
MNAR
0.46
0.68
0.61
0.64
0.63
C.5. Broader Empirical Validation
To further assess the robustness and generalizability of MissScore, we conduct a comprehensive empirical evaluation on 27 real-world tabular datasets covering both classification and regression tasks. The datasets, sourced from the UCI Machine Learning Repository and scikit-learn, encompass a wide range of characteristics, including varying sample sizes, feature dimensions, and class distributions. Following the experimental setup established in ForestDiffusion (Jolicoeu-MMartnneau Martineau
22

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
0.55
9.8
0.80
SDMetrics
SDMetrics
SDMetrics
0.75
0.73
Langevin
0.63
Langevin
0.43
Langevin
Ozaki
0.40
Ozaki
0.63
Ozaki
DSM-mean
meea
0.59
DSM-mean
STaSy-mean
0.53
SlaSy-mean
DSM-mean
0.1
0.3
0.5
0.7
0.9
0.30
Missing Ratio
Missing Ratio
Missing Ratio
(a) MCAR
(b) MAR
(c) MNAR
Figure 5. Fidelity evaluation of MissScore using Langevin and Ozaki samplings, along with other baselines, on the Census dataset varies with missing ratio {0.1, 0.3, 0.5, 0.73 under different missing mechanisms.
al., 2024), we introduce 20% missingness and benchmark MissScore against diverse set of generative baselines including copula-based models (GaussianCopula), VAE-based approaches (TVAE), acore-based methods (STaSy, TabDDPM), and diffusion- or flow-based architectures (Forest-VP, Forest-Flow). Full details on dataset statistics, evaluation metrics, and baseline configurations can be found in Appendices B.1 through B.3 in ForestDiffusion (Jolicocur-Mart) Martineau et al., 2024).
Table 8. Tabular data generation with incomplete data (27 datasets, 3 runs each, 20% missing); values are mean (standard-errow
Method
Wirain
Wiss
COU
F1
F Flaie
Phina |
GaussianCopula
2.60 (0.58)
(0.63)
0.20 (0.04)
0.43 (0.05)
(0.18)
0.48 (0.06)
0.60 (0.04)
(1.00)
0.21 (0.08)
TVAE
2.17 (0.60)
2.40 (0.65)
0.32 (0.04)
0.63 (0.04)
-0.66 (0.95)
0.55 (0.08)
0.45 (0.01)
4.04 (2.30)
0.29 (0.09)
Stasy
3.40 (1.37)
3.67 (1.40)
0.38 (0.05)
0.63 (0.06)
0.27 (0.28)
0.64 (0.06)
0.46 (0.02)
1.09 (0.22)
0.36 (0.10)
TabDDPM
4.36 (1.89)
4.80 (1.90)
0.72 (0.06)
0.71 (0.06)
0.58 (0.11)
0.61 (0.10)
0.42 (0.02)
1.16 (0.35)
0.56 (0.10)
Forest-VP
1.84 (0.51)
2.14 (0.56)
0.53 (0.04)
0.78 (0.03)
0.53 (0.10)
0.71 (0.04)
0.42 (0.01)
1.16 (0.30)
0.43 (0.12)
Forest-Flow
1.82 (0.51)
2.12 (0.56)
0.67 (0.03)
0.84 (0.03)
0.55 (0.11)
0.69 (0.04)
0.43 (0.01)
1.16 (0.32)
0.50 (0.10)
MissScore-Langeee
1.85 (0.52)
2.15 (0.57)
0.65 (0.04)
0.82 (0.03)
0.54 (0.11)
0.68 (0.04)
0.44 (0.02)
1.18 (0.33)
0.49 (0.11)
MissScore-Ozakk
1.79 (0.50)
2.11 (0.55)
(0.03)
0.83 (0.03)
(0.10)
0.70 (0.04)
0.43 (0.02)
1.15 (0.30)
0.52 (0.10)
The full quantitative results are presented in Table 8. MissScore- Ozaior achieves the lowest Wasserstein distance both training and test distributions, indicating that the generated samples closely match the real data in distributional geometry. It also performs competitively on downstream metrics such as fuke and fake, which measure how well synthetic data supports regression and classification tasks. Across most metrics, both LeeSScore-Langevie and MissScore-Ozaki demonstrate performance comparable to or better than strong baselines, including TabDDPM and Forest-Flow, while maintaining robustness under missing data.
every individual metric, we observe that it excels in several recurring settings. First, datasets with well-separated class structures and moderate sample sizes (e.g., seeds), MissScore is able to model the score field accurately and generate high-quality data that supports effective downstream classification. Second, on datasets with balanced but complex class relationships (e.g., parkinsons), it captures nuanced inter-feature dependencies better than methods relying on shallow generative priors. Finally, in moderate-to-hip dimensional datasets (e.g., qsar biodegradation), the use of second-order score matching enables MissScore to model non-linear feature interactions more effectively than first-order methods.
While MissScore does not outperform all baselines
However, the method's limitations emerge in data-scarce regimes. On very small datasets like concrete slump (with only 103 samples), high model flexibility can induce overfitting, undermining score stability. Similarly, for datasets with many classes and few samples per class, such as libras, the variance in score estimation grows substantially, leading to degraded generation quality. These limitations point to directions for future work, including lightweight score architectures and adaptive regularization under low-data conditions. Overall, this extensive evaluation suggests that MissScore offers a strong and principled approach to learning from incomplete tabular data, especially in regimes where capturing higher-order structure crucial.
23

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
D. Additional Information on Causal Discovery
D.I. Related Work
Causal Discovery with Complete data, Causal discovery aims to uncover the underlying causal relationships among variables of interest from purely observational data, specifically identifying causal Directed Acyelic Graph (DAG) for a given dataset. This problem lies at the heart of causal inference, as knowledge of the causal graph enables prediction of the effects of interventions. However, causal discovery from observational data is inherently ill-posed, necessitating additional assumptions, such imposing functional assumptions the data-generating process. We adopt the notion of structural causal model (SCM) to characterize the causal relations among variables. Each SCM M (Z. consists of the exogenous variable set 2- = 171, .-, Za), the endogenous variable set X' = (X1, X Xa). and the function set F= (f1, f2.... :........)......... Here, each function fi computes the variable X; from its parents (or causes) Xps, and an exogenous variable Z, i.e., X; = PA Z.). We focus on a specific class of SCMs, called the additive noise models (ANMs), given by X; = f((P...) + Z;
1,2,..... .d, where Z₁, interpreted as the additive noise variable, is assumed to be independent of variables in XPA, and mutually independent with variables in 2 Zz.
Rolland et al. (2022) proposed an order-based algorithm for this model, further assuming that Si is twice-differentiable nonlinear function and Z; is Gaussian noise. This method enables the identification of leaf nodes based on the diagonal of the Hessian of the log-likelihood. Before proceeding to the method for identifying leaves, we first derive an analytical expression for the score following Lemma in Rolland et al. (2022). The score is written as follows:
d I p( ) |
V 0, logp(x) p(x) = Vz, log 1=1
V 33
log p(a: xpp.))
i=1
log p( - fi)
(where 2, = x; - (i((p)))))
3 p(x, (x(AA))
0f. d llog (1(PA).)
0xj
dz,
Ōx
VECH,
where CH, represents the children of the variable j. As a result, Be, Va, logp(x) a, where a is constant, Consequently, the variance of the diagonal elements of the Hessian is zero (i.e. Varx [H (log p(x)) 0) if and only if node j is a leaf node.
However, existing computational methods for calculating the Hessian struggle to scale efficiently as the number of variables and samples increases, limiting the scalability of Rolland et al. (2022). To address this, Sanchez et (2022) introduced a diffusion-based model that efficiently computes the Hessian, enabling the method scale to larger datasets, both in terms of sample size and number of variables, while maintaining comparable performance to Rolland et al. (2022).
Causal Discovery with Incomplete Data. Several extensions of the PC algorithm have been developed to learn causal graphs from incomplete data (Tu et al.. 2019; Gain & Shpitser, 2018), utilizing only the fully observed samples while mitigating biases in conditional independence tests. Another prominent family of methods relies on Expectation-Ne Maxiiiiiiinaiiiiiiiiiiiion (Dempster et al., 1977), where missing values are iteratively inferred while simultaneously learning the causal structure. Building on the continuous optimization techniques introduced by NOTEARS (Zheng et al., 2018), MissDAG (Gao et al., 2022) extends this approach to continuous identifiable Additive Noise Models (ANMs), using approximate posterior inference via Monte Carlo and rejection sampling when the exact posterior is unavailable. MissOTM (Vo et al., 2024) introduces a score-based method that leverages optimal transport to learn causal structures from incomplete data. A distinct approach is taken by VISL (Morales-Ae Alvarez et al., 2022), which employs amortized variational inference in a Bayesian framework. Unlike MissOTM and MissDAG, VISL assumes a latent low-dimensional factor that captures the essential structure of the data based on observed variables. The latent factors are then used to reconstruct the complete data and discover the underlying causal graph.
24

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
D.2. Evaluation Metries
For each method, we compute the
SHD. Structural Hamming distance between the output and the true causal graph, which counts the number of missing, falsely detected, reversed edges.
Order Divergence. (Rolland et al., 2022) propose this quantity for measuring how well the topological order is estimated. For an ordering т, and a target adjacency matrix A, we define the topological order divergence Doo (T, A) as
Diop (T,, A)
Σ
Σ A ==1 jo
(18)
D.3. Model Architecture
We apply a perturbation of a = 0.1 and jointly optimize Eq. (7) and Eq. (8). The model is a simple MLP with 5 Linear layers, LeakyReLU activation, Layer Normalization, and Dropout rate of 0.2 in the first layer. The learning rate set to 0.001. The first two layers have latent size of maxx1188 3 x d), while the last three Training is conducted with batch size of 128 for 150 epochs. The time efficiency is shown in the figure 4 with ER graph model across various dimensions. All experiments are executed an Intel(R) Xeon(R) Gold 6448H CPU. The algorithm is presented in Algorithm D.3.
latent size of max (1024,
d).
Algorithm 3 MissScore for Causal Discovery with Missing Data
1: Input: Observed data Xobs: score models (1(:;) s2(:;
2: Initialize: π = I: nodes = {1....,d}
3: n, <- shape Xahs
4: for
do
5: Jointly train the score models si (0) and s2 (0) using Xabs with Algorithm I
6: Generate samples Nnew using Algorithm 2 with bootstrapping
7: Estimate the second-order score $2 (Xnew) using $2 (0)
8:
Vj Varx [diag( (s (2xx))))))
9:
i < arg mindes V, The leaf node
10: π + [,, Update topological order
11: nodes < nodes {E) Remove node !
Remove the 8-th column from Xahs
13: end for
14: Obtain the final DAG using CAM pruning associated with the topological order
D.4. Experimental Results
To evaluate MissScore in the context of causal discovery under missing data, we compare it with several satee-of-the-ar methods. These include MissDAG (Gao et al., 2022) and MissOTM (Vo et al., 2024), which are specifically designed for causal structure learning with incomplete data. In contrast, we also include two imputation-basee pipelines that rely on MissForest (Stekhoven & Bühlmann, 2012) complete the data prior to structure learning: DiffAN (Sanchez et al., 2022) (denoted MissDiffAN), diffusion-based score matching approach, and DAGMA (Bello al., 2022) (denoted as MissForest), a continuous optimization method. Together, these baselines represent broad range of causal discovery strategies, including continuous optimization, optimal transport, score-based methods, and imputation-then-leee pipelines.
We begin by conducting experiments using complete data, evaluating our approach alongside various missing mechanisms and different missing ratios. Additionally, we include order divergence in the Table. Our findings reveal that, with complete data, MissScore performs on par with DiffAN, as illustrated in Figure Notably, among all settings, MissScore achieves performance similar current state-of-the-e approaches while offering superior computational efficiency in terms of memory and time. This scalability key advantage where other methods may struggle.
25

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
175
SHD
20
Dimension (d)
Dimension (d)
Figure 6. The data is generated using an ER graph model with different dimensions = {10, 20, 50} and an equal number of edges. Each dataset consists of 1000 samples. Left:
MLP; Right: corresponds MIM.
139
330
24
Dimension (d)
Dimension (d)
VInEAG
Maafora
Muscore
Figure 7. The data is generated under MCAR with missing ratios of 0.1 and 0.3, using ER graph model with different dimensions d= {10, 20, 50} and an equal number of edges. Each dataset consists of 1000 samples. ft corresponds to MLP. Left: SHD with missing ratio 0.1; Right: SHD with missing ratio 0.3.
10
Dimension (d)
Dimension (d)
Figure 8. The data is generated under MCAR with missing ratios of 0.1 and 0.3, using an ER graph model with different dimensions {10, 20, 50} and equal number of edges. Each dataset consists of 1000 samples. f. corresponds MIM. Left: SHD with missing ratio 0.1; Right: SHD with missing ratio 0.3.
26

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
190
23
Dimension
Dimension
Misssoe
MINOTT
HI 200
Wisdfrrst
MINOTM
Wisdhrest
Figure 9. The data is generated under MAR with missing ratios of 0.1 and 0.3, using ER graph model with different dimensions {10, 20, 50} and an equal number of edges. Each dataset consists of 1000 samples. ft corresponds to MLP. Left: SHD with missing ratio 0.1; Right: SHD with missing ratio 0.3.
Dimension (d)
Dimension (d)
Figure 10. The data is generated under MAR with missing ratios of 0.1 and 0.3, using an ER graph model with different dimensions {10, 20, 50} and equal number of edges. Each dataset consists of 1000 samples. fs corresponds MIM. Left: SHD with missing ratio 0.1; Right: SHD with missing ratio 0.3.
23
36
26 Dimension (d)
50
Dimension (d)
MuOOTM
Phaferal
Mus
MuOOTM
Maddad
Hua DAV
Mhoffrrod
Figure
is generated under MNAR with missing ratios of 0.1 0.3, using an ER graph model with different dimensions d = {10, 20, 50] and equal number of edges. Each dataset consists of 1000 samples. ft corresponds to MLP. Left: SHD with missing ratio 0.1; Right: SHD with missing ratio 0.3.
27

---

MissScore: High-Order Score Estimation in the Presence of Missing Data
SHD
72-
60
26 Dimension
Dimension (d)
MisScoree
MESCTM
Minifrrst
Missooee
MISOOTH
145000
Wisherett
Figure is generated under MNAR with missing ratios of 0.1 0.3, using an ER graph model with different dimensions d= = {10, 20, 50} and equal number of edges. Each dataset consists of 1000 samples. fa corresponds to MIM. Left: SHD with missing ratio 0.1; Right: SHD with missing ratio 0.3.
Table Order divergence with missing ratios of {0.1, different missing data mechanisms. The ER graph model considered with varying dimensions d = (10, 20, 50), and an equal number of edges. Each dataset consists of 1000 samples and fi corresponds MLP. Lower order divergence indicating better performance.
MCAR
MAR
MNAR
Dimensions
Methods
a = 0.1
a 0.3
0 0= 0.1
a = 0.3
0: = 0.1
a = 0.3
MissScore
1.62 - 0.99
1.60 + 0.66
1.33 + 0.47
1.89 + 0.75
1.33 + 0.47
1.44 = 0.50
d=10
MissDiffAN
2.00 1.18
2.90 + 1.87
2.75 + 1.30
3.60 + 1.28
2.60 + 1.20
2.30 = 1.35
MissScore
2.10
2.67 2.00
2.70 1.85
1.94 0.29
1.70 + 1.78
0.70 0.46
d=20
MissDiffAN
4.30 2.00
4.22 + 2.35
4.33 + 2.11
3.00 + 1.00
2.00 + 0.89
3.40 = 1.20
MissScore
3.40 2.33
4.00 3.10
3.10
4.60 2.91
2.80 0.98
4.10
3.30
d=50
MissDiffAN
4.50 2.91
3.60 3.32
8.33 3.65
7.40 # 2.84
3.50 + 1.75
3.20 = 3.28
Table 10. Order divergence with missing ratios of a = 0.3) across different missing data mechanisms. The ER graph model considered with varying dimensions d' (10, 20, 50), and equal number of edges. Each dataset consists of 1000 samples and fi corresponds to MIM. Lower order divergence indicating better performance.
MCAR
MAR
MNAR
Dimensions
Methods
0.1
0.3
0.1
0.3
0.1
0.3
MissScore
1.60 -:00 1.20
1.20 0.60
1.35
0.47
1.88 0.93
1.62 + 0.99
1.38 = 0.48
d=10
MissDifFAN
2.22 + 2.10
1.50 0.67
2.50
1.36
3.70 + 1.10
2.40 + 2.20
2.20
1.40
MissScore
1.70 0.90
2.20 1.33
1.82 1.29
1.90 + 1.30
1.56 + 0.68
0.56 = 0.83
d=20
MissDiffAN
2.70 1.49
3.8 1.47
3.70
1.35
3.10 2.21
2.44 + 1.07
2.70 = 1.68
MissScore
4.90 2.12
4.00 + 1.41
4.56 + 2.27
4.60 + 2.24
4.20 + 1.66
4.40 = 1.85
d=50
MissDiffAN
5.56
4.40 1.56
7.90 2.39
7.50 - 2.69
5.00 + 1.61
6.20 = 2.68
28


Reviews:
[{"review_id": "r2qqbMuINg", "reviewer": "Reviewer_ugL7", "summary": "The authors propose a higher-order score estimation method for missing data. They first introduce the correct learning objectives in Theorems 3.3 to 3.5 and then formulate a multi-task objective to learn both the first and second-order scores with variance reduction to improve stability. The authors validate the effectiveness of the proposed method on a Swiss Roll dataset, a simulated Bayesian network dataset, and a Census dataset, comparing it against two baseline methods. The method successfully recovers the first and second-order oracle scores and generates high-quality samples. The authors further investigate the application of the proposed method for causal discovery and conduct experiments on an ER graph, comparing it against four baseline methods.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# To Reviewer ugL7\nWe appreciate your thoughtful comments and the time spent on reviewing our paper. Below, we summarize the concerns and suggestions and provide detailed responses to each.\n\n---\n\n### **Experiment on Real Dataset for Causal Discovery.**\nThe Sachs dataset [1] models a network of cellular signaling pathways with 11 continuous variables and 7,466 samples. It is a widely used benchmark for causal discovery, providing cell signaling data with established causal relationships. The following table presents the Structural Hamming Distance (SHD↓) on the Sachs dataset under MCAR missingness at varying rates (10%, 30%, and 50%):\n\n| Method         | 10%  | 30%  | 50% |\n|----------------|------|------|-----|\n| MissOTM        | 44   | 42   | 43  |\n| MissForest     | 22   | 26   | 26  |\n| MissDAG        | 23   | 21   | 20  |\n| MissDiffAN     | 24   | 23   | 21  |\n| **MissScore**  | 22   | 22   | 19  |\n\n---\n\n### **More Experiments on Real-world Datasets.**\nWe kindly refer to the **Broader Empirical Validation and Additional Comparisons** section in our response to **Reviewer [bpYG]**, where we address similar concerns.\n\n---\n\n### **Whether the theoretical results in Theorems 3.3–3.5 can be extended to the MNAR setting.**\n\nWe thank the reviewer for the insightful question! For Theorems 3.3-3.5, the equivalence between denoising score matching and true score matching holds under MCAR and MAR because missingness is **independent** of unobserved values. This critical independence allows us to estimate $𝔼_{m|x}[1-m]$ from observed data, which ensures unbiased score estimation through the DSM loss. However, under MNAR, the missingness **depends** on the unobserved values themselves, creating fundamental identifiability issues. Even when $\\mathbb{E}_{\\bar{x}|x, m}[1-m]$ is not zero, the DSM objective cannot be shown equivalent to the true score matching objective because the necessary quantities cannot be reliably estimated from observed data alone. Generalizing our theoretical results to MNAR therefore remains an open challenge. That said, we have also conducted experiments under MNAR and observed promising empirical results. We acknowledge this limitation in the conclusion and highlight it as an important direction for future work.\n\n---\n\n### **The data generation process and the definition of ground truth scores used in Fig.1 and Table.1.**\n\nThe data generation process is detailed in Sections 4.1, 4.2, and Appendix B. We generate data from a known multivariate Gaussian distribution $\\mathcal{N}(\\mu,\\Sigma)$, allowing the ground truth scores to be computed analytically as $s_1(x) =-\\Sigma^{-1}(x-\\mu)$ and $s_2(x)=-\\Sigma^{-1}$. Partially observed data is then generated under an MCAR mechanism, and the predicted scores are compared to these analytical ground truth values during evaluation.\n\n---\n\n### **Explanation of the DAG estimation process and its connection to ANM.**\nBriefly, our method leverages a structural property of additive noise models (ANMs), where the diagonal of the log-density Hessian reveals information about causal structure [3]. Specifically, for a leaf node $x_j$, the diagonal Hessian term: $\\frac{\\partial^2 \\log p(x)}{\\partial x_j^2}$ has zero variance across the data distribution. MissScore leverages this by estimating second-order scores directly from incomplete data, avoiding imputation. The DAG proceeds in three steps:\n\n1. Score Estimation: We train first- and second-order score models using our high-order DSM objective with missing data handling.\n2. Node Ordering: At each iteration, the variable with the lowest variance in its diagonal Hessian across samples is selected as a leaf and added to the causal ordering.\n3. DAG Construction: A full topological order is recovered and pruned using the CAM pruning [4] to infer edge directions and finalize the DAG.\n\nWe will revise Section 6 to more clearly summarize the DAG estimation procedure and its connection to the ANM, while referring readers to Appendix D for full details, including the formal ANM formulation, the role of second-order scores in identifying leaf nodes, the complete causal discovery pipeline (Algorithm 3), evaluation metrics, and extensive results under various missingness mechanisms.\n\n---\n\nThank you for reading our rebuttal! We hope the responses provided have fully addressed your concerns. If any questions remain or further clarification is needed, please feel free to let us know. Thank you again for your time, thoughtful feedback, and valuable insights!\n\n---\n\n### **References**\n\n[1] Sachs et al. (2005). *Causal protein-signaling networks from single-cell data.*\n\n[2] Jolicoeur-Martineau et al. (2024). *Generating and imputing tabular data via diffusion and gradient-boosted trees.*\n\n[3] Rolland et al. (2022). *Score matching enables causal discovery of nonlinear additive noise models.*\n\n[4] Bühlmann et al. (2014). *CAM: Causal additive models, high-dimensional order search, and penalized regression.*", "claims_and_evidence": "I believe the theorems are well proved.", "methods_and_evaluation_criteria": "Yes. MSE for whether the score function is learned, fidelity and utility for generation quality, and SHD for causal discovery.", "theoretical_claims": "I check the proof in Appendix A.", "experimental_designs_or_analyses": "See weakness part below.", "relation_to_broader_scientific_literature": "I believe that studying higher-order score estimation for missing data could improve current methods for imputation and generation.", "essential_references_not_discussed": "The literature is discussed comprehensively.", "other_strengths_and_weaknesses": "Strengths: the topic is interesting and the derivation is innovative; the paper is well-written\n\nWeakness: the experimental results are somehow weak to support the effectiveness of the proposed methods. All of the experiments are simulations or small datasets. It would be better for authors to conduct experiments on one real dataset for causal discovery and conduct more real experiments for score estimation.", "other_comments_or_suggestions": "From the background in Section 6, I understand that the proposed method can avoid the two-step estimation process, where data is first imputed and then used for causal discovery. However, it would be clearer if the authors could provide a more detailed explanation of the process for estimating the DAG and its relationship with ANM.", "questions_for_authors": "1. For Theorems 3.3–3.5, do the authors consider it difficult to generalize the conclusions to MNAR? If I understand correctly, all we do is incorporate $\\mathbb{E} _{\\tilde{x} | x, m}$ into the least squares loss. If $\\mathbb{E} _{\\tilde{x} | x, m}[1 - m]$ is not zero, then the denoising score matching objective is also equivalent to the score matching objective. \n\n2. I cannot find the discussion of the data generation process for the experiments used in Figure 1 and Table 1. I am wondering what the ground truth score used in this experiment is. How is the ground truth score for missing data defined?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "p0YTzrDTWX", "reviewer": "Reviewer_Bf4B", "summary": "They introduce a method for estimation high-order scores in the presence of missing data with either MAR or MCAR. For MAR, they need to make use of logistic regression to estimate a needed quantity. They propose stability improvements with variance reductions.\n\nThey show low error in the estimation of the first and second-order score with extra missing data.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# To Reviewer Bf4B\n\nWe appreciate the reviewer’s thoughtful comments and the time spent on reviewing our paper! Below, we summarize the concerns and suggestions and provide detailed responses to each.\n\n---\n\n### **Faster Convergence Rate, Not Faster Wall-Clock Time.**\n\nOzaki sampling achieves **faster convergence rate** in terms of **iteration count**, as supported by theoretical results [1] under certain conditions and validated by our experiments (e.g., Figure 2). However, we acknowledge that this does not necessarily correspond to faster wall-clock time due to the additional computational overhead. To assess this, we measured the average wall-clock training time per iteration and the results are summarized below:\n\n| **Sampling Methods** | **Time (s)** |\n|---------------------|--------------|\n| Langevin | 0.00037 ± 0.00017 |\n| Ozaki | 0.00043 ± 0.00053 |\n\nAs expected, Ozaki sampling incurs a longer wall-clock time compared to Langevin. In our experiments, we consider Ozaki sampling with $s_2$ using only its diagonal, as described in Appendix C. This diagonal approximation reduces the computational overhead, making Ozaki sampling comparable to Langevin in terms of wall-clock time per iteration while retaining its advantage of faster convergence. We will clarify this in the revision that our claims of faster convergence (or speed) refer strictly to iteration count, not wall-clock time. Thank you again for your valuable feedback!\n\n--- \n\n### **Higher Quality of Ozaki Sampling.**\n\nOur wall-clock time analysis shows that Ozaki sampling incurs higher computational cost per iteration compared to Langevin. Despite this, Ozaki sampling often achieves better sample quality. Extensive experiments on 27 real-world datasets demonstrate that Ozaki sampling consistently produces higher-quality samples than Langevin for the same number of iterations. Even when we increased Langevin's sampling time to match that of Ozaki (MissScore-Langevin2), the improvements were modest, with Ozaki still outperforming the longer-running Langevin across all metrics — further emphasizing the advantage of incorporating second-order information. \n\n| Method | W_train ↓  | W_test ↓   | cov_train ↑ | cov_test ↑ | R²_fake ↑  | F1_fake ↑   | F1_disc ↓   | P_bias ↓   | cov_rate ↑  |\n|-----------------------|------------|------------|-------------|------------|------------|-------------|-------------|------------|-------------|\n| MissScore-Langevin|1.85(0.52)|2.15(0.57)|0.65(0.04)|0.82(0.03)|0.54(0.11)|0.68(0.04)|0.44(0.02)|1.18(0.33)| 0.49(0.11)|\n| MissScore-Langevin2|1.83(0.50)|2.13(0.55)|0.67(0.04)|0.83(0.05)|0.55(0.11)|0.69(0.05)|0.44(0.01)|1.19(0.34)| 0.50(0.11)|\n| MissScore-Ozaki|1.79(0.50)|2.11(0.55)|0.70(0.03)|0.83(0.03) | 0.58(0.10) | 0.70(0.04) | 0.43(0.02)|1.15(0.30)|0.52(0.10)|\n\n*Table 1: Tabular data generation with incomplete data across 27 datasets (3 experiments per dataset, 20% missing values). MissScore is evaluated using Langevin and Ozaki with the same number of iterations, and Langevin2 with increased sampling iterations (same sampling time as Ozaki); results show the mean (standard error).*\n\n---\n\n### **A broader benchmark incorporating multiple datasets and advanced metrics.**\n\nWe kindly refer to the **Broader Empirical Validation and Additional Comparisons** section in our response to **Reviewer [bpYG]**, where we address similar concerns.\n\n---\n\n### **Beyond Langevin vs. Ozaki: Broader Applicability of Second-Order Information with Missing Data.**\n\nThank you for your comment. We would like to clarify that:\n- Our primary goal is to address the challenge of missing data through the development of a high-order score estimation framework under partial observability, which avoids the need for full data or expensive automatic differentiation. \n- This framework enables downstream tasks such as causal discovery, which requires second-order information for capturing complex structural dependencies. \n- While our method also leads to better sampling performance (e.g., Ozaki), this is one of the benefits brought by second-order estimation, and not the primary focus of our work.\n\n--- \n\n### **MissForest vs. MissScore.**\n\nMissForest captures dependencies through iterative imputation, making it a strong baseline for causal discovery with missing data. However, it becomes computationally expensive as dimensionality increases (as shown in Fig. 4b). In contrast, our method achieves comparable performance in simulations while being more efficient. Our new experiment on the Sach dataset (response to Reviewer [ugL7]), shows that MissScore maintains stable performance even with more complex dependencies, outperforming MissForest in this regard.\n\n---\n\n### **Other Suggestions.**  \nWe will update the manuscript to use \"Machine Learning (ML) efficiency\" and rasterize the Swiss roll figures to resolve the PDF rendering issue.\n\n---\n\n### **Ref.** \n\n[1] Dalalyan et al. (2017). *User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient.*", "claims_and_evidence": "- novel score-based framework for learning high-order scores in the presence of missing data (true when there is missing data)\n- show efficiently and accurately approximates high-order scores with missing data (true, good evidence)\n- improves both sampling speed and data quality in data generation tasks, with the quality of the generated samples validated across several downstream tasks (Ozaki Sampling does converge in fewer steps than Langevin, but they only show this a 2D toy setting; for quality, the evidence is weak, they need more experiments on different benchmarks)\n\nTo add more details on the faster speed, this seems wrong since Ozaki is probably slower in actually wallclock time. Please add wallclock or FLOPS or similar metrics to actually assess speed, not just in number of steps. Since second-order steps are more expensive.", "methods_and_evaluation_criteria": "For tabular data: Fidelity and utility are okay simple metrics, but they are limited. \n- SDmetrics is mostly 1D and some 2D metrics, its limited. \n- Utility is not the correct name, the Stasy people branded it this way, but this metrics already exist with different name. Its called machine learning (ML) efficiency, see https://arxiv.org/pdf/2209.15421. Please rectify.\n\nYou should consider metrics such as Wasserstein distance, see https://arxiv.org/abs/2309.09968 for a very strong benchmark with many metrics. I recommend using a wide benchmark like this containing many metrics with many datasets, right now only using the Census dataset is not enough. One dataset is very poor evidence because it could be fished for the one that shows the best result on your method.\n\nOnly using swiss roll to say it converge in fewer steps, is not enough.\n\nThe causal discovery experiments that uses MissForest as imputation for methods that don't handle missing data are good.", "theoretical_claims": "I haven't checked the details, but their error rate experiments seems to suggests that their derivations are correct given the low errors.", "experimental_designs_or_analyses": "\"Methods And Evaluation Criteria\"", "relation_to_broader_scientific_literature": "They take existing score estimation and adapt it for when there is missing data (theory and loss function), which is very useful.", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "The main issue right now is the experiments, if the authors can improve the depth of their experiments for both speed and quality (see my suggestions above) it would warrant increasing my score.\n\nI am a bit disappointing that the paper only focus on Langevin versus Ozaki sampling when they could derive a second-order diffusion or something like this. This would make the paper stronger for sure.", "other_comments_or_suggestions": "- please rasterize the figures of swiss rolls, there is too many points, its make the pdf unresponsive.", "questions_for_authors": ".", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ljl587jAYQ", "reviewer": "Reviewer_bpYG", "summary": "The paper introduces MissScore, a novel framework for high-order score estimation in the presence of missing data. Existing high-order score-based generative models assume data completeness, necessitating imputation when data is missing. MissScore extends Denoising Score Matching (DSM) to estimate high-order scores directly, eliminating the need for imputation and addressing a critical gap in the field. The authors establish strong theoretical foundations and derive objective functions tailored to different missing data mechanisms (MCAR, MAR). While the empirical evaluation is somewhat limited in the number of baselines and datasets, the results evidence MissScore’s effectiveness in sampling, data generation, and causal discovery, outperforming existing methods. This work represents a significant theoretical and practical contribution to score-based modeling.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# To Reviwer bpYG\n\nWe sincerely thank you for your valuable feedback and the time spent evaluating our work! Below, we summarize the main concerns and provide detailed responses to each point.\n\n---\n\n### **Related work on VAE- and diffusion-based methods.**\n\nThank you for the suggestions. We agree that these works are highly relevant and will update the related work section to include the provided references and a recent diffusion-based method [8]. The revised paragraph is:\n\n>Deep generative models have emerged as a powerful framework for learning from incomplete data. VAEs are widely used for handling incomplete data due to their flexibility in modeling complex distributions. MIWAE [2] extends the importance-weighted autoencoder to handle missing-at-random (MAR) data, yielding theoretically grounded single and multiple imputations. HI-VAE [1] addresses heterogeneous variables by assigning each data type a tailored likelihood, with a shared variational mechanism that balances learning across different attributes. VAEM [3] similarly supports mixed-type data using a two-stage process—first learning per-feature marginal VAEs, then a dependency VAE to model correlations among the latent codes—while hierarchical variants [4] further increase representational capacity at the cost of more challenging inference. Diffusion models have also gained traction for missing data imputation. TabCSDI [7] adapts conditional score-based diffusion to tabular data, addressing categorical variables and modeling the conditional distribution of missing features. NewImp [6] reinterprets diffusion-based imputation as a gradient flow in probability space, offering theoretical insights, while DIFFPUTER [5] combines the EM algorithm with diffusion models in an iterative framework. ForestDiffusion [8] enhances interpretability and efficiency by replacing neural score estimators with gradient-boosted trees, enabling direct learning from partially observed tabular data. While these methods can be used for imputation along with generation, our approach takes a different direction by directly estimating high-order score functions from incomplete data. This enables downstream tasks such as sampling and causal discovery without explicit imputation.\n\n---\n\n### **Broader Empirical Validation and Additional Comparisons.**\n\nWe evaluate our method on 27 real-world classification and regression datasets for the task of generation, sourced from the UCI Machine Learning Repository and scikit-learn. Details of the datasets, evaluation metrics, and more comparison methods are provided in Appendices B.1, B.2, and B.3 in ForestDiffusion [2]. This comprehensive evaluation benchmarks MissScore against established methods across diverse datasets, ensuring robustness and reliability.\n\n| Method         | W_train ↓ | W_test ↓ | cov_train ↑ | cov_test ↑ | R²_fake ↑ | F1_fake ↑ | F1_disc ↓ | P_bias ↓ | cov_rate ↑ |\n|----------------|--------|--------|----------|----------|--------|--------|--------|--------|---------|\n| Stasy          | 3.40   | 3.67   | 0.38     | 0.63     | 0.27   | 0.64   | 0.46   | **1.09** | 0.36    |\n| TabDDPM        | 4.36   | 4.80   | **0.72** | 0.71     | **0.58** | 0.61 | **0.42** | 1.16   | **0.56** |\n| Forest-VP      | 1.84   | 2.14   | 0.53     | 0.78     | 0.53   | **0.71** | **0.42** | 1.16   | 0.43    |\n| Forest-Flow    | 1.82   | 2.12   | 0.67     | **0.84** | 0.55   | 0.69   | 0.43   | 1.16   | 0.50    |\n| MissScore-Langevin | 1.85   | 2.15   | 0.65     | 0.82     | 0.54   | 0.68   | 0.44   | 1.18   | 0.49    |\n| MissScore-Ozaki    | **1.79** | **2.11** | 0.70 | 0.83     | **0.58** | 0.70 | 0.43   | 1.15   | 0.52    |\n\nTable 2. Tabular data generation with incomplete data (27 datasets, 3 experiments per dataset, 20% missing values); results show the mean.\n\n---\n\n### **Insufficient Support for Scalability Claims.**\n\nThank you for highlighting this point. We agree that the 100-dimensional Gaussian dataset alone may be too simplistic to substantiate a strong claim of high-dimensional scalability. In the causal discovery setting (Section 6), MissScore is evaluated on datasets with up to 50 variables—a scale typically regarded as high-dimensional in this context. We also reported results on 27 real-world tabular datasets, with dimensionalities reaching up to 90. While more extreme high-dimensional scenarios remain unexplored, we will revise the manuscript to more accurately reflect the scalability explored in our current evaluation.\n\n---\n\nThank you for reading our rebuttal! We hope the points discussed above have fully addressed your concerns. Please let us know if any questions remain or if further clarification is needed. We sincerely appreciate your time and thoughtful evaluation!\n\n---\n\n### **References**\n\nRefs [1]-[7] are the same as the list you provided. Apologies for the brevity due to word limits.\n\n[8] Jolicoeur-Martineau et al. (2024). *Generating and imputing tabular data via diffusion and flow-based gradient-boosted trees.*", "claims_and_evidence": "The paper presents strong theoretical and empirical claims regarding the feasibility and effectiveness of high-order score estimation with missing data. It establishes a rigorous mathematical foundation by deriving objective functions for high-order score matching under MCAR and MAR mechanisms, demonstrating that high-order scores can be estimated without imputation.\n\nEmpirically, the authors validate MissScore primarily on synthetic datasets, with only one real-world dataset for synthetic data generation. This is somewhat limited, as generative modeling with missing data is typically evaluated across multiple datasets (e.g., UCI). Additionally, key references and comparisons to deep generative models for missing data are missing, making it difficult to fully assess MissScore’s relative performance. Furthermore, while the paper claims high-dimensional scalability, the use of a 100-dimensional Gaussian dataset may be overly simplistic and insufficient to support this claim.\n\nOverall, expanding the empirical validation to more diverse real-world benchmarks and recent related work would strengthen the evidence.\n\nThe empirical analysis for causal discovery, however, appears reasonable.", "methods_and_evaluation_criteria": "For synthetic data generation, the authors assess MissScore’s effectiveness in estimating high-order scores and improving sampling efficiency. They compare Langevin dynamics and Ozaki sampling to demonstrate the advantages of incorporating second-order score information. However, as noted previously, only one real-world dataset is used in this setting, which limits the generalizability of the results. Expanding the evaluation to multiple real-world datasets, as is common in generative modeling research, would strengthen the empirical validation.\n\nFor causal discovery, the authors evaluate MissScore’s ability to recover causal structures from incomplete data. The method is compared against baselines such as MissForest + DiffAN/DAGMA and MissDAG/MissOTM, using Structural Hamming Distance (SHD) and Order Divergence as evaluation metrics. This part of the evaluation appears more comprehensive and methodologically sound.\n\nOverall, the evaluation methodology is valid and well-structured, but its impact is somewhat constrained by the limited number of real-world datasets used for validation.", "theoretical_claims": "The theoretical contributions are novel, well-structured, and mathematically sound, making this one of the strongest aspects of the paper. The results significantly advance the understanding of high-order score estimation in settings with missing data, and lay the groundwork for further exploration in score-based modeling or more complex missing patters.", "experimental_designs_or_analyses": "The experimental design is well-structured and appropriate for the problem setting, but broader empirical validation would provide stronger support for the method’s real-world applicability.", "relation_to_broader_scientific_literature": "The paper builds on foundational work in score estimation and extends recent advances in high-order score modeling to address missing data. It differentiates itself from imputation-based approaches (e.g., VAEs, GANs) and diffusion-based methods like MissDiff by directly estimating high-order scores from incomplete data. The discussion on VAEs for missing data does not reflect the state of the art. \n\nHowever, the paper lacks direct comparisons to several generative models designed for missing data, both in citations and empirical evaluation. Expanding the related work discussion and empirical benchmarking would better position MissScore within the broader literature.\n\nOverall, the paper makes a significant theoretical contribution, but a more comprehensive comparison to existing deep generative models would strengthen its impact.", "essential_references_not_discussed": "VAEs for missing data have been quite successful in the recent years, and they are briefly commented in the Introduction section as methods that require training additional networks, which is not fully grounded. Many recent work are ignored for assessing the empirical Below I include some unreferenced VAE-based and diffusion-based related works.\n\n\n### References\n\n[1] Nazabal, Alfredo, et al. \"Handling incomplete heterogeneous data using vaes.\" Pattern Recognition 107 (2020): 107501.\n\n[2] Mattei, Pierre-Alexandre, and Jes Frellsen. \"MIWAE: Deep generative modelling and imputation of incomplete data sets.\" International conference on machine learning. PMLR, 2019.\n\n[3] Ma, Chao, et al. \"VAEM: a deep generative model for heterogeneous mixed type data.\" Advances in Neural Information Processing Systems 33 (2020): 11237-11247.\n\n[4] Peis, Ignacio, Chao Ma, and José Miguel Hernández-Lobato. \"Missing data imputation and acquisition with deep hierarchical models and hamiltonian monte carlo.\" Advances in Neural Information Processing Systems 35 (2022): 35839-35851.\n\n[5] Zhang, Hengrui, Liancheng Fang, and S. Yu Philip. \"Unleashing the Potential of Diffusion Models for Incomplete Data Imputation.\" CoRR (2024).\n\n[6] Chen, Zhichao, et al. \"Rethinking the diffusion models for missing data imputation: A gradient flow perspective.\" Advances in Neural Information Processing Systems 37 (2024): 112050-112103.\n\n[7] Zheng, Shuhan, and Nontawat Charoenphakdee. \"Diffusion models for missing value imputation in tabular data.\" NeurIPS 2022 First Table Representation Workshop.", "other_strengths_and_weaknesses": "All strengths and weaknesses are already commented above.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "I do not have any further questions.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The manuscript proposes a framework for high-order score estimation applied to data with missing values within score-based modelling. The reviewers found the manuscript well-written and well-structured, and the theoretical contributions novel and sound. However, the experimental results are weak since they are limited to mainly simulated data or small datasets. The authors addressed this during the rebuttal by adding evaluations of multiple UCI datasets (these results must be included in the final revision of the manuscript). The reviewers also pointed to a lack of methodical and empirical comparison to other generative models designed for missing data. During the rebuttal, the authors proposed a revised version of the related works section, but we still urge the authors to also include these methods in the empirical evaluation. Given that the authors addressed most of the reviewer's concerns during the rebuttal period, we recommend accepting the paper if there is room in the program.
Decision: Accept (poster)