arXiv CSopen access

Generalized Score Matching for Parameter Estimation on Convex Domains

· arxiv_cs

arXiv CS · Papers · License: Open Access

machine learning, deep learning, neural networks

Generalized Score Matching for Parameter Estimation on Convex Domains Saisuchith Mahajan∗ Chandra Sekhar Seelamantula Department of Electrical Engineering Indian Institute of Science, Bengaluru 560012 {nishanths, saisuchithm, css}@iisc.ac.in

arXiv:2609.11521v1 [cs.LG] 10 Sep 2026

Nishanth Shetty∗

Abstract Maximum likelihood (ML) estimation is a principled and statistically efficient approach for learning probabilistic models. However, for unnormalized models, ML estimation requires evaluating the partition function and differentiating through it, which may not always be tractable. Score matching provides a practically viable alternative that circumvents this obstacle by fitting the score in a way that eliminates dependence on the normalizing constant. We derive the generalized score matching objective on a convex subset of Rd constructively starting from Minimum Probability Flow (MPF) learning, and show how classical score matching as well as domain-adapted variants for non-negative data arise naturally within the proposed framework. We show that the resulting objective is a proper local scoring rule of second-order, which provides the theoretical guarantee that the true density is recovered when the objective is minimized. Furthermore, for a model belonging to the exponential family, we establish convexity of the objective together with consistency of the finite-sample estimator under standard regularity conditions. Our derivation sheds new light on the scope and applicability of generalized score matching in various problem settings. We compare generalized score matching-based estimators on constrained domains, where the partition function is analytically intractable. We provide experimental results on parameter estimation for model densities belonging to the exponential family defined over convex subsets of Rd , and a generative modeling use-case to demonstrate broader applicability of the proposed generalized score matching framework.

Introduction

Learning unnormalized probabilistic models is a central challenge in modern machine learning and statistics. Several expressive model architectures, including energy-based models (EBMs) [7, 11– 13, 16, 26, 40], restricted Boltzmann machines (RBMs) [14, 15, 33, 39], and Markov random fields (MRFs) [1, 5, 10, 42], are specified by unnormalized densities. A variety of partition-functionfree estimation methods have been developed over the past several years that allow for parameter estimation without the need for explicit normalization. Score matching [18] is one such technique that circumvents the intractability of estimating the partition function by fitting the score function, defined as the gradient of the log density. Through integration by parts, the objective can be formulated solely in terms of the derivatives of the model density, independent of the normalizing constant and the score function of the target data density. While score matching was originally formulated for densities supported on Rd [18], subsequent work [19, 23, 25, 43–45] has extended the framework to distributions defined on subsets of Rd . ∗ Equal contribution.

Preprint.

A direct application of the original score matching formulation is not straightforward in such settings, as the integration by parts argument may fail in the presence of boundaries where the density or its derivatives are not continuous. An early extension was proposed by Hyvärinen [19], which considered non-negative data and introduced a modified objective designed to ensure validity of the integration by parts argument. This line of work was further generalized by Yu et al. [43], providing a broader framework for densities supported on Rd+ . Lyu [25] presents a general operator-theoretic view of score matching defining a notion of completeness which we refer to as generalized score matching in this work. While the framework is quite general, the choice of operator remains abstract. In parallel, several works have developed complementary extensions of score matching, including formulations for truncated domains [23], extension to ordinal data [41], analyses of the statistical efficiency of the resulting estimators [20, 30], etc. To the best of our knowledge, these prior works explored score matching on subsets of Rd , but a formal derivation of the objective function has not been provided. Our objective is to address this gap. 1.1

Contributions

The genesis of this paper is the following question: Could one derive generalized score matching objectives on subsets of Rd in a constructive fashion that also explains the domain-adapted variants proposed in the prior works ? We answer this question in the affirmative for convex subsets of Rd . The key contributions are are stated below. 1. We provide a constructive method to derive generalized score matching objectives on convex subsets of Rd , including the practically relevant case of Rd+ , starting from Minimum Probability Flow (MPF) [36, 37]. Unlike prior works that define score matching via the Fisher divergence, we construct the generalized score matching objective as the infinitesimal limit of MPF. 2. We show that domain-adapted variants of score matching, such as non-negative score matching [19, 43], arise naturally within our framework through appropriate choice. Our derivation shows that the linear operator [25] in this setting is complete and is tied to the convex function chosen for the domain. 3. We show that the resulting generalized score matching objective defines a proper local scoring rule of second order [29], and that its minimization recovers the true data-generating distribution. 4. For model densities in the exponential family, we prove that the generalized score matching objective is convex in the canonical parameters, and that the corresponding empirical objective yields a consistent estimator under standard regularity conditions. The theoretical developments are complemented by experimental results pertaining to parameter estimation of model densities supported on convex subsets of Rd , and a generative modeling use-case to demonstrate broader applicability of the proposed framework. The key results are stated in the main manuscript and the proofs are provided in the appendix. 1.2

Notation

Scalars are represented using normal font (e.g., x, θ), whereas vectors are represented using boldface d (e.g. x, θ), and matrices are represented qP in uppercase letters (e.g. G, H). For vectors x, y ∈ R , d d×d 2 xk denotes its kth element, ∥x∥ = i=1 xi is the standard Euclidean norm and diag(x) ∈ R constructs a diagonal matrix with diagonal entries as elements of x. The element-wise product is x ◦ y. For a matrix G ∈ Rd×d , Gij denotes the element in the ith row and jth column. diag(G) ∈ Rd extracts its diagonal. For a function f : Rd → R, denote Hf (x) as the Hessian of f evaluated at point x. For a vector field g : Rd → Rd , gk is the kth output and the divergence is given by Pd k (x) (∇·g)(x) = k=1 ∂g∂x . For a matrix field D : Rd → Rd×d , (∇.D)(x) ∈ Rd denotes the columnk Pd ∂Dij (x) wise divergence operation where (∇ · D)(x)j = i=1 ∂x for j ∈ {1, . . . , d}. For any B ⊂ Rd , i 2 ∂B denotes its boundary. For a positive-definite matrix G, ∥v∥G = v ⊤ Gv. We denote by C k the space of functions (real-, vector- or matrix-valued) whose derivatives up to order k are continuous. Throughout this paper, we assume that the distributions under consideration always admit a density function with respect to Lebesgue measure. The true density or data density is denoted as p. We use an 2

1 energy based model to denote the model density as pθ (x) = Z(θ) exp(−Eθ (x)) where Eθ : Rd → R R is the energy function with the partition function given by Z(θ) = exp(−Eθ (x))dx. Rd

Related Work

2.1

Score Matching

Score matching [18] is a parameter estimation method that matches the score function (gradient of log-density) between model and data distributions without computing partition functions. The score function of a probability density pθ (x) is defined as sθ (x) = ∇x log pθ (x) = −∇x Eθ (x) − ∇x log Z(θ) Remarkably, the gradient of the log-partition function vanishes, thus, enabling partition-free estimation. The score matching objective minimizes the expected squared distance between model and data scores: 1 J(θ) = E ∥∇x log pθ (x) − ∇x log p(x)∥2 x∼p 2 Since the data score is unknown, J(θ) can be rewritten using integration by parts (under boundary conditions) as 1 2 J(θ) = E ∥∇x log pθ (x)∥ + Tr(Hlog pθ (x)) + C x∼p 2 where C is a data-dependent constant. This formulation requires only the model score and its divergence, both computable without the partition function. In parallel, Lyu [25] advocates a complementary, operator-theoretic viewpoint: score matching can be defined by replacing the gradient with a suitable linear operator, yielding a broad class of generalized score matching objectives in which normalization constants still cancel. Building on this perspective, Qin and Risteski [30] provide a refined theoretical analysis connecting the statistical efficiency of such generalized objectives to properties of the Markov processes/diffusions naturally associated with the chosen operator. 2.2

Generalized score matching for non-negative data

Hyvärinen [19] proposed a score matching formulation for non-negative data by introducing weights that attenuate boundary effects near 0. This was further generalized by Yu et al. [43] by replacing the elementwise weight x with a more general elementwise weight h(x) = (h1 (x1 ), . . . , hd (xd ))⊤ where hj : R+ → R+ is almost surely positive. The resulting generalized h-score matching loss is " # Jh (θ) := E

x∼p

∇ log pθ (x) ◦ h(x) 2 − ∇ log p(x) ◦ h(x) 2

(1)

Choosing hj (xj ) = x2j recovers the non-negative score matching objective [19]. Then, under certain mild conditions, Yu et al. [43] show that an equivalent tractable form up to an additive constant that does not depend on score of p can be derived and is given by " d # X1 2 ′ hj (xj ) (∇ log pθ (x)j ) + hj (xj )Hlog pθ (x)jj + hj (xj )∇ log pθ (x)j Jh (θ) = E x∼p 2 j=1

Main Results

In this section, we motivate and formally state the key results to show how generalized score matching can be derived from Minimum Probability Flow (MPF) [36, 37]. We consider an energy-based model pθ , a data point x, and a perturbed state y, then the MPF objective is defined as (cf. Equation C − 1 in [37]) Z 1 K(θ) = E g(y, x) exp (Eθ (x) − Eθ (y)) dy , (2) x∼p 2 y 3

(y) 12 where g(y, x) denotes a connectivity function. Observe that exp 12 (Eθ (x) − Eθ (y)) = ppθθ (x) where we refer to the term inside the square root as the likelihood ratio. This ratio quantifies the likelihood assigned to y relative to that assigned to the data point x. To build intuition, consider the case where the connectivity function is chosen such that g(y, x) = 1 whenever y is a sufficiently small perturbation of x, and g(y, x) = 0 otherwise. In this setting, if x ∼ p is sampled from the true data distribution, one expects a well-specified model pθ to assign higher likelihood to x than to nearby perturbed states y. This suggests minimizing the likelihood ratio, averaged over data points and their perturbations, which provides an intuitive interpretation of the MPF objective. Sohl-Dickstein et al. [36, 37] consider a binary connectivity function g, defined as g(y, x) = 1 if y ∈ R(x) and g(y, x) = 0 otherwise, where R(x) denotes a neighborhood of x. In particular, Sohl-Dickstein et al. [37] show that when R(x) is chosen to be a cube centered at x with side length ε (cf. Equation C − 2 in [37]), the MPF objective recovers the score matching objective on Rd in the limit as ε → 0. We refer to this choice of connectivity function as the hard neighborhood case. Beyond this, we also consider a case where g(y, x) defines a valid conditional distribution of y given x. We refer to this setting as the soft neighborhood case. 3.1

Generalized Score Matching

We derive the generalized score matching (GSM) loss [30] from Equation 2 when the true density is supported over Rd . In particular, we consider the soft neighborhood case, where the connectivity function g(y, x) is chosen to be a valid conditional distribution q(y | x). In this case, the objective becomes Z 1 [Eθ (x) − Eθ (y)] dy K(θ) = E q(y|x)exp x∼p 2 y We choose q to be a Gaussian and show in the following that a limiting case of K(θ) yields GSM. Theorem 3.1. Let Eθ : Rd → R be in C 2 , and let D : Rd → Rd×d be a matrix-valued function in C 1 such that D(x) is symmetric and positive definite. Assume that E ∇Eθ (x)⊤ D(x)∇Eθ (x) , E ∇Eθ (x)⊤ (∇ · D)(x) , E [Tr (D(x)HEθ (x))] are fix∼p

x∼p

1 p nite for all θ. Define b(x) = (∇ · D)(x), normalization constant Z = ,ε>0 d/2 (2π) det εD(x) 2 ∥y−x− 2ε b(x)∥D(x)−1 and the Gaussian conditional density qε (y | x) := Z exp − . Then, after 2ε removing the θ-independent terms, dividing K(θ) by ε/4, taking the limit ε → 0 we obtain " # 1 ⊤ L(θ) = E ∇Eθ (x) D(x)∇Eθ (x) − ∇ · (D(x)∇Eθ (x)) x∼p 2 or equivalently, using ∇ log pθ (x) = −∇Eθ (x), we have " # 1 L(θ) = E ∇ log pθ (x)⊤ D(x)∇ log pθ (x) + ∇ · (D(x)∇ log pθ (x)) x∼p 2

(3)

(4)

Equation 4 admits an equivalent formulation in terms of the GSM loss as stated below. Proposition 3.1. Suppose the assumptions of Theorem 3.1 hold, and assume further that p ∈ C 1 and limxk →±∞ p(x)∇ log pθ (x)ℓ D(x)kℓ = 0 for all k, ℓ ∈ {1, 2, . . . , d} and θ. Then, minimizing L(θ) in Equation 4 with respect to θ is equivalent to minimizing i h 1 2 LGSM (θ) = E ∥∇ log pθ (x) − ∇ log p(x)∥D(x) (5) 2 x∼p These results were established under the assumption that the data support is Rd . We now consider the generic case where the support of p is a convex set Ω ⊂ Rd , excluding degenerate cases such as the empty or singleton set. The Gaussian conditional qε considered in Theorem 3.1 is convenient, because its moments are tractable in Rd . However, constructing analogous conditionals with tractable 4

moments on arbitrary domains is challenging. As noted in previous section, MPF recovers score matching on Rd by integrating over a small local neighborhood chosen to be a cube of side ε. Motivated by this, in this case, we consider hard-neighborhood objective. For a neighborhood R(x) ⊂ Ω of x and choosing weight w : Ω × Ω → R+ as g(y, x) in Equation 2, we have   Z 1   [Eθ (x) − Eθ (y)] dy  K(θ) = E  w(y, x)exp x∼p 2 R(x)

We show that, for appropriate choices of R(x) and w, the small-neighborhood limit recovers the GSM objective on Ω. We begin by describing the construction of the neighborhood R(x). To this end, let ϕ : Ω → R be a strictly convex function. The Bregman divergence is defined as Dϕ (y, x) = ϕ(y) − ϕ(x) − ⟨∇ϕ(x), y − x⟩

(6)

Define the local Bregman ball around x to be Crϕ (x) = {y : (y − x)⊤ Hϕ (x)(y − x) ≤ r(x)2 }

(7)

where r : Ω → (0, ∞). Crϕ (x) is an ellipsoid centered at x, defined by the quadratic form induced by Hϕ (x). In this case, the objective with the local Bregman ball is given by Z 1 Krϕ (θ) = E Irϕ (x; θ) where Irϕ (x; θ) = w(y, x) exp [Eθ (x) − Eθ (y)] dy (8) x∼p 2 Crϕ (x) We show that generalized score matching for convex sets can be derived from the setup we introduced considering an appropriately chosen weighting function. Theorem 3.2. Let ϕ be in C 3 and strictly convex on an open convex set Ω ⊂ Rd , and assume that Hϕ (x) ≻ 0 for all x ∈ Ω. Let Eθ : Ω → R be C 2 in x and r : Ω → (0, ∞) be such that Crϕ (x) ⊂ Ω for all x ∈ Ω. Define x̄ = (x + y)/2 and consider the weighting function (y − x)⊤ Hϕ (x̄) (y − x) w(y, x) = exp − 2r(x)2 in the definition of Irϕ (x; θ) in Equation 8. Define

and

assume

that

Gϕ (x) = (det Hϕ (x))−1/2 Hϕ (x)−1 E ∇Eθ (x)⊤ Gϕ (x)∇Eθ (x) , E ∇Eθ (x)⊤ (∇ · Gϕ )(x)

x∼p

(9) and

E [Tr (Gϕ (x)HEθ (x))] are finite for all θ. Then, after removing the θ-independent terms

x∼p

from Irϕ (x; θ), dividing by r(x)d+2 /4, and taking the limit r → 0, we obtain for a positive constant λ 1 Lϕ (θ) = E ∇Eθ (x)⊤ Gϕ (x)∇Eθ (x) − Tr (Gϕ (x)HEθ (x)) − λ∇Eθ (x)⊤ (∇ · Gϕ )(x) x∼p 2 (10) Following Remark D.3, we use the alternative weighting function that yields λ = 1, and define the resulting objective as Lϕ (θ) = E [S ϕ (x; θ)] x∼p

(11)

where S ϕ (x; θ) is redefined as follows: 1 ∇Eθ (x)⊤ Gϕ (x)∇Eθ (x) − ∇.(Gϕ (x)∇Eθ (x)) (12) 2 Equation 11 admits an equivalent formulation in terms of GSM, as stated in the following proposition. Care must be taken in specifying the boundary conditions. We follow the methodology of Liu et al. [23], in particular, Theorem 2, which assumes that the underlying domain has a Lipschitz boundary. Since bounded convex sets admit Lipschitz boundaries (Lemma 1.13 in Chapter 2 of Simon [35]), the proposition below considers bounded convex sets. S ϕ (x; θ) =

Proposition 3.2. Suppose the assumptions of Theorem 3.2 hold, and assume further that Ω is bounded, p ∈ C 1 and for any z ∈ ∂Ω, we have lim p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ nk (z) = 0

x→z

∀ k, ℓ ∈ {1, 2, . . . , d}, ∀ θ

(13)

where x → z takes any sequence in Ω converging to z and (n1 , . . . , nd ) is the unit outward normal vector on ∂Ω. Then minimizing Lϕ (θ) in Equation 11 with respect to θ is equivalent to minimizing h i 1 2 LϕGSM (θ) = E ∥∇ log pθ (x) − ∇ log p(x)∥Gϕ (x) (14) 2 x∼p We address the case of unbounded sets in Appendix E. An immediate consequence of Proposition 3.2 is that the score matching objective in Equation 14 corresponds to the linear operator L defined as 1

(Lg)(x) = Gϕ (x) 2 ∇g(x) 2

(15) 1

