To Reviwer bpYG

We 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.

---

Related work on VAE- and diffusion-based methods.

Thank 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:

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.

---

Broader Empirical Validation and Additional Comparisons.

We 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.

Table 2. Tabular data generation with incomplete data (27 datasets, 3 experiments per dataset, 20% missing values); results show the mean.

---

Insufficient Support for Scalability Claims.

Thank 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.

---

Thank 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!

---

References

Refs [1]-[7] are the same as the list you provided. Apologies for the brevity due to word limits.

[8] Jolicoeur-Martineau et al. (2024). Generating and imputing tabular data via diffusion and flow-based gradient-boosted trees.

To Reviewer Bf4B

We 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.

---

Faster Convergence Rate, Not Faster Wall-Clock Time.

Ozaki 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:

As 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!

---

Higher Quality of Ozaki Sampling.

Our 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.

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).

---

A broader benchmark incorporating multiple datasets and advanced metrics.

We kindly refer to the Broader Empirical Validation and Additional Comparisons section in our response to Reviewer [bpYG], where we address similar concerns.

---

Beyond Langevin vs. Ozaki: Broader Applicability of Second-Order Information with Missing Data.

Thank you for your comment. We would like to clarify that:

• 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.
• This framework enables downstream tasks such as causal discovery, which requires second-order information for capturing complex structural dependencies.
• 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.

---

MissForest vs. MissScore.

MissForest 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.

---

Other Suggestions.

We will update the manuscript to use "Machine Learning (ML) efficiency" and rasterize the Swiss roll figures to resolve the PDF rendering issue.

---

Ref.

[1] Dalalyan et al. (2017). User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient.

To Reviewer ugL7

We 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.

---

Experiment on Real Dataset for Causal Discovery.

The 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%):

---

More Experiments on Real-world Datasets.

We kindly refer to the Broader Empirical Validation and Additional Comparisons section in our response to Reviewer [bpYG], where we address similar concerns.

---

Whether the theoretical results in Theorems 3.3–3.5 can be extended to the MNAR setting.

We 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.

---

The data generation process and the definition of ground truth scores used in Fig.1 and Table.1.

The 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.

---

Explanation of the DAG estimation process and its connection to ANM.

Briefly, 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:

1. Score Estimation: We train first- and second-order score models using our high-order DSM objective with missing data handling.
2. 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.
3. DAG Construction: A full topological order is recovered and pruned using the CAM pruning [4] to infer edge directions and finalize the DAG.

We 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.

---

Thank 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!

---

References

[1] Sachs et al. (2005). Causal protein-signaling networks from single-cell data.

[2] Jolicoeur-Martineau et al. (2024). Generating and imputing tabular data via diffusion and gradient-boosted trees.

[3] Rolland et al. (2022). Score matching enables causal discovery of nonlinear additive noise models.

[4] Bühlmann et al. (2014). CAM: Causal additive models, high-dimensional order search, and penalized regression.