Score Matching with Missing Data
Authors: Josh Givens
Song Liu
Henry Reeve
Venue: ICML 2025 oral (2025)
Keywords: Score Matching
Missing Data
Variational Inference
Primary area: general_machine_learning->unsupervised_and_semisupervised_learning
TL;DR: We adapt score matching to missing data and demonstrate its general applicability by applying it to graphical model estimation.

Abstract: Score matching is a vital tool for learning the distribution of data with applications across many areas including diffusion processes, energy based modelling, and graphical model estimation. Despite all these applications, little work explores its use when data is incomplete. We address this by adapting score matching (and its major extensions) to work with missing data in a flexible setting where data can be partially missing over any subset of the coordinates. We provide two separate score matching variations for general use, an importance weighting (IW) approach, and a variational approach. We provide finite sample bounds for our IW approach in finite domain settings and show it to have especially strong performance in small sample lower dimensional cases. Complementing this, we show our variational approach to be strongest in more complex high-dimensional settings which we demonstrate on graphical model estimation tasks on both real and simulated data.

Score Matching with Missing Data
Josh Givens Song Liu
Henry W Reeve 2
Abstract
density itself. Furthermore, the validity of the score match- ing objective itself requires only very mild assumptions of the family of proposed scores further ensuring this flexibil- ity. Alongside the classical method (Hyvärinen (005). var- ious adaptations of score matching have arisen to improve performance, decrease computational cost, and extend the approach to a wider range of settings (Hyvärinen, 2007
Score matching is a vital tool for learning the dis- tribution of data with applications across many areas including diffusion processes, energy based modelling, and graphical model estimation. De- spite all these applications, little work explores its use when data is incomplete. We address this by adapting score matching (and its major exten- sions) work with missing data in flexible set- ting where data can be partially missing over any subset of the coordinates. We provide two sep- arate score matching variations for general use, an importance weighting (IW) approach, and a variational approach. We provide finite sample bounds for our IW approach in finite domain set- tings and show it to have especially strong perfor- mance in small sample lower dimensional cases. Complementing this, we show our variational ap- proach to be strongest in more complex high- dimensional settings which we demonstrate on graphical model estimation tasks on both real and simulated data.
Vincent 2011 Song et al. 2020 Liu et el. 2022).
In this work, we extend the score matching framework to handle missing data training time. Specifically, we learn the full score function from partially missing multidimen- sional input data, a paradigm term missing score match- ing. Crucially, our approach is compatible with any pa- rameterised score model, enabling its application to both explicit score formulations and more general approaches such as neural networks (NNs). We propose two methods to adapt the original score matching method as well as its popular adaptations, truncated, sliced, and denoising score matching (Hyvärinen, 2007] Vincent] 2011 Liu & Wangl 2017 Song et al. 2020). These two distinct but closely re- lated methods complement each other allowing for wide range of problems to be tackled. The first method is a sim- pler importance weighting (IW) approach which we refer to as marginal IW score matching. For this method we obtain finite sample bounds in the bounded domain set- ting under certain conditions. We also provide experimen- tal results demonstrating its efficacy in lower dimensional settings and where less data is available. Our second ap- proach is a more computationally sophisticated variational approach which we refer to as marginal variational score matching. We demonstrate the efficacy of this approach in more complex, high dimensional settings by applying it to the problem of graphical model estimation with both real and synthetic datasets.
1. Introduction
Over the last decade, score matching has established itself as a powerful tool with downstream use in many areas of machine learning. Examples include: energy based mod- elling (Swersky et al. 2011: Bao et al. 2020 Li et al.,
2019b), mode-seeking clustering (Sasaki al. 2014), and
perhaps most prominently of all Diffusion processes (Song & Ermon 2019 Song et al.) 2021b Tashiro et al. 2021; Song et al. 202 Huang et al., 2021). Score matching aims to learn the score of a distribution which is the gra- dient of the log of the probability density function (PDF) (s(x)= = V log ((z)) In contrast to modelling the density directly, the score does not need to integrate to one meaning there is need to calculate a normalising constant. This allows it to be much more more flexibly modelled than the
In section 2 we discuss relevant works for score matching and related fields. In section 3 we will introduce our prob- lem more formally including score matching and any nota- tion used. Section[4wil be used to introduce our methods. Section 5 will present results on some real and simulated datasets. In Section 6 we give our conclusion.
'School of Mathematics, University of Bristol, Bristol, UK "School of Artificial Intelligence, Nanjing University, China, Cor- respondence to: Josh Givens <josh. givens @ briktol.ac.ukz.
2. Related Works
Proceedings of the 1nnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
While there has been some work which utilises score matching with missing data, these approaches mostly do so

---

Score Matching with Missing Data
exclusively through the lens of diffusion models. Specif- ically works such as MissDiff (Quyang et al. 2023) and Ambient Diffusson (Daras et al. 2023) require the score function itself to take the form of a neural network (NN) which learns the scores of the fully-observed and corrupted scores simultaneously. This prohibits their use in situations where our model for the score is some explicit parameteri- sation whose parameters we want to learn as is the case in settings such as energy based modelling Li et al.] (2023);
timisation procedure (as acknowledged in their follow on work), and |Bao et al. (2021) is only given for denoising score matching, not for classical or sliced score matching.
3. Setting
3.1. Notation
For n E N let [n] == {l....,, n]. For random variable Z we use supp(2) Z) for the support of Z. For f : Rd - R we write 0z == of
[ao et al. (2020): [Salimans & Ho (202 1) and Gaussian graphical models (Lin et al., 2016 Yu et al. 2018). Ambi- ent Diffusion also requires the data to be further artificially learning paradigm making both Ambient Diffusion and MissDiff subject to various levels of out of sample learning without specific adjustments for this phenomenon.
where a (x1,...,,ra) and ăzz Vaf(a) (0, fz))............................ ((f())., the gradient of f. For f : Rd - Rd take f(x), as the jth component of f(x) and write Va f(x) == f.)(a) +...... aaf(ar) Finally for a. b E R. take a o b to be the Hadamard product.
corrupted in order to create a
We now introduce some indexing notation which we will be using for RVs and functions throughout. This will prove useful when identifying the missing non-missing compo- nents of our data. Let Z be random variable taking values in Rdd We use Zj refer to the j component Z and for A C [d] take Zx = {Zz)jex- We use negation in indexing to mean the complementing coordinates. More precisely we let -j denote [d] \ {o} and let -À denote [d] A. We typically use 7 to denote an independent copy of Z. For a function f X + y and æx E Xx, x'xx E X_ we take [ax, xx²x) to be [(z) where
Looking more generally at distribution estimation with missing data, multiple works in the field of generative mod- elling have looked to tackle the problem of providing a gen- erative model for a distribution given corrupted samples from it. Prominent among these are MisGAN (Li et al., 2019a), which presents a marginalised GAN framework and MCFlow (Richardson al. 2020)- which presents EM like normalising flow framework. Neither of these approaches allow for flexible specification of parametric density estimate however with MCFlow requiring the den- sity to be a normalising flow and MisGAN having no model for the density whatsoever.
Tj ifj Z if
2j =
To our knowledge, the only approach which seems to adapt score matching to missing data in a parameter preserving manner is presented in (Uehara al., 2020) using an iter- ative EM-like procedure. However they themselves admit that there is little intuitive understanding of when this ap- proach will converge. Additionally, due to the nature of the score matching objective, the expectation step cannot be directly approximated using Monte Carlo estimation and instead requires fractional importance weighting, a method which employs nested Monte Carlo estimates introducing bias into the training objective.
We will take X to be RV taking values in X C Rd repre- senting our original dataset and X' to be a RV representing some
weighting distribu- tion. i.e., the "artificial distributions" we will utilise in our method. Similarly. we take E, E' to be expectations with respect to (W.r.t) X, X' respectively.
Throughout we take p to be the pdf of the RV, X, and pa to be a model therein. We let represent an unnormalised density (i.e. N-1 9 = p for some normalising constant N > 0.) We will write marginalisation//ondiiinnng for both true and model densities implicitly with p(xx) = fx p(x )dæ and p(=x|==>> being the conditional den- sity of Xx|X = - for example.
Parallel to this, some papers have looked to extend score matching to the latent variable setting, an area with much commonality to missing data (Vértes& & Sahani, 2016 Bao al. 2020 202T). Latent variable modelling differs in two crucial aspects from missing score matching. Firstly the components which are unobserved (the latent variables) re- main constant between samples, and secondly there is not necessarily notion of ground truth for the unobserved components in when data is corrupted. Additionally each of these works has limitations; Vértes & Sahani (2016) only applies to exponential families, [ao et al. (2020) re- quires a gradient unrolling step in its optimisation which is computationally expensive and can lead to errors in the op-
Now that we have introduced our notation we can move onto the key area of focus for our work, score matching.
3.2. Score Matching
First proposed by (Hyvärinen 2005). score matching aims to learn the gradient of the log-density (score). The advan- tage of this framework over full density approaches such as maximum likelihood estimation (MLE) is that we are not restricted to parametric models which integrate to 1. This allows us to be much more flexible in how we param-
2

---

Score Matching with Missing Data
eterise in turn making high dimensional distribution mod- elling more feasible. We now introduce the approach.
integration by parts) giving us that under certain regularity conditions on g, 8, and X,
Let X be RV over R with PDF We say that is the unnormalised density of X if - g(x) = p(z) where p is the PDF of X and N is the normalising constant of q- Define the score, of X to be
(T( = Ex(X), (20, + (xx
jeed
s a = V logp(x)) V log (xx).
E
Σ (x)(( Fri(6) C. jed
The aim of score matching is to learn s from collection of IID copies of X which we denote D = (x()) J Fol- lowing Hyvärinen (2005), we introduce generic parame- terised proposal score Sy for A E θ S IR" and aim to min- imise the Fisher Divergence between the true distribution and our proposal distribution which is given by
This can again be approximated via data using standard Monte Carlo approximation. Full details on the conditions required for this approach alongside the proof can be found in (Liu et al., 2022). Two other key extensions of score matching are sliced score matching (Song et al. [0200] and
denoising score matching (Vincent 2011). We introduce these extensions in Appendix D with our corresponding adaptations to missing data given in Appendix A.1 Now,
((0)= = (X) so {=)))
The key result from Hyvärinen (2005) which enables us to practically implement score matching is that under certain (fairly minimal) regularity conditions, which we provide in Appendix D.I
we give our missing data scenario.
3.3. Missing Data Scenario
we have
Instead of observing samples from X we assume that we observe samples from the corrupted version of the RV given by X. To define X we introduce mask RV M over {0, 1]d and then define X by
L(0)= = E 2v 80 (x) |80 (X) -(()) -
(1)
where here and throughout, we take C' to represent any con- stant which does not depend upon 0. Crucially, L(0) is now an expectation of observable random variables. Hence we can now approximate this with our data and take as
if M, Ø if M 0
X,
X,
=
where X, Ø represents that coordinate being miss- ing. We will be focussing on the missing completely at random scenario where M 1 X. However, we do pro- vide an extension to missing not random data in Ap- pendix A.1.4 We introduce the RV A on P([d) defined by A = fie E [d] M, 1} so that A gives the non-corrupted coordinates of and take À to be a sample of A. Crucially given samples from X, we also have samples from XA-
0 argmin
1
Σ
2Vx x so(x())+ + [1=
2V
0 n
TRUNCATED SCORE MATCHING
A limitation of standard score matching is that it requires lim,,,-oe p(x)= = 0 for all x, E R. Thus it cannot be used for many distributions with compact support if the density does not converge to zero at the (topological) boundary. Initial work to adapt score matching to truncated distribu- tions was presented in (Hyvärinen] 2007) for distributions on [0. oo) then further expanded in (Liu et al. 2022 Yu ett al. 2022) to general compact spaces X. For our com- pact space X' C Rd we use O.'' to denote the (topological) boundary. We now minimise some weighted version of the Fisher divergence whose weights go to zero at the bound- ary. Specifically let g X + IR be a function satisfying lim g(x),= for any x' E aX, j E [d]. Our objec- tive is then
Our aim is to adapt the score matching objective to estimate the full score s by a parameterised score Sy using samples from the corrupted data D = 1 {x(()
(x) 12 12
4. Marginal Score Matching
To motivate our approach we look at how we might use MLE in the case where the normalising constant and con- ditional normalising constants were calculable. For pa our parametric model of the density, we would choose 0 to be
PT(0)= = E (x) - x(x)
0 argmax log Po (x
Just as in classical score matching we obtain an equivalence (though this time via Green's theorem rather than simple
where Po is the associated corrupted data density when X ~ Po. As our data is missing completely at random
3

---

Score Matching with Missing Data
this is actually equivalent to maximising
If (e) also holds and there exists some 0" such that so- (X) s(X) a.s. Then if 0 is minimiser of Loa (0) we have that (x Np(X) a.s. for some constant N, i.e. the minimiser is the true unnormalised densiry.
log Pe:A, (x where (ax) == fxc
po(x)dd=_
X-
For notational simplicity we will thus reframe our problem as working with marginal samples {XX, Ji=1-
Through this result we have shown, much like with stan- dard score matching, that under certain regularity condi- tions our objective uniquely minimised by the true un- normalised density. We then approximate this objective by
4.1. Marginal Score Matching
Our approach is to directly alter the score matching objec- tive similarly. Just as densities have associated marginal densities so do scores have associated marginal scores.
LM.. (0)
x 8A. + 1| (X x, {1
Definition 4.1 (Marginal Score function). Let be a score function with s(x) > log ((x) for q an unnormalised PDF. Then the associated marginal score function is
and choose 0 argming LM.. (0).
Unfortunately this approach in its current state is practi- cally infeasible as the integrals involved in deriving the marginal scores for any non-trivial problem will be in- tractable. Hence, we must devise a way to estimate the marginal scores without having to compute the integrals. We tackle this issue in Section 4.2 but first provide similar result for the case of truncated score matching.
8x(xx) = Vax log
q( q(x)dx-xx
Rd-IA
This definition of marginal scores restricts 8 to a genuine score function. For this reason we will also want Se to always be a genuine score function or at least to have an anti-derivative. The simplest way to achieve this is to work with qo X + (0, oc) as our baseline and de- fine sa(x) = v log (zx). We will also take pe pa(x) (Ex go (x) qo (z) which we assume to be unknown.
4.1.1. TRUNCATED SCORE MATCHING
Truncated score matching can be adapted similarly to stan- dard score matching by simply having marginal weighting functions 9x X - [0, oo) for each subset À E supp(A) and taking the marginal truncated Fisher divergence to be
With this notion of a marginal score we can define our marginal Fisher divergence to be
FTM (0)= E (xx (Bx(XN) (xA(X))-8so(XA))
(0) E|| (XA) (Xx)(1)
(3)
where 8x:0 is defined analogously 8)- As with normal score matching can relate this objective to one involving terms of x. We first need the following assumptions.
using integration by parts gives the following equivalence
LTM (0)=0)
ga (XA)) 0no o
Assumption 4.2. For any 0 > 0, supp(A))
JEA
(a) pe is well defined, i.c. fx go (z)dd < oc;
+
20,gN (XA)
(5)
(b) (X o < oo;
JEA
(c) px (xx) is differentiable and exoo is twice differentiable;
FFrM (6)
(d) PA (xx) (0:0))000 - +0 as ||>||| | >0:;
Proof given in Appendix C.1.1 We then take LTM, as the Monte-Carlo estimate of LTM. We also construct simi- lar objectives from sliced and denoising score matching as well as a similar result for missing not at random data in Appendix A.I We now move to the task of estimating the marginal scores in these objectives.
(c) Px (Xx) px(Xx) almost surely (a.s.) for all À E SA) implies that (X) = p(X) a.s..
Assumption (a) ensures that our proposal unnormalised density always genuine unnormalised density. Assump- tions (b)-(d) are similar to the standard assumptions given for standard score matching. Assumption (e) is identifi- ability assumption which is required to be feasibly able to learn the true data distribution from our corrupted data.
4.2. Importance Weighting
Our first proposal is an importance weighting approach. Let p' be a density over Nd-I which we can both eval- uate and sample from then
Proposition 4.3. Given Assumptions 4.2 (d-(d) hold
(+) =E[2Vx X (X,)) + |SA.0 (
(4)
f (o(x))d = Bxx
(xx,xXX' p'(X",)
(6)
(0) C.
4

---

Score Matching with Missing Data
Algorithm 1 Marginal IW Score Matching
Our E ② be the minimiser of LTM. (0). 1 C R
with diam(@)- A then for sufficiently large n.r
Input: ( Jieja), qu. p', 2,0077 r E N.
Set 0 8o-
plog(dnrA/0) min (on m.
repeat
P (TM((OO))) 2811
< 8.
for i=1 to n do
Sample {X JRer from p' (.)
Use X (x -A; ]xE|r] to get Monte-Carlo esti- mates, AA (x(()) (X A, of the marginal scores by (7)
(x.().))
Note that r is the number of importance weighting samples for each data sample and therefore is something we can choose ourself. This means that with this approach we can achieve approximately v convergence rates. A downside however is that to achieve this we need be of order least which would lead to (n²²) computational cost. In practice we find relatively strong performance choosing small. Setting it at r 10 in our experiments.
end for
Use SA..:: (xX) to obtain LM/TMM.. (0) ) by
Compute V T VoLM//MM.r (0) and update the value of 0. until Maximum iteration reached.
This allows to define our marginal score estimate.
Remark 4.7. The error presented is measured with respect to our Marginal Fisher Divergence, rather than the full Fisher Divergence (which would be the preferred accuracy metric). Relating these two quantities requires connecting the fully observed distribution its marginals, a task that depends the specific form of the distribution. Inves- tigating the assumptions and conditions under which this connection can be made offers an interesting and valuable direction for future research.
Definition 4.4 (Marginal Score Estimate). For a given À E supp(A) xx E Xx, score model, 89, and E N we take our estimate of 80:2,r (xxx) to be
2a ( x()))
8x,1:0 = Vas log
(7)
p'(x (K)
where X - x are IID copies of X" p'.
4.3. Gradient First Approach
4.2.1. IW SAMPLE OBJECTIVE
A key limitation with an IW approach is that it will struggle in higher dimensional scenarios. Additionally the impor- tance weighting is embedded inside other functions which leads to the same nested expectation issue as the EM ap- proach of Uehara et al. (2020), causing bias in our esti- mator. As an alternative to this we build upon a variational approach initially discussed in the context of latent variable models in Vértes & Sahani (016)) Bao et al. (2020) 2021).
We can now plug these marginal score estimates into our sample objective for either normal or truncated score matching. We use M/TM to denote analogous definitions and results for both marginal and truncated marginal score matching. Let (x Ji=1 be our samples from XA. We then take our IW sample objective to be as (M/TM.((0) but with & a (x((1) replacing 8 (Xx The full ob- jective is given in Appendix [E.III We refer to this sample objective as (M/TTM..r()) and take estimate be
The core idea is to start with LM as before and then take gradients w.r.t. our parameters before then writing our ob- jective in terms of expectations over X-XAA;-- As we don't 't then need to take gradients of these expectations w.r.t. 6, we can estimate them with any black-box method we de- sire, opening the door for variational approximation to be used. This approach has been explored for exponential family distributions (Vértes & Sahani 2016) and for de- noising score matching (Bao et al. 22221) however we pro- vide the most general version of this result which can be ap- plied to any of the score matching methods and any model class. We first introduce the following key Lemma.
Ó = argmin LAA/TM,,.
Algorithm | gives our high level estimation algorithm.
Remark 4.5. Algorithm [ can directly be applied to both sliced and denoised score matching by replacing equation by equations (33 and (15) respectively.
4.2.2. FINITE SAMPLE BOUNDS
A benefit of truncated score matching is that it allows to work distributions with densities bounded below which enables to give finite sample bounds for the error of our estimated score w.r.t. our marginal objective. We briefly present these now with more detail given in Appendix A.2
Lemma 4.8. Fix X S [d], .x Xy. We have that for any function hg : X IR
( (I), X²x)
(8)
VE ho =E (I), X'11
(9)
Theorem 4.6. Suppose assumption 4.2 alongside assump- tions A./
+Cov (s0 (2),, X-x),
A.II
A.13 from the Appendix hold and let
5

---

Score Matching with Missing Data
where V represents the gradient w.r.l. either x or 0 and here E', Cov' are w.r.t X'=xXxx 2 po(-22>)-
the dependence on our current parameter estimate however we will omit this for brevities sake. The following propo- sition from Bao et al. 2020 shows how to train ø.
This results allows us to obtain our alternative objective.
Corollary 4.9. Let Lsa be defined as in
(4) We have that
Proposition 4.12 (Bao et al. 2020 For distributions p' PPPP let F(y' (p) and L (p' (p) be the Fisher and KL. divergences
between p' and p. We have that for any X E [d], x E Xx
(0) =I
4 (s)( 0,))(()))
(10)
jEA
P/ (X ((o)))
E' (Xx, X (XA,
X=
KL (I EE' log
40 g0 X"
+B
(ola)( | log (1's (1(20))))
where for any function ho R.d + R. x _ [d].
M (ho) (Xs,
x= X-xx)
+ Cov' (Vo log go (XA, X ho x_x)
and E', Cov' are w.r.t X X'A (o(() with E being w.r.l. XA p.
(xxx) where expectations are w.r.l. X-A ~ ((xx) and B is a constant not depending upon & (but will depend on 0.) In other words we can fit to the conditional density ( xx)) given only the unconditional unnormalised den- sity go (x)..) or full score so(xx..-)
Proofs for both results are given in Appendix
Crucially E,, Co can be estimated freely. This allows us to use variational inference to approximate Po( Po(x=x1>>) l(xx) and in turn the expectations and covariances in (10).
This allows us to train p'g (-|xx) to approximate the condi- tional density, po (x))) In our case we won't be learning this variational model for a fixed æ or even fixed observed coordinates A. Hence we take our objective to be one
Remark 4.10. We provide additional implementation de-
tails for computing this gradient estimate in Appendix A.5
We also discuss equivalences between this objective and our marginal IW objective in A.3
I'' ('S(X' (
JKL(Ø. 0) = E log
go (X(A
We explore estimation of E', Cov'in Section 4.3.2 but first we provide a similar result for truncated score matching.
Jr(Ø, = E |/T.C. log (6(X''XXA) 20 (X X"x
4.3.1. 1. TRUNCATED SCORE MATCHING We define a similar objective for truncated score matching.
with (LA, ~ P'((X_,Xx))(X)) We then take and Jr, Jr to be the Monte-Carle approximations with from the same distribution.
Corollary 4.11. With LTM defined as in we have that
samples
VOlTTM( =R 2
9x(XN)s (W +.)80())))
Remark 4.13. Jr has the advantage of not needing to know the normalising constant of d' = No r/ either.
JEA
Remark 4.14. As depends upon 6, we need update it each time we update A. In practice we find taking 10 gradient steps of o for each gradient step of 0 to work well.
+ 0,g0 (x),(()),
(11)
with V and E' defined as in Corollary 4.9
With this, we define Monte-Carlo estimate of ( with samples {(x())
VOLM TM
be the
Proof given in Appendix [.I.I Similar results for sliced and denoising score matching are given in Appendix A.
A. x ) )iL.)eeiixii where X(1 are our orig- inal corrupted data samples from p and (iA) are our variational samples from P's (.I. We can now state our full variational approach which is given in Algorithm
4.3.2. VARIATIONAL APPROXIMATION
We can now use variational approximation to estimate the expectations and covariances in Corollaries 4.9 & 4.1l
Specifically, let p=)(z=\|=>> be some generative condi-
Remark 4.15. Algorithm 2 can directly be applied to both sliced and denoised score matching by replacing equation (T)) by equations (14) and (16) respectively.
tional distribution dependent upon parameter é. We want to train/ P/g to approximate pg. We may write (0) to highlight
6

---

Score Matching with Missing Data
Algorithm 2 Marginal Variational Score Matching
Fisher Divergence
Marg- IW (Ours)
Input: {Xx, A; Jieja qu. p Go, do, L e N.E.N Set f 8g. Φ On-
Marg-Var (Ours)
Zeroed
repeat
2
EM
for I 1= 1 to L do
For i E [x] sample ( from p/s x.)
A.
1-
Use {(XA) Ag ".0,.xx) x"x) to get Monte-Carlo ap-
proximates of JKL/F(0, (0,0) given by JKL,
Compute VoJKL/F(Ø0 0) and update Φ.
200
400
600
800
1000
Sample Size
end for
For E [nx] sample (X())) hRE from (-
Figure 1: Average Fisher Divergence for Gaussian score estimates alongside 95% C.I.s Lower is better.
Use {(XX) A; x (aNEEMMI) to get our Monte- Carlo estimate, VoLM/TM(()) using equation (TT)
Use this gradient estimate update 0. until Maximum iterations reached.
not converge, a phenomenon we discuss more in Appendix In Appendix [I..I we present the average mean and precision estimation error for this experiment. In Appendix B.1.2 present the untruncated results and illustrate how the naive marginalisation poorly models strong relationship between dimensions and 10.
D.2
5. Results
Here we go through simulated results comparing our IW approach (Marg-IW) in Algorithm T and our variational ap- proach (Marg-Var) in Algorithm 2 to the EM approach of Uehara al. (2020). We also compare to naive marginal- isation approach involving zeroing out the missing dimen- sions and only taking the observed output dimensions of the score, which we call Zeroed Score Matching. This ap- proach is the natural adaptation of MissDiff from Ouyang et al. (2023) away from NN to explicitly parameterised models. We describe Zeroed Score Matching and its re- lation to MissDiff Appendix D.2] In our experiments, highlight a unique strength of our methods by applying them to explicitly parameterised score models. We could however, equally apply them to more complex, noninter- pretable models such as NNs. More implementation details can be found in Appendix E.3
5. 1.2. NON-GAUSSIAN MODEL
For this experiment we tested our parameter estimation for a an ICA inspired unnormalisable model of the form
p(æ) exp
Σ hd
,
Here we parameterise our model identically with the aim of estimating 6*. We vary the dimension of X and plot the estimation error with sample size of 1,000 and each coordinate missing independently with probability 0.5, The results are presented in Figure 2
Theta Error
5.1. Parameter Estimation
3-
1. TRUNCATED GAUSSIAN MODEL
In this experiment 10-dim normal distribution is set up with fixed mean and random covariance before being trun- cated on the first 3 dimensions. 1000 samples are taken and corrupted independently on each coordinate with prob- ability 0.2. For each of our methods a Gaussian score fit and the Fisher divergence between this score and the truth computed. This repeated 200 times with the mean Fisher divergence alongside 95% C.I.s then presented in figure More details in Appendix E.3. Marg-IW and EM per- form best with Marg- Var approaching asymptotically. We see the effect of Zeroed's ss naive marginalisation as it does
Marg-IW (Ours)
Marg-Var (Ours)
Zeroed
EM
0
10
20
30
40
50
Dimension
Figure 2: Average Fisher Divergence for Gaussian score estimates alongside 95% C.I.s. Lower better.
Our variational method (Marg-Var consistently yields the lowest error. Moreover, as the dimensionality increases, the performance gap between Marg- Var and the other meth- ods widens. This supports the notion that our approach
All code and data for the experiments presented can also be found at htttps://qithut com/ joshg/jshiivve vens/ ScoreMat.ch www.:aatthiiiiiiiiicm whwiswitiiisiiiiiiiiii
7

---

Score Matching with Missing Data
more accurately able to capture complex marginalisations than the competing approaches which fail as the dimen- sion grows. We note that all other methods perform com- parably with the performance of EM and Marg-IW being indistinguishable, a pattern we observe throughout our ex- periments. This similarity is unsurprising both approaches use self normalised importance weighting to approximate conditional expectations with respect to our current score estimate while being broadly motivated by fitting the marginal scores. Nevertheless, the precise mechanism for this similarity remains unclear and warrants further explo- ration. Additional experiments exploring the effect of sam- ple size and missingness probability on estimation accuracy are given in appendix B.1.3
1.0
0.9
AUC
0.8
Marg-IW (Ours)
0.7
Marg-Var (Ours)
Zeroed
0.6
EM
0.2
0.4
0.6
0.8
Missingness Probability
Figure 3: Mean AUC of star graph edge detection with varying missingness alongside 95% C.I.s. Higher is better.
5.2.2. UNSTRUCTURED DENSE GRAPHICAL MODEL
5.2. Gaussian Graphical Model Estimation
Here we create GGM by making each edge occur inde- pendently with probability 0.5. The rest of the experiment was constructed as before. Results are given in Figure 4 Again we can see that our variational approach performs
Gaussian graphical models (GGM) are a popular way of modelling dependence between dimensions of data. Let assume that the underlying data follows a Gaussian distri- bution with mean f E Rd and precision P E Rddd In this setting, Bayesian network (BN) can represent the depen- dencies between the dimensions of X with the (undirected) edges of the BN exactly being the non-zero off-diagonal entries of the precision, P. Hence estimating the precision matrix P gives the BN. Score matching has been shown to be an effective way of achieving this with L.1 -eereuullariattinn on the off-diagonal of P to push terms to 0 (Lin et al.) 2016 Yu et al. 2018). Decreasing the level of Llreguularisttion then gives range of classifiers with increasing True and False positive rates (TPR/FPR) as the level of regularisa- tion decreases. Score matching can also be applied to trun- cated GGMs where aim to learn the original BN but only observe the samples inside some truncated region.
0.85
0.80
AUC
0.75
Marg-IW (Ours)
0.70
Marg-Var (Ours)
0.65
Zeroed
EM
0.60
0.2
0.4
0.6
0.8
Missingness Probability
Figure Mean AUC of dense graph edge detection with varying missingness alongside 95% C.I.s. Higher is better.
We apply our methods to learn GGMs and truncated GGMs with missing data as well. We use varying levels of Ll reg- ularisation on our objective via proximal stochastic gradi- ent descent in our optimisation [Beck] 2017).
best though not as clearly as in the previous example. We believe this to be because for more unstructured problems, naive marginalisation performs moderately well.
5.2.3. INCREASING NUMBER OF STARS
5.2.1. 1. STAR SHAPED TRUNCATED GRAPHICAL MODEL Here we create star shaped GGM in which one node has a high probability of being connected with each other node independently and all other connections have probability We truncate the data along a random hyperplane such that 20% of the distribution lies outside of the truncation boundary. Each coordinate is then MCAR independently with the same probability. We run multiple experiments with this probability ranging from 0.2 to 0.9 and present the results in figure 3 As we can see here Marg-Var performs best with all other approaches performing comparably. For illustrative purposes, we plot individual ROC curves from this experiment Appendix B.2.3
To explore this further, we construct and experiment where we vary the number of star centres (high degree nodes) while keeping the edge density constant. We present the re- sults in Figure As we increase the number of star centres, Marg- Var no longer noticeably outperforms the other ap- proaches. This is because as the number of stars increases, (i.e. the structure of the graph decreases) naive marginali- sation is a better approximation. This is illustrated on the marginal precisions themselves in Appendix B.2.
5.2.4. S&P 100
Here we took closing price data over 5 years for the 100 stocks in the S&P 100 with each stock being dimension
8

---

Score Matching with Missing Data
0.85
Marg-IW (Ours)
Marg-Var (Ours)
0.9
Zeroed
AUC
0.80
EM
AUC
0.8
Marg-IW (Ours)
0.7
Marg-Var (Ours)
0.75
Zeroed
0.6
0.70
0.2
0.4
0.6
0.8
Missing Probability
Number of Star Centres
Figure Mean AUC of various methods when compared to non-corrupted score matching with 95% confidence in- tervals on genetic yeast data. Higher is better.
Figure 5: Mean AUC with 95% C.I.s for edge detection as number of star centres in graph varies. Higher is better.
and each day being sample. Gaussian graphical models with various levels of connectivity were then constructed using standard score matching on the fully observed data. The data was then artificially corrupted and each missing score matching approach applied. The AUC was then cal- culated for each method taking the GGM from fully ob- served score matching as the ground truth. More details given in appendix E.3.3 The results shown in figure
6. Conclusion
conclude, score matching is a versatile method whose applications the heart of modern machine learning prob- lems. In this work we have tackled the problem of adapt- ing score matching partially missing data. We have pre- sented two separate but related approaches to this method, one using importance weighting and another using varia- tional approximation. We have also provided extensions of these methods to truncated score matching, sliced and de- noising score matching. For truncated score matching with our IW approach we have provided finite sample bounds on the accuracy of the estimated score in terms of the marginal truncated Fisher divergence.
0.9
AUC
0.8
Marg-IW (Ours)
Marg-Va Var (Ours)
0.7
Zeroed
EM
We have provided several simulated and real world experi- ments demonstrating methods' efficacy for both param- eter estimation and downstream GGM edge detection. We have shown the benefits and drawbacks of each approach with IW performing best in lower dimensional settings with less data and the variational approach performing best in more complicated higher dimensional settings.
0.4
0.6
0.8
Missing Probability
Figure 6: Mean AUC of various methods when compared to non-corrupted score matching with 95% confidence in- tervals on stocks in S&P 100. Higher is better.
There is. however still much work to be done in this area. From theoretical perspective, while have finite sam- ple bound the error of our loss, marginal nature of the loss makes it unclear exactly how this translates to param- eter or general score model accuracy, leaving room for fur- ther theoretical exploration. From an implementation per- spective, variational inference the presence of missing data requires accounting for the randomness of "latent" and "observed" variables. The standard variational inference technique be further refined to accommodate this set- ting. Finally, since our method is compatible with denoised score matching, it can naturally be extended to diffusion- based model. This paves the way for future work on ap- plying our approach to generative modelling with diffusion processes in the presence of missing data.
As we can see Marg- Var clearly out performs all the other approaches which appear to perform equivalently.
5.2.5. YEAST DATA
Here data first introduced in Brem & Kruglyak (005)) used consisting of readings expression for 7086 genes/ORRs across 262 yeast segregants. Each gene repre- sents dimension with each segregant representing sam- ple. We subset the data to take the 106 genes present in at least 95% of the samples with the aim of learning the rela- tionship between them. The same approach as the previous section applied with the results shown in f Again Marg- Var clearly outperforms the other approaches which all perform comparably.

---

Score Matching with Missing Data
Acknowledgements
Daras, G., Shah, K., Dagan, Y., Gollakota, A., Dimakis, A., and Klivans, A. Ambient diffusion: Learning clean distributions from corrupted data. In Thirty-seventh conference on neural information processing systems, 2023. WRLLittps::/opeerevve net forum??
Josh Givens was supported by a PhD studentship from the EPSRC Centre for Doctoral Training in Computational Statistics and Data Science (COMPASS).
id Id=wBBBLyyKBY
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Fattahi, S. and Sojoudi, S.
Graphical lasso and thresholding: Equivalence and closed-form solutions. Journal of Machine Learning Research, 20(10):1-44, 2019. URL http:/ http:://mmm.rrrrom http://all...oopppppe/////0/0
17-501 .html
References
Huang. C.-W,, Lim, J. H., and Courville, A. C. A varia- tional perspective on diffusion-based generative models and score matching. Advances in Neural Information Processing Systems, 34:22863-22876 2021.
Bao, F. LI, C., Xu, K., Su, H., Zhu, J., and Zhang, B. Bi-level score matching for learning energy-based latent variable models. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.). Advances in neural information processing systems, volume 33, pp. 181 18-11222. Curran Associates, Inc., neurips 2020. URL http://eeeeeapppe///////////
Hyvärinen, A. Estimation of non-normalized statistical
models by score matching. Journal of Machine Learning Research, 6(24):695-709, 2005.
Cc cc/paper_fiiieceeeeeeee fi les
d25a (0053333322811300000711111111113 Pape:
Hyvärinen, A. Some extensions of score match- ing. Computational Statistics & Data Analy- sis, 51(5):2499-2512 2007. ISSN 0167-9473,
pdf
Bao, F., Xu, K., Li, C., Hong, L., Zhu, J., and Zhang, B. Variational (gradient) estimate of the score function in energy-based latent variable models. In Meila, M. and Zhang, T. (eds.), Proceedings of the 38th interna- tional conference on machine learning, volume 139 of Proceedings of machine learning research, pp. 651-661. PMLR, July 2021. URL https: mlr. wwwpsppseeee.......... https//v19//aaae http::::s///pprreeeinn.g html
doi: http://bb..rr////.00...sst2222.. URL https///ww..sceeneeeiieeee cclenne/articee///////8777777770
Li, S. C.-X., Jiang, B., and Marlin, B. MisGAN: Learning from Incomplete Data with Generative Adversarial Net- works. In International Conference on Learning Rep- resentations, 2019a. URL TRRL:atps:///poccc btttps://pperreieee. net/forun d= DV33cKm
Beck, A. The proximal gradient method. In First- order methods in optimization, pp. 269-329. So- ciety for Industrial and Applied Mathematics, URL
Li, Z., Chen, Y., and Sommer, F. T. Learning energy-based models high-dimensional spaces with multi-scale de- noising score matching, 2019b. arXiv: 1910.07762 [stat.ML].
2017. doi: http://wwww...aaaaooorooooaaaa//
https:
ch10
Brem, R. B. and Kruglyak, L. The landscape of genetic complexity across 5,700 gene expression traits in yeast. Proceedings of the National Academy of Sciences of the United States of America, 102(5):1572-1577) February 2005. ISSN 0027-8424 1091 16490. doi: 10.1073/pnas. 0408709102. Place: United States.
Li, Z., Chen, Y., and Sommer, F. T. Learning energy- based models in high-dimensions spaces with multi- scale denoising- score matching. Entropy. An Interna- tional and Interdisciplinary Journal of Entropy and In- formation Studies, 25(10), 2023. ISSN 1099-4300. doi:
10.3390/e251011667 URL https://wwwt mdpii.
com/ 00w11/99990///11111 tex tex.pubmedid: 37895489.
Number: 1367
Burda, Y., Grosse, R. B., and Salakhutdinov, R. Impor- tance weighted autoencoders. In Bengio, Y. and LeCun, Y. (eds.), 4th international conference on learning rep- resentations, ICLR 2016, san juan, puerto rico, may 2-4, 2016, conference track proceedings, 2016. URL 7/arxiv. wwwii.orr//bb.... aaaaaba/11599 00519 tixbbibso dblp computer science bibliography, https:/dblpporg tex timeetmst Thu, 25 Jul 2019 14:25:37 +0200.
Lin, L., Drton, M., and Shojaie, A. Estimation of high- dimensional graphical models using regularized score matching. Electronic Journal of Statistics, 10(1):806 854, 2016. doi: 10. 10.1214/16-EJS1122 doi 11.4..21141661B1
URL Publisher: In- stitute of Mathematical Statistics and Bernoulli Society.
10

---

Score Matching with Missing Data
Liu, Q. and Wang. D. Learning Deep Energy Mod- els: Contrastive Divergence vs. Amortized MLE, July 2017. http::/aaxiii..a/// a ab777 00797 uIXiv1100.007997 [es, stat].
Swersky, K., Ranzato, M., Buchman, D., Freitas, N. D., and Marlin, B. M. On autoencoders and score matching for energy based models. In Proceedings of the 28th in- ternational conference on machine learning (ICML-11), pp. 1201-1208, 2011.
Liu, S., Kanamori, T., and Williams, D. J. Estimat- ing density models with truncation boundaries using score matching. Journal of Machine Learning Re- search, 23(186):1-38, 2022. URL http:///mir. ore/ htt/pappers/v23/221 0218.
Tashrro, Y., Song, J., Song, Y., and Ermon, S. CSDI: Con- ditional score-based diffusion models for probabilistic time series imputation. In Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), Advances in Neural Information Processing Systems, 2021.
Ouyang, Y., Xie, L., Li, C., and Cheng, G. MissDiff: Train- ing diffusion models on tabular data with missing values. In ICML 2023 workshop on structured probabilistic in- ference & generative modeling, 2023. URL https: 7/openzeview ftree/foouueee wmmmmmmmddsseeedid
Uchara, M., Matsuda, T., and Kim, J. K. Imputation timators for unnormalized models with missing data. In Chiappa, S. and Calandra, R. (eds.), Proceedings of the twenty third international conference on artifi- cial intelligence and statistics, volume 108 of Proceed- ings of machine learning research, pp. 831-841. PMLR, ngs. mlr.
Richardson, T. W., Wu, W., Lin, L., Xu, B., and Bernal, E. A. MCFlow: Monte Carlo Flow Models for Data Im- putation. In 2020 IEEE/CVF Conference Computer Vision and Pattern Recognition (CVPR), pp. 14193- 14202, June 2020. doi: 10.1009/VPP44260.22222 01421. ISSN: 2575-7075
August 2020. URL Uttpttts:://proeeee press /vess/vv08// wwt.DD0/uuhaaaa... hhtml
Vincent, P. A connection between score matching and de- noising autoencoders. Neural Computation, 2367::1661 1674, 2011. doi: 10.
Salimans, T. and Ho, J. Should EBMs model the en- ergy or the score? In Energy based models workshop ICLR 2021, 2021. URL httpp::/ / net/ htt/://orummdieeeeee
Vértes, E. and Sahani, M. Learning doubly intractable la- tent variable models via score matching. InSymposium on advances in approximate Bayesian inference (AABI), 2016.
TF2 jRNb
Sasaki, H., Hyvärinen, A., and Sugiyama, M. Clustering via mode seeking by direct estimation of the gradient of log-density. In Calders, T., Esposito, F., Hüllermeier, E., and Meo, R. (eds.), Machine learning and knowl- edge discovery in databases, pp. 19-34, Berlin, Heidel- berg, 2014. Springer Berlin Heidelberg, ISBN 978-3- 662-44845-8.
Yang, E. and Lozano, A. C. Robust gaussian graphi- cal modeling with the trimmed graphical lasso. In Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., and Garnett, R. (eds.), Advances in neural information processing systems, volume 28. Curran Associates, Inc.,
2015. URL http:://pprocceeeings.sssst neurss.
Song, Y. and Ermon, S. Generative Modeling by Esti- mating Gradients of the Data Distribution. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
cc/ /taeppprfff1ee//ppppe////11//11 Ttbbbicc22990aal13555555888
Paper
pdf
Yoon, J., Jordon, J., and Schaar, M. Gain: Missing data imputation using generative adversarial nets. In Interna- tional conference on machine learning, pp. 5689-5698, 2018. tex tex.organization: PMLR.
Song, Y., Garg, S., Shi, J., and Ermon, S. Sliced score matching: A scalable approach to density and score es- timation. In Adams, R. P. and Gogate, V. (eds.), Pro- ceedings of the 35th Uncertainty in Artificial intelli- gence Conference, volume 115 of Proceedings of Ma- chine Learning Research, pp. 574-584. PMLR, jul 2020.
Yu, S., Drton, M., and Shojaie, A. Graphical models for non-negative data using generalized score matching. In Storkey, A. and Perez-Cruz, F. (eds.). Proceedings of the twenty-first international conference on artificial intelligence and statistics, volume 84 of Proceedings of machine learning research, pp. 1781-1790. PMLR,
Song, Y., Durkan, C., Murray, I., and Ermon, S. Maximum likelihood training of score-based diffusion models. Ad- vances in Neural Information Processing Systems, 34: 1428, 2021a.
April 2018 URL pttps::/ppoceedng.ss mmr.. press / v84/ ww4/v//4lll00 8b. html
Song, Y., Sohl-Dickstein, J.. Kingma, D. P., Kumar, A., Ermon, S., and Poole, B. Score-based generative model- ing through stochastic differential equations. In Interna- tional Conference on Learning Representations, 2021b.
Yu, S., Drton, M., and Shojaie, A. Generalized score matching for general domains. Information and Infer- ence: A Journal of the IMA, 11(2):739-780, 2022.

---

Score Matching with Missing Data
A. Additional Theoretical Results
Here we present some interesting results which we feel help further build up the landscape of our method but were unable to fit within the main body of the paper.
A.1. Additional Methods
Firstly some additional adaptations of score matching. Most of these are relatively immediate adaptations following our framework for missing score matching although there are some important aspects and caveats which make them worth officially documenting. Missing and sliced score matching are introduced in detail in Appendix D
A. TRUNCATED SCORE MATCIIING
We have already presented truncated score matching in the paper however we present in more details alongside its assumptions here.
Assumption A.1. For any E supp(A), e e:
X is connected, open and Lipschitz;
Px Ix, 1x.0 E H (Xx):
Px: Ix are continuously differentiable and go:a is twice continuously differentiable;
for any x's E 0xX,, and j E x we have
lim
(x 0.
where v((x)) is the normal vector to the boundary x
This now leads us to our proposition on the validity of our population objective.
Proposition A.2. Suppose that assumptions 4.2 & A. hold. Then we have
(T) = (x) -xX (TM
(12)
where C' is does not depend upon 0. As direct result for 0 minimiser of LTa( we huve that 8g(X) = s(X) a.s.
Proof. Proof given in Appendix C.I.I
1.1.2. MISSING SLICED SCORE MATCHING
For readers who are unfamiliar with sliced score matching we provide brief introduction Appendix D.3 For sliced
score matching the only adaptations we need make is use our marginal scores and now alter our projection vectors be over the appropriate subspace. Thus our objective becomes
(f) = Fl2 XXx (V $A (Xx))) T FSM (0) C
VA + VA (A,((A))
where for any À E supp(A), VA is RV on RI satisfying C positive definite and B|[[V[||
< co.
To write this and it's gradient in terms of the full score, Sy we can again use Lemma 4.8
This gives the following results Proposition A.3.
(0) -2E E' Vxa (V. (0(A vN) + (xXXxx E'((V ET X'
(13)
VoLsM (0) -2E VA (xx (V (o() Vi)) + VnA (VT so())) so X2 (V so())
(14)
12

---

Score Matching with Missing Data
where for any function hg R.a R, (ha)a =E Voho XA, + Cov' (Vo log qo (Xa, (xx and E', Cov' are w.E.t. X AA ~ po vo((|Xx) XA)
with E being w.r.t. XA ~ p.
Proof. We first have that
(0) = 2
V
NA
(X,,
VA + X=x))))))
2
Viy V (X XLMMVV (VAE
JEA
V V x, X C... Tv VA
JEN
+ (V EXA,
21E' Vxa (V1 60 va)] E'((V (V
+ E' E'VV 3 (X x N'aaa
2E
Xxx V s0 (XA.X" x) VN) +EV((Tsso X E'((T sa(XA, x
The second results directly from applying Lemma 4.8 again alongside the chain rule.
1.1.3. MISSING DENOISED SCORE MATCHING
As with sliced score matching the adaptation is relatively immediate however we do first need to make some further restrictions on our noising process. Specifically we require that for any t E [0. and j. E [d] we have X(t), T x(t), X (0)5-
In most practical implementations each coordinate independently noised therefore satisfying this condition. We require this allow us to easily write the marginal transition kernel for any supp(A) given by (x)) We then make population objective
(DM = E (t)) ((x))))))(((((())))))))+(+ (x() (x)))((x(((xxxx()))
We can again write this in terms of sy as do in the following proposition Proposition A.4.
}
(0)
v(t)
Σ
E'
(X,n + Vxx.ep)ogxx(XX(()))xx000
(15)
JEN
VOLDM(()
v(t) Σ
E' (x((Xa(),, (F10)s((Xx1())xxx2))))
(16)
jea
+ (Xor(so(Xx,())XX (X,
Vxx()) (xxx)))))
13

---

Score Matching with Missing Data
where for uny function hg: R.d R, (ha)a) =E' (X X'-NI++ Cov' (Vo log go (x) and E', Cov' are w.F. X (X)) po())Xa) with E being w.r.l. (1( ~ Pt-
Proof. Using Lemma 4.8
we have that
|
(0)
R
v(()
Σ
AA (X)()(()) + (x (0))
JEA
= E
v(t)
Σ
E'
s0 xxx + (x))
JEA
A second application of the lemma then gives,
VOLDM = E (tt)
Σ
E' X-,, t) (xa),
JEA
+
+ Pope
A. 1.4. MISSING NOT AT RANDOM DATA
So far we have assumed that our data is missing completely at random so that X^/A = ~~ Xx. In other words, we could treat corrupted samples as though they were simply marginal samples and still perform valid inference. Often however, such an assumption is unrealistic and the probability parts sample begin missing depends upon the sample itself. Generally this split into two cases. Missing at Random (MAR) and Missing not at Random (MNAR). MAR data occurs when the probability of coordinate being missing depends only upon other coordinates of the sample. This means that M, X5 X-j. In MNAR data we allow Mj to depend upon X, as well meaning that an observations value determines its own probability of being missing.
Here we will focus in the MNAR scenario and treat the MAR scenario as special case of this. The core idea of this approach will be to work with "joint" score rather than a marginal score. Before we do this we need to set-up our MNAR case. Specifically for E supp(A) define the event
Ex = (X
and define (X) X)) Throughout will assume each fpx to be known. This is often unrealistic assumption however this allows us the flexibility of having a method which is independent of how the 4a are learned.
To work with this MNAR data we need to define some adaptations of densities and score functions.
Definition A.5. X with PDF and event E define p(:; E)) to be the "joint" density satisfying
J.
(la:E)de P([x E B)UE)
B
for all B Bx-
From this and with our particular events redefine the missing score as,
sx (xxx))= = V œx log (zx: Ex)
(17)
(f
Vay log
((a:E))a E)dda->>
(18)
( ((aa))exaa)ddd
= Vay log
14

---

Score Matching with Missing Data
Remark A.6. this missing score is not the same as the marginal score. We slightly abuse notation here using the same notation did for the marginal score. This is however reasonable as for the MCAR case the marginal score and the missing score are identical.
With this newly defined score, we can proceed similarly to the MCAR case and use the objective (m) defined in 4 or ( but with our new defined score. We now show a provide a similar justification for this approach as the MCAR case but first need to introduce an additional assumption.
Assumption A.7. For each À supp(A). ((EX(Xx) 0 a.s..
Remark A.8. We do not require every missingness pattern have positive probability just that if missingness pattern does have positive probability, it has positive probability for every possible underlying sample.
This then leads us to our desired result.
Proposition A.9. Suppose with are in our MNAR set-up and assume that assumptions 4.2 & A.7 hold and that there exists 8t with 80- (X) = 8x(X) a.s.. Then if 0 is minimiser of (0a(0) where the missing scores are defined by then sg(X)= s(X) a.s.
The proof for this is similar to the MCAR case and is given in Appendix C.1.2
Now we have our objective we need to see how we can derive (xxx) Again we can do this similarly to the MCAR case. Let Qo be our estimate of the unnormalised density then
(x:0(x>> Va log PAio (x:: (xx:
Vax log /
pa(:; Ex)dæ-A
Vay log
X-1
Va, log
90(2)(2(Z)dzz
X-
(0 X' (EX)
Vax log DEp
v(X_x))
4o 40 x (k) )o (Ix x (1k)
1
V log
p' p(
k=1
As a result we can approximate our objective analogously to our approach for MCAR data.
A.2. Finite Sample Bounds for Truncated Importance Weighted Score Matching To be able to prove finite sample bound results we first need to present some key definitions.
Definition A.10 (Approximate Truncated Marginal Score Matching Objective). For n,r E N, take our sample objective to be
(TM 11
(A + 2V xx T T x
gn.
with $x,r;0 (xxx) being estimated marginal score from Definition 4.4
Additionally define
(0(x:)
(=)) f1,(2,,)
Vacdot Vacdot(Vaa0001 (z)o(x)) p' y2(x__)
fo, pl A1
f2,(2,
p(xx-A)
We now set-up the following assumptions Assumption A.11. There exists a > s.t. for all E X, X E supp(A), 1,,
15

---

Score Matching with Missing Data
(ax, 0)1i, (x(xx), | ( a.
1 < fo,((1200)
Remark A.12. is this assumptions which restrict us from obtaining a similar result the non-truncated case as it is unrealistic to have both (ax) and p(ax)) -+ ||æx|| oo.
Assumption A.13. For each À supp(A), e {0, 2, we have that for any 8, 0' E G:
[1xx - - fo.x Mc()))(000)
where M& _x),, Mx X",, are sub-Gaussian with parameters dT oi. respectively for all æ E &
Remark A.14. This assumption is immediately satisfied if @ is compact and ft, 0) are pointwise Lipschitz w.r.t 0.
Hence this assumption slightly weaker than a uniformly lipschitz assumption
We can state our theorem
Theorem A.15. Assume that assumptions
4.2
A.1
A.11
hold and let Our E θ C RP be the minimisers of A.13
LTM.. (0). Ifo I2? then for sufficiently large n,
plog dnr diam (0)/8)
plog(n diam diamm())/))) (9)/5)
attrr
(C+
log(n/0)
P FrM((nn)) 2 B
32
+ Ba
6.
nr
where 31,3. 32 depend upon a, Ba depends upon a, xampp())) and C depends upon a,
[T[Ma X=)0
Proof. The proof for this alongside multiple intermediary results can be found in C.2
Here we have shown convergence of our sample/approximate objective to the population objective. This combined with proposition A.2 which states that our population objective minimised by the true score suggests that our approach does give valid inference learning the score. A key limitation of our result is that to obtain convergence, we require r-000. Furthermore, obtain log(n)/n rate convergence we need r be of the same order as n. As the computational complexity of our algorithm in (nr)), this suggests that to obtain desired convergence to the population objective will have (n²) computational complexity.
Remark A.16. Our dependency on our Lipschitz constants only enters into the term with the associate sub-Gaussian parameters entering only into the a.
Remark A.17. Dependence upon g simply requires g and Vg bounded. This is achieved compact X by g(x)) mix d(æxxe o²) and on non-compact space by g(x) min d(x, x') A1
A.3. Relationship between IW and Variational objectives
Despite being derived quite differently from the marginal score matching objective. We show below that the two objectives are actually identical in some cases. Specifically, when the IW density p' doesn't depend upon the observed data A we can treat the importance weighted approach as an importance weighting approximation of the gradient estimate in (10) We demonstrate this through the two results below
16

---

Score Matching with Missing Data
Lemma A.18. For some density p' which generates samples (X-X beer let
wk (x (
((( )
Wk wk Σ Wkk
k'=l
1
(xxx = V xa log
W/k
R=1
Eiv (ggo(xx())) ==
Wh (x>, X (KK)
(f =
Tkga xx x )
X11)) x((k)
Trge x x
(+
wu onf (2x (ax. x
k=1
k=l
Then
0.: (xx) Eiv ( xx x²1
(19)
100 x" w = Biwl ge (a), X ,)1 + (s X- w: go x²x)).
(20)
where V represents the gradient w.r.t. x 0. In other words, can take importance weights first then gradients (L.HS) or gradients and then importance weights (RHS).
Proof. Proof given in Section C.4
Corollary A.19. We have that
Xo = 28, (2 (x), x" Ne) Bie Veso(zx, + Coviu (xx, s0(2))
2(Ei Vosa X- Coviw (o X_. se(xx, X(x
(21)
+ 2(E. Vod, (1) x Cov (Vol
Proof. Proof given in Section C.4
A.4. Exploring the Marginal Fisher Divergence for Normal Distributions While intuitively the Fisher divergences of the marginal distributions should act effective proxies for the Fisher diver- gence for the fully observed distributions, we would like to be able to examine the relationship between the two more explicitly. We do know that marginal Fisher divergences will be zero when then fully observed distributions equivalent however here we give a more detailed examination in the case of normal distributions.
Suppose that X ~ N(u, P-1)
p(x) - expl
-x) -))TTT( 14) +C
with C' a constant w.r.t. æ. We then have that
s(x) - P(æ 11)
If suppose that our unnormalised density/score model is of the form
(o(x) exp { (x (o0 Toox mo0
80 Po(x -40)
17

---

Score Matching with Missing Data
Then with the marginal Fisher taken to be
(0) = EA.x, [XA
where here for each supp(A), 8,, S0: are the true marginal scores. Using properties of the normal distribution and the Schur complement we know that the precision of X is given by
{(P))) xa))
PA. - Pa.-2 -
AAJ
Plugging this in we get
(0) =E PPAA- PAAAA)A (PA-AP-A-aP-A Ph PAIA-APo -
P.-AA-PP--AAAPP--A - (PosA, Po:r, Poj A; A,- Po: A. A.)
This shows why naive marginalisation by zeroing out missing coordinates of our score would not work as in this case the Fisher divergence would be given by
(0) -E ((PAA - - P Po:A,A)XA ((PA.A PA. PA--PIA-AP-AAAA PAAA)AA PossAMo:))
2
P-xx for all A E supp(A) meaning it will not give us the which encourages Po:r, to be close to PA. PA.-AP - true density.
A.5. Variational Pseudo-loss
When using (TO) it is helpful be able view it the gradient of some pseudo-loss allowing it to plug into standard ML framework where calculate the loss, take the gradient w.r.t. parameter using auto-differentiation. and update our parameter estimate based on this. The below result show how we can do this by creating loss with certain instances of our parameter detached from the computational graph.
Proposition A.20. Let
J(0,0' xx x_x)== = 2E ss (xx, X"x (E X + Cov' (log 20 (ax,, X-x); so (xx,, X-x)
2(E s0 ( X' x² + Cov' (log 1o (ax, X" A So (xA X x))
+ 2(E'(0,,8 (0,88(20) (x X" + Cov' (loggo(0x, X" 0,00 x)))
where E,, Cov' are w.r.t. X 1Xx xxi 0¹ Then
② 00
L 2x))
J(0 00, xx)) 0=0"
Proof. This just follows directly from the exchangeability expectations and gradients (when the gradient w.r.t some- thing independent of the expectation distribution.)
this loss (by replacing all instances of H' with and then detaching them from the computation graph) to treat our problem as a standard gradient descent problem.
Hence
Note that while we can treat this like a loss for our optimisation, our intent is not actually to minimise it. The estimated form of the loss given in the proof of Corollary 4.9 which is given in C.3 but we state again explicitly here for convenience.
(l = E
Σ jEA
E' [sa
2E [s (XA, X'
(XA,
20' (XA,
18

---

Score Matching with Missing Data
Precision Distance (Froben
Marg-IW (Ours)
Mean Distance (L2)
2.0
2.0
Marg-IW (Ours)
Marg-Var (Ours)
Marg-Var (Ours)
1.5
Zeroed
1.5
Zeroed
EM
EM
1.0
0.5
0.5
0.0
0.0
200
400
600
800
1000
200
400
600
800
1000
Sample Size
Sample Size
(b) Average error in estimation of the Precision (Frobenius Norm).
(a) Average in estimation of the Mean (L2 norm).
Figure 8: Average parameter estimation enror for truncated Gaussian score estimates alongside 95% confidence intervals under various methods.
B. Additional Experimental Results
Here we present some additional experimental results not in the main body of the paper.
B.1. Parameter Estimation
B.1.1. TRUNCATED GAUSSIAN MODEL
Here we present the accompanying mean and precision error results for Gaussian model estimation experiment presented in Section 5.1.1 These results are presented in Figure
B. 1.2. UNTRUNCATED GAUSSIAN MODEL.
Here we present the untruncated version of the experiment presented the main paper. Details of the distribution are the same as presented in Appendix
but without the truncation.
Fisher Divergence
2.5
Marg-IW (Ours)
Marg-Var (Ours)
2.0
Zeroed
1.5
EM
1.0
0.5
0.0
200
400
600
800
1000
Sample Size
Figure Average Fisher Divergence with 95% C.I.s for various approaches. Lower is better.
As we can see obtains similar results here as in the truncated case.
We also illustrate what the true covariance and precision matrix look like for this example alongside the naive marginali- sation in order to highlight where Zeroed Score Matching goes wrong.
In Figure 10
we can see the covariance and precision of a sample distribution where we can clearly see the strong depen- dence of dimensions and 10 relative to the others.
In Figure we can see the naive and true marginal precisions when dimension missing. For this plot the values have
19

---

Score Matching with Missing Data
1
1.0
1
3
2
3
3
0.5
4
4
1
0.0
0
6
6
7
7
-1
-0.5
8
8
-2
9
9
10
-1.0
10
-3
1
2
3
4
5
6
7
8
9
10
2
4
5
7
8
10
(a) Covariance
(b) Precision
Figure 10: Covariance and precision from sample distribution from our normal experiment.
been cube-rooted in order to emphasize the difference between zero and non-zero entries. Here we can see that the naive marginalisation wouldn't capture the dependence between dimension 10 and the other dimensions that gets introduced when dimension I is removed. This means that a naive marginalisation would assume that dimension 10 must have a direct dependence on dimensions 2-9 even when that is true. Interestingly, the rest of the marginalisation seems very similar suggesting that in some potentially less structured cases, naive marginalisation approximation. This supports the results we see in our GGM estimation where highly structured graphs like star graphs are much more affected naive marginalisation than unstructured graphs.
provide semi-reasonable
B. 1.3. NON-GAUSSIAN ESTIMATION
Here we present further experiments exploring the non-Gaussian model presented in Section 5.1.2 Here we fix the dimen- sion as 10. In Figure 12a we fix the missing probability as 0.5 and vary the sample size. In Figure 12b we fix the sample size as 1000 and vary the missing probability.
From Figure 12a we see that both EM and Marg-IW have the smallest estimation error for larger sample sizes. Zeroed Score Matching has the largest estimation error due to its inability to appropriately marginalise the distribution. In Figure 12b we observe that Marg- Var has the smallest estimation error with performance convergence to that of Marg-IW and EM the missing probability increases.
B.2. GGM Estimation
B.2.1. VARYING NUMBER OF STAR CENTRES
Here we present illustrations of the marginalisations for our star-shaped graphs with node and then 5 nodes both with the same edge density.
In Figure 13 we show the covariance, precision, marginal covariance, and marginal precision for star graph with 1 centre where the marginal terms are with dimension removed. As we can see clearly the only meaningful structure left in the graph after marginalisation are in the negative precision terms which the model naive marginalisation fails to capture.
In Figure 14 we show the same thing for the case a star graph with centres. As we can see in the centre case, the naive marginalisation picks up most of the structure as there are fewer negative terms which ignores and also lots of additional positive terms which does successfully pick up.
B.2.2. VARYING NUMBER OF DIMENSIONS
Here we use our same star-shaped GGM in the main paper but with a varying number of dimensions. Throughout 1,000 samples are used and each coordinate is missing independently with probability 0.7. Results are presented in Figure [5
As we can see, for higher dimensions the variational approach clearly performs best however at lower dimensions the other
20

---

Score Matching with Missing Data
1.5
1.5
1
1.0
1.0
3
3
4
0.5
4
5
0.0
0.0
6
6
7
-055
7
-0.5
8
8
9
1.0
9
11.
10
10
-15
1.5
1
2
3
4
5
6
7
8
9
10
2
3
4
5
7
8
10
(a) Naive Marginalisation
(b) True Marginalisation
Figure 11: Marginalisation of the precision to remove dimension by the naive approach (i.e. subsetting the precision) and the correct approach. All values cube-rooted for contrastive purposes.
approach catch up and even overtake This because lower dimensions, IW can effectively model the marginalisation of the score and so the more complicated variational approach is not required.
B.2.3. INDIVIDUAL ROC CURVES
Here we present individual ROC curves from Section 5.2.1 with missing probability of 0.5. Here we specifically present the ROC curves from the first runs out of the 50 performed for the experiment. These ROC curves are displayed in Figure We observe that Marg- Var has consistently the best ROC curves of any of the methods.
16
C. Additional Proofs
C.1. Marginal Score Matching Objectives proof of Proposition
4.3
For the first claim have that
E
1:.0 (Xx)) (XAa)
X
| |X||s 2E[SA, sa(XA)
PPAA A) / pa(ar))sn (daa)d
C
E sA:0 X-| -2 Σ
AEsupp(A)
P(A X) / VeSPA(E>)" (x00)))da
C' + Bl[l
Σ
AEsupP(A)
C' + +X
Σ
P(A = A) Σ /
j=1 Xr
AEsupp(A)
Hence by integration by parts we have
(X))));0((A)) SA (Xx))
E
(Xx,))
2 Σ
P(A=
lim +8x
Px (a( / (x
AEsup((A)
JEA
E| (X+)]
V XA
8 (xx)) +
justifying first claim. Hence if 9 minimised objective then it minimises
E ((()))((XN)
(X.
21

---

Score Matching with Missing Data
3.0
3.0
2.5
2.5
Theta Error
Theta Error
2.0
2.0
1.5
Marg-IW (Ours)
Marg-IW (Ours)
1.0
1.0
Marg-Var (Ours)
Marg-Var (Ours)
0.5
Zeroed
Zeroed
EM
EM
0.0
0.0
600
800
1000
1200
1400
0.2
0.3
0.4
0.5
0.6
0.7
Sample Size
Missing Probability
Varying sample size with missing probability 0.5 and di- (b) Varying missing probability with sample size 1000 and mension 10.
dimension 10.
Figure 12: Average parameter estimation error (L.2 norm) for ICA inspired model with 95% C.I.s for various methods. Lower is better
As there exists a "true" 0 we know that this objective is minimised at 0 and so we must have 8x,0(X>) = 8x(x)) a.s. for all E supp(A). By our assumption this then gives (X) = p(X) a.s.
C.I.I. TRUNCATED SCORE MATCHING
Proof of Proposition A.2 We mostly use (Liu al. 2022). We firstly have that
Elga(Xa)llsA,00X)
(XA)l
(X)x((X)))xx(XXx1])) EXo. (X))) 2Riea (X. (X1))A,(XN)) (xx)))
Now have that
E[ga (XA)8A (X(a)
Σ PPA = (XWSAAAX) (XWWAA A]
AEsupp(A)
P(A 2) Σ / (x(E))S0NN(EN) θ
Σ
px(xx)dxx
dx,
(Esupp(A)
JEN
Σ
PPA = )) Σ
ps pm (m>
JEA x
0x.
(Esupp(A)
Σ
P(A = 2) Σ
/ Ser, (ox (oa
axx,
(Esupp(A)
JEA
P(A = X) Σ /
(x() θ
Σ
dxe 50.2 (ZN); + (x(x)))x ps px
JEA x
dor
(d)
= VA T A (X.) + (X(())sA
where (a) is given by Green's Theorem and (b) is given by our limiting condition. Plugging this back into our original result gives
E (XN)) = (1ax)(x)))) sx.o(xa)) + VX,) (X,)(N))]
LTm((())
From our conditions we also know that 8" is a minimiser of J(8) hence by our conditions the unique minimiser. Therefore 8" is the unique minimiser of (f)
22

---

Score Matching with Missing Data
1.0
5
5
2
10
10
0.5
15
15
1
20
20
25
0
25
0.0
30
30
35
35
-0.5
40
-2
45
45
50
50
-1.0
1
10
15
20
25
30 35
40
45
50
1
5
10
15
20
25
30
35
40
45
50
(a) Covariance
(b) Precision
1.0
1.0
5
5
10
10
0.5
0.5
15
15
20
20
25
0.0
25
0.0
30
30
35
35
-055
-05
40
40
45
45
50
1..
50
-1.0
1
10 15
20 25
30 35 40 45 50
1
5
10
15
20
25
30 35 40 45 50
(c) Naive Marginalisation
(d) True Marginalisation
Figure 13: Covariance, precision, and marginalisations of the precisions to remove dimension 1 by the naive approach (i.e. subsetting the precision) and the correct approach. All values cube-rooted for contrastive purposes.
C.1.2. MNAR PROOFS
Proof of Proposition A.9 Firstly we have that
E SA..E (XA)= SA,E.:0 (X.)) = E [X] SALE) +C
where C does not depend upon 0. Examining the second term closer see that
Σ
sx(x))
#A:0 (xx FA
xxs
(Esupp(A) ()
Σ
Σ
8x(=5)) Ex)daa
AEsupp(A) Xr
Σ
Σ
VarPA(Ex; Ex);
XEsupp(A) jea
XA
Σ
Σ
Px(Ex; FA)) 0j
$2.0 (xx) as result of integration by parts
(EsupP(A) JEX
X
E.A [Vx, (AN)
23

---

Score Matching with Missing Data
5
1.0
5
1.0
10
10
15
15
0.5
20
20
25
0.0
25
0.0
30
30
35
-0.5
35
-0.5
40
40
110
45
11.0
45
50
50
1
10
15
20
25
30
35
40
45
50
1
5
10
15
20
25
30
35
40
45
50
(a) Covariance
(b) Precision
5
5
1.0
1.0
10
10
15
15
0.5
0.5
20
20
25
0.0
25
0.0
30
30
35
-055
35
-0.5
40
40
-1.0
1.0
45
45
50
50
1
10
15
20 25
30
35 40 45 50
1
5
10
15
20
25
30 35 40 45 50
(c) Naive Marginalisation
(d) True Marginalisation
Figure 14: Covariance, precision, and marginalisations of the precisions to remove dimension 1 by the naive approach (i.e. subsetting the precision) and the correct approach. All values cube-rooted for contrastive purposes.
Hence we have
[1esss(x)]]+++2xx sAo((Na]] +
Hence, just as in Proposition 4.3 this is minimised when
8x(Xx) = (x=0(xa)
a.s. for all À E sApp(A). We then have that for any À E supp(A)
(x = (x,0(()) for all 2y E XA
PA (Ix) PEX:O(EX) for all σλ Xx
P(EA|KA)PA (TT) Ph,((E))((AA(EA)) for all as E X'
4 (X(A)) (1:O(EN)) for xA X'
24

---

Score Matching with Missing Data
0.95
0.90
AUC
0.85
Marg-IW (Ours)
0.80
Marg-Var (Ours)
Zeroed
0.75
EM
1
5
10
20
50
100
Dimension
Figure 15: Mean AUC of various methods for edge detection of star-shaped GGM increase the dimension presented alongside 95% confidence intervals.
C.2. Finite Sample Bound Proofs
Here we give some key results alongside their proof allow to obtains finite sample bounds for truncated score match- ing. This first result is what really underpins approach.
Proposition C.1. Let (0,d) d) be a compact metric space for any denote the ri-covering number of by ( θ). Now define random functions Z0), Z, (0) IR for all (with deterministic) such that for any 0 E θ
P (1ZZ(00) - Z(0) 8(E,r)
we additionally assume that Ze, Lipschitz with constants C,,, C respectively then we have then
(,) - > ((0))))((0)))++)((() n (C+++ S N(n,
sup
6EB
Proof. Let 0,...., (N((,) be an cover of θ then we have that
sup (0) 2(0)
sup
min ()
(lZ,(0)
Z(0))+ 1(00)+|Z,(00)) 2(01)13
Ic0
(0, - Z }
< sup AE0
min jE(N(n.))
( + 2(00)}
max
(jE(N()))
< sup
min
{C,, + C)10 -010--111 jE.N(n
max
(1Z,(0) - 2(0)))} by our Lipschitz condition
#E0
[
< (C, C'n++
max N()
(0,) - 2(0,))) by definition of 0,,...,0((n..
Therefore have the following relationship between events:
sup 7,(0) Z(0) c 2M1 C 0EB
U
12,(01) Z(0)) >8
jE(N(n,9))
25

---

Score Matching with Missing Data
First 4 ROC curves for missing probability 0.5
0.8
0.8
0.6
0.4
Marg-IW (Ours)
Marg-IW IW (Ours)
Marg-Var (Ours)
Marg-Var Var (Ours)
0.2
0.2
Zeroed
Zeroed
EM
EM
0.0
0.8
0.2
0.8
1.0
0.8
0.8
0.4
0.4
Marg-IW (Ours)
Marg-IW (Ours)
Marg-Var (Ours)
Marg-Var (Ours)
0.2
Zeroed
Zeroed
EM
EM
0.0
1.0
Figure 16: ROC Curves in 4 separate runs for GGM estimation of a model truncated normal distribution with star-shaped Precision matrix.
26

---

Score Matching with Missing Data
This therefore gives
sup - -(xx -) 03) > AEE
>E+m(C, (CCC, +c))
n
(2,(0)) - 2(01)12 e) )
<P
(E(N(n,)
N)n
P ((IZ,(01) 2(0))}>))
==1
S NN, O)((Em).
Remark C.2. This result does not require Cr,, C' too be deterministic. A feature that we will be exploiting later on.
To be able to say meaningful statements about our Lipschitz bounds for our proof, it will also be helpful to make subgaus- sian statements about nested sums. We give to results to enable this now
Lemma C.3. Let X,Y be independent RVs and define function : X' x y 7 R. Now suppose that for any x E X, g(x, Y) is nau-Gaassian with parameter a. We then have that EXg(X, Y)[Y] is sub-Gaussian with parameter a
Proof. As g(x, sub-Gaussian we have that for any À 0
(22)
E Texp (g(xx, Y) Big(x, Y)D) & exp
02
Our aim is to then use this to bound
E (exp (A -l,() Y)
We first have that E[g( = = which turn gives
E (exp (A (Blg YY]
E E Texp (> El,((X [У]
by Jensen's inequality
IE [E (exp () ()) X))()) [X]]
S sup E [exp {) (g(x, Y) - -[o(z,))])) by independence
22 0²
exp
Lemma C.A. Let X, be independent RVs on spaces X, y and g : X x y R.s.t. for any x X, y g(x, Y) and g(X,v) v) are sub-Gaussian with parameters ox respectively. Let (x(0) Jim and [y(k)]r be IID copies of X, X,Y Y respectively then
e'
exp {
P
9(x()) (())) Eig(X,YY) >E s nexp
Proof. Again let W 1 ((1)) n (+)) V) -E(
x)] We aim bound W well as
27

---

Score Matching with Missing Data
For W() we have that
P(W()) >
>)
=[(((W)) x[X(1)]
)
E
Σ (( -
- -Eg(( x)(x(1)
< x(()
c2orr } = =exp { sSoor
E
exp
2
Therefore we have that
P
W
IP
Ü (
1=1
P(W()) =)
1=1
s
S <nexp
Additionally from Lemma
C.3 we have that [X, Y) is sub-Gaussian with parameter 01 giving that
8 n
P
(XX -lo( -
E F[g(X
>e) sexp
Hence combining these get
P
9(x E(a(X, >5 S nnexp
exp
8
To proceed define the intermediary step between the population objective and the sample objective this is
1 1 11 (x
LTM:n (0) =
(x +
X
Proposition C.5. For E let (X be IID copies of X' ~ p' and assume that supp(X) S supp(X'). For e and qu : X + [0, oo), with [[gal] ||| < oc define RVs Yo, Ya,r by
Z(0) == LTM.1 = (X((x))
(xx, (x,(( + (1(x)(()1122) (Xnga(XN)) sA,o(XA)
+ (11 + I (XW
Z,)0)= == LTMILI(0) = -A
Suppose that the following hold all x E X, X E supp(A)
0 < <o < (0x 0) < bo
[x, (1)x(1,0)11== by
(2,,((x,0) bn
(x(2))) co
28

---

Score Matching with Missing Data
V wx (+x) 9 Then have that
m
(() Z(0) e) (4+ (4+ exp
a²
with a depending upon oo, bo, b1, b2, co.
Proof. First define
(1) = 1 -
xa x(00).
(0) =
R=1
Using the definition of the marginal and estimated marginal score we then have that
Y2( + [0)1 1
T Yi(0)
Z = (XA)))
Yo( 2 Yo(0)²
2
Yo(0)
(0)
1 (0)1
V (1(0)
(i(0)
Z,( = (XWN
+ 2
(0)
2 (o)
Yo(0)
We can therefore write write [0) Z(0 as
(0)
1 (0)1,
Y20) 1|2(0))
(Xa)
(0) Z(0) a
+
(0)
2 (0)2
Yo(0) Y2(0)²
Yo(
(0)
12,r (0) Y2(0)] (0) - -Y2(0)|(2((0)
Yo(
Yo(0)oo,-
1
[0)]
-|]|i~-|XXM (0)1)
(6()))
+
+ 2
2
Yo(0)²
(0)))
V (Y), (1))
(0) - Yo(o)] V (XA (1)
+
Yo(0)
(o.))o..-
Y2, (0) Y2(0)
(0) - 12(0))b2
ag
1 (0))
[0}]²
(0)2 Yo(0)2
+
2
2
ad
c1 + (0) Yo(0)) cbb + (0) (0)
oB
60
12, (0)
Y2)))
(0) Yo(0)) b2
00
ag
1 (0) Y1 (0) 01
1
Yo,- (0) - Yo(0) bob2
+
+
2
of
2
of
(0) (0)
(o) - Yo(0)) cib1
+
40
aB
Now if we define the events
a
ag
a
}
Eo( =
(0)- - Yo(0)>
6bz
3bn b
6ez by
(0) [1V Y1(0)|>>
ag
}
3b1
6c1
E0 (2- Y2(0)>>
86
29

---

Score Matching with Missing Data
then have that
((()))) (12,(0) Z(0) a) S UE
1-1
Using standard Hoeffding bounds for Eo(O), E2(0) and union bounds in conjunction with Hoeffding bounds for ( we get that
ad
aj
((E)(O)) 2 1- 2exp
-62r
Digb
36671908
(0) 2 2dexp
-
d0
aj 36c
9b
ag 2
(()) 2 2 exp
36b2
Hence
m2m
(() Z(0))> 8) (4+ 2d)
exp
a
with
6bgb2
31
60g by by
312 o
60g bj
6b2
a max
0g
ap
au
dg
proof of Theorem A.15
Our strategy will be follows:
Use our bound on !Z) 7(0 from Proposition C.5 to bound [LTM, I..(()) (0)
Use covering number argument alongside Lemma C.4 to bound suppl TTMM;,, (0) - I.n (0))-
Use similar approach to bound suppl (.) (f)).
Combine these bound [TTM;,,- [f)l
For the first step we have that
PP P( (f)) SP
Ü
( (0) [1) E)
1,r
"
P([L,,, (0) - (0)
P((Z,(()) Z(0) 6)
i=1
emm
n(4 + 2d) exp
a²
30

---

Score Matching with Missing Data
where
: Ex= (X,, X 0312
(X 00 Ei fu.a (1) (k)
(0) 94, (x
+
2
Eim fon( ix' 00)
x (K)
9A,
+ fo(XX )
E ( X 0) xx
E
f (1) X X 0. xxx
(0) 9r. (0))==(xx(()) (
1
fo A x
1..0) .0) xxxx)
2
[oaaxs . 0) 1)) xx((
E
Elgs, (X (X ( X x
( X
We now try and derive a Lipschitz constant for both Inr and T.n. Define 9 X - R by
1 M2((a) +
62
b2
a1
+ 9
g(x)
+
Mo(x) +
Mi (a)
Ap
2ag
og
2ng
Op
Then we have that For La, we have
1 n
g g ,xx")) (8)
LTMm (0) La,r (f') S p(0,0'
nr
i=1 k=1
similarly we have
1 2
(0) s p(0,0) 0*)
Using an identical argument to Proposition C.1 can get the following bound for any m > 0:
E a
P(supILTM(00) LTM. (0)>> e + mm (Cn,r + C.)) N(m. ()n(4 + 2d) exp
de@
Hence we now need bound C, and Cn,y. We know that both of these terms converge to == Efg(XA, X,). To obtain rates on this convergence we make sub-Gaussian statements about g. Specifically for x E X' have that ((X,, x-x) is GabbGansiin with parameter
b2 ag
bf
cib1
a1
02, +
0oA +
ag
2a
2a
ag
We can therefore immediately bound Cn - C using Lemma
C.3 and Hoeffding bounds to get
PCC, - C' > s)
= - C > c|A)]
e2n }
E
exp
80 ]
e2n
< exp
80
31

---

Score Matching with Missing Data
where a max Acsupp() AA- To bound Cn,r C we can use Lemma C.4 to get
( (1, r(k)) - -X) > E [A)
E [P C =)
T21
k=1
E2m
=²n
}
I E
71 cxp
+ exp
80 -A
800
e'm
=²n
= n2 exp
exp
80 2
80
with o', d'ا define identically to a, a with LLA replaced with o/ and σr replaced with d'_x- As a result we get that
P(sup/LTMm..r (0) (n(0))> + 27 (C(+))) AE
Nm,, 0)n(4 + 2d) exp {
=2m a²
neexp {
52m
{
E2n
2 exp
80'2
802
Now we have the bound supø LTM. (0) (0)l We aim to bound supo LTM: (0) LrM( To that end we
have
b2
1 bj
cbb a'
ag
2 ag
ao
a.s. and so we can use Hoeffding bounds to get
E2n
P LTL LTM(())>>>>> 8) S exp
2a
Again using argument similar to Proposition
C.1
we have for any 1/2 > 0
P(sup (0) - LTM(0))> E2 + m2 + - N(2, ex
{
E2n
+ P(C,, -0> E3)
2a" 2
Combining these two results gives that for any & >> 0
(())))LTM(0) - LTMm. +2 + 7/1 (2C + 25a) + m2( + E3)) 0E6
< 0EE
- =1 (C 53) U sup
LTM: (0) LTMM (0) 02 m²(2C +s))
< Psupp (TMSupLLMM(0 (0) - L.. (0)) + 2/) (C + 2a)) + PP(ssp| lE8
LTM, (0)
LTM., (0) 82 + m2 +
BEB
E2n
< <N(1, ()n(4+ 2d) exp { -==maa? 2| N(12 θ)
exp
PCC - C'> 53) + P(Cn,r -C E3)
2a
$ăn
))
s3m
{
< N(m, () 2d) exp { m=² N2
exp
+ nexp
+ 2 2xxp
(22)
20x
80
802
Take 7n 1/r and 1/2 1/n so that for sufficiently large n, N(m, 0) S exp (plog (33/2) diam(O)r)) and N Q)
exp log (33/2) diam Plugging this into (2)
gives
P suppLrm((0) LTMM (0) 5 + + 1/r 25s) + 1/n(2C++E) 000
3 S nn4 + 2d) exp p plog , diam(O)r eina } exp {
3 plog tiam(()n
+ (2 + n) exp
2a'
8(0 V
32

---

Score Matching with Missing Data
Now if we take each of these terms be equal 5/3 gives for sufficiently large n,r
log(dnr diam( + $2
p plog(n (iann()am())))
P suppLrm((0) LTM:n, -0)..
V
+ Bg nttr)
log(n/6) <6. 6.
c++
n
where 81 82 = =9a', 10(a V o'). As there exists 8 O minimiser of LTM(() now have that [L(O.) [(0*) < [L.))) L.,
L(0*)]-
Finally we know that Jt) = LTM + C' and under a correctly specified model L(0*) = O. Therefore we have that
((]J(On,))(> 2=) S P(supl LTM (0) LTM(0))> s)
lEe
and so have our result simply replacing 3g with 2
C.3. Gradient First Proofs
Proof of Lemma 4.8 First we have that for any X, xh.
(0:)( Væy log po(ar)da-x
J po (oxddxx_>
(x)
p(x-->)
po(
Vax log (o)
po (x=
A (s0((())))
Taking expectations on both side w.r.t. XA) proves equation (
For we have again for any X, 22
/ V(po(x= (x)
V (e X'xx]
/ Vole=sle>)V + / /Volt=
po Vg + / /ooo==/oa/Voogrrzzz
E'(Vgo X" + E log Lo X "
- E'IV log po (x)gg
E'
log X' )) X2
40,1 (m))) [go ( X'xx)
E'(V X" + E log Pe ( X" "_)) X"x)]
Pe (2) - EE' le
X X2
X'
( ( X' w + Cov (s (I X' ,), 9o (I x X' ))
Where again here V represents the gradient w.r.t. either xx or We can again take expectations on both sides w.r.t. A, XA get our result.
33

---

Score Matching with Missing Data
Proof of Corollary 4.9 First define
(m()) = Σ 20 + sx,0
JEA
so that [0 XA)] = (0) Then using our two score identities,
8
(9) &
we can re-write L(0;
Lm (0: =
Σ
X x)1 by 8
JEX
Σ
E'' x' - X x)1 + Var(so(
by 2
Σ
E' (5 (E), X" X->))| + xx x 2E
JEX
Where all expectations w.r.t X A = =x; 0 Now using the fact that
Vo (E'Iso(ex,
( V''' x" X"1]
and using again on each term of the above gives followed by taking expectations w.r.t. A, XA gives the result.
C.3.1. TRUNCATED SCORE MATCHING
Proof of Corollary 4.11 Now using Lemma 4.8 we know that this can be re-written as
Σ
91 (0 X + 2
X
JEX
2E +1
XLx)]
JEA
Now taking gradient w.r.t. we get
(oL(0; = Σ (xx())) - 2E (E(ss((a)))(((V(
(0)))))(())))))))))))))))))))))) sz)));))
JEA
2 (V p8g Cov'
Cov 2e
20,g(0x) E' (
C.4. IW and Gradient First Relationships Proof of Lemma 1. We first have that,
$0.2 (xx)) = Væ, log
Σ W1
1=1
2i=1
Tax" W/
E i=1
W;
ú V ma log qo ( x²())
Eiw V loggo (I), x x1(21
Ein (so x((()))
34

---

Score Matching with Missing Data
Where the penultimate result uses the fact that Vay wi w; log x )
Additionally,
V Ei=1 ( (((6))
(Vax Σx km w%)
EE Rm Wago(ax) xx(())
Vax Eiolg X_xx) For the second term we can again use Vx, Wk
Ex=1 Wk
Ek=1
W/ V logpo x²x to write this as
Vay, logpo X" X²x Bii 900 (2 X'2x]
The first term we can write
Eラ=1 we log po ( X (1) E wk (2EA, X'x)
Σi=1 Wi
Ei=1
EE [lo pe (æx, X x100 (2 (x=xx''xx)) + Bivl V., V (2 X'xxl
Combining these 2 gives our desired results.
Proof of Corollary A.79 The exact same arguments give us the second result but with the importance weighting identities &
19
and 20
replacing
D. Additional Score Matching Details & Extensions
In this section give some additional details score matching and introduce some pre-existing score matching exten- sions and methods. For some of these approach can adapt our method to work with them while others act comparison points for our approach.
D.1. Classical Score Matching
The assumptions for the classical score matching result presented in
Assumption D.1. For
(a) The pdf p(æ) is differentiable w.r.t. a;;
(b) Our score estimating function 8g is differentiable w.r.t. æ;;
(c) so (X) [[]] [X)][]] oo;
(d) (x)) O whenever ||æ||-+ co.
D.2. Zeroed Score Matching and MissDiff
In this section we take x to be the vector in Rd with (a:); = Iz if i E X and 0 otherwise and then take X be the RV equivalent.
MissDiff is generative modelling technique that aims to learn generative models from corrupted tabular data using diffu- sion models with denoising score matching. Throughout their method is assumed to be neural network. The core idea is to replace the standard score matching objective with
(x,xA(0)xx.1) (x())))((((()))))) Vx.0 log (Xxx))()))
Essentially, missing values are zero imputed and then the score is tested only on output dimensions whose input dimensions are non-missing. The idea of this approach that trains sæ ( (t), f)x to approximate the true marginal score 8 (æx ()), which in turn encourages so (=(()), to approximate the full score se t).
For comparison purposes, we adapt this in two ways to create Zeroed Score Matching. Firstly we change the objective to standard score matching and secondly we no longer require 8 to be a neural network. Thus our adapted version of the MissDiff objective is
Lieeoee (0) V (Xx)a (x(X)) [{5o(X])
(23)
35

---

Score Matching with Missing Data
The key issue with this approach is that with 89 no longer necessarily a neural network it is not reasonable (or indeed even possible) for so model both the joint and marginal scores. In other words we cannot expect both so (ax)) s (21x) and sa(x) s(x) making it naive marginalisation method for the score.
We can give a brief example with multidimensional normal distribution. Supposed that X ~ (O) with P = -- Then the score is s(x) = - Px which we would model with so(x) -oox. Under the Zeroed scheme, we would take the marginal score to be
(),0(x) (-Po2))
However we know that if X ~ N(0, Ex)
with Eg Po then X2 N(, Ex,2:0) and should take
(x,0(x)) -E 1ax = - )xa)))>)) Txx x
and crucially, PA,A # ((P-))x.x) unless P=
We illustrate this for our simulated experimental settings in Appendices B.1.2 & B.2.1
We also explore the implications of this for the marginal Fisher divergence for the normal in Appendix A.4
D.3. Sliced Score Matching
One issue with score matching is that (xcdotss is computationally expensive to compute for large d. solution to this proposed in [Song al. (220) where, instead of testing the full score 1-dimensional projections/slics are tested instead. This is done by introducing another RV V on (X,Bx) satisfying [E[]] 0 and E[VVT] positive semi-definite.
The original objective is then
Fs(O) = (VV ((V))=)((((((()) so(x)))²]
Which then leads to the following equivalence
(0) = R2 [xx so (x)))) V + (V sx(())
Fs C.
which is less computationally expensive w.r.t. and can be approximated with samples of X and V. For the proof of this results and the precise conditions see [Sng al. (2020).
D.4. Denoised Score matching
Denoised score matching is another adaptation which removes the need to takes derivatives of the score all together Vin- cent 2011). As we will see later on it is also the method used most prominently in diffusion processes.
In denoising score matching we construct a collection of RVs (x()) with X(0) = X and X ((m(() Ip)) We assume that the noise (t) grows sufficiently so that X(T) is approximately an Isotropic Gaussian, The aim now to estimate 8(x, t) == V log pr (æ) where pt is the PDF of X X(1). The denoising score matching objective is then.
(x)()) - Vx( Vx (xg)(X()) = E 5x 1) (X((())X()
where here t is treated as random and uniformly distributed on (l...., T}. p(x()))=(0)) is the transition kernel from (0) to (() and X(L) is a self-specified weighting function over time. Due to our choice of noising process, V Va(e) ()) = t (t) (0))
Remark D.2. Convention is to up-weii weight larger values of as earlier parts of the reverse diffusion process (hence later parts of the original diffusion process) are seen as the most complex and where most of the data's 's structure learned.
36

---

Score Matching with Missing Data
Our estimate is thus
i = argmin #ERR
Σ Y1
(sa)))ls((x())) Vxw (oxx(x()))xx)))))
where ,,,..., tn are sampled uniformly from [0,
Remark D.3. Originally denoising score matching was proposed for a single fixed however for the purpose of generative modelling and annealed Langevin dynamics, multiple noise levels or even a continuous noising processes are used.
E. Additional Details
E.1. Objectives
E.I.I. MARGINAL IW SCORE MATCHING OBJECTIVE
be IID copies of Xx, with known PDF p' (æA,). Our full sample objective given by Let (X Tk=1
xx xx0 x x
E
E
La (0) =
Σ 24 x(:
+
(24)
(x(x D' F' (1i07)
(o(x(())
E
p'(X (11))
E. 1.2. MARGINAL TRUNCATED IW SCORE MATCHING
Our full sample objective for truncated IW missing score matching is given by
3 1 X x2
3g ( (0
3
1 =k=1 ((X(
12
k=1
((
L L.m... (0)= Σ Σ g(X g(XA(0) 71 JEN
20,
(25)
20(x)))
(0(x())
E
1 Σ
p'(X (ill)
// 1'(X' (i
0x xx ((
# Σ
+ 20, 20,g(Xx,)0
2x(x)))
Σ
v'(X (ri))
E.2. Variational Modelling Details
For our purposes our variational model 1/S. has to able able to model px(xx_)/xx for any value of À C d. For our model take X' x (()()))),
Hence require yal take in any subset of coordinates and output the complementing coordinates We achieve this by creating a'g. d-dim in to d-dim out Neural Network with 2 hidden layers of 200 nodes as per Burda et al. 2016). For the input we replace æy with æ the zero filled version inline with the approach taken in MissDiff (Ouyang al. 2023). That (a*), 0 if j & X and is 2, otherwise. As output then simply take the appropriate coordinates. We can write this more succinetly as
me(xx) = (1)))
We also experimented with making 16, a 2d dim input NN and taking
M (mx) = p/s (x22. m)
where m is the d-dimensional binary mask for the corruption similarly to GAIN YYoon et al., 2018) however this did not seem to provide any advantage. We also experimented with making as depend upon Φ however we found this made training much more unstable.
37

---

Score Matching with Missing Data
E.3. Experiment Implementation Details
E. NORMAL DISTRIBUTION ESTIMATION
The mean is taken fixed at μl =
0.5)T. ( We randomly construct the covariance by first sampling eigenvalues uni- formly on (0.5, 1.5] and then sampling choosing eigenvectors uniformly on the unit hypersphere for the first dimensions. We then construct the the 10th dimension strong dependence on only the first dimension by taking X 2 X + 1 Z with Var(Z) Var(X1)) The data is then truncated to be above the 10% quantile or each of the first three dimensions.
In each case batches of 100 samples were taken and learning rate of 0.01 was used with Adam used as the optimisation algorithm. Our score model was parameterised in terms of the Cholesky decomposition of the precision matrix in order to ensure the Precision estimate stayed positive definite. For our Importance weighting and the EM approach of Uehara et al. 2020), an isotropic Gaussian with mean O and coordinatewise variance of 16 was used.
E.3.2. GAUSSIAN GRAPHICAL MODEL ESTIMATION
For Gaussian graphical model estimation we add 1.1 regularisation thereby modifying our objective be
(0) +7
Σ |Pall
12j3''d
where 9 = (J, with P being our precision estimate. We minimise the objective by performing proximal gradient descent on IsTM((), Specifically for a learning rate v, a current estimate of 0 given by 0, and an estimate of the gradient given by T/. We take our estimate to be (yp - vn) where
B 7 for 3 > 77
(33) =
0
for - ym S 3 2 mm
8 m for 3 < -111
In our experiments we start with a precision estimate being P / and with large value o and, run 200 iterations, and then decrease the value of every 10 subsequent iterations for 100 sequentially smaller values of y. At the end of each block of 10 iterations the precision matrix is taken and an adjacency matrix is produced by thresholding the entry's 's values at some small value. The TPR and FPR are then calculated for each of these increasingly dense matrices and then an ROC curves plotted using these values. The AUC of this ROC is then computed which the statistic reported the plots.
We took LI regularisation ensure that at the highest level the graph had no edges and the lowest level the graph had all possible edges. For this experiment, this was achieved with y E (10-1.7, 1000 10-4). Throughout we took the threshold for edge presence to be 0.002.
E.3.3. REAL WORLD DATA
For these experiments we chose the range of L1 regularisation similarly to ensure a full range of edge densities (here this meant 7 E (10', 10 1) and then constructed a semi-automated procedure for choosing the threshold for edge presence. There is precedent for choosing the detection threshold after LI regularisation as per Fattahi & Sojoudi 2019). We did this by choosing the non-zero threshold at the smallest value that gave a sufficiently smooth increase in edge density between the snapshots where we sample our estimated adjacency matrices.
Specifically the smoothness we were trying to achieve was avoid sudden decreases in the edge density as our regularisation level decreased. Our rough measure of this was to sum up all the negative jumps between sequential adjacency matrix estimates where the previous jump had been positive. This sum was then taken as our measure of "jumpiness" with larger values representing a larger level. Visually inspecting the change in positive level over time we find that a level of 0.01 for high-dimensional cases and level 0.05 for low dimensional cases represented a relatively smooth change in edge density. The smallest non-zero threshold that satisfied this was then chosen by iterative shrinking grid search.
To test the performance of our adjacency matrix estimates, estimated the adjacency matrix using standard score match- ing on the non-corrupted data. We estimated these adjacency matrices at 5 pre-determined values of edge densities given specifically, 0.05, 0.1, 0.15, 0.20.25. This lead to 5 different "ground truth" adjacency matrices. For each method, the AUC was calculated for each of these "ground truth" adjacency matrices and these 5 AUCS averaged. For each missing probability, 25 random samples of the corruption were produced and this AUC metric calculated. The average of these was then plotted alongside 95% confidence intervals.
38

---

Score Matching with Missing Data
The S&P 100
taken from the S&P 500 data between 2013 and 2018 given in https://ww. kagg e com/
dataset. http///aaangnnaaannaaeeeeee sandtt500 with the 100 stocks that made up the S&P 100 taken from roughly the mid-point of the time period which we obtained from http::s//ee...iiipeed....o.
bldic=656413597
The yeast data was obtained from / was accessed pontttss//wwwe http://wwww/..i....//o///......t
nobi. n nlm. nih. http://////e///a/////c..////////
and other subsets have previously been studied in the context of GGM estimation in
Yang & Lozano (0115).
All data can also be found in the
GitHub
repository
at
http:://giiaaub..mm///giiigvvee
39


Reviews:
[{"review_id": "1MQwPApIJa", "reviewer": "Reviewer_iocR", "summary": "The given paper presents novel framework for Score Matching (SM) in missing data framework. Particularly, the setting assumes that the (multi-dimensional) random variable under consideration has missing coordinates. The authors solve this problem by using an auxiliary variable that denotes the masking random vector. Using this, the authors introduce 'Marginal Score Matching' where the score matching is performed only on the visible coordinates of the random variable. The authors then provide several (tractable) variants of marginal score matching - truncated SM, denoising SM,  sliced SM. However, these methods are still intractable due to marginalization integral present in marginal score definition. To this end, the authors provide two ways to estimate these loss objective - (a) using importance weighting and, (b) using variational approximation. The authors present empirical results on synthetic datasets, PGMs and real world tabular datasets like Stock prices and Yeast data. They compare their method against MissDiff and EM to show the efficacy of their method.", "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": "Thank you for your useful comments and feedback, we address some specific points below.\n\n>The proposed method makes sense, however, I am not sure about the evaluations. Primarily\nbecause of the limited metrics (AUC) used to demonstrate the efficacy of proposed method.\n\nThanks for pointing this out! \nSince we focus on the parameter estimation problem, \nwe have now added new metrics specifically designed for the parameter estimation problem in standard and truncated normal estimation: https://ibb.co/kk8LKNV and here: https://ibb.co/pjDzCpGQ respectively.\nAdditionally we have included a new ICA inspired parameter estimation experiment (more details in 2nd paragraph of response to reviewer 2txV) which can be found here: https://ibb.co/7JTJCc6b .\n\nIn Gaussian Graphical Model tasks, we focussed on AUC as we were treating it as an edge detection problem and wanted to demonstrate the efficacy over a variety of sparsity settings. The use of classification metric in graphical model learning tasks through the False Positive Rate (FPR) and True Positive Rate (TPR) was used in previous works [2] and [3], thus, we are following the convention, using the AUC to summarize the overall performance of the graphical model learning. \nFor further illustration, we have included a plot of the ROC curves for the first 4 runs of the star based precision with a missing probability of 0.5 here: https://ibb.co/35dgkVnK and will be sure to include them in the appendix of the final paper. In these we can see that the proposed Marg-Var seems to consistently outperform the other approaches.\n\n>Marg-IW seems to have similar performance as EM. Is there any insight into this?\n\nThanks for raising this interesting question.\nWe believe this is because the both methods rely on importance weighting to approximate marginal/conditional expectations. Additionally they both aim to maximize the marginal score (EM can be seen as maximizing the marginal likelihood iteratively.) \n\nHowever, the detailed reason of their performance similarity is definitely an interesting observation which is worthy of further exploration and we will highlight this future direction in the revision. \n\n>I am confused as to why these objectives cannot be extended to neural networks? Can't one optimize neural network parameters using these objectives?\n\nYes, You are correct that this approach can be extended to estimating neural network models (such as energy-based model) from missing observations, similar to how its variant for handling non-missing observations [4].  \n\nSince most neural network models are also unnormalizable, \nthis extension to our method would further expand the usage of such models. We see this flexibility of handling various unnormalizable models as one of the core strengths of our method. \nWe decided to focus on analysing the performance of score matching through simpler and more interpretable models such Gaussian Graphical Models as this was something which we felt was under-explored in existing missing data literature. Moreover, estimating truncated graphical model using score matching has been higlighted as an application of score matching [2,3]. \n\nApplication of our methods to various neural network architectures and its use in downstream tasks such as learning a diffusion model from missing data is definitely a strong direction for future research, and we will make sure to highlight this in our conclusion. \n\n## References\n[2] Lin, L., Drton, M., and Shojaie, A. (2016). Estimation of high-dimensional graphical models using regularized\nscore matching. Electronic Journal of Statistics, 10(1):806 – 854. Publisher: Institute of Mathematical\nStatistics and Bernoulli Society.\n\n[3] Shiqing Yu, Mathias Drton, Ali Shojaie (2022). Generalized score matching for general domains, Information and Inference: A Journal of the IMA. \n\n[4] Song, Yang, and Diederik P. Kingma. \"How to train your energy-based models.\" arXiv preprint arXiv:2101.03288 (2021).", "claims_and_evidence": "Most of the claims made in the paper are Theoretical which are rigorously justified via mathematical proofs. \n\nThe empirical claims seem a little unclear. While the authors point out problems with EM-based methods in Related works, it seems that their method doesn't provide significant advantage performance wise (E.g., Fig. 1).", "methods_and_evaluation_criteria": "The proposed method makes sense, however, I am not sure about the evaluations. Primarily because of the limited metrics (AUC) used to demonstrate the efficacy of proposed method.", "theoretical_claims": "I went through the proofs of most of the results. The authors have provided solid mathematical grounding to the proposed method. I am happy and satisfied with the claims and proofs.", "experimental_designs_or_analyses": "The experiments are primarily performed on synthetic datasets. While it helps in understanding the validity of the method in simple settings, it is not clear how it would perform in complex settings (such as image generation/inpainting, etc). However, since the paper has theoretical inclination, I won't count it as a negative. \n\nPerformance wise, it looks like the proposed method, particularly Marg-Var has very slight advantage over EM whereas Marg-IW seems to be performing as well as EM. Can the authors provide an explanation for this?", "relation_to_broader_scientific_literature": "The work essentially extends the various score matching objectives to the missing data setting. While there are methods like Ambient Diffusion or Cold Diffusion that operate in more generic setting, the current work has theoretical flavor to it which extends the existing approaches to missing data settings.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "#### Strengths\n1. The paper is well written and presented. Although heavily theoretical, I was able to most of the claims and proofs.\n2. The proposed setting has practical relevance, specially for tabular data.\n3. Going through the proofs, it seems most of the SM objective can be adapted to missing data scenario.\n\n#### Weaknesses\n1. The experimental section looks a little weak. Perhaps including tabular data benchmarks could help with this.\n2. Marg-IW seems to have similar performance as EM. Is there any insight into this?\n3. I am confused as to why these objectives cannot be extended to neural networks? Can't one optimize neural network parameters using these objectives?", "other_comments_or_suggestions": "1. Line 102: \\del_j -> \\del_1 ?", "questions_for_authors": "N/A", "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": "BWSJZM7cZT", "reviewer": "Reviewer_rz1W", "summary": "The paper explores score matching for missing data, proposing two distinct approaches: importance-weighted score matching and a variational approach. The importance-weighted method relies on reweighting score matching objectives using an auxiliary distribution over missing variables, while the variational method frames score estimation as an optimization problem using a parametric variational model. The authors provide theoretical justifications for both methods and evaluate them empirically on synthetic Gaussian data, as well as two real-world datasets: S&P 100 stock data and yeast gene expression data. The results indicate that their methods improve Fisher divergence estimates and AUC scores for missing data imputation, showing modest but consistent improvements over baselines such as Expectation Maximization (EM) and a modified version of MissDiff, a recent diffusion-based model for 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": "Thank you for your useful comments and feedback as well as the additional references which we will make sure to include. We address some specific points below\n\n>However, the scope of the experiments is limited. The reliance on Gaussian data raises concerns about the generalizability of the method to more complex distributions ... Expanding the evaluation to non-Gaussian data and stronger baselines would make the results more conclusive.\n\nThanks for your comments! To further enhance our experimental section, we have included a new experiment on parameter estimation for an intractable ICA inspired model where $p(x)\\propto\\exp\\{-\\sum_{ij}\\theta_{ij}x_i^2x_j^2\\}$. Links to experimental results on estimate $\\theta$ under varying sample size, missingness proportions, and dimensions here: https://ibb.co/7JTJCc6b. As we can see in higher dimensional settings, Marg-Var clearly outperforms other methods while in lower dimensional settings, both Marg-Var and Marg-IW perform best (with EM performing similarly to Marg-IW).\n\n> The MissDiff baseline is simplified, making it unclear whether the proposed methods would outperform a full diffusion-based model.\n\nIn this paper, we primarily focus on the parameter estimation tasks. MissDiff was originally designed for estimating the score function in a diffusion process, not the parameterized density model. As such, Missdiff, in its original formulation could not be directly applied to our tasks as it does not provide parameter estimates. This also holds for any other diffusion based imputation approaches such as CoRR which again do not give us a direct parameter estimation procedure.\n\nWe include the MissDiff baseline more as an example of what naive marginalisation of the scores would look like and to demonstrate that the idea behind MissDiff can't naturally adapted to the parameter estimation problem. We will make sure to clarify this further in the paper and change the name of the MissDiff-Param approach in our results to make it clear that it is a variant of MissDiff that has been adapted to our problem.", "claims_and_evidence": "The paper presents two main methodological contributions—importance-weighted score matching and a variational approach for marginal score estimation—and claims that these methods improve score estimation and missing data imputation. The theoretical justifications for both methods are sound, with well-formulated derivations that clearly support the proposed objectives.\n\nHowever, the empirical evidence is somewhat limited in scope. The experiments primarily focus on Gaussian data, and only two real-world datasets (S&P 100 and yeast gene expression) are used for evaluation. While the results demonstrate modest improvements over baselines in terms of Fisher divergence and AUC scores, the gains are not substantial. Additionally, the choice of baseline for comparison raises concerns—MissDiff is used in a modified form without its original neural network component, making it unclear how the proposed method would compare against a full diffusion-based model. Furthermore, the importance-weighted approach assumes knowledge of an auxiliary distribution  $p{\\prime}$  over missing variables, which may not be available in practical applications, potentially limiting its real-world applicability.\n\nOverall, while the theoretical claims are well-supported, the empirical evaluation could be more comprehensive, particularly in terms of dataset diversity and baseline fairness.", "methods_and_evaluation_criteria": "The proposed methods—importance-weighted score matching and a variational approach for marginal score estimation—are well-motivated and theoretically sound. The evaluation metrics, including Fisher divergence and AUC scores, are appropriate for assessing score estimation and imputation performance.\n\nHowever, the scope of the experiments is limited. The reliance on Gaussian data raises concerns about the generalizability of the method to more complex distributions. Additionally, the MissDiff baseline is simplified, making it unclear whether the proposed methods would outperform a full diffusion-based model. Expanding the evaluation to non-Gaussian data and stronger baselines would make the results more conclusive.", "theoretical_claims": "The paper provides clear and strong theoretical justifications that support the proposed importance-weighted score matching and variational marginal score estimation approaches. The derivations are well-structured, and the key results appear correct.", "experimental_designs_or_analyses": "The empirical soundness is limited due to the concerns previously discussed. The experiments primarily focus on Gaussian data, and while two real-world datasets are included, they may not fully capture the complexity of missing-data scenarios. Additionally, the MissDiff baseline is simplified, making the performance comparisons less conclusive. A broader evaluation with non-Gaussian data and stronger baselines would improve the robustness of the findings.", "relation_to_broader_scientific_literature": "The paper builds on prior work in score-based generative modeling and missing data imputation, particularly leveraging score matching techniques for handling incomplete data. The importance-weighted approach follows principles from importance sampling and prior work on weighted score matching, while the variational approach aligns with existing variational methods used in generative modeling. Additionally, the paper compares its methods to MissDiff, a recent diffusion-based approach for missing data, though the baseline is simplified in the experiments. While the paper is well-positioned within the literature, a more thorough discussion of how the proposed methods compare to alternative generative approaches could strengthen its positioning.", "essential_references_not_discussed": "I would suggest that the authors discuss the following related works, which explore diffusion-based approaches for missing data imputation. Their discussion could help better contextualize the contributions of the current paper.\n\n[1] Zhang, Hengrui, Liancheng Fang, and S. Yu Philip. \"Unleashing the Potential of Diffusion Models for Incomplete Data Imputation.\" CoRR (2024).\n\n[2] 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[3] 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 relevant strengths and weaknesses have been pointed out in the review.", "other_comments_or_suggestions": "I have only found one potential typo:\n- In line 263 (left column): \"Setting it at $r=10$ in our experiments.\" could be better phrased as \"We set it at $r=10$ in our experiments.\"", "questions_for_authors": "I don't have further questions at this point.", "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": "ub1vpPeu9w", "reviewer": "Reviewer_ftZ3", "summary": "The paper proposes a principled score-based method for learning jointly-specified probabilistic models in the presence of missing data in the training set. The key idea is to match the scores of marginal distributions by marginalising the missing variables first. Since marginalisation is often computationally intractable, the paper introduces two approaches based on importance sampling and variational inference. The method is validated for truncated Gaussian parameter estimation and Gaussian graphical model estimation.", "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": "Thank you for your helpful comments and all the errata found, we will make sure to correct them. We respond to specific questions below.\n\n>In Assumption 4.2 you set $\\theta>0$, is that a typo? I am not sure why the parameters should be positive.\n\nYes, sorry thank you for spotting this, it is a typo it should be for any $\\theta$ in our parameter space.\n\n>For the VI approach you parametrise the conditional variational distribution using a neural network, which is likely the reason why it performs worse than the other approaches in Figure 1. I wonder if for this simple setting you could have specified the variational distribution jointly as, fit it via and then sample the conditionals\nin equations 10/11.\n\nThis is a good point, and you're correct that for the truncated normal distribution we can indeed approximate with a normal distribution. We have used a more general neural network model in experiments as in many cases, there may not be a joint distribution from which we can easily draw conditional samples (see our ICA inspired example in our response to reviewer 2txV). We wanted to showcase our method in the most generalised way possible.\n\nThanks for raising this excellent point. We will highlight this possible variant method for the truncated normal regime in the revision.\n\n>The experimental setting focuses on fairly simple Gaussian models, and explores parameter estimation and graphical model estimation problems. While the paper could benefit from at least one experiment on a more \"interesting\" model, such as an EBM\n\nWe appreciate this feedback and have ran experiments on an non-Gaussian, ICA inspired model with an intractable likelihood of $\\left(p(x)\\propto\\exp{-\\sum_{i,j}\\theta_{ij}x_i^2x_j^2}\\right)$. The task was estimating $\\theta_{i,j}$. We ran experiments with varying sample sizes, missingness proportions, and dimensions and results are presented here: https://ibb.co/7JTJCc6b. The the accuracy is reported in parameter estimation error (Frobenius norm of $\\theta$ treating it as a matrix).\n\n\n>The first equation in Section 4, regarding model estimation from incomplete data via MLE holds not only for MCAR but also for MAR data, even though you only consider MCAR. I believe that the proposed score-based methods would similarly hold for MAR data, without the need for the additional MNAR results given in the appendices. This would make the paper stronger in general, as MCAR is a fairly strong assumption, rarely true in real scenarios.\n\nThanks for your comments! \nOur approach may not hold for the MAR case. The MLE approach holding for MAR data is specific to the MLE objective rather than the marginal framework in general. However, one can use our MNAR approach for MAR data and in the MAR case the conditional probabilities of being missing should be identifiable. \n\nOne can see that MAR adjusted score splits into the sum of the original score and a term involving the missingness pattern (i.e. for MAR data we have $s_{\\lambda}(x_\\lambda)=\\nabla\\log\\varphi_\\lambda(x_\\lambda)+\\nabla\\log p(x_\\lambda)$). However, we cannot apply  score matching since we cannot ignore the term involving the missingness probability. Thanks for pointing this out and we will add this discussion in the revision. \n\nAs a small side note, we believe the MAR setting to often be not much less restrictive that MCAR, (especially in cases where each coordinate can be missing). Specifically, any case where any two coordinates can be missing at once is likely to be either MCAR or MNAR at least according to the definition in Little \\& Rubin (2020). As an example say $X_2$ being missing depends on $X_1$ and each coordinate being missing is independent while this seems like MAR as $X_2$ being observed directly depends upon only $X_1$, for the case were $X_1,X_2$ are both missing (i.e. $\\Lambda=[d]\\setminus\\\\{1,2\\\\}$), we clearly have $\\mathbb{P}(\\Lambda=[d]\\setminus\\\\{1,2\\\\}|X)\\neq\\mathbb{P}(\\Lambda=[d]\\setminus\\\\{1,2\\\\}|X_{[d]\\setminus\\\\{1,2\\\\}})$ making it not MAR.\n\n>In appendix A.1.4 you assume that is known to have a method that is independent of how is estimated. However, it would be important to note that when the data is MNAR, and, consequently, the density of interest may not be identifiable at all, e.g. [1].\n\nThis is a good point and we should have better highlighted that learning the $\\varphi_{\\lambda}$ is difficult and often non-identifiable without further assumptions or supplementary data. Indeed learning the joint distribution with MNAR data is in general a non-identifiable task. We will highlight this fact in the appendix and also include the provided citation.\n\nWe thank the reviewer for the other citations and will include them in the revision, making sure to highlight that multi-level score matching was introduced specifically for diffusion processes.", "claims_and_evidence": "The proposed methods and their proofs seem to be sound. In fact, the paper is very well-written and the accurately-placed remarks provide well-appreciated observations and caveats.", "methods_and_evaluation_criteria": "To my knowledge, this is first method score-based method in the context of missing data that may be unbiased subject to the typical conditions for importance sampling and variational inference. The potential unbiasedness will be very appreciated by the practitioners working with missing data.", "theoretical_claims": "The claims are correct and capture the important nuances as far as I can tell. I have checked the proofs in Appendices A and D.4.1, which seem correct.\n\nOne short question:\n\n* In Assumption 4.2 you set $\\theta > 0$, is that a typo? I am not sure why the parameters should be positive.", "experimental_designs_or_analyses": "The experimental setting focuses on fairly simple Gaussian models, and explores parameter estimation and graphical model estimation problems. While the paper could benefit from at least one experiment on a more \"interesting\" model, such as an EBM, I believe the number of theoretical results and extensions to various score-matching flavors make up for it. \n\nI've got one small question still. For the VI approach you parametrise the conditional variational distribution using a neural network, which is likely the reason why it performs worse than the other approaches in Figure 1. I wonder if for this simple setting you could have specified the variational distribution _jointly_ as $p_{\\phi}'(X) = \\text{TruncNorm}(X, \\phi)$, fit it via $J_F$ and then sample the conditionals $p_{\\phi}'(X'\\_{-\\Lambda} \\mid X\\_{\\Lambda})$ in equations 10/11.", "relation_to_broader_scientific_literature": "The paper correctly compares to the two existing methods in the literature and highlights their weaknesses. The paper also proposes extensions of the proposed method to various flavors of score matching, namely, truncated, sliced, and denoising.", "essential_references_not_discussed": "I think the exposition of the missing data scenario is fairly succint but sufficient. However, for the interested readers, references to the classical literature that formalises missing data problems could be useful, e.g. [1].\n\nReferences:\n\n[1] Little and Rubin (2020). Statistical Analysis with Missing Data: Third Edition. (Section 1.3)", "other_strengths_and_weaknesses": "As I have expressed above, the paper solves an important problem, is well-written, and I believe will be found useful by practitioners dealing with missing data.", "other_comments_or_suggestions": "Suggestions:\n\n* The first equation in Section 4, regarding model estimation from incomplete data via MLE holds not only for MCAR but also for MAR data, even though you only consider MCAR. I believe that the proposed score-based methods would similarly hold for MAR data, without the need for the additional MNAR results given in the appendices. This would make the paper stronger in general, as MCAR is a fairly strong assumption, rarely true in real scenarios. \n* In appendix A.1.4 you assume that $\\varphi_{\\lambda}$ is known to have a method that is independent of how $\\varphi_{\\lambda}$ is estimated. However, it would be important to note that when the data is MNAR, $\\varphi_{\\lambda}$ and, consequently, the density of interest may not be identifiable at all, e.g. [1]. \n* In Appendix D.4 the denoising score matching is introduced from the diffusion perspective, where instead of learning a score of a single distribution, we learn a score of multiple time-dependent distributions. As such, references to [2] or [3] may be relevant here as [4], that is referenced in the text, did not explicitly cover this scenario (although this is correctly highlighted in the remark at the end of the section).\n\nTypos:\n* Line 056: \"component [of] Z\"\n* Line 109: \"missing [and] non-missing\"\n* Line 176: \"score matching, [we] can relate\"\n* Line 242: \"results\"->\"result\"\n* Line 327: \"$\\hat J_F, \\hat J_F$\" -> \"$\\hat J_{KL}, \\hat J_F$\"\n* Line 437: Last sentence of Section 6 has repeated \"score matching\"\n* Line 687: \"begin\" -> \"being\"\n* Line 720: \"show a provide\" -> \"provide\"\n* Line 727: \"with\" -> \"we\"\n* Line 1761: \"\\cdot\"\n* Line 1780: \"takes\" -> \"take\"\n* Subsection D.4.1 should probably be its own section?\n\nReferences:\n\n[1] Nabi et al (2020). Full Law Identification In Graphical Models Of Missing Data: Completeness Results\n\n[2] Song and Ermon (2020). Generative Modeling by Estimating Gradients of the Data Distribution\n\n[3] Song et al (2021). Score-Based Generative Modeling through Stochastic Differential Equations\n\n[4] Vincent (2011). A Connection Between Score Matching and Denoising Autoencoders", "questions_for_authors": "I have no 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": {}}, {"review_id": "AYX60hRfHA", "reviewer": "Reviewer_2txV", "summary": "The paper addresses the problem of parameter estimation from missing data with score matching objective. The authors introduce two general frameworks for estimating the (gradient of) marginal score and demonstrate strong empirical performance on synthetic and real 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": "Thank you for your insightful comments and feedback as well as the additional references, we will make sure to include them in our paper.\n\n>Since the goal is parameter estimation, it would be more convincing if the authors could conduct experiments on parameter recovery tasks and compare the discrepancy between the true and estimated parameters in various settings of missing data.\n\nWe agree that parameter estimation/recovery experiments would be more aligned with our goals. We link to plots on parameter estimation accuracy for the normal case here: https://ibb.co/kk8LKNV and the truncated case here: https://ibb.co/pjDzCpGQ. \nThey show largely the same trends as with the fisher divergence and we will make sure to include them in the paper. These results will be added to the main body of the paper. \n\nOur experiments did initially estimate the parameters however we just reported the performance through the Fisher divergence as this allowed us to present a singular performance metric rather than two numbers (one for the mean and the other for the precision matrix.) Parameter estimation error comparisons will be added to the main paper in the revision. \n\n> It would be convincing to see how the methods can be applied to various models or data types as the proposed framework is claim to be general.\n\nWe agree with the reviewer's comments and have ran experiments on an non-Gaussian, ICA inspired model with an intractable likelihood of $p(x)\\propto \\left(\\exp{-\\sum_{i,j}\\theta_{ij}x_i^2x_j^2}\\right)$. The task was estimating $\\theta_{i,j}$. We ran experiments with varying sample sizes, missingness proportions, and dimensions and results are presented here: https://ibb.co/7JTJCc6b. The the accuracy is reported in parameter estimation error (Frobenius norm of $\\theta$ treating it as a matrix).\n\n> So far it remains unclear to me regarding the theoretical advantages of the score matching estimates compared to MLE and other alternatives namely the Wasserstein estimates, for instance in terms of properties of the estimate, complexity or convergence rate. This is in fact important to understand how the proposed method applies in practice.\n\nThe key advantage of score matching estimates over MLE or Wasserstein distance minimizer is that we can work with unnormalisable models such as truncated Gaussians and our aforementioned ICA inspired model. Additionally, even for some normalisable models, if we cannot compute a tractable conditional probability model from joint model, we cannot apply EM (which maximizes the likelihood) easily. \n\nAs a direct theoretical comparison between score matching and MLE, [1] was able to show cases where score matching is both computationally efficient and obtains the same optimal rate of convergence as MLE, despite score matching handles unnormalizable models. \n\n## Questions\n> For Marg-IW approach, how is p' chosen in practice? \n\nIn practice kept $p'$ simple, choosing to be an isotropic gaussian with zero mean and standard deviation proportional to the standard deviation of the non-missing data.\n\n>How is the function in Corollary 4.7 chosen? Algorithm 2 only discusses how to estimate the conditional densities, without mentioning the function .\n\n$h$ is just a generic function used to define the functional $\\psi$ so $h$ takes the values given in equation (10). In other words this is just notation to try and keep equation 10 as succinct as possible\n\n>Line 706: What do the notations $B$ and $\\mathcal{B}_\\mathcal{X}$ refer to?\n\nApologies, did not properly introduce this notation. $\\mathcal{B}_{\\mathcal{X}}$ is supposed to represent all the possible events on the RV X (i.e. the sigma algebra on our covariate space $\\mathcal{X}$) so that $B$ represent any arbitrary event on $\\cal X$. We will make sure to properly introduce it in the final paper.\n\n> Line 259: Is it supposed to be $\\log(n)/n$ convergence rate? Furthermore, in the main text, please define the notations used for Theorem 4.5 such as $\\beta$, $\\delta$ for completeness.\n\nYes sorry it is meant to be $\\sqrt{\\log(n)/n}$ our aim was to highlight that this a $1/\\sqrt{n}$ rate up to log factors which diminish rapidly. Also we will make sure to introduce $\\beta$ which is a constant >0 and $\\delta$ which is any sufficiently small probability. \n\n## References\n[1] Pabbaraju, C., Rohatgi, D., Sevekari, A. P., Lee, H., Moitra, A., and Risteski, A. (2023). Provable benefits\nof score matching. In Oh, A., Naumann, T., Globerson, A., Saenko, K., Hardt, M., and Levine, S., editors,\nAdvances in neural information processing systems, volume 36, pages 61306–61326. Curran Associates,\nInc.", "claims_and_evidence": "The work provides solid theoretical results to demonstrate the applicability of the proposed frameworks. However, the presented empirical evidence is not sufficiently convincing to me. \n\nIn Section 5.1, the authors conduct an experiment on synthetic data where they only report the Fisher divergence between the true and estimate score functions. Since the goal is parameter estimation, it would be more convincing if the authors could conduct experiments on parameter recovery tasks and compare the discrepancy between the true and estimated parameters in various settings of missing data. Then we would know how effectively the score matching estimates perform compared to MLE.", "methods_and_evaluation_criteria": "The paper learns the parameters by minimizing the explicit score matching objective over the marginal score function. The objective naturally follows the derivation of the original score matching objective from complete data and achieves convergence under the standard assumptions.", "theoretical_claims": "There is no issue with the correctness of the theoretical claims that I am aware of.", "experimental_designs_or_analyses": "As discussed above, there is a lack of experiments on how well the proposed method can recover the ground-true parameters, for example the mean and/or variances in the Gaussian case. Furthermore, the current results are limited to Gaussian models. It would be convincing to see how the methods can be applied to various models or data types as the proposed framework is claim to be general.", "relation_to_broader_scientific_literature": "The paper introduces a score matching framework for estimating model parameters from missing data. At a high level, score matching aims to find parameters that minimize the Fisher divergence between the observed data distribution and the model distribution. This approach is part of the broader family of parameter estimation methods based on minimum distance estimation.", "essential_references_not_discussed": "Parameter estimation from incomplete data is a fundamental statistical problem, to which maximum likelihood estimation via the EM algorithm has been the dominant approach. MisGAN (Li et al., ICLR'19) is another method that uses GAN for estimation. A recent alternative to MLE is the Wasserstein estimate from the optimal transport framework, wherein one seeks the parameters that minimize the Wasserstein distance between the two distributions. Vo et al. (ICML'24) proposes a framework called OTP for learning parameters from missing data.\n\n*Li, S. C. X., Jiang, B., & Marlin, B. MisGAN: Learning from Incomplete Data with Generative Adversarial Networks. In International Conference on Learning Representations.* \n\n*Vo, V., Le, T., Vuong, L. T., Zhao, H., Bonilla, E. V., & Phung, D. Parameter Estimation in DAGs from Incomplete Data via Optimal Transport. In Forty-first International Conference on Machine Learning.*", "other_strengths_and_weaknesses": "**Strengths:** \n\nThe paper is generally well-written and solid. The extension of explicit score matching to dealing with missing data is a straightfoward approach, yet interesting to me. I also appreciate the authors' attempt to consider MNAR cases, though the assumption that the probability of missingness pattern is given is rather unrealistic, rendering the proposal in Proposition A.9 of little use to me. However, it is acceptable to limit the scope of the current work to M(N)AR cases since it is generally known that without making very specific assumptions, the parameters are non-identifiable from MNAR data. \n\n**Weaknesses:** \n\nSimilar to EM, the method requires that the unnormalized density $q$ is tractable for estimation and further developments to compute the marginal scores. Perhaps this is a reason why the experiments are limited to Gaussian models. \n\nAt the high level, the proposed algorithm somewhat mimics EM where the Marg-IW and Marg-Var variants can be viewed as two ways to perform inference in the E step: the former assumes the density $p'$ is known, similar to the fact the the posterior is assumed tractable; the latter is similar to doing variational inference. Meanwhile, MisGAN and OTP frameworks can sidestep the intractabilty of the density to be learned, though algorithmically, these frameworks also follow an EM-like style: an inference step for missing value imputation, followed by a step of learning parameters from the pseudo-complete data. \n\nSo far it remains unclear to me regarding the theoretical advantages of the score matching estimates compared to MLE and other alternatives namely the Wasserstein estimates, for instance in terms of properties of the estimate, complexity or convergence rate. This is in fact important to understand how the proposed method applies in practice.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "* For Marg-IW approach, how is p' chosen in practice? \n* How is the function $h$ in Corollary 4.7 chosen? Algorithm 2 only discusses how to estimate the conditional densities, without mentioning the function $h$.\n* Line 706: What do the notations $B$ and $\\mathcal{B}_{\\mathcal{X}}$ in  refer to? \n* Line 259: Is it supposed to be $\\log(n)/n$ convergence rate? Furthermore, in the main text, please define the notations used for Theorem 4.5 such as $\\beta_1$, $\\delta$ for completeness.", "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 authors consider statistical estimation with score matching, in the presence of missing data.

They focus on the case where the model of interest is a model with intractable normalising constant (the case of score-based models, where the score is directly modelled is less emphasised in the paper). They propose an objective to maximise, which is the expected value of the norm between the observed scores (the missing values are marginalised) and modelled ones ($F_M(\theta)$, defined in Eqn. (3)). When no values are missing, this reduces to the standard expected score-matching loss. They then use the standard clever tricks for the score-matching literature to obtain empirical versions of $F_M(\theta)$, both when the data is compactly supported and not. They then attack this empirical loss using two approaches to deal with the missing values: one based on importance sampling, and one based on variational inference. They perform experiments on toy datasets, Gaussian graphical models, and added new experiments on independent component analysis (ICA) during the discussion period.

All in all, this is an excellent paper, that nicely extends several blends of score matching to the missing data setting. During the discussion period, **experiments on ICA significantly strengthened the paper, and it is crucial that these experiments are included in the final version of the paper**.

Among the other things that should be clarified in the final version are the advantages/drawbacks of the proposed methods against the EM approach and other methods to estimate models with intractable normalising constants with missing values, and whether or not the work on denoising score matching could be applied to diffusion models.

Finally, two things that I wondered while reading the paper, and that were not really discussed during the discussion: what are the theoretical properties of the obtained parameters estimates (consistency, asymptotic normality...)? Could this be extended to score matching for data that lives on a manifold? These two points could be optionally discussed in the final version.
Decision: Accept (oral)