in the framework of Lyu [25]. This holds because ∥v∥G = v ⊤ Gv = v ⊤ G 2 G 2 v = G 2 v for a positive definite matrix G. The following proposition shows that the operator L is complete. Proposition 3.3. Let L be the operator defined in Equation 15, then L is complete, that is, for 1 )(x) 2 )(x) densities p1 (x) and p2 (x) with support Ω, suppose (Lp = (Lp p1 (x) p2 (x) almost everywhere, then p1 (x) = p2 (x) almost everywhere. Remark 3.1. It immediately follows from Proposition 3.3 that if LϕGSM (θ) = 0, then pθ = p a.e. By choosing appropriate convex functions ϕ, we retrieve score matching objectives in different settings. 2

Corollary 3.3. For ϕ(x) = 12 ∥x∥ defined on Ω = Rd , we get Hyvärinen’s score matching objective h i 1 2 LϕGSM (θ) = E ∥∇ log pθ (x) − ∇ log p(x)∥ 2 x∼p The proof follows directly from Proposition 3.2 since Hϕ (x) = I =⇒ Gϕ (x) = I. Pd d Remark 3.2. For ϕ(x) = − Q i=1 log xi defined on Ω = R+ , we have Gϕ (x) ⊤ d x21 , . . . , x2d and the objective function turns out to be i=1 xi diag " d # Y 2 1 ϕ xi ∇ log pθ (x) ◦ x − ∇ log p(x) ◦ x LGSM (θ) = E 2 x∼p i=1

This objective is similar to the non-negative score matching objective proposed in [19], but with Qd the extra term i=1 xi , which is due to the presence of the determinant in the definition of Gϕ (x). Pd Extending further, we can work with a strictly convex function ϕ(x) = i=1 ϕi (xi ) where ϕi : R+ → R is strictly convex. In this case, ⊤ ! 1 1 1 , . . . , ′′ Gϕ (x) = Qd p diag ϕ′′1 (x1 ) ϕd (xd ) ϕ′′i (xi ) i=1 and the objective in Equation 14 " # 1 1 2 1 1 ϕ LGSM (θ) = E Qd p ∇ log pθ (x) ◦ h(x) 2 − ∇ log p(x) ◦ h(x) 2 2 x∼p ϕ′′i (xi ) i=1 h i⊤ 1 1 where h(x) = ϕ′′ (x , . . . , . This is similar to what Yu et al. [43] consider but with the ′′ ϕd (xd ) 1) i 1 extra factor Qd √ ′′ . i=1

ϕi (xi )

Remark 3.3. Existing works, including [19, 23, 30, 41, 43], take generalized score matching objectives as the starting point. In contrast, through Theorem 3.1, Proposition 3.1, Theorem 3.2, and Proposition 3.2, we provide a principled derivation of these objectives from a unified formulation. To the best of our knowledge, this is a novel development and constitutes one of the main contributions of our paper. 6

Remark 3.4. The objective Equation 2 considered in the manuscript is motivated by its connection to Minimum Probability Flow (MPF) and serves as the primary starting point of our derivation. As shown in Appendix G, the same proof strategy extends to a broader class of objectives, yielding corresponding generalizations of the main results. This highlights the flexibility of derivation beyond the specific MPF formulation.

Generalized Score Matching is a Proper Scoring Rule

Given a set of probability distributions P with each element having support X , a scoring rule [29] is defined as the loss S(x, Q) incurred when a sample x ∼ P ∈ P is realised and the model distribution choice was Q ∈ P. The expectation of S(x, Q) denoted by S(P, Q) is given by S(P, Q) = Ex∼P [S(x, Q)]. A scoring rule is considered proper if S(P, Q) ≥ S(P, P ) for all P, Q ∈ P, and strictly proper if the inequality is strict whenever Q ̸= P . For a simply connected domain X ⊂ Rn and a twice-differentiable strictly positive model density q(x) on X , Parry [28] showed that all previously known multidimensional scoring rules can be generated by n X 1 −1 Gij (x) qi (x) qj (x) ϕ[q](x) = − q(x) 2 i,j=1

where ϕ[y] := ϕ(x1 , . . . , xn , y, y1 , . . . , yn ) is differentiable in x, and twice differentiable, jointly ∂q (strictly) concave and 1-homogeneous in (y, y1 , . . . , yn ), where qi := ∂x , G(x) = [Gij (x)] is a i symmetric positive definite matrix. The associated (strictly) proper local scoring rule generated is ! ! n X ∂Gij (x) qj (x) qij (x) 1 qi (x)qj (x) + S(x, Q) = Gij (x) − (16) q(x) 2 q(x)2 ∂xi q(x) i,j=1 2

q where qij := ∂x∂i ∂x . This is a fairly general framework and encompasses several previously known j loss functions such as maximum likelihood and score matching. We discuss the relevance of proper scoring rules of the second order in Appendix H. Additionally, the following proposition formally establishes that the proposed objective in Equation 12 constitutes a proper scoring rule. Proposition 4.1. S ϕ (x; θ) as defined in Equation 12 is a proper scoring rule with G(x) = Gϕ (x).

Analysis for Exponential Family

For a model density pθ from the exponential family, log pθ (x) = θ ⊤ t(x) − ψ(θ) + b(x)

(17)

where θ ∈ Θ ⊂ Rr , t : Rd → Rr represents the sufficient statistics, ψ(θ) is the normalizing constant, and b(x) is the base measure with t and b being almost surely differentiable. Let {xi }N i=1 drawn i.i.d. from the density p and the finite sample version of the objective in Equation 4 is given by " # N 1 X 1 ⊤ ∇ log pθ (xi ) D(xi )∇ log pθ (xi ) + ∇ · (D(xi )∇ log pθ (xi )) L̂(θ) = N i=1 2 Proposition 5.1. L̂(θ) can be expressed as quadratic 1 ⊤ ⊤ θ ΓN θ + g N θ+C (18) 2 where C ∈ R is a constant independent of θ, Jt (x) denotes the Jacobian of t evaluated at x, and L̂(θ) =

1 X ΓN = Jt (xi )D(xi )Jt (xi )⊤ N i=1 gN =

N i 1 Xh E(xi ) + Jt (xi )(∇ · D(xi )) + Jt (xi )D(xi )∇b(xi ) N i=1

where, in turn, E(x) ∈ Rr is a vector with entries d X

El (x) =

j,k=1

Djk (x)

∂ 2 tl (x) , ∂xj ∂xk

l = 1, 2, . . . , r.

Assuming that ΓN is positive semidefinite almost surely, L̂(θ) is a convex function of θ. We now show that estimator that minimizes L̂(θ) is consistent. Theorem 5.1. Let θ0 = arg minθ L(θ) be the true parameter that minimizes the objective in Equation 4. Define Γ0 = E[Γ1 ], g0 = E[g1 ], Σ0 = E[(Γ1 θ0 + g1 )(Γ1 θ0 + g1 )⊤ ] Further, assume that ΓN is a.s. positive definite, Γ0 , Γ−1 0 , g0 and Σ0 exist and are entry-wise finite. Then, the unconstrained minimizer of L̂(θ) is a.s. unique with closed-form solution θ̂N = −Γ−1 N gN . Moreover, the estimator is consistent and asymptotically normal, i.e., √ d a.s. −1 θ̂N −−→ θ0 and N (θ̂N − θ0 ) − → N (0, Γ−1 0 Σ0 Γ0 ) as N → ∞

Experiments

In this section, we evaluate the proposed generalized score matching estimators against existing methods for parameter estimation on constrained domains. Specifically, we focus on distributions supported on positive orthant Rd+ , the (d − 1)-simplex, and the standard simplex polytope S d = {x ∈ Pd Rd : xi > 0, i=1 xi < 1}, where MLE is intractable. We also consider application to generative modeling which pertains to training an implicit VAE [38] on MNIST [21] and CelebA [24] using the proposed GSM loss, demonstrating that the framework extends beyond parameter estimation. 6.1

Truncated Gaussian Model

We consider a truncated Gaussian density of the form 1 2 pµ,K (x) ∝ exp − ∥x − µ∥K 1Ω (x), 2 where 1Ω denotes the indicator function restricting the support to Ω, µ ∈ Rd is the location parameter, and K ∈ Rd×d is the symmetric positive definite precision matrix. We examine two choices of constrained support: the simplex polytope Ω = S d and the positive orthant Ω = Rd+ . For each domain, experiments are evaluated across multiple sample sizes N ∈ {200, 500, 800}, with performance aggregated over 50 independent trials. For S 10 , we compare our proposed estimators against the baselines of Truncated Score Matching [23] and Yu et al. [43]. We denote our choices of ϕ as P 4/3 P P P P 4/3 ϕ1 (x) = 49 (P i xi + (1 − i xP ), ϕ2 (x) = i xi log xi + (1 − i xi ) log(1 − i xi ), and i) ϕ3 (x) = − i log xi − log(1 − i xi ), alongside baseline choice h1 (x) = x. Figure 1 illustrates Mean Squared Error (MSE) for µ and K across sample sizes, and Table 1 reports quantitative MSE at N = 800. We observe that ϕ1 achieves the lowest median MSE across all sample sizes N for both parameters, yielding a median MSE of 0.0072 for µ at N = 800, compared to 0.0227 for ϕ2 , 0.0828 for ϕ3 , 0.0829 for h1 , and 0.1050 for Truncated SM. For precision matrix estimation, the MSE for h1 remains around 105 –106 across all N , roughly an order of magnitude higher than all other evaluated estimators. Furthermore, while ϕ3 , Truncated SM, and h1 exhibit extreme estimation error outliers extending up to 102 –104 the error spread for ϕ1 contracts consistently as sample size increases. Experimental results for the positive orthant (Ω = Rd+ ) are deferred to Section K.2. 6.2

Discussion and Practical Considerations

The boundary conditions in Proposition 3.2 (Equation 13) naturally motivate choosing a generator that vanishes at the boundary, i.e., Gϕ (x) → 0 as x → ∂Ω. While this requirement mirrors the boundary attenuation ideas introduced by Hyvärinen [19] and Yu et al. [43] for non-negative data, our framework induces these decay dynamics intrinsically from the Hessian of a chosen convex ϕ 8

Estimation Error for µ vs Sample Size

Estimation Error for K vs Sample Size

104

106

102

MSE (log scale)

103 101

105

100

10 1 10 2

104

10 3

N = 200

N = 500

N = 800

Number of Samples

φ1 (x) = 94 φ2 (x) =

µX

xi4/3 + 1 − µ

X ¶4 ¶

xi log(xi ) + 1 −

N = 200

xi 3

X ¶

xi log 1 −

φ3 (x) = −

X ¶

Truncated SM

N = 500 Number of Samples µ

log(xi ) − log 1 −

X ¶

N = 800

h1 (x) = x

Figure 1: Comparison of parameter estimation error (MSE) for µ and K on the simplex polytope S 10 across sample sizes N ∈ {200, 500, 800}, evaluated over 50 independent trials. Baselines include Truncated Score Matching (Truncated SM) from Liu et al. [23] and h(x) = x from Yu et al. [43]. via the generator Gϕ (x). A general convex polytope Ω = {x ∈ Rd : a⊤ k x < bk , 1 ≤ k ≤ m} encapsulates our primary experiments (Section 6.1,K.2,K.3), where the distance to each facet is given by the affine slack sk (x) = bk − a⊤ we instantiated this k x. In our empirical evaluations, Pm 9 4/3 geometry through three barrier choices: the power barrier ϕ (x) = , entropic 1 k=1P 4 sk (x) Pm m barrier ϕ2 (x) = k=1 sk (x) log sk (x), and logarithmic barrier ϕ3 (x) = − k=1 log sk (x). The exponent 4/3 in ϕ1 is chosen because it recovers the MLE for a 1D exponential (cf. Section K.1). Notice that all three choices have attenuation behavior near the boundary. On bounded domains such as S d and ∆d−1 , estimators from Yu et al. [43] yield degraded performance because their underlying generator does not attenuate near boundary. More importantly, our results show that different attenuation rates lead to different empirical performances. While ϕ1 , ϕ2 , and ϕ3 all attenuate at the boundary, their specific decay rates differ, which suggests that the rate of boundary attenuation is a factor in finite sample estimator performance. While these empirical findings help rule out ineffective choices of ϕ, a principled theoretical procedure for selecting the optimal ϕ remains an open problem.

Conclusions and Outlook

Starting from the MPF objective we showed that, under small perturbations and an appropriately defined neighborhood, a limiting analysis gives rise to the generalized score matching objective for convex subsets of Rd . Classical score matching and variants, such as non-negative score matching, arise as special cases within this framework. In particular, extending the analysis to unbounded convex subsets requires a different proof strategy from the bounded setting. Overall, the proposed formalism provides a unified treatment of several existing score matching formulations. We proved that the generalized score matching objective defines a proper local scoring rule of second order. For densities from the exponential family, we proved that the generalized score matching objective is convex in the canonical parameters, and that the corresponding empirical objective yields a consistent estimator under standard regularity conditions. A limitation of the proposed framework is that computing the generator Gϕ is computationally expensive in very high dimensions, which places a practical constraint on the choice of ϕ (cf. Remark 3.2). An important direction for future work is to develop a systematic approach for choosing ϕ. Characterizing the MSE-optimal choice of ϕ, studying the existence of ϕ that recovers the maximum likelihood estimate, establishing non-asymptotic guarantees for the corresponding estimators, etc. are all interesting directions for future work.

References [1] Julian Besag. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society: Series B (Methodological), 36(2):192–225, 1974. doi: https: //doi.org/10.1111/j.2517-6161.1974.tb00999.x. URL https://rss.onlinelibrary.wiley. com/doi/abs/10.1111/j.2517-6161.1974.tb00999.x. 1 [2] Patrick Billingsley. Convergence of probability measures. Wiley Series in Probability and Statistics: Probability and Statistics. John Wiley & Sons Inc., New York, second edition, 1999. ISBN 0-471-19745-9. A Wiley-Interscience Publication. J [3] Christopher M. Bishop. Pattern Recognition and Machine Learning (Information Science and Statistics). Springer, 1 edition, 2007. ISBN 0387310738. D.2 [4] Valentin De Bortoli, Alexandre Galashov, J Swaroop Guntupalli, Guangyao Zhou, Kevin Patrick Murphy, Arthur Gretton, and Arnaud Doucet. Distributional diffusion models with scoring rules. In Forty-second International Conference on Machine Learning, 2025. URL https: //openreview.net/forum?id=N82967FcVK. H.1 [5] Y. Boykov, O. Veksler, and R. Zabih. Fast approximate energy minimization via graph cuts. IEEE Transactions on Pattern Analysis and Machine Intelligence, 23(11):1222–1239, 2001. doi: 10.1109/34.969114. 1 [6] Ciwan Ceylan and Michael U Gutmann. Conditional noise-contrastive estimation of unnormalised models. ArXiv, abs/1806.03664, 2018. URL https://api.semanticscholar.org/ CorpusID:47019615. A.2 [7] Yilun Du and Igor Mordatch. Implicit generation and modeling with energy based models. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. URL https://proceedings.neurips.cc/paper_files/paper/2019/file/ 378a063b8fdb1db941e34f4bde584c7d-Paper.pdf. 1 [8] Werner Ehm and Tilmann Gneiting. Local proper scoring rules of order two. The Annals of Statistics, 40(1):609 – 637, 2012. doi: 10.1214/12-AOS973. URL https://doi.org/10. 1214/12-AOS973. H.1 [9] Ethan N. Epperly, Joel A. Tropp, and Robert J. Webber. XTrace: Making the most of every sample in stochastic trace estimation. SIAM Journal on Matrix Analysis and Applications, 45 (1):1–23, 2024. doi: 10.1137/23M1548323. K.6 [10] Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6 (6):721–741, 1984. doi: 10.1109/TPAMI.1984.4767596. 1 [11] Will Grathwohl, Kuan-Chieh Wang, Joern-Henrik Jacobsen, David Duvenaud, Mohammad Norouzi, and Kevin Swersky. Your classifier is secretly an energy based model and you should treat it like one. In International Conference on Learning Representations, 2020. URL https://openreview.net/forum?id=Hkxzx0NtDB. 1 [12] Michael Gutmann and Aapo Hyvärinen. Noise-contrastive estimation: A new estimation principle for unnormalized statistical models. In Yee Whye Teh and Mike Titterington, editors, Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, volume 9 of Proceedings of Machine Learning Research, pages 297–304, Chia Laguna Resort, Sardinia, Italy, 13–15 May 2010. PMLR. URL https://proceedings.mlr.press/v9/ gutmann10a.html. A.2 [13] Geoffrey E. Hinton. Training products of experts by minimizing contrastive divergence. Neural Computation, 14(8):1771–1800, 2002. doi: 10.1162/089976602760128018. 1 [14] Geoffrey E. Hinton. A Practical Guide to Training Restricted Boltzmann Machines. Springer Berlin Heidelberg, Berlin, Heidelberg, 2012. ISBN 978-3-642-35289-8. doi: 10.1007/ 978-3-642-35289-8_32. URL https://doi.org/10.1007/978-3-642-35289-8_32. 1 10

[15] Geoffrey E. Hinton, Simon Osindero, and Yee-Whye Teh. A fast learning algorithm for deep belief nets. Neural Computation, 18(7):1527–1554, 2006. doi: 10.1162/neco.2006.18.7.1527. 1 [16] J J Hopfield. Neural networks and physical systems with emergent collective computational abilities. Proceedings of the National Academy of Sciences, 79(8):2554–2558, 1982. doi: 10. 1073/pnas.79.8.2554. URL https://www.pnas.org/doi/abs/10.1073/pnas.79.8.2554. 1 [17] M.F. Hutchinson. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in Statistics - Simulation and Computation, 19(2): 433–450, 1990. doi: 10.1080/03610919008812866. URL https://doi.org/10.1080/ 03610919008812866. K.6 [18] A. Hyvärinen. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(24), 2005. URL http://jmlr.org/papers/v6/ hyvarinen05a.html. 1, 2.1, K.1 [19] Aapo Hyvärinen. Some extensions of score matching. Computational Statistics & Data Analysis, 51(5):2499–2512, 2007. ISSN 0167-9473. doi: https://doi.org/10.1016/j.csda.2006.09.003. URL https://www.sciencedirect.com/science/article/pii/S0167947306003264. 1, 2, 2.2, 2.2, 3.2, 3.3, 6.2, K.1 [20] Frederic Koehler, Alexander Heckett, and Andrej Risteski. Statistical efficiency of score matching: The view from isoperimetry. In The Eleventh International Conference on Learning Representations, 2023. URL https://openreview.net/forum?id=TD7AnQjNzR6. 1 [21] Yann LeCun, Corinna Cortes, and CJ Burges. Mnist handwritten digit database. ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist, 2, 2010. 6 [22] Sebastian Lerch, Thordis L. Thorarinsdottir, Francesco Ravazzolo, and Tilmann Gneiting. Forecaster’s dilemma: Extreme events and forecast evaluation. Statistical Science, 32(1):106– 127, 2017. ISSN 08834237, 21688745. URL http://www.jstor.org/stable/26408123. H.1 [23] Song Liu, Takafumi Kanamori, and Daniel J. Williams. Estimating density models with truncation boundaries using score matching. J. Mach. Learn. Res., 23(1), January 2022. ISSN 1532-4435. 1, 3.1, 3.3, 6.1, 1, K.2, 1, K.3, 4, 3, 4, 2, 6, K.4.3, K.4.3 [24] Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang. Deep learning face attributes in the wild. In Proceedings of International Conference on Computer Vision (ICCV), December 2015. 6 [25] S. Lyu. Interpretation and generalization of score matching. In Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial Intelligence, 2009. ISBN 9780974903958. 1, 2, 2.1, 3.1 [26] Erik Nijkamp, Mitch Hill, Tian Han, Song-Chun Zhu, and Ying Nian Wu. On the anatomy of mcmc-based maximum likelihood learning of energy-based models. Proceedings of the AAAI Conference on Artificial Intelligence, 34(04):5272–5280, Apr. 2020. doi: 10.1609/aaai.v34i04. 5973. URL https://ojs.aaai.org/index.php/AAAI/article/view/5973. 1 [27] Tianyu Pang, Kun Xu, Chongxuan Li, Yang Song, Stefano Ermon, and Jun Zhu. Efficient learning of generative models via finite-difference score matching. In Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS ’20, Red Hook, NY, USA, 2020. Curran Associates Inc. ISBN 9781713829546. K.6 [28] Matthew Parry. Extensive scoring rules. Electronic Journal of Statistics, 10(1):1098 – 1108, 2016. doi: 10.1214/16-EJS1132. URL https://doi.org/10.1214/16-EJS1132. 4 [29] Matthew Parry, A. Philip Dawid, and Steffen Lauritzen. Proper local scoring rules. The Annals of Statistics, 40(1):561 – 592, 2012. doi: 10.1214/12-AOS971. URL https://doi.org/10. 1214/12-AOS971. 3, 4, H.1 11

[30] Yilong Qin and Andrej Risteski. Fit like you sample: Sample-efficient generalized score matching from fast mixing diffusions. In Shipra Agrawal and Aaron Roth, editors, Proceedings of Thirty Seventh Conference on Learning Theory, volume 247 of Proceedings of Machine Learning Research, pages 4413–4457. PMLR, 30 Jun–03 Jul 2024. URL https://proceedings.mlr.press/v247/qin24a.html. 1, 2.1, 3.1, 3.3 [31] Jongha Jon Ryu, Abhin Shah, and Gregory W. Wornell. A unified view on learning unnormalized distributions via noise-contrastive estimation. In Forty-second International Conference on Machine Learning, 2025. URL https://openreview.net/forum?id=Wwj6jjxZet. A.2, A.2, G.1 [32] Tobias Schröder, Zijing Ou, Jen Lim, Yingzhen Li, Sebastian Vollmer, and Andrew Duncan. Energy discrepancies: A score-independent loss for energy-based models. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 45300–45338. Curran Associates, Inc., 2023. A.3 [33] David Sherrington and Scott Kirkpatrick. Solvable model of a spin-glass. Phys. Rev. Lett., 35: 1792–1796, Dec 1975. doi: 10.1103/PhysRevLett.35.1792. URL https://link.aps.org/ doi/10.1103/PhysRevLett.35.1792. 1 [34] Nishanth Shetty and Chandra Sekhar Seelamantula. Monte carlo score matching for image generation. In ICASSP 2025 - 2025 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 1–5, 2025. doi: 10.1109/ICASSP49660.2025.10889041. K.6 [35] Leon Simon. Introduction to geometric measure theory. Tsinghua Lectures, 2(2):3–1, 2014. URL https://math.stanford.edu/~lms/ntu-gmt-text.pdf. 3.1 [36] Jascha Sohl-Dickstein, Peter Battaglino, and Michael R. DeWeese. Minimum probability flow learning. In Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML’11, page 905–912, Madison, WI, USA, 2011. Omnipress. ISBN 9781450306195. 1, 3, 3, A.1 [37] Jascha Sohl-Dickstein, Peter B. Battaglino, and Michael R. DeWeese. New method for parameter estimation in probabilistic models: Minimum probability flow. Phys. Rev. Lett., 107:220601, Nov 2011. doi: 10.1103/PhysRevLett.107.220601. URL https://link.aps.org/doi/10. 1103/PhysRevLett.107.220601. 1, 3, 3, A.1 [38] Y. Song, S. Garg, J. Shi, and S. Ermon. Sliced score matching: A scalable approach to density and score estimation. In Proceedings of the Thirty-Fifth Conference on Uncertainty in Artificial Intelligence, UAI, 2019. URL http://auai.org/uai2019/proceedings/papers/204. pdf. 6, K.6, K.6, K.3, 7, 8 [39] Tijmen Tieleman. Training restricted boltzmann machines using approximations to the likelihood gradient. In Proceedings of the 25th International Conference on Machine Learning, ICML ’08, page 1064–1071, New York, NY, USA, 2008. Association for Computing Machinery. ISBN 9781605582054. doi: 10.1145/1390156.1390290. URL https: //doi.org/10.1145/1390156.1390290. 1 [40] Max Welling and Yee Whye Teh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML’11, page 681–688, Madison, WI, USA, 2011. Omnipress. ISBN 9781450306195. 1 [41] Jiazhen Xu, Janice L. Scealy, Andrew T.A. Wood, and Tao Zou. Generalized score matching. Journal of Multivariate Analysis, 210:105473, 2025. ISSN 0047-259X. doi: https://doi.org/10. 1016/j.jmva.2025.105473. URL https://www.sciencedirect.com/science/article/ pii/S0047259X25000685. 1, 3.3 [42] Jonathan S. Yedidia, William T. Freeman, and Yair Weiss. Understanding belief propagation and its generalizations, page 239–269. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 2003. ISBN 1558608117. 1 12

[43] Shiqing Yu, Mathias Drton, and Ali Shojaie. Generalized score matching for non-negative data. Journal of Machine Learning Research, 20(76):1–70, 2019. URL http://jmlr.org/ papers/v20/18-278.html. 1, 2, 2.2, 2.2, 3.2, 3.3, 6.1, 6.2, 1, E.1, K.1, K.2, 1, 2, 3, K.3, 4, 3, 4, 5, 2, 6, 5, K.4.2 [44] Shiqing Yu, Mathias Drton, and Ali Shojaie. Generalized score matching for general domains. Information and Inference: A Journal of the IMA, 11(2):739–780, 01 2021. ISSN 2049-8772. doi: 10.1093/imaiai/iaaa041. URL https://doi.org/10.1093/imaiai/iaaa041. [45] Shiqing Yu, Mathias Drton, and Ali Shojaie. Interaction models and generalized score matching for compositional data. In Soledad Villar and Benjamin Chamberlain, editors, Proceedings of the Second Learning on Graphs Conference, volume 231 of Proceedings of Machine Learning Research, pages 20:1–20:25. PMLR, 27–30 Nov 2024. URL https://proceedings.mlr. press/v231/yu24a.html. 1, K.3, 4, 3, 2, 6

Connections to Related Objectives

In this section, we discuss the connections between the initial objective in Equation 2 and related formulations that have appeared in the literature. Since we started with Minimum Probability Flow (MPF), we begin with a brief overview of MPF and discuss its connection to other related objectives. A.1

Minimum Probability Flow Learning

Minimum Probability Flow (MPF) [36, 37] is a parameter estimation framework that avoids the computation of intractable partition function. MPF introduced the continuous time Markov process with dynamics that transport probability mass between states. Given transition rates Γ(y, x) and the connectivity function g(y, x) between states y and x, the probability flow from state x to state y is 1 Γθ (y, x) = g(y, x) exp [Eθ (x) − Eθ (y)] 2 The MPF objective function measures the expected probability flow from the empirical data distribution p to non-data states over infinitesimal time ϵ. MPF starts from a KL divergence; a first-order Taylor expansion yields the objective in Equation 2 (Equation C − 1 from Sohl-Dickstein et al. [37]). Intuitively, we are interested the find parameter θ that minimizes the probability flow from data states to non data states. The consistency of the MPF estimator under suitable regularity conditions was established by Sohl-Dickstein et al. [36], to which we refer the reader for a detailed treatment. A.2

Conditional Noise Contrastive Estimation

The central idea of Noise Contrastive Estimation (NCE) [12] is to learn a classifier that distinguishes samples from the data distribution p from samples of noise distribution pn . It is known that the noise distribution pn must be carefully chosen to guarantee good convergence of the resulting estimator, generally considered hard in practice. To address this limitation, Ceylan and Gutmann [6] introduced conditional NCE (CondNCE) where noisy samples are generated conditionally on the observed data samples. This framework was further generalized by Ryu et al. [31] through the introduction of the f −CondNCE framework, based on general convex function f . We consider the objective proposed by Ryu et al. [31] (cf. Equation 4) p(x)q(y | x) pθ (x)q(y | x) p(y)q(x | y) Hf (θ) = y∼p, E Df , − x∼p, E f p(y)q(x | y) pθ (y)q(x | y) p(x)q(y | x) x∼q(·|y)

x∼p, y∼q(·|x)

y∼q(·|x)

′

[−f (ρθ (x, y)) + ρθ (y, x)f (ρθ (y, x)) − f (ρθ (y, x))]

(19)

where f : R≥0 → R is strictly convex function, Df denotes the Bregman divergence defined as in Equation 6 and ρθ (x, y) = ppθθ (x)q(y|x) (y)q(x|y) . Following Ryu et al. [31], we consider symmetric channel q(y | x) = q(x | y), then the objective reduces to pθ (x) pθ (y) ′ pθ (y) pθ (y) HfS (θ) = x∼p, E −f ′ + f −f (20) pθ (y) pθ (x) pθ (x) pθ (x) y∼q(·|x)

√ Consider the above objective for f (x) = − x, then we have "s # pθ (y) S Hf (θ) = x∼p, E pθ (x) y∼q(·|x)

which coincides with Equation 2 when the connectivity function g(y, x) is chosen as conditional density q(y | x). A.3

Energy Discrepancy

Schröder et al. [32] proposed Energy Discrepancy (ED), a framework that learns the energy function directly without relying on score derivatives. ED constructs a loss based on the energy difference 14

between data points and perturbed samples. Formally, let x be a data point and y be a perturbed sample sampled from conditional density q(· | x). The contrastive potential Eq (y) is defined as (cf. Equation 4 in [32]) Z Eq (y) = − log q(y | x) exp(−Eθ (x))dx x

The Energy Discrepancy objective is then defined as the expected difference between the model energy at the data point and the contrastive potential at the perturbed point   1 EDq (p, Eθ ) = E [Eθ (x)] − x∼p, E [Eq (y)] 2 x∼p y∼q(·|x)

We propose a simplifying modeling assumption by replacing the energy function corresponding to the contrastive potential induced by q, i.e., Eq (y), with the model energy, Eθ (y). While this substitution is not exact, it is intuitively justified when the perturbation density q(y | x) is sufficiently localized and Eθ varies smoothly. In this regime, the log-sum-exp integral defining Eq (y) is dominated by contributions from x in the immediate vicinity of y. Consequently, Eq (y) can be locally approximated by Eθ (y), i.e., Eq (y) ≈ Eθ (y). We define 1 ˆ EDq (p, Eθ ) = x∼p, E Eθ (x) − Eθ (y) 2 y∼q(·|x)

By invoking the convexity of the exponential function and Jensen’s inequality, we obtain an upper ˆ q , that is bound on ED 1 ˆ exp EDq (p, Eθ ) ≤ x∼p, E exp [Eθ (x) − Eθ (y)] 2 y∼q(·|x)

Rather than minimizing the left-hand side directly, we consider minimizing this upper bound as a surrogate objective. Observe that this surrogate coincides with Equation 2 when the connectivity function g(y, x) is chosen as the conditional density q(y | x), corresponding to the soft neighborhood case. A factor of 12 was introduced in the definition of energy discrepancy to show the connection to MPF objective explicit.

Proof of Theorem 3.1

For ease of reading, we first present a sketch of the proof, followed by the complete proof. Proof Sketch. The proof proceeds via Taylor expansion of the exponential term of the integrand of Equation 21. The key steps are as follows: (i) we define δ = y − x, which follows δ ∼ N (εb(x), 2εD(x)); (ii) We expand exp 21 [Eθ (x) − Eθ (y)] to second order in δ using Taylor expansion, which introduces terms involving ∇Eθ (x) and the Hessian HEθ (x); (iii) We evaluate Iε (x; θ) by computing the Gaussian moments, where E[δ] = ε(∇ · D)(x) and E[δδ ⊤ ] = εD(x) + O(ε2 ). Using the circulant property of trace, we obtain an expression in terms of ε; (iv) after subtracting the θ-independent constant, normalizing by ε/4, and taking ε → 0, we use the divergence identity ∇ · (D(x)∇Eθ ) = ∇Eθ (x)⊤ (∇ · D)(x) + tr(D(x)HEθ (x)) to recover Equation 3. Proof. Define the integral

Z Iε (x; θ) =

qε (y | x) exp

1 [Eθ (x) − Eθ (y)] dy 2

Recall that the conditional density is qε (y | x) = Z exp −

∥y − x − 2ε b(x)∥2D(x)−1 2ε

(21)

where Z =

(2π)d/2

√1

det εD(x)

is the normalization constant and b(x) = (∇ · D)(x). Due to the

construction, it’s easy to see that y | x ∼ N (x + 2ε b(x), εD(x)). Define δ = y − x, then δ | x ∼ N 2ε b(x), εD(x) . So Equation 21 will be Z 1 Iε (x; θ) = wε (δ | x)exp [Eθ (x) − Eθ (x + δ)] dδ 2 where ⊤ 1 ε ε 1 −1 p exp − δ − b(x) (εD(x)) δ − b(x) wε (δ | x) = 2 2 2 (2π)d/2 det εD(x) Since Eθ is C 2 and assuming ε to be small, we can expand Eθ (x + δ) around x to second order. 1 ⊤ Eθ (x + δ) = Eθ (x) + ∇Eθ (x) δ + δ ⊤ HEθ (x)δ + O(∥δ∥3 ) 2 where HEθ (x) denotes the Hessian of Eθ at x. Therefore, 1 ⊤ Eθ (x) − Eθ (x + δ) = −∇Eθ (x) δ − δ ⊤ HEθ (x)δ + O(∥δ∥3 ) 2 Applying the exponential function and expanding to second order in δ: 1 exp [Eθ (x) − Eθ (x + δ)] 2 1 1 ⊤ = exp − ∇Eθ (x) δ − δ ⊤ HEθ (x)δ + O(∥δ∥3 ) 2 4 2 1 1 1 ⊤ ⊤ ∇Eθ (x) δ + O(∥δ∥3 ) = 1 − ∇Eθ (x) δ − δ ⊤ HEθ (x)δ + 2 4 8 Since δ has variance of order ε, we have ∥δ∥ = O(ε1/2 ), and thus the error terms satisfy O(∥δ∥3 ) = O(ε3/2 ). Using the fact that for any scalar a = v ⊤ w, we have a2 = Tr(v ⊤ wv ⊤ w) = Tr(wv ⊤ wv ⊤ ), and for any matrix A, δ ⊤ Aδ = Tr(Aδδ ⊤ ), we can rewrite: 1 1 1 ⊤ [Eθ (x) − Eθ (x + δ)] = 1 − ∇Eθ (x) δ − Tr HEθ (x)δδ ⊤ exp 2 2 4 1 ⊤ + Tr ∇Eθ (x)∇Eθ (x) δδ ⊤ + O(ε3/2 ) 8 Substituting into Iε (x; θ) and using the linearity of integration and trace: Z 1 [Eθ (x) − Eθ (y)] dy Iε (x; θ) = wε (δ | x) exp 2 1 1 ⊤ = 1 − ∇Eθ (x) E[δ] − Tr HEθ (x)E[δδ ⊤ ] 2 4 1 ⊤ + Tr ∇Eθ (x)∇Eθ (x) E[δδ ⊤ ] + O(ε3/2 ) 8 where the expectation is taken with respect to the Gaussian distribution N ( 2ε b(x), εD(x)). Also ε ε b(x) = (∇ · D)(x) 2 2 E[δδ ⊤ ] = εD(x) + E[δ]E[δ]⊤ = εD(x) + O(ε2 ) E[δ] =

Substituting these moments: ε ε ⊤ Iε (x; θ) = 1 − ∇Eθ (x) (∇ · D)(x) − Tr (D(x)HEθ (x)) 4 4 ε ⊤ + Tr D(x)∇Eθ (x)∇Eθ (x) + O(ε3/2 ) 8 16

⊤

Using the cyclic property of trace, Tr(D(x)∇Eθ (x)∇Eθ (x) ) = ∇Eθ (x) D(x)∇Eθ (x): ε ε ε ⊤ ⊤ Iε (x; θ) = 1 − ∇Eθ (x) (∇ · D)(x) − Tr (D(x)HEθ (x)) + ∇Eθ (x) D(x)∇Eθ (x) + O(ε3/2 ) 4 4 8 ⊤

Using the relation ∇ · (D(x)∇Eθ (x)) = ∇Eθ (x) (∇ · D)(x) + Tr (D(x)HEθ (x)), the objective becomes h i ε ε ⊤ K(θ) = Ex∼p 1 + ∇Eθ (x) D(x)∇Eθ (x) − ∇ · (D(x)∇Eθ (x)) + O(ε3/2 ) 8 4 Removing the θ independent terms, taking the dividing by 4ε and taking the limit ε → 0 yields the desired result in Equation 3.

Proof of Proposition 3.1

The proof follows by applying integration by parts in Equation 4. Proof. Consider Equation 4 " # d X d X 1 ∂ ⊤ ∇ log pθ (x) D(x)∇ log pθ (x) + Tr (D(x)Hlog pθ (x)) + D(x)kℓ L(θ) = Ex∼p ∇ log pθ (x)ℓ 2 ∂xk k=1 ℓ=1 X d X d Z 1 ⊤ D(x)kℓ Hlog pθ (x)kℓ p(x)dx = Ex∼p ∇ log pθ (x) D(x)∇ log pθ (x) + 2 k=1 ℓ=1 x d X d Z X ∂ + p(x)∇ log pθ (x)ℓ D(x)kℓ dx (22) ∂x k k=1 ℓ=1 x | {z } C(x)

Now consider d X d Z X ∂ C(x) = p(x)∇ log pθ (x)ℓ D(x)kℓ dx ∂xk x k=1 ℓ=1

d X d Z X

Z ···

k=1 ℓ=1

=−

d Z d X X k=1 ℓ=1

=−

k=1 ℓ=1

=−

d X d Z X

D(x)kℓ

∂ (p(x)∇ log pθ (x)ℓ ) dx ∂xk

∇p(x)k D(x)kℓ ∇ log pθ (x)ℓ dx −

d X d Z X k=1 ℓ=1

:0 Z ∂ xk =+∞ ] [p(x)∇ log pθ (x) D(x) dx − D(x)kℓ (p(x)∇ log pθ (x))dx ℓ kℓ /k xk =−∞ ∂xk x

d X d Z X k=1 ℓ=1

∇ log p(x)k D(x)kℓ ∇ log pθ (x)ℓ p(x)dx −

D(x)kℓ Hlog pθ (x)kℓ p(x)dx

x d X d Z X k=1 ℓ=1

d X d X = −Ex∼p ∇ log pθ (x)⊤ D(x)∇ log p(x) − k=1 ℓ=1

D(x)kℓ Hlog pθ (x)kℓ p(x)dx

Z D(x)kℓ Hlog pθ (x)kℓ p(x)dx x

where in the third step we have used integration by parts along with the regularity condition (mentioned in Proposition 3.1). Substituting this in Equation 22 yields 1 ⊤ ⊤ L(θ) = Ex∼p ∇ log pθ (x) D(x)∇ log pθ (x) − ∇ log pθ (x) D(x)∇ log p(x) 2 Observe that minimizing L(θ) with respect to θ is equivalent to minimizing 1 L(θ) + Ex∼p ∇ log p(x)⊤ D(x)∇ log p(x) 2 17

since the added term does not depend on θ. Therefore, we have 1 LGSM (θ) = L(θ) + Ex∼p ∇ log p(x)⊤ D(x)∇ log p(x) 2 h i 1 2 = Ex∼p ∥∇ log pθ (x) − ∇ log p(x)∥D(x) 2 where the last equality follows from expanding the weighted norm. This completes the proof.

Proof of Theorem 3.2

To prove Theorem 3.2, we first establish the following auxiliary lemmas. Since the proofs involve repeated use of multivariate integrals, we introduce the following shorthand notation. For an integral over Rd of the form Z f (x) dx1 · · · dxd , we write, for indices i < j, dxi:j = dxi · · · dxj . D.1

Supporting Lemmas

Lemma D.1. Let x ∈ Rd and let B ⊂ Rd denote the unit Euclidean ball centered at the origin. For any coordinate index i ∈ {1, 2, . . . , d}, define Z 1 I= xi exp − x⊤ x dx 2 B Then I = 0. Proof. Without loss of generality, assume i = 1. By symmetry of the integrand and the domain B, the result holds for any choice of i. So, we have ! Z 1 Z √1−Pdk=2 x2k 1 2 1 2 I= exp − xd · · · x1 exp − x1 dx1 dx2:d √ P 2 2 2 xd =−1 x1 =− 1− d k=2 xk Since the integrand x1 exp(− 12 x21 ) is an odd function of x1 , we can conclude that I = 0. Lemma D.2. Let x ∈ Rd with d ≥ 2, and let B ⊂ Rd be the unit Euclidean ball centered at the origin. For distinct indices i ̸= j where i, j ∈ {1, 2, . . . , d}, define Z 1 ⊤ I= xi xj exp − x x dx 2 B Then I = 0. Proof. Without loss of generality, take i = 1 and j = 2; by symmetry the value is the same for any distinct pair. So ! Z 1 Z √1−Pdk=2 x2k 1 2 1 2 I= x1 exp − x1 dx1 dx2:d exp − xd · · · √ P 2 2 2 x1 =− 1− d xd =−1 k=2 xk The innermost integral is an odd function of x1 , Remark D.1. Following the definitions of Lemma D.2, consider the matrix valued intergral Z 1 I= exp − x⊤ x xx⊤ dx 2 B 18

where the integration is done element-wise. Then Iij = 0 for i ̸= j from Lemma D.2. For any i ∈ {1, 2, . . . , d}, the integral Z 1 Iii = exp − x⊤ x x2i dx 2 B does not depend on i. Assuming that Iii = C3 where C3 is constant, we get I = C3 I where I is identity matrix of size d × d. Lemma D.3. Let x ∈ Rd and let B ⊂ Rd be the unit Euclidean ball centered at the origin. For indices i, j, k ∈ {1, . . . , d} (not necessarily distinct), define Z 1 I= xi xj xk exp − x⊤ x dx 2 B Then I = 0. Proof. Observe that product xi xj xk contains an odd power of at least one dimension. Following the steps as proof of Lemma D.2, we conclude that I = 0. Lemma D.4. Let x ∈ Rd with d ≥ 2, and let B ⊂ Rd be the unit Euclidean ball centered at the origin. For distinct indices i ̸= j, where i, j ∈ {1, 2, . . . , d}, define Z Z 1 ⊤ 1 ⊤ 4 2 2 I1 = xi exp − x x dx and I2 = xi xj exp − x x dx 2 2 B B Then I1 = 3I2 . Note that by symmetry, both integrals are independent of choice of indices i and j. Proof. Without loss of generality, take i = 1, j = 2, by symmetry, the results holds for any distinct pair. Consider the transformation for I1    y1√+y2  x1 2 −y√ 1 +y2   x2    2  x    3 =  y x=  3 .   ..   ..  .  xd yd which in matrix form is x = Qy =

Q2 0

0 Id−2

where

Q2 =

√1 2 − √12

√1 2 √1 2

As Q is orthogonal, x⊤ x = (Qy)⊤ (Qy) = y ⊤ y. Because of this, the integration region remains B and the absolute value of Jacobian determinant is |det Q| = 1. Applying this change of variable 4 Z y1 + y2 1 ⊤ √ I1 = exp − y y dy 2 2 B Z 1 1 ⊤ 4 3 2 2 3 4 = (y + 4y1 y2 + 6y1 y2 + 4y1 y2 + y2 ) exp − y y dy 4 B 1 2 Z Z Z 1 1 ⊤ 3 1 ⊤ 1 1 ⊤ 4 2 2 4 y exp − y y dy + y y exp − y y dy + y exp − y y dy = 4 B 1 2 2 B 1 2 2 4 B 2 2 | | | {z } {z } {z } I1

where the terms with odd powers vanish by Lemma D.3. Rearranging the above equation will give us I1 = 3I2 . 19

Lemma D.5. Let x ∈ Rd , and let B ⊆ Rd be the unit Euclidean ball centered at the origin. For indices i, j, k, ℓ ∈ {1, . . . , d}, define Z 1 I= xi xj xk xℓ exp − x⊤ x dx 2 B Then I = C6 (δij δkℓ + δik δjℓ + δiℓ δjk ) where δpq is the Kronecker delta (equal to 1 if p = q and 0 otherwise), and C6 is a constant. Proof. If the product xi , xj , xk , xℓ contains any index with odd total power, then I = 0 by following similar arguments from Lemma D.1, Lemma D.2, Lemma D.3. Therefore, we need indices to appear as an even power and occurs only when four indices can be paired into two distinct pairs or when all four indices are equal. R  x4i exp − 12 x⊤ x dx if i = j = k = ℓ  B R   1 ⊤ 2 2  x x dx if two equal pairs are (i, j) and (k, ℓ) with i ̸= k x x exp − RB i k 2 1 ⊤ 2 2 I= x x exp − 2 x x dx if two equal pairs are (i, k) and (j, ℓ) with i ̸= j RB i2 2j    x x exp − 12 x⊤ x dx if two equal pairs are (i, ℓ) and (j, k) with i ̸= j  B i j  0 otherwise Using Lemma D.4, we can unify all the cases which is given by I = C6 (δij δkℓ + δik δjℓ + δiℓ δjk ) where C6 =

1 3

1 x4i exp − x⊤ x dx 2 B

By symmetry, this constant is independent of the choice of i. The reason we choose C6 with x4i and not x2i x2j (i ̸= j) is to incorporate the case when d = 1. D.2

Proof

For ease of reading, we first present a sketch of the proof followed by the complete proof. Proof Sketch. The proof proceeds via Taylor expansion of objective in small radius limit. The key steps are as follows: (i) we apply an invertible linear transformation that maps the Bregman ball Crϕ (x) to a Euclidean ball. This coordinates transformation eases the integration in subsequent steps; (ii) since r(x) is small, we expand exp 12 [Eθ (x) − Eθ (y)] to second order; (iii) similarly, we (y−x)⊤ Hϕ (x̄) (y−x) expand the soft connectivity function exp − to first order by first expanding 2r(x)2 Hϕ (x̄) around Hϕ (x); (iv) by symmetry, all odd moments vanish, while the even moments contribute to terms in Lϕ (θ). After normalizing by r(x)d+2 and taking r(x) → 0, we recover Equation 10. Proof. Observe that r exists because Ω is open set. The objective with the weighting function is given by    Z (y − x)⊤ Hϕ (x̄)(y − x) 1   [E (x) − E (y)] dy  Lϕr (θ) = Ex∼p  exp − exp θ θ 2r(x)2 2  y∈Crϕ (x)  | {z } w(y,x)

We approximate the exponentials assuming r(x) is small i.e. y is close to x. Here, we have Hϕ : Rd → Rd×d . For y near to x, we approximate Hϕ (x̄) by Hϕ (x). The i-th partial derivative is ∂ d×d . By first-order Taylor expansion, ∂xi Hϕ (x) ∈ R d y−x 1X ∂ Hϕ (x̄) = Hϕ x + ≈ Hϕ (x) + Hϕ (x)(y − x)k 2 2 ∂xk k=1

where (y − x)k denotes the k-th coordinate of y − x. So (y − x)⊤ Hϕ (x̄) (y − x) = (y − x)⊤ Hϕ (x)(y − x) +

d 1X 4 (y − x)⊤ Hϕ (x)(y − x) (y − x)k + O ∥y − x∥ 2 k=1

Substituting into the weighting function, we have (y − x)⊤ Hϕ (x)(y − x) × w(y, x) = exp − 2r(x)2 ! d 1 X ⊤ ∂ exp − (y − x) Hϕ (x)(y − x) (y − x)k × 4r(x)2 ∂xk k=1 1 4 O ∥y − x∥ exp r(x)2 (y − x)⊤ Hϕ (x)(y − x) × = exp − 2r(x)2 ! d 1 X 1 6 ⊤ ∂ 1− O ∥y − x∥ (y − x) Hϕ (x)(y − x) (y − x)k + 4r(x)2 ∂xk r(x)4 k=1 (23) where we performed a first-order Taylor expansion of exp z. Using the fact that ∥y − x∥ is of order r(x), the second term in first equality is of order r(x) and third term is of the order r(x)2 . As we are doing taylor expansion of first order (in terms of r(x)), we can safely assume that third term is 0 and exponential of that will be 1. Following the same steps as proof of Theorem 3.1, we perform a second-order taylor expansion of energy exponential term. So we have 1 1 1 [Eθ (x) − Eθ (y)] = 1 − ∇Eθ (x)⊤ (y − x) − (y − x)⊤ HEθ (x)(y − x)+ exp 2 2 4 1 2 3 ∇Eθ (x)⊤ (y − x) + O ∥y − x∥ (24) 8 Using the approximations from Equation 23 and Equation 24 in Equation 8, we obtain Z (y − x)⊤ Hϕ (x)(y − x) Irϕ (x; θ) = exp − × 2r(x)2 y∈Crϕ (x) ! d 1 X ∂ 1 6 1− (y − x)⊤ Hϕ (x)(y − x) (y − x)k + O ∥y − x∥ × 4r(x)2 ∂xk r(x)4 k=1

1 1 1 − ∇Eθ (x)⊤ (y − x) − (y − x)⊤ HEθ (x)(y − x)+ 2 4 ! 1 2 3 ∇Eθ (x)⊤ (y − x) + O ∥y − x∥ dy 8 We will evaluate the integral term by term. The order terms will be handled at the last. Applying the following change of variables u=

1 1 Hϕ (x) 2 (y − x) r(x)

or equivalently

y − x = r(x)Hϕ (x)− 2 u

Under this transformation, we have 1

(y − x)⊤ Hϕ (x)(y − x) = r(x)2 u⊤ Hϕ (x)− 2 Hϕ (x)Hϕ (x)− 2 u = r(x)2 u⊤ u. 21

Thus, the Bregman ball Crϕ (x) transforms to unit Euclidean ball B = {u : u⊤ u ≤ 1}. Also, we have ⊤ 1 du dy = det r(x) Hϕ (x)− 2 1

= r(x)d det(Hϕ (x))− 2 du where we have used the fact that Hϕ (x) is symmetric and positive definite. We now evaluate the integral term by term. 1. We have Z (y − x)⊤ Hϕ (x)(y − x) dy I1 = exp − 2r(x)2 y∈Crϕ (x) 1

= C1 r(x)d det(Hϕ (x))− 2 R where C1 = u∈B exp − 21 u⊤ u du is a constant. 2. We have (y − x)⊤ Hϕ (x)(y − x) exp − ∇Eθ (x)⊤ (y − x)dy 2 ϕ 2r(x) y∈Cr (x) Z 1 1 ⊤ 1 d+1 − 12 det(Hϕ (x)) ∇Eθ (x)⊤ Hϕ (x)− 2 udu exp − u u = − r(x) 2 2 u∈B Z d X 1 ⊤ d+1 − 12 ∇Eθ (u) Hϕ (x) i ui du det(Hϕ (x)) = − r(x) 2 u∈B i=1

I2 = −

1 2

=0 The last steps uses Lemma D.1. 3. We have Z 1 (y − x)⊤ Hϕ (x)(y − x) I3 = − exp − (y − x)⊤ HEθ (x)(y − x)dy 4 y∈Crϕ (x) 2r(x)2 Z 1 1 r(x)d+2 1 ⊤ − 21 det(Hϕ (x)) =− exp − u u u⊤ Hϕ (x)− 2 HEθ (x)Hϕ (x)− 2 udu 4 2 u∈B Z 1 1 1 1 ⊤ r(x)d+2 −2 det(Hϕ (x)) exp − u u Tr u⊤ Hϕ (x)− 2 HEθ (x)Hϕ (x)− 2 u du =− 4 2 u∈B Z 1 1 1 r(x)d+2 1 =− det(Hϕ (x))− 2 exp − u⊤ u Tr Hϕ (x)− 2 uu⊤ Hϕ (x)− 2 HEθ (x) du 4 2 u∈B Z d+2 1 1 1 r(x) 1 =− det(Hϕ (x))− 2 Tr Hϕ (x)− 2 exp − u⊤ u uu⊤ du Hϕ (x)− 2 HEθ (x) 4 2 u∈B d+2 1 C3 r(x) =− det(Hϕ (x))− 2 Tr Hϕ (x)−1 HEθ (x) 4 R where last step follows from Remark D.1, with C3 = u∈B exp − 12 u⊤ u u2i du. 4. We have Z 2 1 (y − x)⊤ Hϕ (x)(y − x) I4 = exp − ∇Eθ (x)⊤ (y − x) dy 2 8 y∈Crϕ (x) 2r(x) Z d+2 1 1 1 ⊤ r(x) − 12 exp − u u ∇Eθ (x)⊤ Hϕ (x)− 2 uu⊤ Hϕ (x)− 2 ∇Eθ (x)du = det(Hϕ (x)) 8 2 u∈B =

1 C3 r(x)d+2 det(Hϕ (x))− 2 ∇Eθ (x)⊤ Hϕ (x)−1 ∇Eθ (x) 8

where we apply Remark D.1 as in the previous step. 22

5. We have 1 I5 = − 4r(x)2 =−

1 4r(x)2

! X d (y − x)⊤ Hϕ (x)(y − x) ⊤ ∂ (y − x) Hϕ (x)(y − x) (y − x)k dy exp − 2r(x)2 ∂xk y∈Crϕ (x) k=1 Z d X ∂ (y − x)⊤ Hϕ (x)(y − x) (y − x)k (y − x)i (y − x)j dy Hϕ (x)ij exp − ∂xk 2r(x)2 y∈Crϕ (x)

k,i,j=1

d X

= −C5

k,i,j=1

Z d X 1 ∂ −1 −1 −1 exp − u⊤ u up uq us du Hϕ (x)ij Hϕ (x)kp2 Hϕ (x)iq 2 Hϕ (x)js2 ∂xk 2 u∈B p,q,s=1

=0 1

where C5 =

r(x)d+1 det(Hϕ (x))− 2 . The last step follows from Lemma D.3. 4

6. Denote ∂x∂ k Hϕ (x)ij = ∂k Hϕ (x)ij . Then Z 1 (y − x)⊤ Hϕ (x)(y − x) I6 = exp − × 8r(x)2 y∈Crϕ (x) 2r(x)2 ! d X ∂ (y − x)k (y − x)⊤ Hϕ (x)(y − x) ∇Eθ (x)⊤ (y − x) dy ∂xk k=1

1 8r(x)2

d X

∇Eθ (x)ℓ ∂k Hϕ (x)ij ×

k,i,j,ℓ=1

(y − x)⊤ Hϕ (x)(y − x) exp − 2r(x)2 y∈Crϕ (x)

= C6′

d X

(y − x)k (y − x)i (y − x)j (y − x)ℓ dy

∇Eθ (x)ℓ ∂k Hϕ (x)ij ×

k,i,j,ℓ=1 d X

−1 −1 −1 −1 Hϕ (x)kp2 Hϕ (x)iq 2 Hϕ (x)jt 2 Hϕ (x)ℓs 2

p,q,t,s=1

T (x)

! 1 ⊤ exp − u u up uq ut us du 2 u∈B }

r(x)d+2 det(H (x))− 2

ϕ where C6′ = . 8 R 1 1 ⊤ 4 3 u∈B exp − 2 u u ui du, we have

d X

T (x) = C6

Using the result of Lemma D.5, where C6

−1

Hϕ (x)kp2 Hϕ (x)iq 2 Hϕ (x)jt 2 Hϕ (x)ℓs 2 (δpq δts + δpt δqs + δps δqt )

p,q,t,s=1

= C6

d X

−1

p=1

+ C6

−1

Hϕ (x)kp2 Hϕ (x)ip 2

d X

d X t=1

−1

Hϕ (x)kp2 Hϕ (x)jp2

p=1

+ C6

d X p=1

Hϕ (x)jt 2 Hϕ (x)ℓt 2

d X

−1

Hϕ (x)iq 2 Hϕ (x)ℓq 2

q=1 −1

−1

Hϕ (x)kp2 Hϕ (x)ℓp2

d X

Hϕ (x)iq 2 Hϕ (x)jq2

q=1 −1

−1

Since Hϕ (x) is a symmetric matrix, we have p=1 Hϕ (x)kp2 Hϕ (x)ip 2 = Hϕ (x)−1 ki . Applying this identity to all index pairs yields −1 −1 −1 −1 −1 T (x) = C6 Hϕ (x)−1 H (x) + H (x) H (x) + H (x) H (x) ϕ ϕ ϕ ϕ ϕ ij jℓ kj iℓ kℓ ki

Substituting T (x) back into I6 , we obtain

I6 = C6′ C6

d X

−1 −1 −1 −1 −1 ∇Eθ (x)ℓ ∂k Hϕ (x)ij Hϕ (x)−1 ki Hϕ (x)jℓ + Hϕ (x)kj Hϕ (x)iℓ + Hϕ (x)kℓ Hϕ (x)ij

k,i,j,ℓ=1

(25) Using the definition given in Equation 9 and matrix derivative identities (cf. equations C.21, C.22 from Appendix C in Bishop [3]), we have 1

∂ ∂ 1 ∂e− 2 log det(Hϕ (x)) Gϕ (x) = Hϕ (x)−1 + Hϕ (x)−1 1 ∂xk ∂xk det(Hϕ (x)) 2 ∂xk 1 1 −1 −1 ∂k Hϕ (x)Hϕ (x)−1 − ∂k Hϕ (x) Hϕ (x)−1 =− 1 Hϕ (x) 1 Tr Hϕ (x) det(Hϕ (x)) 2 2 det(Hϕ (x)) 2 Define

d X ℓ=1

d X ∂ ∇Eθ (x)ℓ Gϕ (x)kℓ ∂xk k=1

d X

−1 1

det(Hϕ (x)) 2 ℓ=1

∇Eθ (x)ℓ

k=1

d X

−1 2 det(Hϕ (x))

1 2

−1 2 det(Hϕ (x))

1 2

d X

∇Eθ (x)ℓ

d X

1 Hϕ (x)−1 ∂k Hϕ (x)Hϕ (x)−1 kℓ + Tr Hϕ (x)−1 ∂k Hϕ (x) Hϕ (x)−1 kℓ 2

2 Hϕ (x)−1 ∂k Hϕ (x)Hϕ (x)−1 kℓ + Tr Hϕ (x)−1 ∂k Hϕ (x) Hϕ (x)−1 kℓ

ℓ=1

k=1

d X

∇Eθ (x)ℓ

ℓ=1

d X

−1 2 det(Hϕ (x))

1 2

Hϕ (x)−1 ∂k Hϕ (x)Hϕ (x)−1 kℓ + Hϕ (x)−1 ∂k Hϕ (x)Hϕ (x)−1 ℓk

k=1

+ Tr Hϕ (x)−1 ∂k Hϕ (x) Hϕ (x)−1 kℓ =

−1 ∇Eθ (x)ℓ ∂k Hϕ (x)ij Hϕ−1 (x)ki Hϕ−1 (x)jℓ + Hϕ (x)−1 ℓi Hϕ (x)jk

ℓ,k,i,j=1

! −1 + Hϕ (x)−1 kℓ Hϕ (x)ji

Comparing with Equation 25, we obtain d

X X ∂ 1 I6 R=− =⇒ I6 = −2C6 C6′ det(Hϕ (x)) 2 ∇Eθ (x)ℓ 1 ′ ∂xk 2 det(Hϕ (x)) 2 C6 C6 ℓ=1

k=1

Substituting C6′ =

r(x)d+2 det(Hϕ (x))− 2 , we have 8

X ∂ C6 r(x)d+2 X I6 = − ∇Eθ (x)ℓ 4 ∂xk ℓ=1

k=1

1 1

det(Hϕ (x)) 2

Hϕ (x)

−1 kℓ

−1

det(Hϕ (x)) 2

Hϕ (x)

kℓ

7. We have Z (y − x)⊤ Hϕ (x)(y − x) 1 × exp − I7 = 16r(x)2 y∈Crϕ (x) 2r(x)2 ! d X ⊤ ∂ ⊤ (y − x) Hϕ (x)(y − x) (y − x)k (y − x) HEθ (y − x) dy ∂xk k=1

d X 1 HEθ (x)ab Hϕ (x)ij 16r(x)2 k,i,j,a,b=1 Z (y − x)⊤ Hϕ (x)(y − x) (y − x)k (y − x)i (y − x)j (y − x)a (y − x)b dy exp − 2r(x)2 y∈Crϕ (x) | {z } Q(x)

Applying the usual change of variables, we have Z (y − x)⊤ Hϕ (x)(y − x) Q(x) = exp − (y − x)k (y − x)i (y − x)j (y − x)a (y − x)b dy 2r(x)2 y∈Crϕ (x)

= C7′

d X

−1

1 exp − u⊤ u up uq ur us ut du 2 u∈B

Hϕ (x)kp2 Hϕ (x)iq 2 Hϕ (x)jr2 Hϕ (x)as2 Hϕ (x)bt 2

p,q,r,s,t=1 1

where C7′ = r(x)d+5 det(Hϕ (x))− 2 . Regardless of the values of p, q, r, s, t, an odd power of at least one u coordinate will remain, concluding that Q(x) = 0, hence I7 = 0. 8. The result for I8 is analogous to I7 (similar to the relationship between steps 3 and 4). We have Z 1 (y − x)⊤ Hϕ (x)(y − x) I8 = − exp − 32r(x)2 y∈Crϕ (x) 2r(x)2 ! d X ∂ 2 (y − x)⊤ Hϕ (x)(y − x) (y − x)k ∇Eθ (x)⊤ (y − x) dy ∂xk k=1

=0 Since this calculation is very similar to step 7, we omit the detailed steps. So our integral without the order terms is C3 r(x)d+2 C3 r(x)d+2 1 d −1 ⊤ −1 C r(x) − Tr H (x) H (x) + ∇E (x) H (x) ∇E (x) 1 ϕ Eθ θ ϕ θ 1 4 8 det(Hϕ (x)) 2 ! d d X C6 r(x)d+2 X ∂ 1 −1 − ∇Eθ (x)ℓ Hϕ (x) 4 ∂xk det(Hϕ (x)) 12 ℓ=1

k=1

kℓ

The order term involves the following integral Z (y − x)⊤ Hϕ (x)(y − x) 1 6 exp − × O ∥y − x∥ × 2r(x)2 r(x)4 y∈Crϕ (x) ! 2 1 1 1 3 ⊤ ⊤ ⊤ 1 − ∇Eθ (x) (y − x) − (y − x) HEθ (x)(y − x) + ∇Eθ (x) (y − x) + O ∥y − x∥ dy 2 4 8 The odd powers will go to 0, so the above integrals yields the terms 1 det(Hϕ (x))

1 2

D1 r(x)d+2 + D3 (θ)r(x)d+4 + D4 (θ)r(x)d+4

where D1 is constant and D3 , D4 are functions of θ. Since our expansion retains terms only up to order r(x)d+2 , we neglect the higher-order terms involving D3 and D4 , and retain only the 25

contribution of D1 . So we have Irϕ (x; θ) =

C3 r(x)d+2 Tr Hϕ (x)−1 HEθ (x) + 4 det(Hϕ (x)) ! r(x)d+2 D1 C3 r(x)d+2 ⊤ −1 ∇Eθ (x) Hϕ (x) ∇Eθ (x) + 1 8 det(Hϕ (x)) 2 ! d d X C6 r(x)d+2 X 1 ∂ −1 − Hϕ (x) ∇Eθ (x)ℓ 4 ∂xk det(Hϕ (x)) 12 1

1 2

C1 r(x)d −

ℓ=1

k=1

kℓ

Removing terms that are not dependent on θ (terms involving C1 , D1 ), the lowest order is r(x)d+2 . So we consider the following limit ϕ Ir (x; θ) Er→0 (r(x)d+2 /4) One can construct function r such that r(x) is bounded for all x ∈ Ω (see the Remark D.2 for more details). Consequently, the integrand is dominated by an integrable function, where integrability follows from the assumptions of the theorem. Therefore, by the dominated convergence theorem, the limit and expectation can be interchanged. The objective then becomes # " C3 ⊤ ⊤ ∇Eθ (x) Gϕ (x)∇Eθ (x) − C3 Tr (Gϕ (x)HEθ (x)) − C6 ∇Eθ (x) (∇ · Gϕ )(x) E x∼p 2 Minimizing the above objective is equivalent to minimizing the objective in Equation 10 with 6 λ= C C3 . Remark D.2. Since Ω is an open set, for every x ∈ Ω, there exists a radius function r′ such that Crϕ′ (x) ⊂ Ω. We now show that one can choose such a radius function to be uniformly bounded. Let B > 0 be any fixed constant, and define r(x) = min{r′ (x), B}. By construction, r(x) ≤ B for all x ∈ Ω. Moreover, since reducing the radius can only shrink the corresponding set, we have Crϕ (x) ⊂ Crϕ′ (x) ⊂ Ω. Thus, r is a uniformly bounded radius function satisfying the required containment condition. Remark D.3. The constant λ in Equation 10 depends on the choice of the weighting function used in the local construction. In particular, consider the alternative weighting (y − x)⊤ Hϕ (x)(y − x) w(y, x) = exp − × 2r(x)2 ! d X C3 ⊤ ∂ 1− (y − x) Hϕ (x)(y − x) (y − x)k . 4r(x)2 C6 ∂xk k=1

For this alternative weighting, the weighting function itself does not need to be expanded. Instead, we expand the integrand appearing in the MPF objective and proceed with the remainder of the proof. The resulting in this case gives a coefficient of 1 for ∇Eθ (x)⊤ (∇ · Gϕ )(x), and hence λ = 1. Moreover, λ = 1 is precisely the value for which the resulting objective is a proper second-order scoring rule, as established in Proposition 4.1.

Proof of Proposition 3.2 and Its Extension to Unbounded Sets

E.1

Proof

The proof follows by applying integration by parts on Equation 11. Proof. For the energy-based model, we have Eθ (x) = − log pθ (x) − log Z(θ), which yields ∇Eθ (x) = −∇ log pθ (x)

and 26

HEθ (x) = −Hlog pθ (x)

Rewriting Equation 11 in terms of log pθ , we obtain # " d X d X 1 ∂ ϕ ⊤ L (θ) = Ex∼p ∇ log pθ (x) Gϕ (x)∇ log pθ (x) + Tr (Gϕ (x)Hlog pθ (x)) + Gϕ (x)kℓ ∇ log pθ (x)ℓ 2 ∂xk k=1 ℓ=1 X d X d Z 1 = Ex∼p ∇ log pθ (x)⊤ Gϕ (x)∇ log pθ (x) + Gϕ (x)kℓ Hlog pθ (x)kℓ p(x)dx 2 k=1 ℓ=1 x∈Ω d X d Z X ∂ + Gϕ (x)kℓ dx (26) p(x)∇ log pθ (x)ℓ ∂xk k=1 ℓ=1 x∈Ω | {z } C(x)

Denote ds is surface element on ∂Ω, consider d X d Z X ∂ C(x) = Gϕ (x)kℓ dx p(x)∇ log pθ (x)ℓ ∂x k k=1 ℓ=1 x∈Ω Z d X d Z X = p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ nk (x)ds − x∈∂Ω

k=1 ℓ=1

=−

d X d Z X k=1 ℓ=1

=−

k=1 ℓ=1

=−

Gϕ (x)kℓ

∂ (p(x)∇ log pθ (x)ℓ )dx ∂xk

∂ (p(x)∇ log pθ (x)ℓ ) dx ∂xk

∇p(x)k Gϕ (x)kℓ ∇ log pθ (x)ℓ dx −

x∈Ω

d X d Z X k=1 ℓ=1

x∈Ω

d Z d X X

Gϕ (x)kℓ

d Z d X X k=1 ℓ=1

∇ log p(x)k Gϕ (x)kℓ ∇ log pθ (x)ℓ p(x)dx −

x∈Ω

Gϕ (x)kℓ Hlog pθ (x)kℓ p(x)dx

x∈Ω d X d Z X k=1 ℓ=1

⊤

= −Ex∼p ∇ log pθ (x) Gϕ (x)∇ log p(x) −

d X d Z X k=1 ℓ=1

Gϕ (x)kℓ Hlog pθ (x)kℓ p(x)dx

x∈Ω

Gϕ (x)kℓ Hlog pθ (x)kℓ p(x)dx

x∈Ω

where in the second and third step we have used integration by parts along with the regularity condition given in Equation 13. Substituting this in Equation 26 yields 1 Lϕ (θ) = Ex∼p ∇ log pθ (x)⊤ Gϕ (x)∇ log pθ (x) − ∇ log pθ (x)⊤ Gϕ (x)∇ log p(x) 2 Observe that minimizing Lϕ (θ) with respect to θ is equivalent to minimizing 1 Lϕ (θ) + Ex∼p ∇ log p(x)⊤ Gϕ (x)∇ log p(x) 2 since the added term does not depend on θ. Therefore, we have h i 1 1 2 LϕGSM (θ) = Lϕ (θ) + Ex∼p ∇ log p(x)⊤ Gϕ (x)∇ log p(x) = Ex∼p ∥∇ log pθ (x) − ∇ log p(x)∥Gϕ (x) 2 2 where the last equality follows from expanding the weighted norm. This completes the proof. E.2

Treatment of Unbounded sets

In this section, we address how to handle the case where the set Ω is unbounded. The difficulty is that the integration by parts argument used above may no longer be valid in this setting. To overcome this issue, we impose additional regularity conditions beyond those stated in Proposition 3.2. Consider a function f : [0, ∞) → R+ defined as   1 if x ≤ 1   1 exp(− x−1 ) f (x) = 1 − exp(− 1 )+exp − 1 if 1 < x < 2 ( 2−x )  x−1  0 if x ≥ 2 27

The function f is flat in the regions ≤ 1 and ≥ 2 and infinitely smooth in (1, 2). Using this, define ρn : Rd → R+ as ρn (x) =

1 if ∥x∥ ≤ n f (∥x∥ − n + 1) if ∥x∥ > n

where n ∈ N. Observe that ρn is equal to 1 inside the ball of radius n, and equal to 0 outside the ball of radius n + 1. In the intermediate region, it transitions smoothly from 1 to 0. For any k, ℓ ∈ {1, 2, . . . , d}, assume that the following quantities 1.

R x∈Ω

∇ log pθ (x)ℓ ∂x∂ k Gϕ (x)kℓ p(x)dx

|∇ log pθ (x)ℓ Gϕ (x)kℓ nk (x)| p(x)ds

x∈∂Ω

x∈Ω

|∇p(x)k Gϕ (x)kℓ ∇ log pθ (x)ℓ | dx

x∈Ω

|Gϕ (x)kℓ Hlog pθ (x)kℓ | p(x)dx

x∈Ω

|Gϕ (x)kℓ ∇ log pθ (x)ℓ | p(x)dx

are finite and Equation 13 holds. Using assumptions 3, 4, we infer that Z Gϕ (x)kℓ x∈Ω

∂ (p(x)∇ log pθ (x)ℓ ) dx ≤ ∂xk

Z |∇p(x)k Gϕ (x)kℓ ∇ log pθ (x)ℓ | dx x∈Ω

Z +

|Gϕ (x)kℓ Hlog pθ (x)kℓ | p(x)dx x∈Ω

<∞ Then for any k, ℓ and n ∈ N, we have Ex∼p

∂ ρn (x)∇ log pθ (x)ℓ Gϕ (x)kℓ ∂xk

≤ Ex∼p

∂ ∇ log pθ (x)ℓ Gϕ (x)kℓ ∂xk

<∞ Consider Z lim

n→∞

ρn (x)p(x)∇ log pθ (x)ℓ x∈Ω

∂ Gϕ (x)kℓ dx = ∂xk

Z lim ρn (x)p(x)∇ log pθ (x)ℓ

x∈Ω n→∞

Z =

p(x)∇ log pθ (x)ℓ x∈Ω

∂ Gϕ (x)kℓ dx ∂xk

which is exactly equal C(x) from the proof of Proposition 3.2. The first step is due to dominated convergence theorem and second step is using the fact that limn→∞ ρn (x) = 1 for any x ∈ Rd . Each integral on left hand side can be written ∂ ρn (x)p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ dx = ∂x k x∈Ω

ρn (x)p(x)∇ log pθ (x)ℓ x∈Ω∩B̄n+1

∂ Gϕ (x)kℓ dx ∂xk

where B̄n+1 is closure of open ball centered at origin of radius n + 1. As Ω ∩ B̄n+1 is convex and bounded, we can apply integration by parts for each individual integral. Knowing this, we have Z Z ∂ ∂ p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ dx = lim ρn (x)p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ dx n→∞ ∂x ∂x k k x∈Ω Zx∈Ω = lim ρn (x)p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ nk (x)ds n→∞ x∈∂Ω Z ∂ (ρn (x)p(x)∇ log pθ (x)ℓ )dx − lim Gϕ (x)kℓ n→∞ x∈Ω ∂xk Z = lim ρn (x)p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ nk (x)ds n→∞ x∈∂Ω Z ∂ − lim ρn (x)Gϕ (x)kℓ (p(x)∇ log pθ (x)ℓ )dx n→∞ x∈Ω ∂xk Z ∂ ρn (x)dx − lim p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ n→∞ x∈Ω ∂xk Z = lim ρn (x)p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ nk (x)ds x∈∂Ω n→∞ Z ∂ − lim ρn (x)Gϕ (x)kℓ (p(x)∇ log pθ (x)ℓ )dx n→∞ ∂x k x∈Ω Z ∂ ρn (x)dx − lim p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ ∂xk x∈Ω n→∞ where all the interchanging of limits and integral is valid due to dominated convergence theorem. For last term which has ∂x∂ k ρn (x), interchange is possible because ρ ∈ C ∞ and it’s first derivative is continuous over a compact set implying that it is bounded. As limn→∞ ∂x∂ k ρn (x) = 0 and first term is 0 due to assumption of Equation 13, we have Z Z ∂ ∂ p(x)∇ log pθ (x)ℓ Gϕ (x)kℓ dx = − lim ρn (x)Gϕ (x)kℓ (p(x)∇ log pθ (x)ℓ )dx n→∞ ∂x ∂x k k x∈Ω Zx∈Ω ∂ =− Gϕ (x)kℓ (p(x)∇ log pθ (x)ℓ )dx ∂x k x∈Ω One can now follow same proof as bounded case to arrive at completing the squares argument. Remark E.1. The treatment of unbounded domains requires stronger technical assumptions and may not be directly applicable in all practical settings. For structured spaces such as Ω = Rd+ , one can instead follow the approach of Yu et al. [43], where suitable conditions are imposed to justify the use of Fubini-Tonelli theorem. Adopting analogous assumptions in our setting would likewise allow us to carry out the completing-the-squares argument. With this remark, we aim to convey that when additional structure is available on set Ω, the assumptions imposed here can be relaxed.

Proof of Proposition 3.3 1

Proof. This follows as Gϕ (x) 2 ∇ log p1 (x) = Gϕ (x) 2 ∇ log p2 (x) implies ∇ log p1 (x) = (x) ∇ log p2 (x) almost everywhere. This is same as ∇ log pp12 (x) = 0. As Ω is convex and convex (x) subsets are connected in Rd , we conclude that pp12 (x) = c almost everywhere for some constant c. As both p1 , p2 are densities, we conclude that p1 (x) = p2 (x) almost everywhere.

Generalization of Theorem 3.1 and Theorem 3.2

In this section, we derive the GSM objective from a more general starting objective than the one considered in the main text (Equation 2). Motivated by the discussion in Section A.2, we propose the following objective as the starting point Z pθ (x) pθ (y) ′ pθ (y) pθ (y) Lf (θ) = E g(y, x) −f ′ + f −f dy (27) x∼p pθ (y) pθ (x) pθ (x) pθ (x) 29

√ where f : R≥0 → R is strictly convex function. Choosing f (z) = − z recovers the MPF objective (Equation 2). Further more, if the connectivity function g(y, x) is chosen as conditional density q(y | x) and satisfies the symmetry condition q(y | x) = q(x | y), then the above objective coincides with Equation 20. We first establish a supporting lemma and then present the corresponding generalizations. Lemma G.1. Consider the term pθ (x) pθ (y) ′ pθ (y) pθ (y) h(θ) = −f ′ + f −f pθ (y) pθ (x) pθ (x) pθ (x) For y in a sufficiently small neighborhood of x, the second order Taylor expansion of h around y is given by 2 1 3 ′′ ⊤ ⊤ ⊤ −f (1) + f (1) ∇Eθ (x) δ − δ HEθ (x)δ − 2∇Eθ (x) δ + O ∥δ∥ 2 where δ = y − x. Proof. Define ρθ (x, y) = ppθθ (x) (y) = exp (Eθ (y) − Eθ (x)). Then h(θ) = −f ′ (ρθ (x, y)) + ρθ (y, x)f ′ (ρθ (y, x)) − f (ρθ (y, x)) We expand each term separately around y. First, 2 1 ′ (f (1) + f ′′ (1)) ∇Eθ (x)⊤ δ 2 1 ′ 3 ⊤ − f (1)δ HEθ (x)δ + O ∥δ∥ 2

f (ρθ (y, x)) = f (1) − f ′ (1)∇Eθ (x)⊤ δ +

Similarly 2 1 f ′ (ρθ (y, x)) = f ′ (1) − f ′′ (1)∇Eθ (x)⊤ δ + (f ′′ (1) + f ′′′ (1)) ∇Eθ (x)⊤ δ 2 1 ′′ 3 ⊤ − f (1)δ HEθ (x)δ + O ∥δ∥ 2 and 2 1 f ′ (ρθ (x, y)) = f ′ (1) + f ′′ (1)∇Eθ (x)⊤ δ + (f ′′ (1) + f ′′′ (1)) ∇Eθ (x)⊤ δ 2 1 ′′ 3 ⊤ + f (1)δ HEθ (x)δ + O ∥δ∥ 2 Substituting the above expansions into the definition of h(θ) and collecting terms up to second order yields the desired expression. We now state the following generalization of Theorem 3.1. Theorem G.1. Suppose the assumption of Theorem 3.1 hold. Choosing g(y, x) to be conditional density qε (y | x) in Equation 27 defined as   2 ε y − x − b(x) 1 2 D(x)−1  p qε (y | x) = exp − d/2 2ε (2π) det εD(x) where ε > 0. Then for very small ε, the objective can be written as " 1 ′′ Lf (θ) = −f (1) + εf (1) E ∇Eθ (x)⊤ D(x)∇Eθ (x) − Tr (D(x)HEθ (x)) x∼p 2 # 3 − ∇Eθ (x)⊤ (∇ · D)(x) + O ε 2

For Theorem 3.2, we consider a hard neighborhood case same as in main text. In this case, the objective would be "Z # pθ (x) pθ (y) ′ pθ (y) pθ (y) ϕ ′ Lf,r (θ) = E w(y, x) −f + f −f dy x∼p pθ (y) pθ (x) pθ (x) pθ (x) Crϕ (x) (28) where Crϕ (x) is defined as in Equation 7. We now state the following generalization Theorem G.2. Suppose the assumption of Theorem 3.2 hold. Define x̄ = (x + y)/2 and consider the weighting function (y − x)⊤ Hϕ (x̄) (y − x) w(y, x) = exp − 2r(x)2 in Equation 28. Then for very small r, the objective can be written as " C 3 ′′ Lf (θ) = −Cf (1) + f (1) E r(x)d+2 ∇Eθ (x)⊤ Gϕ (x)∇Eθ (x) x∼p 2 # ⊤ d+4 − C3 Tr (Gϕ (x)HEθ (x)) − C6 ∇Eθ (x) (∇ · Gϕ (x)) + O r(x) where Gϕ is defined as in Equation 9 and C, C3 , C6 are positive constants. Remark G.1. Both of the above generalizations closely parallel Theorem 3.3 in Ryu et al. [31] and were inspired by that result. In particular, Theorem 3.3 derives the score matching, whereas our generalizations extend the same perspective to the generalized score matching. Remark G.2. We omit qthe proofs of Theorem G.1 and Theorem G.2, as they follow by replacing the (y) Taylor expansions of ppθθ (x) used in the proofs of Theorem 3.1 and Theorem 3.2 with the expansion given in Lemma G.1. The similarity of the arguments is further reflected in the fact that the constants C3 and C6 appearing in Theorem G.2 are identical to those arising in the proof of Theorem 3.2.

Proper Scoring Rules of the Second Order

H.1

Relevance

Proper scoring rules ensure that minimizing the score matching loss would lead to the model density matching the true data density almost everywhere. This property is important in several applications, for instance, forecasting, where the use of an improper scoring rule can lead to an incorrect assessment of extreme events, which can result in inaccurate predictions [22], and generative modeling, where Bortoli et al. [4] demonstrate that a carefully constructed scoring rule can be leveraged to accelerate sampling. Specifically, in Section K.1, we consider the problem of estimating the rate parameter λ of an exponential random variable defined over the domain R+ , and show that the original score matching objective yields the estimate λ̂ = 0. This is clearly inaccurate because the original score matching loss is not a proper scoring rule in this setting. We proceed to show that convex functions that are appropriately defined on this domain lead to better estimates of λ. Moreover, the relevance of second-order locality, beyond properness, is three-fold. Parry et al. [29] show that proper local scoring rules of odd order do not exist, and that the log-likelihood is the unique proper local rule of order 0. Therefore, order 2 is the minimal non-trivial order for a proper scoring rule. They also establish that all proper local scoring rules of order ≥ 2 can be evaluated without the normalising constant. This explains why score matching circumvents the partition function. Finally, Parry et al. [29] and [8] provide a complete characterisation of all second-order local proper scoring rules via a generating matrix G(x). H.2

Proof of Proposition 4.1

Proof. This proposition holds when one compares Equation 12 and Equation 16, followed by substituting G(x) = Gϕ (x). Expanding the divergence term in Equation 12, we have d X d X 1 ∂ S ϕ (x; θ) = ∇Eθ (x)⊤ Gϕ (x)∇Eθ (x) − Tr (Gϕ (x)HEθ (x)) − ∇Eθ (x)j Gϕ (x)ij 2 ∂xi i=1 j=1 31

Following the same as proof of Proposition 3.2, writing the above equation in terms of log pθ , we have S ϕ (x; θ) =

d d d X d X 1 XX Gϕ (x)ij ∇ log pθ (x)i ∇ log pθ (x)j + Gϕ (x)ij Hlog pθ (x)ji 2 i=1 j=1 i=1 j=1 d X d X

∇ log pθ (x)j

i=1 j=1

∂ Gϕ (x)ij ∂xi

Using the relation Hlog pθ (x) =

1 1 ∇pθ (x)∇pθ (x)⊤ Hpθ (x) − pθ (x) pθ (x)2

Substituting it back, we have S ϕ (x; θ) =

d X d X

d d 1 1 XX 1 Gϕ (x)ij Hpθ (x)ij − Gϕ (x)ij ∇pθ (x)i ∇pθ (x)j p (x) 2 i=1 j=1 pθ (x)2 i=1 j=1 θ

d X d X

∇ log pθ (x)j

i=1 j=1

∂ Gϕ (x)ij ∂xi

Comparing it with Equation 16, we conclude that S ϕ (x; θ) is a proper scoring rule with G(x) = Gϕ (x) and q = pθ .

Convexity for Exponential Family

Proof. We compute the first order and second order partial derivatives of log pθ (x) from Equation 17 d

l=1

∂ log pθ (x) X ∂tl (x) ∂b(x) ∂ 2 log pθ (x) X ∂ 2 tl (x) ∂ 2 b(x) = θl + , = θl + ∂xj ∂xj ∂xj ∂xj ∂xk ∂xj ∂xk ∂xj ∂xk We rewrite L̂(θ) as N

L̂(θ) =

1 X X Djk (xi ) N i=1 j,k=1

∂ 2 log pθ (xi ) 1 ∂ log pθ (xi ) ∂ log pθ (xi ) + ∂xj ∂xk 2 ∂xj ∂xk

∂Djk (xi ) ∂ log pθ (xi ) ∂xj ∂xk (29)

Substituting the partial derivatives from above in Equation 29, we get N

1 X X L̂(θ) = N i=1

j,k=1

d X

d ∂ 2 tl (x) ∂ 2 b(xi ) ∂Djk (xi ) X ∂tl (xi ) ∂b(xi ) Djk (xi ) θl + Djk (xi ) + θl + ∂xj ∂xk ∂xj ∂xk ∂xj ∂xj ∂xj l=1 l=1 ! !! d d X X 1 ∂tl (xi ) ∂b(xi ) ∂tm (x) ∂b(xi ) + Djk (xi ) θl + θm + 2 ∂xj ∂xj ∂xk ∂xk m=1 l=1

Collecting all the terms linear in θ, quadratic in θ and constant with respect to θ separately, we get ! N 1 X L̂(θ) = L(xi ; θ) + Q(xi ; θ) + C(xi ; θ) N i=1 32

(30)

where L(xi ; θ) =

d X

θl Djk (xi )

∂ 2 tl (xi ) ∂Djk (xi ) ∂tl (xi ) ∂tl (xi ) ∂b(xi ) + θl + θl Djk (xi ) ∂xj ∂xk ∂xj ∂xk ∂xj ∂xk

θl Djk (xi )

X ∂Djk (xi ) X X X ∂ 2 tl (xi ) ∂b(xi ) + [Jt⊤ (xi )]kl θl + Djk (xi ) [Jt⊤ (xi )]jl θl ∂xj ∂xk ∂xj ∂xk

j,k,l=1

d X j,k,l=1

d X

θl Djk (xi )Htl (x)kj +

j,k,l=1

d X

j,k=1

l=1

d X j,k=1

El (xi )θl +

l=1

d X

∂Djk (xi ) ⊤ [Jt (xi )θ]k + ∂xj

j,k=1

l=1

d X k,j=1

∂b(xi ) Djk (xi )[Jt⊤ (xi )θ]j ∂xk

(∇ · D(xi ))k [Jt⊤ (xi )θ]k + ∇b(xi )⊤ D(xi )Jt⊤ (xi )θ

k=1

=E(xi )⊤ θ + (∇ · D(xi ))⊤ Jt⊤ (xi )θ + ∇b(xi )⊤ D(xi )Jt⊤ (xi )θ !⊤ = E(xi ) + Jt (xi )(∇ · D(xi )) + Jt (xi )D(xi )∇b(xi ) where El (xi ) =

d P

Djk (xi )Htl (x)kj for l ∈ {1, . . . , d}.

j,k=1

Q(xi ; θ) =

1 2

d X

θl θm Djk (xi )

j,k,l,m=1

∂tm (xi ) ∂tl (xi ) ∂xj ∂xk

d d d d X X ∂tm (xi ) ∂tl (xi ) 1 XX Djk (xi ) θm θl 2 j=1 ∂x ∂xk j m=1 k=1

d d 1 XX

2 j=1

k=1

d d 1 XX

2 j=1

l=1

Djk (xi )

d X

[Jt (xi )⊤ ]jm θm

m=1

d X

[Jt (xi )⊤ ]kl θl

l=1

[Jt (xi )⊤ θ]j Djk (xi )[Jt (xi )⊤ θ]k

k=1

⊤ 1 = Jt (xi )⊤ θ D(xi )Jt (xi )⊤ θ 2 1 = θ ⊤ Jt (xi )D(xi )Jt (xi )⊤ θ. 2 C(xi ; θ) =

d d d d d d ∂b(xi ) ∂b(xi ) X X ∂ 2 b(xi ) X X ∂Djk (xi ) ∂b(xi ) 1 XX Djk (xi ) Djk (xi ) + + 2 j=1 ∂xj ∂xk ∂xj ∂xk j=1 ∂xj xj ∂xk j=1 k=1

k=1

Then Equation 30 can be written as ! N 1 ⊤1 X ⊤ L̂(θ) = θ Jt (xi )D(xi )Jt (xi ) θ 2 N i=1 !⊤ N N 1 X 1 X + E(xi ) + Jt (xi )(∇ · D(xi )) + Jt (xi )D(xi )∇b(xi ) θ + C(xi ; θ) N i=1 N i=1 ! N N 1 P 1 P and since we know that ΓN = Jt (xi )D(xi )Jt (xi )⊤ , C = C(xi ; θ) and gN = N i=1 N i=1 ! N 1 P E(xi ) + Jt (xi )(∇ · D(xi )) + Jt (xi )D(xi )∇b(xi ) , and retrieve the loss in Equation 18. N i=1

Proof of Theorem 5.1

Proof. From Proposition 5.1, we have L̂(θ) =

1 ⊤ ⊤ θ ΓN θ + g N θ+C 2

Existence and Uniqueness of the Minimizer The gradient of L̂(θ) with respect to θ is ∇θ L̂(θ) = ΓN θ + gN . Setting the gradient to zero yields the first-order condition ΓN θ = −gN . By assumption, we know that ΓN is positive definite a.s.. Therefore, the minimizer is a.s. unique and has the closed-form solution θ̂N = −Γ−1 N gN . Almost Sure Convergence Since ΓN and gN are sample averages of i.i.d. random variables {xi }N i=1 drawn from the true distribution p, we apply the strong law of large numbers. We know that Γ0 and g0 exist and are entry-wise finite. Therefore, by the strong law of large numbers, a.s. a.s. ΓN −−→ Γ0 and gN −−→ g0 as N → ∞. The true parameter θ0 minimizes the population loss Ex∼p [L̂(θ)]. By taking expectations in L̂(θ) and applying the first-order condition, Γ0 θ0 = −g0 . a.s. a.s. Since ΓN −−→ Γ0 and Γ0 is invertible (which guarantees Γ−1 0 exists), and gN −−→ g0 . From the −1 continuous mapping theorem for the matrix-valued function w(A) = A [2], we get a.s.

−1 θ̂N = −Γ−1 N gN −−→ θ0 = −Γ0 g0

Asymptotic Normality. Since ΓN θ̂N = −gN and Γ0 θ0 = −g0 , subtracting these equations gives ΓN (θ̂N − θ0 ) = −(gN − g0 ) − (ΓN − Γ0 )θ0 , −1 Multiplying both sides by Γ−1 N yields θ̂N − θ0 = −ΓN [(gN − g0 ) + (ΓN − Γ0 )θ0 ]. Multiplying √ by N , we obtain √ √ N (θ̂N − θ0 ) = −Γ−1 N [(gN − g0 ) + (ΓN − Γ0 )θ0 ]. N √ PN Now observe that N [(gN − g0 ) + (ΓN − Γ0 )θ0 ] = √1N i=1 Zi where " #

Zi = Jt (xi )D(xi )∇b(xi ) + Jt (xi )(∇ · D)(xi ) + E(xi ) − g0 + (Jt (xi )D(xi )Jt (xi )⊤ − Γ0 )θ0 By construction, E[Zi ] = 0 since we defined Γ0 and g0 as Γ0 = E[Jt (xi )D(xi )Jt (xi )⊤ ],

g0 = E[Jt (xi )D(xi )∇b(xi ) + Jt (xi )(∇ · D)(xi ) + E(xi )]. PN d Leveraging the central limit theorem, we know that √1N i=1 Zi − → N (0, Σ0 ). Furthermore, since a.s.

a.s.

−1 ΓN −−→ Γ0 and Γ0 is invertible, we have Γ−1 N −−→ Γ0 . Finally, from Slutsky’s theorem, we get

√

1 X d −1 N (θ̂N − θ0 ) = −Γ−1 Zi − → N (0, Γ−1 0 Σ0 Γ0 ). N · √ N i=1

Experiments

In the subsequent sections, we provide additional experimental results. K.1

1D data - Exponential Distribution

Suppose we have N samples namely x1 , . . . , xN sampled from density p and let the model density be exponential distribution that is pθ (x) = θexp(−θx) where θ > 0. The maximum-likelihood estimate (MLE) is given by N θ̂M LE = PN i=1 xi 34

The objective function for original score matching is # " 2 1 2 d2 1 d log pθ (x) + 2 log pθ (x) = Ex∼p θ E x∼p 2 dx dx 2 The estimate from score matching is θ̂SM = 0. Observe that this estimate is not useful because it doesn’t depend on data and always gives 0. In fact, in this setting the original score matching is not a proper scoring rule. Considering the objective proposed in Yu et al. [43], for any positive function h : (0, ∞) → R+ , we have " # 2 2 1 d d d 1 ′ 2 ′ E h(x) log pθ (x) + h (x) log pθ (x) + h(x) 2 log pθ (x) = E h(x)θ − h (x)θ x∼p 2 x∼p 2 dx dx dx In this case, the estimator will be PN

θ̂GSM −N N = Pi=1 N

h′ (xi )

i=1 h(xi )

If we choose h(x) = x , then we get score matching proposed in Hyvärinen [19]. The estimate will be PN 2 i=1 xi θ̂N SM = PN 2 i=1 xi As the data x lies on R+ , we define convex function ϕ(x) : (0, ∞) → R. Then gϕ (x) =

1 3

ϕ′′ (x) 2

. In

this case, then objective (Equation 11) becomes " # 2 2 1 d d d 1 ′ 2 ′ E gϕ (x) log pθ (x) + gϕ (x) log pθ (x) + gϕ (x) 2 log pθ (x) = E gϕ (x)θ − gϕ (x)θ x∼p 2 x∼p 2 dx dx dx The estimate in this case will be PN

′ i=1 gϕ (xi )

θ̂OU RS = PN

i=1 gϕ (xi )

We have 3 2

1. If we choose ϕ(x) = x log x, then gϕ (x) = x and the estimator will be

1 2 i=1 xi 3 P 2 2 N i=1 xi

i=1

2. If we choose ϕ(x) = − log x, then gϕ (x) = x3 and the estimator will be PNi=1x3i . 4 3

If we choose ϕ(x) = 49 x , then gϕ (x) = x and our estimator will be exactly equal to MLE. Interestingly Hyvärinen [18] showed that for gaussian case (where the support is Rd ) original score matching gives the same estimator as the MLE. Note that if we choose h(x) = x, then also we get MLE. To compare the estimators derived above, we additionally consider data generated from an exponential distribution with parameter equal to 2 i.e. p(x) = 2 exp(−2x), x > 0. For each estimator, we evaluate the empirical mean and standard deviation of the estimates as functions of the sample size N . For a fixed N , these statistics are computed over 50 independent runs. We exclude the original score matching estimator from the comparison, as it is identically zero in this setting. The results are summarized in Figure 2. K.2

Truncated Gaussian Model

Following the experimental setup described in main text, we report results on unbounded positive orthant Ω = Rd+ with d = 10. In this setting, the estimators from Liu et al. [23] are omitted as their formulation requires densities with bounded support. Consequently, we compare our proposed estimators against the estimators proposed in Yu et al. [43] which explicitly designed weight functions 35

Empirical Mean vs Sample Size 3.0

Empirical Std vs Sample Size

MLE

h(x) = x 2 φ(x) = xlog x φ(x) = −log x

Empirical Mean

2.8

h(x) = x 2 φ(x) = xlog x φ(x) = −log x

0.8

Empirical Std

True Rate

2.6

MLE

1.0

0.6

2.4

0.4

2.2

0.2

2.0 25

100

125

150

Number of samples

175

200

100

125

Number of samples

150

175

200

Figure 2: Comparison of various estimators of θ when model density is exponential distribution

Table 1: Performance comparison of parameter estimators (µ and K) for the truncated Gaussian = S 10 at N = 800 across 50 independent trials, where ϕ1 (x) = Pmodel on Ω P P P P 4/3 4/3 9 + (1 − i xi ) , ϕ2 (x) = i xi log(xi ) + (1 − i xi ) log (1 − i xi ), ϕ3 (x) = i xi 4 P P − i log(xi ) − log (1 − i xi ), and h1 (x) = x. We report the Mean, Median, and Standard Deviation of the MSE. MSEµ M ETHODOLOGY OURS OURS OURS T RUNCATED SM [23] Y U ET AL . [43]

MSEK

ϕ(x); h(x)

M EAN

M EDIAN

S TD .

M EAN

M EDIAN

S TD .

ϕ1 ϕ2 ϕ3 — h1

0.077 0.254 739.031 170.848 103.371

0.007 0.023 0.083 0.105 0.083

0.234 0.728 5051.283 1192.049 715.077

2.48 × 104 3.58 × 104 1.50 × 105 5.19 × 104 6.05 × 105

2.17 × 104 2.91 × 104 1.26 × 105 4.77 × 104 5.81 × 105

1.51 × 104 2.25 × 104 7.79 × 104 2.35 × 104 1.05 × 105

P 4/3 h(x) for the positive orthant. We denote our choices of ϕ as ϕ1 (x) = 94 i xi , ϕ2 (x) = P P 2 i log xi , alongside the baseline choices h1 (x) = x and h2 (x) = x i xi log xi , and ϕ3 (x) = − from Yu et al. [43]. Figure 3 illustrates MSE for µ and K across sample sizes, and Table 2 reports quantitative MSE at N = 800. We observe that the baseline h1 achieves the lowest median MSE across all sample sizes N , reaching 0.0016 for µ and 4.70 × 103 for K at N = 800. Among the proposed estimators, ϕ1 achieves the lowest median MSE across all N , tracking closely with the quadratic choice h2 on K (8.47 × 103 versus 7.93 × 103 at N = 800). While ϕ2 and h2 display extreme outliers at N = 200 (with mean MSEs of 17.37 and 65.53 on µ), their error distributions concentrate tightly at N = 500 and N = 800. In contrast, ϕ3 maintains the highest median MSE and widest spread of outliers across all sample sizes, with its median MSE for K exceeding 105 .

Table 2: Performance comparison of parameter estimators (µ and K) for the truncated Gaussian P 4/3 model on Ω = R10 at N = 800 across 50 independent trials, where ϕ1 (x) = 94 i xi , ϕ2 (x) = + P P 2 i xi log xi , ϕ3 (x) = − i log xi , h1 (x) = x, and h2 (x) = x . We report the Mean, Median, and Standard Deviation of the MSE. MSEµ

MSEK

M ETHODOLOGY

ϕ(x); h(x)

M EAN

M EDIAN

S TD .

M EAN

M EDIAN

S TD .

OURS OURS OURS Y U ET AL . [43] Y U ET AL . [43]

ϕ1 ϕ2 ϕ3 h1 h2

0.0038 0.0272 1.9524 0.0019 0.0060

0.0032 0.0108 0.0960 0.0016 0.0043

0.0023 0.0587 5.5063 0.0011 0.0056

8.67 × 103 1.81 × 104 2.04 × 105 4.89 × 103 8.44 × 103

8.47 × 103 1.67 × 104 1.90 × 105 4.70 × 103 7.93 × 103

2.44 × 103 6.00 × 103 7.12 × 104 1.18 × 103 2.33 × 103

Estimation Error for µ vs Sample Size 102

106

101

MSE (log scale)

Estimation Error for K vs Sample Size

107

103

100

105

10 1 10 2

104

10 3

N = 200

N = 500 Number of Samples

N = 800

N = 200

φ1 (x) = 94 xi4/3

φ3 (x) = −

φ2 (x) =

h1 (x) = x

xi log xi

N = 500 Number of Samples

N = 800

h2 (x) = x 2

log xi

Figure 3: Comparison of parameter estimation error (MSE) for µ and K on the positive orthant R10 + across sample sizes N ∈ {200, 500, 800}, evaluated over 50 independent trials. Baselines include 2 hi (x) = xi and hi (x) = xi from Yu et al. [43]. K.3

Dirichlet Model

We consider a Dirichlet distribution with concentration parameters α in Rd+ , whose density is given by pα (x) ∝

d Y

i −1 xα 1x∈∆d−1 i

i=1

where ∆d−1 = {x ∈ Rd : xi > 0, i=1 xi = 1} denotes the d − 1 simplex. Since ∆d−1 has empty interior in Rd , Theorem 3.2 cannot be applied directly. However, we can parameterize the simplex using d − 1 coordinates and formulate the score matching problem on a subset of Rd−1 . Specifically, by dropping the last coordinate, we obtain the parameterization y = (x1 , . . . , xd−1 ) ∈ S d−1

xd = 1 −

d−1 X

i=1

The corresponding density in the new coordinates is given by " #⊤  d−1 X p̃α (y) = pα  y1 , . . . , yd−1 , 1 − yi  i=1

∝

1−

d−1 X

!αd −1 d−1 Y

i=1

yiαi −1 1y∈S d−1

i=1

We compare our proposed estimators against Truncated Score Matching (Truncated SM) [23], Generalized Score Matching for Compositional Data (GSM CD) [45] configured with h(x) = x, and the estimator h(x) = x from Yu et al. [43]. We denote our choices of ϕ as ϕ1 (x) = P 4/3 P P P P 9 4/3 + (1 − i xP ), ϕ2 (x) = i xi log xi + (1 − i xi ) log(1 − i xi ), and ϕ3 (x) = i) 4 (P i xi − i log xi − log(1 − i xi ). Figure 4 illustrates MSE for α across sample sizes, and Table 3 reports quantitative MSE at N = 800. We observe that the proposed estimator with ϕ2 achieves the lowest mean (0.2423) and median (0.1799) MSE across all sample sizes, followed by ϕ1 and ϕ3 . All three proposed choices achieve lower median MSE than the baselines across every evaluated sample size. Among the baselines, Truncated SM attains a median MSE of 0.7314 at N = 800, while GSM CD with h(x) = x yields a median MSE of 1.3099. The coordinate weighting h(x) = x from Yu et al. [43] displays the highest error throughout, maintaining a median MSE around 6.8–8.6 across all sample sizes. As N increases from 200 to 800, the error spreads of ϕ1 , ϕ2 , and ϕ3 decrease monotonically. 37

Estimation Error for α vs Sample Size

MSE (log scale)

101

100

10 1

N = 200

φ1 (x) = 94

µX

Truncated SM

xi4/3 + 1 −

N = 500 Number of Samples X ¶4 ¶

xi 3

φ2 (x) =

GSM CD

xi log(xi ) + 1 −

X ¶

N = 800

xi log 1 −

X ¶

φ3 (x) = − h(x) = x

log(xi ) − log 1 −

X ¶

Figure 4: Comparison of parameter estimation error (MSE) for α of the Dirichlet distribution on the simplex ∆9 (d = 10) across sample sizes N ∈ {200, 500, 800}, evaluated over 50 independent trials. Baselines include Truncated Score Matching (Truncated SM) [23], GSM CD with h(x) = x from Yu et al. [45], and h(x) = x from Yu et al. [43]. Table 3: Performance comparison of parameter estimators for the Dirichlet distribution concentration parameter ∆9 (d= 10) at N = 800 across 50 independent trials, where Pα on the simplex P P P P 4/3 4/3 9 ϕ1 (x) = 4 + (1 − i xi ) , ϕ2 (x) = i xi log(xi ) + (1 − i xi ) log (1 − i xi ), i xi P P and ϕ3 (x) = − i log(xi ) − log (1 − i xi ). We report the Mean, Median, and Standard Deviation of the MSE.

M ETHODOLOGY

ϕ(x); h(x)

M EAN MSE

M EDIAN MSE

S TD . MSE

OURS OURS OURS T RUNCATED SM [23] GSM CD [45] Y U ET AL . [43]

ϕ1 ϕ2 ϕ3 — h(x) = x h(x) = x

0.5066 0.2423 0.6311 1.1192 1.5547 7.9691

0.3871 0.1799 0.4879 0.7314 1.3099 7.6004

0.4244 0.1806 0.5178 1.0901 0.9400 2.1477

K.4 K.4.1

Ablation Studies Robustness Across Diverse Ground Truth Parameter Initializations

The primary experiments performed for Truncated Gaussian model on S d , Rd+ and Dirichlet model on ∆d−1 were for fixed choices of ground truth parameters. To verify that the performance advantages of the proposed estimators are not artifacts of a specific parameter initialization, we systematically evaluate estimator robustness across 50 distinct ground truth parameter configurations. To aggregate performance across diverse parameter initializations, we evaluate estimators using three complementary summary metrics: Win Rate, Average Rank, and the median MSE across configurations. For each ground-truth configuration c ∈ {1, . . . , C}, we compute the median MSE across 50 independent Monte Carlo trials for every candidate estimator. An estimator is assigned a win for configuration c if it achieves the strictly lowest trial-median MSE, yielding a cumulative win PC rate C1 c=1 1rankc =1 . Similarly, candidate methods are ranked from 1 (best) to M (worst) on each PC configuration based on trial-median MSE, from which we report the average rank C1 c=1 rankc . 38

Table 4: Ablation study evaluating the robustness of parameter estimators for the truncated Gaussian model on the standard simplex polytope S 10 across 50 diverse ground-truth configurations (N = 800, 50 trials per configuration). We report the Win Rate, Average Rank, and Median MSE across all configurations. K

µ ϕ(x); h(x)

M ETHODOLOGY OURS OURS OURS T RUNCATED SM [23] Y U ET AL . [43]

ϕ1 ϕ2 ϕ3 — h1

W IN R ATE

AVG . R ANK

M EDIAN MSE

W IN R ATE

AVG . R ANK

M EDIAN MSE

0.90 0.02 0.06 0.00 0.02

1.16 2.20 3.68 3.04 4.92

0.0118 0.0209 0.1005 0.0400 2.2088

1.00 0.00 0.00 0.00 0.00

1.00 2.02 4.00 2.98 5.00

4.06 × 104 6.20 × 104 2.80 × 105 7.95 × 104 1.52 × 106

Table 5: Ablation study evaluating the robustness of parameter estimators for the truncated Gaussian model on the positive orthant R10 + across 50 diverse ground-truth configurations (N = 800, 50 trials per configuration). We report the Win Rate, Average Rank, and Median MSE across all configurations. K

µ M ETHODOLOGY

ϕ(x); h(x)

OURS OURS OURS Y U ET AL . [43] Y U ET AL . [43]

ϕ1 ϕ2 ϕ3 h1 h2

W IN R ATE

AVG . R ANK

M EDIAN MSE

W IN R ATE

AVG . R ANK

M EDIAN MSE

0.06 0.04 0.04 0.82 0.04

2.94 3.78 4.18 1.24 2.86

0.1066 0.2346 0.2588 0.0268 0.0653

0.06 0.00 0.00 0.94 0.00

2.66 3.92 5.00 1.06 2.36

4.02 × 103 8.68 × 103 1.83 × 105 2.39 × 103 3.54 × 103

Evaluating average rank alongside win rate prevents misleading conclusions where an estimator wins narrowly in select regimes but suffers catastrophic numerical degradation in others. Table 4, Table 5 and Table 6 shows the quantitative results for Truncated Gaussian model on S 10 , R10 + and Dirichlet model on ∆9 respectively. We observe the following 1. For the truncated Gaussian on S 10 , ϕ1 achieves the lowest error across configurations, securing a 100% win rate on precision matrix estimation (K) with an average rank of 1.00 and median MSE of 4.06 × 104 . On µ, it attains a 90% win rate with an average rank of 1.16 and median MSE of 0.0118, with ϕ2 ranking second (average rank 2.20). In contrast, h1 ranks last on estimation of K across all 50 configurations with a median MSE of 1.52 × 106 . 2. For the Dirichlet model on ∆9 , ϕ2 achieves the lowest error across configurations, obtaining a 96% win rate, an average rank of 1.04, and a median MSE of 0.0912. ϕ3 and ϕ1 follow with average ranks of 2.28 and 3.16, respectively. All three proposed estimators achieve lower average ranks than Truncated SM [23] (4.22), GSM CD with h(x) = x [45] (4.66), and Yu et al. [43] (5.64). 3. On R10 + , the baseline h1 achieves the lowest error, obtaining a 94% win rate on K (average rank 1.06, median MSE 2.39 × 103 ) and an 82% win rate on µ (average rank 1.24, median MSE 0.0268). Among our proposed estimators, ϕ1 achieves average ranks of 2.94 on µ and 2.66 on K, trailing h2 (x) = x2 (average rank 2.36 on K). ϕ3 ranks lowest on estimation of K with an average rank of 5.00 and a median MSE of 1.83 × 105 . K.4.2

Effect of Power Barrier Exponent

As noted in Section 6.2, the rate of attenuation at the boundary influences finite sample estimator performance. To isolate the impact of this decay rate, we study the power barrier family by systematically varying its exponent. For a general convex polytope Ω = {x ∈ Rd : a⊤ k x < bk , 1 ≤ k ≤ m}, let sk (x) = bk − a⊤ x denote the slack. We consider ϕ of the form k m

ϕ(x) =

X p 1 sk p(p − 1) k=1

Table 6: Ablation study evaluating the robustness of parameter estimators for the Dirichlet distribution concentration parameter α on the simplex ∆9 (d = 10) across 50 diverse ground-truth configurations (N = 800, 50 trials per configuration). We report the Win Rate, Average Rank, and Median MSE across all configurations.

M ETHODOLOGY

ϕ(x); h(x)

W IN R ATE

AVG . R ANK

M EDIAN MSE

OURS OURS OURS T RUNCATED SM [23] GSM CD [45] Y U ET AL . [43]

ϕ1 ϕ2 ϕ3 — h(x) = x h(x) = x

0.02 0.96 0.02 0.00 0.00 0.00

3.16 1.04 2.28 4.22 4.66 5.64

0.7411 0.0912 0.2661 1.5181 3.8273 8.9018

Estimation Error for µ vs p

Estimation Error for K vs p h1 (x) = x

MSE (log scale)

10 1

10 2

10 3

h1 (x) = x

105

104

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

Figure 5: Estimation error for µ (left) and K (right) as a function of the power barrier exponent p on R10 + (d = 10, N = 800). Solid curves and shaded regions display the median MSE and inter-quartile range (IQR) across 50 trials, respectively. The dashed horizontal line with shaded band indicates the median and IQR of the linear baseline h1 (x) = x [43]. where p ∈ (0, 2) \ {1}. The corresponding Hessian is Hϕ (x) =

m X sp−2 ak a⊤ k k k=1

Because the generator matrix scales inversely with the Hessian, enforcing boundary attenuation requires elements of Hessian to diverge near boundary, which implies p < 2. We evaluate this sweep on the 10 dimensional truncated Gaussian model on R10 + at sample size N = 800, across exponents p ∈ {0.2, 0.5, 0.8, 1.2, 4/3, 1.6, 1.8} over 50 independent trials. Figure 5 summarizes the results, reporting the median MSE alongside the inter-quartile range (IQR). As shown in the Figure 5, for smaller exponents (p ∈ {0.2, 0.5}), both µ and K suffer from substantial error degradation, with median MSE reaching 0.132 for µ and exceeding 105 for K. These exponents induce rapid generator decay near boundaries, resulting in broad IQR bands. Performance improves steadily once p > 1, where the chosen power barrier in experiments ϕ1 (p = 4/3) recovers reliable parameter estimates (median MSE of 0.0038 on µ and 8.43 × 103 on K), and p = 1.6 achieves the lowest error among all evaluated powers (median MSE of 0.0021 on µ and 5.68 × 103 on K), approaching the performance of the linear baseline h1 (x) = x [43]. However, as p approaches 2 (p = 1.8), error rebounds sharply on both parameters (median MSE rises to 0.0149 for µ and 1.17 × 104 for K) with widening of the IQR, indicating that under attenuating boundary samples reintroduces boundary noise into the GSM objective. K.4.3

Evaluation Under Non-Convex Support

In practical applications, observed data may be supported on non-convex domains. Although our theoretical framework assumes a convex support, we investigate the empirical behavior of our 40

estimator when applied to non-convex geometries via a convex relaxation. Specifically, we consider the non-convex domain Ω = {x ∈ Rd : ∥x∥ < a} ∪ {x ∈ Rd : b < ∥x∥ < c},

0 < a < b < c.

Because our methodology requires a convex domain, a natural heuristic is to evaluate the estimator on the convex hull of Ω, which corresponds to the open ball Ω′ = {x ∈ Rd : ∥x∥ < c}. We construct a logarithmic barrier potential directly on Ω′ : 2 ϕ(x) = − log c2 − ∥x∥ . We generate samples from a standard Gaussian density (µ0 = 0, K0 = Id ) truncated to Ω with radii a = 0.2, b = 0.5, and c = 1.0 in dimension d = 3. We compare our convex-hull estimator against Truncated Score Matching [23], which accounts for domain boundaries directly by weighting the objective. Estimation Error for µ vs Sample Size

Estimation Error for K vs Sample Size

103 102

101

MSE (log scale)

102

100

101

10 1 10 2 10 3

N = 50

N = 100

100

N = 200

Number of Samples

φ1 (x) = − log(c 2 − kxk 2 )

N = 50

N = 100 Number of Samples

N = 200

Truncated SM

Figure 6: Parameter estimation error for µ (left) and K (right) on the non-convex domain Ω with a = 0.2, b = 0.5, and c = 1.0 across sample sizes N ∈ {50, 100, 200} over 50 independent trials. Because the barrier ϕ is constructed with respect to the relaxed domain Ω′ , its generator Gϕ (x) does not vanish on the internal boundaries (∥x∥ = a and ∥x∥ = b), causing standard boundary conditions to fail. Consequently, we expect estimation quality to degrade. As shown in Figure 6, precision matrix estimation (K) degrades substantially, with the MSE plateauing around 150 across all sample sizes. Surprisingly, location parameter estimation (µ) remains accurate: the median MSE decreases from 0.0959 at N = 50 to 0.0079 at N = 200, outperforming Truncated Score Matching. To explain this asymmetric performance, we examine the integration by parts step that connects the tractable score matching loss to the weighted Fisher divergence: h i 1 2 LGSM (θ) = Ex∼p ∥∇ log pθ (x) − ∇ log p(x)∥Gϕ (x) ϕ 2 1 = Ex∼p ∇ log pθ (x)⊤ Gϕ (x)∇ log pθ (x) + Ex∼p [∇ · (Gϕ (x)∇ log pθ (x))] 2 − B(θ) + C, where C is a parameter-independent constant and B(θ) denotes the boundary integral over ∂Ω: Z B(θ) = p(x) n(x)⊤ Gϕ (x)∇ log pθ (x) ds. ∂Ω

Here ∂Ω comprises three concentric spheres: ∥x∥ = a, ∥x∥ = b, and ∥x∥ = c. On the outer boundary ∥x∥ = c, Gϕ (x) → 0 by construction, so its contribution vanishes. The remaining boundary integral evaluates over the internal surfaces: Z Z x ⊤ x ⊤ B(θ) = p(x) Gϕ (x)∇ log pθ (x) ds − p(x) Gϕ (x)∇ log pθ (x) ds, a b ∥x∥=a ∥x∥=b 41

where the sign inversion on ∥x∥ = b reflects the inward-pointing normal relative to the outer annulus. For a standard Gaussian centered at the origin, the density p(x) ∝ exp − 21 ∥x∥ ⊤

is constant on

⊤

any sphere ∥x∥ = r. Furthermore, x Gϕ (x) = γ(r)x for a scalar constant γ(r). Substituting the score ∇ log pθ (x) = −K(x − µ), the integrand on each spherical shell becomes proportional to x⊤ K(x − µ) = x⊤ Kx − x⊤ Kµ. By symmetry (x 7→ −x), the linear term integrates to zero: Z x⊤ Kµ ds = 0. ∥x∥=r ⊤

In contrast, the quadratic term x Kx is strictly positive and non-zero: Z r2 Tr(K) x⊤ Kx ds = Area(Sd−1 ) ̸= 0. r d ∥x∥=r Consequently, the boundary term simplifies to B(θ) = −CK Tr(K), for a positive constant CK independent of µ. This yields an important insight — the boundary discrepancy does not depend on µ, meaning ∇µ B(θ) = 0. Therefore, the gradient of our sample objective with respect to µ remains an asymptotically unbiased estimator of the true Fisher divergence gradient, allowing µ̂N → µ0 as N → ∞. Conversely, ∇K B(θ) ̸= 0, which introduces an irreducible asymptotic bias into the loss for K. This explains why K exhibits non-convergent error, while µ converges stably. In comparison, Truncated Score Matching [23] uses the exact boundary distance weighting to ensure boundary vanishing across all bounding surfaces, which preserves theoretical consistency for both µ and K. While this consistency is reflected in the steady decrease of its precision error with sample size, Truncated Score Matching exhibits noticeably higher estimation error and wider spread for µ at small sample sizes compared to our approach. We hypothesize that this could be due to finite-sample gradient fluctuations arising from the non-smooth Euclidean distance function; however, a thorough investigation of this behavior is beyond the scope of this work. Remark K.1. This decoupling between µ and K relies on the joint spherical symmetry of the domain cuts and the centered density (µ0 = 0). If the true mean were translated away from the origin, p(x) would vary across the internal boundaries, breaking the cancellation and coupling µ to the non-zero boundary integral. Nonetheless, this experiment highlights that boundary attenuation is a sufficient, rather than strictly necessary, condition for recovering subset parameters in constrained settings. K.5

Data Driven Ranking of Candidate Choices of ϕ

In practical setting, a practitioner must select a suitable potential function ϕ without access to the ground-truth parameter. To address this, we propose a diagnostic procedure to rank candidate choices of ϕ purely based on the observed data. A natural starting point is Theorem 5.1, which establishes −1 that the asymptotic variance of our estimator is given by Γ−1 0 Σ0 Γ0 . Intuitively, a potential ϕ that induces a higher asymptotic variance is generally less desirable. Computing this variance directly presents two issues: it relies on population expectations, and it requires the unknown true parameter θ0 . We can circumvent these issues by replacing the population expectations with their finite sample empirical estimates, and substituting the true parameter with the estimated parameter θ̂N . While this yields a practical heuristic rather than an exact finite sample bound, it remains firmly grounded in the asymptotic theory of our framework. b denote this sample-estimated asymptotic covariance matrix. We propose scoring candidate Let Σ b which we refer to as the estimated total variance. A potentials using the trace of this matrix, Tr(Σ), lower estimated total variance suggests a more statistically efficient choice of ϕ. As an alternative heuristic, we consider the geometric curvature of the objective. Because the empirical loss is quadratic in θ (true in case of exponential family, cf. Equation 18), its Hessian is exactly ΓN . We propose ranking the candidates using the condition number of this Hessian, denoted as κ(ΓN ). A lower condition number implies a better conditioned optimization landscape. To evaluate these heuristics, 42

we ran ranking experiments on our primary settings: the truncated Gaussian model on S 10 and R10 +, and the Dirichlet model on the simplex ∆9 . We fixed the sample size to N = 800 for these diagnostic tests; because our proposed methodology evaluates gradients at the empirical estimate θ̂N , ranking criteria computed on very small sample sizes may be excessively noisy and unreliable. To quantify the quality of our proposed scores, we compare the data-driven rankings against the true rankings, which are determined by the actual parameter Mean Squared Error (MSE) computed using the ground truth. We report Top-1 Accuracy (the proportion of trials where the heuristic successfully identifies the candidate with the lowest MSE) alongside Worst Regret and Mean Regret to measure the relative error penalty incurred when a suboptimal candidate is selected. We also report both Spearman (ρ) and Kendall (τ ) rank correlation coefficients between the predicted and true ranks. For the truncated Gaussian model, we separate the ranking evaluation for µ and K. Table 7: Evaluation of data-driven candidate ranking heuristics across 50 independent trials at N = 800. Top-1 Accuracy measures how often the heuristic selects the MSE-optimal potential. Worst and Mean Regret quantify the relative sub optimality when a misselection occurs. Rank correlations (ρ and τ ) evaluate monotonic agreement with ground-truth error ordering. D OMAIN

PARAMETER

S CORE

T OP -1 ACC .

W ORST R EGRET

M EAN R EGRET

S 10

b T R(Σ) κ(ΓN )

0.980 0.980

0.2714 0.2714

0.0054 0.0054

0.9900 0.9900

0.9867 0.9867

b T R(Σ) κ(ΓN )

0.780 0.780

245.8676 245.8676

5.5259 5.5259

0.5700 0.5700

0.5600 0.5600

b T R(Σ) κ(ΓN )

1.000 1.000

0.0000 0.0000

1.0000 1.0000

b T R(Σ) κ(ΓN )

0.980 0.980

0.6111 0.6111

0.0122 0.0122

0.9500 0.9500

0.9467 0.9467

b T R(Σ) κ(ΓN )

0.580 0.080

10.5635 30.1927

0.7870 3.9781

0.3800 −0.1100

0.3200 −0.1200

R10 +

∆9

b Table 7 summarizes the results of this evaluation. In general, the estimated total variance metric Tr(Σ) performs well, successfully identifying the most stable choices of ϕ with high accuracy and positive rank correlations across all evaluated domains. The condition number κ(ΓN ) performs comparably well on the truncated Gaussian settings but fails on the Dirichlet model (yielding only 8% Top-1 accuracy and negative rank correlation). This indicates that while the condition number effectively captures the curvature of the loss function, it is not a universally reliable metric for ranking ϕ across diverse geometries. For the Dirichlet model, the estimated total variance yields a moderate Top-1 Accuracy of 58%. Upon inspecting the individual trial logs, we observed that the performance gap between ϕ1 and ϕ2 in this specific setting is minimal. Consequently, while the diagnostic occasionally selects the second best potential which dampens the Top-1 accuracy score — the actual performance penalty for doing so is marginal, resulting in relatively low overall regret. Most importantly, the estimated total variance reliably rejected highly unstable choices across the evaluated trials. Remark K.2. It should be noted that this diagnostic approach becomes computationally infeasible in b requires computing very high-dimensional settings, as computing the estimated total variance Tr Σ the inverse of ΓN .

K.6

Implicit Variational Autoencoders with Score Matching

To demonstrate that the proposed framework is applicable to generative modeling, beyond parameter estimation problems, we consider implicit variational autoencoders (VAEs) [38], where score matching naturally arises in the training procedure. We emphasize the novelty of this evaluation: to the best of our knowledge, we are the first to apply the generalized score matching objective, rigorously derived in Theorem 3.1, to this specific generative problem. We demonstrate how the 43

proposed objective can be scaled to high-dimensional settings. Consider the objective in Equation 4: " # 1 ⊤ ⊤ L(θ) = E ∇ log pθ (x) D(x)∇ log pθ (x) + Tr (D(x)Hlog pθ (x)) + ∇ log pθ (x) (∇ · D)(x) x∼p 2 While this formulation circumvents the computation of the partition function Z(θ), it introduces an additional computational challenge, namely computation of the trace of D(x)Hlog pθ (x). A practical alternative is to consider an unbiased estimator of the trace term. This is typically done using Hutchinson-style trace estimator [9, 17]. Such estimators have been used gainfully in several prior works [27, 34]. The key idea is to estimate the trace by considering projections onto random directions and computing the corresponding quadratic form. In particular, we consider " 1 L̂(θ) = x∼p, E ∇ log pθ (x)⊤ D(x)∇ log pθ (x) + v ⊤ D(x)Hlog pθ (x)v 2 v∼pv # + ∇ log pθ (x)⊤ (∇ · D)(x)

(31)

where pv denotes the sampling distribution of the random direction vectors v, chosen independently of x and satisfying E[vv ⊤ ] = I. If we choose D(x) = I in Equation 31, then we recover the original score matching objective, with the trace term replaced by corresponding Hutchinson estimator. This resulting objective coincides with the variance-reduced version of sliced score matching (SSM-VR) considered by Song et al. [38]. For implicit VAEs, we follow the training methodology of Song et al. [38]. We replace the SSM-VR objective used to train the score network with the GSM objective given by Equation 31, which we refer to as Hutchinson-GSM (H-GSM). The model architecture, hyperparameter choices, and training procedure are identical to those followed by Song et al. [38]. In the implementation, the prior of the latent is chosen to be a standard Gaussian supported on Rd . Therefore, we work with the unconstrained GSM objective derived in Theorem 3.1. This allows us to directly explore the effect of different choices of D. Specifically, we consider three choices of D with D1 (x) = diag(σ(x)), D2 (x) = diag(exp(− |x|)), and D3 (x) = diag(exp(−x2 )), where σ(x) denotes the element-wise sigmoid function, and |x|, x2 and exp(x) denote element-wise absolute value, square and exponential operations, respectively. We compare implicit VAEs trained using SSM-VR and the H-GSM objective on MNIST and CelebA datasets. Following the evaluation protocol of Song et al. [38], we report the negative log-likelihood (NLL) for MNIST and Fréchet Inception Distance (FID) for CelebA over various choices of D. The results for MNIST are summarized in Table 8. As observed, H-GSM performs competitively with SSM-VR, in particular NLL for D1 (x) = diag(σ(x)) with latent dimension 8, and NLL for D2 (x) = diag(exp(− |x|)) with latent dimension 32. In both cases, the H-GSM objective exhibits higher variance across runs than SSM-VR. The results for CelebA are summarized in Table 9. We observe that H-GSM with D3 (x) = diag(exp(−x2 )) performs competitively with the SSM-VR baseline, with its FID approaching that of SSM-VR as training progresses. The D2 (x) = diag(exp(− |x|)) variant is closer to the baseline SSM-VR than D1 (x) = diag(σ(x)). In terms of variability across runs, D3 exhibits a variance comparable to or lower than that of SSM-VR at later checkpoints, whereas D1 and D2 generally exhibit higher variance. These results further highlight the sensitivity of H-GSM performance to the choice of the diagonal matrix D(x). We provide uncurated generated samples for CelebA and MNIST in Figure 7 and Figure 8, respectively. Remark K.3. The SSM-VR results reported in Table 8 and Table 9 were obtained by re-running the original implementation provided by Song et al. [38] under the same experimental settings. Remark K.4. This experiment illustrates an important distinction between Theorem 3.1 and Theorem 3.2. In the unconstrained setting of Theorem 3.1, one can specify the weighting matrix D directly, without requiring an underlying convex function ϕ. By contrast, in the convex domain setting of Theorem 3.2, the generator is induced by the choice of ϕ. In this sense, the unconstrained formulation offers greater flexibility in practice but requires stronger assumptions. 44

Table 8: NLL comparison of implicit VAE training objectives across different latent dimensions averaged over 5 runs. The expressions σ(x), exp(− |x|), and exp(−x2 ) are evaluated element-wise. SSM-VR 8 32

H-GSM

D1 (x) = DIAG(σ(x))

D2 (x) = DIAG(exp(− |x|))

D3 (x) = DIAG(exp(−x2 ))

96.17 ± 0.10 89.31 ± 0.12

96.72 ± 0.39 89.81 ± 0.20

99.61 ± 1.95 90.37 ± 0.65

130.02 ± 6.05 95.91 ± 2.34

L ATENT D IMENSION

Table 9: FID comparison of implicit VAE training objectives on CelebA dataset at different training checkpoints averaged over 5 runs. The latent dimension is fixed at 32. The expressions σ(x), exp(− |x|), and exp(−x2 ) are evaluated element-wise. SSM-VR I TERATION 10K 20K 30K 40K 50K 60K 70K 80K 90K 100K

H-GSM

D1 (x) = DIAG(σ(x))

D2 (x) = DIAG(exp(− |x|))

D3 (x) = DIAG(exp(−x2 ))

101.07 ± 1.54 83.10 ± 1.10 76.34 ± 1.45 72.12 ± 1.38 68.89 ± 0.69 67.80 ± 1.70 66.94 ± 2.06 65.70 ± 1.97 64.49 ± 1.48 64.10 ± 1.82

140.41 ± 10.13 119.66 ± 11.49 106.44 ± 10.34 100.20 ± 11.57 93.31 ± 11.75 87.77 ± 10.27 85.85 ± 10.72 82.20 ± 9.48 79.80 ± 7.57 78.11 ± 7.15

122.01 ± 6.03 101.79 ± 5.88 87.57 ± 5.91 80.66 ± 5.96 75.77 ± 6.27 73.22 ± 5.90 71.21 ± 5.74 69.48 ± 6.55 67.29 ± 4.67 66.60 ± 5.06

121.51 ± 2.43 99.75 ± 5.68 86.45 ± 4.99 79.82 ± 3.94 74.31 ± 1.50 70.90 ± 2.14 68.87 ± 1.71 66.97 ± 1.13 65.37 ± 1.31 64.15 ± 0.53

Remark K.5. Although the experiments in this section are based on the unconstrained objective over Rd , the scalability discussion based on the Hutchinson estimator applies equally to the generalized score matching objective for densities supported on convex domains (Equation 11).

SSM-VR

H-GSM (D1 )

H-GSM (D2 )

H-GSM (D3 )

<latexit sha1_base64="npLeoRJF130/nAc6FkDMFsrSol0=">AAAB9HicbVDLTgJBEOzFF+IL9ehlIpjgQbLLAT0SNZGLCUZ5JLAhs8MsTJidXWdmSQjhO7x40Bivfow3/8YB9qBoJZ1UqrrT3eVFnClt219WamV1bX0jvZnZ2t7Z3cvuHzRUGEtC6yTkoWx5WFHOBK1rpjltRZLiwOO06Q2vZn5zRKVioXjQ44i6Ae4L5jOCtZHc6tnN/S0q5K+7Tv60m83ZRXsO9Jc4CclBglo3+9nphSQOqNCEY6Xajh1pd4KlZoTTaaYTKxphMsR92jZU4IAqdzI/eopOjNJDfihNCY3m6s+JCQ6UGgee6QywHqhlbyb+57Vj7V+4EyaiWFNBFov8mCMdolkCqMckJZqPDcFEMnMrIgMsMdEmp4wJwVl++S9plIpOuVi+K+Uql0kcaTiCYyiAA+dQgSrUoA4EHuEJXuDVGlnP1pv1vmhNWcnMIfyC9fENapyP7g==</latexit>

<latexit sha1_base64="oyiiQCr00O+ZtH+bhYUbkjcsNIU=">AAAB7XicbVA9TwJBEJ3zE/ELtbTZSExsJHcUaEm0sTFBkYMELmRv2YOVvd3L7p4JIfwHGwuNsfX/2PlvXOAKBV8yyct7M5mZFyacaeO6387K6tr6xmZuK7+9s7u3Xzg49LVMFaENIrlUrRBrypmgDcMMp61EURyHnDbD4fXUbz5RpZkUD2aU0CDGfcEiRrCxkl+v3577991C0S25M6Bl4mWkCBlq3cJXpydJGlNhCMdatz03McEYK8MIp5N8J9U0wWSI+7RtqcAx1cF4du0EnVqlhyKpbAmDZurviTGOtR7Foe2MsRnoRW8q/ue1UxNdBmMmktRQQeaLopQjI9H0ddRjihLDR5Zgopi9FZEBVpgYG1DehuAtvrxM/HLJq5Qqd+Vi9SqLIwfHcAJn4MEFVOEGatAAAo/wDK/w5kjnxXl3PuatK042cwR/4Hz+AKOijoc=</latexit>

<latexit sha1_base64="LzNOoApJyKuhk/ngPZTEb9zQJB4=">AAAB9HicbVDLTgJBEOzFF+IL9ehlIpjgQbLLAT0SNZGLCUZ5JLAhs8MsTJidXWdmSQjhO7x40Bivfow3/8YB9qBoJZ1UqrrT3eVFnClt219WamV1bX0jvZnZ2t7Z3cvuHzRUGEtC6yTkoWx5WFHOBK1rpjltRZLiwOO06Q2vZn5zRKVioXjQ44i6Ae4L5jOCtZHc6tnN/S0q5K+7pfxpN5uzi/Yc6C9xEpKDBLVu9rPTC0kcUKEJx0q1HTvS7gRLzQin00wnVjTCZIj7tG2owAFV7mR+9BSdGKWH/FCaEhrN1Z8TExwoNQ480xlgPVDL3kz8z2vH2r9wJ0xEsaaCLBb5MUc6RLMEUI9JSjQfG4KJZOZWRAZYYqJNThkTgrP88l/SKBWdcrF8V8pVLpM40nAEx1AAB86hAlWoQR0IPMITvMCrNbKerTfrfdGaspKZQ/gF6+MbbCKP7w==</latexit>

<latexit sha1_base64="DMiw4Hj8EnnHJxSIpS+4cF0fhzo=">AAAB9HicbVDLTgJBEOzFF+IL9ehlIpjgQbKLCXokaiIXE4zySGBDZodZmDD7cGaWhGz4Di8eNMarH+PNv3GAPShYSSeVqu50dzkhZ1KZ5reRWlldW99Ib2a2tnd297L7Bw0ZRILQOgl4IFoOlpQzn9YVU5y2QkGx53DadIbXU785okKywH9U45DaHu77zGUEKy3Z1bPbhztUyN90z/On3WzOLJozoGViJSQHCWrd7FenF5DIo74iHEvZtsxQ2TEWihFOJ5lOJGmIyRD3aVtTH3tU2vHs6Ak60UoPuYHQ5Ss0U39PxNiTcuw5utPDaiAXvan4n9eOlHtpx8wPI0V9Ml/kRhypAE0TQD0mKFF8rAkmgulbERlggYnSOWV0CNbiy8ukUSpa5WL5vpSrXCVxpOEIjqEAFlxABapQgzoQeIJneIU3Y2S8GO/Gx7w1ZSQzh/AHxucPbaiP8A==</latexit>

Figure 7: For the CelebA dataset, we compare the visual sample quality of implicit VAEs trained using the proposed H-GSM objective for various choices of D against the SSM-VR baseline [38]. The generated samples are of comparable or superior visual quality with respect to the baseline (SSM-VR) across the choices of D(x) with D1 (x) = diag(σ(x)), D2 (x) = diag(exp(− |x|)), and D3 (x) = diag(exp(−x2 )). The latent dimension is fixed at 32 for all settings shown here. A quantitative evaluation of these models via Fréchet Inception Distance (FID) is provided in Table 9.

Latent Dimension 8 <latexit sha1_base64="An/Iziti1Q/gVZwF1LhdcI1Y+Xo=">AAAB/XicbVC7TgJBFL2LL8QH66OzmQgmVmSXAimJWlhYYCKPBDZkdpiFCbOPzMya4Ib4KzYWGmPrf9j5N87CFgqeZJKTcx9z7nEjzqSyrG8jt7a+sbmV3y7s7O7tF82Dw7YMY0Foi4Q8FF0XS8pZQFuKKU67kaDYdzntuJOrtN55oEKyMLhX04g6Ph4FzGMEKy0NzONbrGig0DXzaZB2oXK9PDBLVsWaA60SOyMlyNAcmF/9YUhivUIRjqXs2VaknAQLxQins0I/ljTCZIJHtKdpgH0qnWTufobOtDJEXij0007m6u+JBPtSTn1Xd/pYjeVyLRX/q/Vi5dWdhAVRrE8ki4+8mCMVojQKNGSCEsWnmmAimPaKyBgLTJQOrKBDsJdPXiXtasWuVWp31VLjMosjDydwCudgwwU04Aaa0AICj/AMr/BmPBkvxrvxsWjNGdnMEfyB8fkDARuUSg==</latexit>

Latent Dimension 32 <latexit sha1_base64="x2Uj9xYdV7cAlDWcVaN6IOqRsTk=">AAAB/nicbVDLSsNAFL2pr1pfUXHlZrAVXJUkQnVZ1IULFxXsA9pQJtNJO3TyYGYilFDwV9y4UMSt3+HOv3HSZqGtBwYO5z7m3OPFnEllWd9GYWV1bX2juFna2t7Z3TP3D1oySgShTRLxSHQ8LClnIW0qpjjtxILiwOO07Y2vs3r7kQrJovBBTWLqBngYMp8RrLTUN4/usKKhQjcsoGHWhSrnTqVvlq2qNQNaJnZOypCj0Te/eoOIJHqHIhxL2bWtWLkpFooRTqelXiJpjMkYD2lX0xAHVLrpzP4UnWplgPxI6KetzNTfEykOpJwEnu4MsBrJxVom/lfrJsq/dFMWxom+kcw/8hOOVISyLNCACUoUn2iCiWDaKyIjLDBROrGSDsFePHmZtJyqXavW7p1y/SqPowjHcAJnYMMF1OEWGtAEAik8wyu8GU/Gi/FufMxbC0Y+cwh/YHz+AHFKlIE=</latexit>

SSM-VR <latexit sha1_base64="oyiiQCr00O+ZtH+bhYUbkjcsNIU=">AAAB7XicbVA9TwJBEJ3zE/ELtbTZSExsJHcUaEm0sTFBkYMELmRv2YOVvd3L7p4JIfwHGwuNsfX/2PlvXOAKBV8yyct7M5mZFyacaeO6387K6tr6xmZuK7+9s7u3Xzg49LVMFaENIrlUrRBrypmgDcMMp61EURyHnDbD4fXUbz5RpZkUD2aU0CDGfcEiRrCxkl+v3577991C0S25M6Bl4mWkCBlq3cJXpydJGlNhCMdatz03McEYK8MIp5N8J9U0wWSI+7RtqcAx1cF4du0EnVqlhyKpbAmDZurviTGOtR7Foe2MsRnoRW8q/ue1UxNdBmMmktRQQeaLopQjI9H0ddRjihLDR5Zgopi9FZEBVpgYG1DehuAtvrxM/HLJq5Qqd+Vi9SqLIwfHcAJn4MEFVOEGatAAAo/wDK/w5kjnxXl3PuatK042cwR/4Hz+AKOijoc=</latexit>

H-GSM (D1 ) <latexit sha1_base64="npLeoRJF130/nAc6FkDMFsrSol0=">AAAB9HicbVDLTgJBEOzFF+IL9ehlIpjgQbLLAT0SNZGLCUZ5JLAhs8MsTJidXWdmSQjhO7x40Bivfow3/8YB9qBoJZ1UqrrT3eVFnClt219WamV1bX0jvZnZ2t7Z3cvuHzRUGEtC6yTkoWx5WFHOBK1rpjltRZLiwOO06Q2vZn5zRKVioXjQ44i6Ae4L5jOCtZHc6tnN/S0q5K+7Tv60m83ZRXsO9Jc4CclBglo3+9nphSQOqNCEY6Xajh1pd4KlZoTTaaYTKxphMsR92jZU4IAqdzI/eopOjNJDfihNCY3m6s+JCQ6UGgee6QywHqhlbyb+57Vj7V+4EyaiWFNBFov8mCMdolkCqMckJZqPDcFEMnMrIgMsMdEmp4wJwVl++S9plIpOuVi+K+Uql0kcaTiCYyiAA+dQgSrUoA4EHuEJXuDVGlnP1pv1vmhNWcnMIfyC9fENapyP7g==</latexit>

H-GSM (D2 ) <latexit sha1_base64="LzNOoApJyKuhk/ngPZTEb9zQJB4=">AAAB9HicbVDLTgJBEOzFF+IL9ehlIpjgQbLLAT0SNZGLCUZ5JLAhs8MsTJidXWdmSQjhO7x40Bivfow3/8YB9qBoJZ1UqrrT3eVFnClt219WamV1bX0jvZnZ2t7Z3cvuHzRUGEtC6yTkoWx5WFHOBK1rpjltRZLiwOO06Q2vZn5zRKVioXjQ44i6Ae4L5jOCtZHc6tnN/S0q5K+7pfxpN5uzi/Yc6C9xEpKDBLVu9rPTC0kcUKEJx0q1HTvS7gRLzQin00wnVjTCZIj7tG2owAFV7mR+9BSdGKWH/FCaEhrN1Z8TExwoNQ480xlgPVDL3kz8z2vH2r9wJ0xEsaaCLBb5MUc6RLMEUI9JSjQfG4KJZOZWRAZYYqJNThkTgrP88l/SKBWdcrF8V8pVLpM40nAEx1AAB86hAlWoQR0IPMITvMCrNbKerTfrfdGaspKZQ/gF6+MbbCKP7w==</latexit>

H-GSM (D3 ) <latexit sha1_base64="DMiw4Hj8EnnHJxSIpS+4cF0fhzo=">AAAB9HicbVDLTgJBEOzFF+IL9ehlIpjgQbKLCXokaiIXE4zySGBDZodZmDD7cGaWhGz4Di8eNMarH+PNv3GAPShYSSeVqu50dzkhZ1KZ5reRWlldW99Ib2a2tnd297L7Bw0ZRILQOgl4IFoOlpQzn9YVU5y2QkGx53DadIbXU785okKywH9U45DaHu77zGUEKy3Z1bPbhztUyN90z/On3WzOLJozoGViJSQHCWrd7FenF5DIo74iHEvZtsxQ2TEWihFOJ5lOJGmIyRD3aVtTH3tU2vHs6Ak60UoPuYHQ5Ss0U39PxNiTcuw5utPDaiAXvan4n9eOlHtpx8wPI0V9Ml/kRhypAE0TQD0mKFF8rAkmgulbERlggYnSOWV0CNbiy8ukUSpa5WL5vpSrXCVxpOEIjqEAFlxABapQgzoQeIJneIU3Y2S8GO/Gx7w1ZSQzh/AHxucPbaiP8A==</latexit>

Figure 8: For the MNIST dataset, we compare the visual quality of samples generated by implicit VAEs trained using the proposed H-GSM objective for various choices of D against the SSM-VR baseline [38]. The proposed generalized score matching-based estimator consistently generates visual samples of comparable fidelity to the baseline across various latent dimensions and choices of D(x) with D1 (x) = diag(σ(x)), D2 (x) = diag(exp(− |x|)), and D3 (x) = diag(exp(−x2 )). Across all methods, latent dimension 8 outperforms latent dimension 32 and this is consistent with the results reported by Song et al. [38]. 47

Record · ID 673532 · SHA-256 f6aacbdb9b2e4dc1
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.