Direct Fisher Score Estimation for Likelihood Maximization
Authors: Sherman Khoo
Yakun Wang
Song Liu
Mark Beaumont
Venue: NeurIPS 2025 spotlight ()
Keywords: Simulation-based inference
Likelihood-free inference
Score Matching
Gradient-Based Optimization
Likelihood Maximization
Primary area: probabilistic_methods

Abstract: We study the problem of likelihood maximization when the likelihood function is intractable but model simulations are readily available. We propose a sequential, gradient-based optimization method that directly models the Fisher score based on a local score matching technique which uses simulations from a localized region around each parameter iterate. By employing a linear parameterization for the surrogate score model, our technique admits a closed-form, least-squares solution. This approach yields a fast, flexible, and efficient approximation to the Fisher score, effectively smoothing the likelihood objective and mitigating the challenges posed by complex likelihood landscapes. We provide theoretical guarantees for our score estimator, including bounds on the bias introduced by the smoothing. Empirical results on a range of synthetic and real-world problems demonstrate the superior performance of our method compared to existing benchmarks.

Direct Fisher Score Estimation for Likelihood Maximization
Sherman Khoo
Yakun Wang
Song Liu"
Mark Beaumont
Abstract
We study the problem of likelihood maximization when the likelihood function is intractable but model simulations are readily available. We propose sequential, gradient- based optimization method that directly models the Fisher score based on local score matching technique which uses simulations from localized region around each parameter iterate. By employing linear parameterization for the surrogate score model, our technique admits closed-form, least-squares solution. This approach yields fast, flexible, and efficient approximation to the Fisher score, effectively smoothing the likelihood objective and mitigating the challenges posed by complex likelihood landscapes. We provide theoretical guarantees for our score estimator, including bounds on the bias introduced by the smoothing. Empirical results on a range of synthetic and real-world problems demonstrate the superior performance of our method compared to existing benchmarks.
1 Introduction
Implicit simulator- models are now routine in many scientific fields, such as biology CCsillryy al. 2010]. cosmology [Schafer and Freeman, 2012], neuroscience (Sterratt
2011]. engineering [Bharti et al. 2021]. and other scientific applications [Toni et al., 2009]] In traditional statistical models, there a prescribed probabilistic model, which provides an explicit parameterization for the data distribution, allowing development of the likelihood function for further inference. In contrast, simulation-bal based models define the distribution implicitly through a computational simulator. Thus, while simulation of the data for various parameter settings is possible, the probability density function of the data or the likelihood function is often unavailable in closed form. This problem setting is known likelihood-fiee inference
iimulation-basee inference (SBI) [Cranmer al. 2020].
Traditionally, methods this setting have focused on a Bayesian inference technique known as approximate Bayesian computation (ABC) [Beaumont et al. 22022. Fundamentally. the ABC method builds an approximation to the Bayesian posterior distribution by drawing parameter samples from the prior distribution, generating datasets from the drawn parameter values, and filtering parameter values through a rejection algorithm based on the distance of the generated summary statistic of the dataset from the observations. In recent years, there has been rise of generative modeling unsupervised learning in machine learning, which aims to recover a data distribution given a set of samples. Such generative models are often built from neural networks [Mohamed and Lakshminarayanan 2017]. and given the fundamental similarity with the SBI problem setting, this has led to significant cross-pollination across the two fields, with the development of many SBI inference methods using generative neural networks [Durkan et al. 2018 Papamakarios et al. 2017].
While significant progress in SBI has come from Bayesian approaches, such methods are often computationally demanding; furthermore, many scientific disciplines retain preference towards
*School Mathematics, University of Bristol
School of Biological Sciences. University of Bristol
Correspondence to: Sherman Khoo bineerman.hooo @ bristol..c.uk>>
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

maximum likelihood approaches. In contrast, practitioners who favor faster point estimation through maximum likelihood lack comparably mature tools. To close this gap. we propose fast, simulation- efficient, and robust gradient-basee base technique for estimating the maximum likelihood for SBI. We build on a popular technique in generative modeling known as score matching (Hyvärinen and Dayan 2005], which has seen significant use in score- score-hased generative models (Song and Ermon| 2019]. now a cornerstone for many state-of-the-art approaches in generative modeling (Yang et al. 2023]. We adapt score matching to estimate the Fisher score, that is, the gradient of the log-likelihood function with respect to the parameters, within a localized region. This estimated gradient can then be used in any first-order gradient-basee -aseed stochastic optimization algorithm such stochastic gradient descent (SGD) to obtain an approximate maximum likelihood estimator (MLE) and serves as a potential avenue for uncertainty quantification for the MLE through the empirical Fisher information matrix.
Our contributions in this work are as follows.
We propose lightweight, simulation-efficn estimation of simulator models based on a novel local Fisher score matching technique
and robust method for maximum likelihood
We derive theory for our local Fisher score matching technique and establish a connection with the Gaussian smoothing gradient estimator, offering a unifying perspective for zeroth- order optimization techniques and likelihood optimization for SBI
We demonstrate the effectiveness of our method in real-world experiments for applied machine learning and cosmology problems, showcasing both its efficiency and robust performance compared to existing approaches
2 Background
2.1 Score Matching
Density estimation is the problem of learning a data distribution px(x) using only an observed dataset, x PD. An approach to this problem is learning the density with energy-based model (EBM), which parameterizes the model through its scalar-valued energy function Eg I + R where x E X' C RR, givingg the model density (x(x) exp(-Eo (xp(-Ex(x))/3) (x))/Zo.
Since the energy function is an unnormalized density function, it can be flexibly parameterized, usually through a neural network. However, note that the normalization constant Zo- = cxp( Ea (x))xxxxx is still function of 0, and therefore will still need to be computed training the EBM through standard likelihood maximization. Since this multidimensional integral is often intractable and requires costly approximation method, score matching [Hyvärinen and Dayan 2005] often used to bypass the computation of the normalization constant. This is done by considering
alternate training objective instead of MLE, based on the score function $9 _² + R where Sp = Vx log pa = V Eg In fact, since equivalence in the score amounts to equivalence in the distribution, matching the scores is equivalent to performing density estimation. One starting point is the explicit score matching objective (ESM). Defining the gradient operator on a scalar-valued function as Vx == (22, Zagg where Bar
9 is the partial derivative operator for (I1...., IR), and the Jacobian operator vector- li.f we have:
valued function f : Rm Rn as J.,
0f Loz
(f) = Ex Ex>po(x), 2 [15o(x) - Vx loggpD(Xx)]
However, this objective is not tractable due to the need to evaluate Vx (x) Hence, this objective is transformed to:
LESM(O) (0) = (xx
2o2 (x)]l + tr (Jxss(x)) + (constants w.r.t. 8)
Although this objective can be directly estimated, and thus optimized and used in the training of an EBM, it is computationally expensive due to the presence of the Jacobian term, motivating further extensions to the standard score matching objective, such as the denoising score matching objective Vincent 2011 and the sliced score matching objective (Song al. 2020
2

---

2.2 Maximum Likelihood Estimation and Fisher Score
Maximum likelihood estimation (MLE) is foundational tool in statistical inference, under standard regularity conditions. it is consistent and asymptotically efficient [Casella and Berger [024 Section 10J. Central to the MLE is the Fisher score, defined as the gradient of the log-likeldoood with respect to the parameters, Vo log p(x 0). From an optimization point of view, the score provides the direction of steepest ascent of the log-likelihood MLE approaches. From inferential point of view, the covariance of the Fisher score is equal to the Fisher information matrix (FIM), which, through the Cramér-Rao lower bound [Rao 1992]. lower bounds the variance of any unbiased estimator. Furthermore, the distribution of the MLE is asymptotically normal with covariance equal to the inverse the FIM, which underpins Wald-type confidence intervals and hypothesis tests (Van der Vaart 2000 Section 5].
parameter space, and thus drives gradient-based
2.3 Notation and Problem Setup
We consider a statistical model where the data x E X C Rdx are generated from a distribution Pa parameterized by 8 Ω C Rd In the simulation-base implicitly defined, so we can draw samples from this model for any choice of & but the closed-form expression for the probability density function, and hence the likelihood function is not known.
inference setting, this statistical model is
Given set of independent and identically distributed observations, D = {x, J drawn from the true data-generating process X, Pe-, where 0" denotes the true parameter, the maximum likelihood
estimator is MMLE arg max p(D |0)-
As the likelihood function L(0 D) I p(x; 0) is not available for SBI models, typical likelihood maximization cannot be applied directly. We thus propose a Fisher score matching estimator, ÓFSM. Our method is fundamentally a first-order optimization approach, and our main focus on the direct estimation of the gradient of the log-likelihood function at each parameter iteration, which is done with a novel local Fisher score matching objective. We first discuss our Fisher score estimation technique in Section before proceeding with the MLE procedure in Section
Likelihood-free Fisher Score Estimation
Score matching (Hyvärinen and Dayan 2005 is classical method in density estimation, but is not directly applicable in likelihood gradient maximization, as it typically targets the Stein score, i.e., the gradient with respect the data Vx log (x() instead of the Fisher score, which is the gradient with respect to the parameters V log py (x). Hence, we propose to adapt score matching into a novel local Fisher score estimation technique which estimates the gradient of the log-likelihood for fixed
parameter point O, at any data sample x. V Vol(0; x) x0
V log log
0=0.
3.1 Local Fisher Score Matching Objective
Around the target parameter point 6,, we introduce local proposal distribution q(O 0t), which we typically take as an isotropic Gaussian distribution, q(O 0.) = ((0,,o2)) When combined with the statistical model Po. we induce a joint distribution in both the data and parameter space that has probability density p(x ()q(O 0z). Note that by drawing parameter samples from the local proposal distribution and then drawing corresponding data samples for the parameter samples,
can easily draw samples from this joint distribution.
To estimate the score function, we use score model Sw R2 + Rdo where Sw (x) has parameters W. Our starting point is the adapted, localized score matching least-squares loss for the Fisher score.
;.) (xx-p(x10),0001100 V. logp(x | Sw (x)))]
As we are within the simulation-based inference framework, we do not have closed form expression for the Fisher score Vo logp(x () and hence this objective function not tractable, We first expand the square of Equation (L) which allows us to rewrite J (W; 0x) as:
3

---

(WW:8)) 00)0
||| Sw S Sw (x))
Vo logp(x 0) (constants w.r.t. W)
We focus on the cross-term, lEx (1xe=p(((000000(((00) [Sw (x) Ve (x 0)] Using an integration-by-
parts trick, this term can be transformed to - -xxpp((()) Sw (x) Vo log q(0
(t)) Note
that we have eliminated the dependence on the intractable likelihood function p(x 0). Thus, this allows us rewrite J (W; 0e) as follows.
Theorem 3.1 (Local Fisher Score Matching (FSM)). Let J(W) be defined in Equation ① Under suitable boundary conditions, can be rewritten (up to an additive constant w.r.f. W) as
;(); [x--p((00,00-01000))
[>)(]]] + 2SS (x) V log (()) (.)
(2)
The complete details for Theorem 3.1 are provided in Appendix A.1 Given that we can draw proposal samples {0(2) j=ر where (63 q(O 0e) and corresponding data samples (x yn where X
jk-1
(3) p(x ((1)), the objective J (W;00, 0z) can be approximated by Monte Carlo estimation.
3(W:0) 0.) m
Sw (x)") (xx Ve log 0:)le=00))
(3)
Next, we show the optimal solution for the local FSM objective, J(W; 8,). The proof of Theorem 3.2 Appendix A.2
Theorem 3.2 (Bayes-optimal Local Fisher Score). The optimal score model for the FSM objective J(W: 0,), is
S (x: 0,) _0-p(((x,0.) Ve log p(x 0)]
As the score matching objective Equation ① is taken as expectation over the parameter proposal distribution q() 0c), the Bayes-optimal score model for this objective is generally biased and instead of being the true score at the point O., it is instead an average of the score over the posterior induced from the proposal distribution and the statistical model, that is, p(0 x, 0.). Thus, this score matching objective targets smoothed likelihood around 9,. We elaborate this in more detail in Section 5.1
3.2 Score Model Parameterization
A key aspect of the Fisher score matching technique is the choice of parameterization for the surrogate score model, Sw (x; 0z), which approximates the Fisher score the target parameter iterate 0. V log pe(x) p=) For computational tractability, we propose using lightweight linear surrogate score model based on the following derivation.
Let the surrogate score model be defined as Sw (x: 0,) = WT x, where W E Rd Xdy is the weight matrix for our model. Recall that we first draw a set of parameters { from the proposal distribution q( 0z). Then, define the j-th data matrix as X, E RRaxd constructed from training samples {x 1i-1 drawn from the model p(x () at g(3), and the j-th Gram matrix as G, = X X Using the linear score model for the local Fisher score matching objective function in Equation
and solving for the first-order conditions, we obtain the normal equation,
G =
Vo log g(0 (,) (
(4)
1=1
We can thus obtain closed-form solution for the linear Fisher score matching estimator as:
m
W
Va (0 1 | [=)
(5)

---

Once W is obtained, we can use this to construct our Fisher score estimator S(x; 8t) W x. We provide complete derivation and further discussion in the Appendix A.3 Although the local Fisher score matching objective is general framework that is agnostic to the choice of parameterization of the model, using a linear model essentially recasts the model estimation procedure as multivariate linear regression, benefiting from well-understood theory and efficient implementations. Although linear model might not be sufficient to fully capture the full data-parameter relationship, it provides a strong baseline that find works well empirically compared to more flexible neural network-based model, which incurs significant computational costs the form of an inner optimization loop and increased variance. We provide empirical comparisons with the neural network-based score model in the relevant experimental sections of the appendix, and details of the implementation in Appendix A.4 A..4
4 Likelihood-free MLE with Approximate Fisher Score
Using our local Fisher score matching (FSM) method as described in Section 3 we describe how maximum likelihood estimation (MLE) can be performed in the likelihood-fif --reee setting. Unlike many
SBI methods that attempt to estimate the likelihood globally, our method is inherently sequential by focusing only on local Fisher score estimation at the parameter point Oz.
Given set of N independent and identically distributed observations, D {x fixed parameter point 0,, we obtain an estimated FSM model S(x; 0t) using training samples (05)) 1- fx* (5)
(5) p(x @())). As the FSM model is a function of x, we can evaluate at any observation x,, providing us with an approximate gradient of the log-likelihood evaluated at 0,, (ee D) = Σ S(x;; 0r). This can then be used directly in any iterative stochastic gradient-based algorithm such stochastic gradient descent (SGD) (Robbins) and Monro 1951]. Adam [Kingma and Ba 2015], or R.SPProp [Tieleman] 2012]. where at each parameter iteration 0,, a new FSM model S(x; 0z) is estimated. The FSM-MLE algorithm with SGD presented in Algorithm 1
Ti=1. drawn from 863) q( fe). X
Algorithm 1 FSM-MLE Algorithm (SGD)
Input: N independent and identically distributed observations D = {x, ) initial parameter do. step size 7) and proposal distribution q(O 0%)
Initialize 4 0 while t < T do
1. For current iterate 6,, sample {0(i) yy=1 from proposal distribution 0 ~ g(0) (.). and then
sample corresponding data samples {x[]] (3)
=k= from xxx) ~ p(x ((:)).
2. Estimate Fisher score model Sxx; 0t) using training samples (106) 3==1 {x ( p0 JR-1)
3. Set 0r+1 Oz + nS(:: where S(D; 0.) Ei S(x;0,))
4. t 1 ++1
end while
4.1 Fisher Score Proposal Distribution
Our local FSM approach crucially uses a proposal distribution q(0 ) in the parameter space, defining a local region for the estimation of our Fisher score model. Although most distributions with differentiable and unbounded density can be used, we use an isotropic Gaussian distribution q (.) N( 0e, fe.a21), I). which has a simple, closed-form solution and direct theoretical interpretation as discussed in Section This introduces a single scalar hyperparameter or that controls the width of the proposal distribution. While further extensions such as diagonal covariance matrix or adaptive covariance could be explored, we keep to the isotropic Gaussian proposal distribution as provides simple and effective baseline. We provide further discussion on the choice of proposal distribution and calibration scheme for in Appendix A.5
5

---

Theoretical Analysis
In this section, provide theoretical analysis of our proposed local Fisher score matching technique and the stochastic gradient optimization based on this technique.
5.1 Connection Gaussian Smoothing
Gaussian smoothing is a popular zeroth-order optimization technique that estimates gradients using only function evaluations when the gradient function known [Nesterov and Spokoiny 2017 Duchi et al. 2012]. As the Gaussian smoothing gradient estimator targets a smoothed function, it
widely applicable even for non-smooth functions, which would not be amenable with standard gradient estimation, and has been shown to be robust to local optima [Starnes et al. 2023 and
applicable for many challenging machine learning problems [Salimans 2017].
Although standard Gaussian smoothing is straightforward for black-box optimization problems, note that it directly applicable in the simulation-basee inference setting the intractable likelihood L(0) - p(x 0) is not explicitly accessible. Nonetheless, we show here that our proposed local Fisher score matching technique can be directly cast likelihood-free analogue of Gaussian smoothing. Specifically, under a Gaussian proposal distribution, (0 0,) ) = No 6,, I) the Bayes-optimal Fisher score exactly the gradient of smoothed likelihood. We provide full proof of Theorem [.1 in Appendix AA66
Theorem 5.1 (Equivalence Gaussian Smoothing). Under isotropic Gaussian proposal, q(O
= NNO Oz, 0²,,, the optimal FSM estimator is equivalent to the gradient of the smoothed likelihood
0+)=
Vo, ((01;);)) = -Egmp(0(x.00) Vo logr(x
where E(0;;x) log. p(x g 0²) do and p(0 x, 0e) x p(x ())(( (.) is the induced posterior from the proposal distribution q(O 0e)
Observe that the smoothed likelihood can be further rewritten as
1.0
x(0xx) log E22.000 L(0, + zz
0,1
o= 0=5
a 1.0
where L(0;x)= = p(( 0). This is exactly the Gautsian-smoothed likelihood function, except importantly that explicit evaluations of the like- lihood L(0) were not used. Instead, our Fisher score matching technique only obtains samples from the model p(x 0) for the FSM estimation. Hence, our method directly inherits many of the robustness benefits of Gaussian smoothing while still being applicable in the SBI setting.
0=
True
x
Parameter iterates
0.6
0.4
Figure demonstrates the effects of smoothing in one-dimensional, shifted exponential likeli- hood model with single observation. The true likelihood is zero for 0 Z 0LLE, and hence any gradient- based optimization which is initialized beyond the boundary will be stuck that region. However, using our smoothed likelihood (de- picted with differing values of proposal variance a²), we are able to obtain non-zero gradient even outside the nominal support, allowing us successfully optimize the likelihood function.
0.2
0.0
Figure I: Optimizing a non-smooth, exponential likelihood with FSM estimator 0.5) for 10 parameter iterates from initial point 0g =
We can also further view the FSM procedure as a form of Empirical Bayes (EB) MMorris 1983), by interpreting the proposal distribution q( 0,) as local prior centered at 0,, which, together with the simulator model, defines an EB marginal likelihood function E(0,; x). Theorem 5.1 then shows that our Bayes-optimal FSM estimator is exactly the hyperparameter gradient of the EB marginal likelihood. Hence, this provides complementary Bayesian interpretation of our FSM method in addition to the optimization viewpoint of Gaussian smoothing.

---

5.2 Properties of the FSM estimator
We now provide theoretical guarantees for our FSM estimator under a Gaussian proposal distri- bution by characterizing its bias. In particular, by establishing the bias in terms of the smoothing hyperparameter 6, we highlight fundamental trade-off in the FSM estimation procedure.
Theorem 5.2 (Bias characterization of the FSM estimator). Let 0* be the true parameter, and denote Xo Pg. as random observations sampled from the true model. Suppose there exists unique maximum likelihood estimator for this model, and that the log-likelihood is L-smooth. Recall that g(Xx: 8,) Vg p(xo 0) 8==, is the true Fisher score, S* (xx: (xx::r) En-pp((xx06)) Vo log p(x 0)
is the optimal FSM estimator. For fixed parameter point tr.
The bias 0g is bounded by
Exol| S 0c) g(xo; S vdo Exo Exx [R(x))]
where R(x) pix 0()
is likelihood ratio term and d is the dimension of the parameter space
We provide full proof in Appendix A.7 From
Theorem 552, we can see that, increasing o, we increase the bias of the FSM estimator. Intu-
y= (Strong curvature)
T= 2.0 (Weak curvature)
Standard Devistion (Stdl Dev)
Mosn squared Error (MSE)
NGE
10
wiii
1100
itively, this is because a governs the degree of smoothing, which induces "smearing" effect of the FSM gradient estimates. On the other hand, for the linear FSM estimator, note that in the estimator W, we have the proposal gradient term V log q( 0e)lo-p() - (0C) 0z). and hence taking - 0 inflates the variance of W. Thus, there is a fundamental bias-variance trade-off the choice of a.
10a
Optinal ot 511 20
optiral of Sei Dev
10-4
10-1
10-2
10a
10:
10
10
10
10f
102
10-2
10-1
107
101
Proposal Vartance
Proposal Variance
Figure 2: Mean-squared error and standard devia-
tion of the FSM estimator S with varying proposal hyperparameter o², for two Gaussian likelihoods (O007 with differing curvature.
FigureZeemirrin empiricall illustrates the bias-variance trade-off of the linear FSM estimator with an isotropic Gaussian proposal distribution for two Gaussian likelihood models with differing curvature. In particular, the figure also shows the effect the log-likelln likelihood curvature, or the gradient- Litsschtt constant L from Theorem 5.2, on the MSE-optimal choice of the proposal scale a. When the curvature stronger (larger L),, smoothing tends to introduce more bias, and the optimal is smaller control the bias. Conversely, when curvature is weaker (smaller L), a larger is optimal to reduce the variance the score estimator.
p(x p(x0,) encodes the estimation error from using training samples around the parameter iterate points 6, to estimate an FSM estimator that evaluated at observations Xg ~ Pg-. Hence, for parameter iterates 0, that are far from the true parameter 0**, we are likely to get subpar estimation of the true gradient, while as we approach the true parameter, our estimation is likely to improve. However, increasing a, we can sample from a wider parameter space and are therefore more likely to obtain parameter samples that cover d". Thus, o also encodes an inherent exploration-exploitatinn trade-off.
Furthermore, note that the likelihood ratio term, [(x)= p(x(6")
5.3 Convergence Guarantees
As we have shown that our FSM gradient estimator closely relates to the Gaussian smoothing gradient estimator in Section 5.1 we can leverage established results showing the asymptotic convergence of stochastic gradient- based optimization methods with such biased gradient estimators. In particular, instead of using the final parameter iterate of the gradient- based optimization procedure as the approx- imated MLE OT Rs ÔMLE, we instead propose using an averaged SGD estimator 0y = Σi= 12.10 0, based on Polyak-Ruppet averaging [Polyak and Juditsky 1992 Ruppert 1988). which enjoys stronger theoretical guarantees. We provide the relevant convergence arguments in Appendix (A.8
A further benefit is that since we can obtain the quantification of the algorithmic uncertainty using the averaged SGD, 0r ÔMLLE from Appendix A.8 and the statistical uncertainty of the MLE 0MI. 0" from standard statistical theory, we can provide result showing the quantification of the joint uncertainty using the averaged SGD 0r an approximate MLE.
7

---

Theorem 5.3. Let OMLENN be the MLE,fo N i..dd. samples. Suppose that the number of iterations in the optimization algorithm T dominates the number of observations N such that v + 0 as N, - cc. Then, assuming that (O 0MLE.N) = =p((1) uniformly over both N, T and that the standard regularity conditions for the MLE are met, we have as N, - 00,
(O 0²) - ((0)II(0))))
where I(0)) is the Fisher information matrix evaluated at the true parameter
We provide the proof in Appendix A.9 Given that the Fisher information matrix I(0*) can be approximated using the Fisher score by drawing samples x; Pg and evaluating Z(0*) R
N Eim1 Vo loggo((x, p(x: OMLE) Ve log p(x; (MML)", we can also estimate this with our FSM method and take advantage of this result to obtain uncertainty quantification based Theorem 5.3
6 Related Work
The method closest to ours is the approximate MLE approach of Bertl et al. (2017), which first estimates the likelihood through kernel density estimation (KDE) before applying simultancous perturbation stochastic approximation (SPSA) [Spall] 1992] algorithm, which amounts to using a finte-ddifferences gradient estimator on the likelihood function estimated using KDE. In contrast, our FSM method directly estimates the Fisher score, merging density and gradient estimation into one step and thereby reducing both model complexity and computational overhead. After posting the first version of this manuscript on arXiv, we became aware of related independent work by [Sui et al.) [025]. J. which proposes similar Fisher score matching estimator. Their focus is on Fisher score estimation more broadly, whereas our work targets simulation
MI.E specifically.
The use of MLE in the wawuulatiionbaseeee model setting was first addressed in the seminal work on SBI Diggle and Gratton [1984), although the inference of SBI is more typically addressed within the Bayesian framework, as exemplified by the ABC algorithm. Naturally, since the maximum posteriori estimate (MAP) of the posterior distribution under a uniform prior corresponds to the MLE within the prior support, |ubiio and Johansen (2133) suggested leveraging the ABC algorithm and using KDE to obtain the MLE, and more recent neural surrogate SBI methods, such as SNLE [Papamakarios et al. 2019], while not specifically targeted for MLE, could be used in the same way. Another similar line of research is the work lonides et al. [2017] and Park 2023). which develop the MLE methodology in the SBI setting for partially observed Markov models. Research focused on developing SBI methods using score matching is growing field [Geffner et al. 2023 Sharrock et al.
2024 Jiang et al. 2025). J. however, this has been limited to amortized Bayesian inference, and, to Knowledge, we are the first work that has adapted score matching for the purpose of direct Fisher score estimation and MLE in the SBI setting.
7 Experimental Results
We evaluate our local Fisher score matching (FSM) technique on both controlled numerical studies and challenging real- woald SBI problems distribution ( | 0.) = (O, o²1) with a linear FSM estimator, and an empirical comparison with the neural network-based FSM estimator is provided in the relevant experiment sections in the Appendix.
2l For all experiments, we use an isotropic Gaussian proposal
As primary baseline, we compare against the approximate MLE method of Bertl et al. 2017 here referred to as KDE-SP. which estimates a likeliiid via kernel density estimation (KDE) and then uses a simultancous perturbation (SP) estimator to compute gradients:
(((() {((+) Ê((-)
where 5 is Rademacher random vector with i.i.d. entries, 0= + co, is perturbation constant, and ε(0 Xobs) = log p Xobs 0) is the log-likelihood estimated from the KDE by simulating data samples around the target parameter 0 and evaluating at the observations Xobs We provide further details about the implementation Appendix A.10
Code is available htttp::/giite wttp:/uw...ooo/mm////oaaa
8

---

Gradient 1lā -9ll
0= C= 0.0001
0.01
0= c= 10
0= 100
FSM Gradient
FSK Gradient
FSM Cradient
FSM Gradient
3000
30
0.10
1000
0.00
0.000
107
107
10F
10f
10
10a
107
10f
10F
10
16a
102
10*
10
10f
102
30
10
10"
10
Simulation budget
Simulation budget
Simulation budget
Simulation budget
Figure 3: Comparison of FSM and KDE-SP in 2D Gaussian model under varying hyperparameters (proposal variance or perturbation constants). Error bars show 95% CIs over 100 repeated gradient approximations.
7.1 Numerical Studies
To investigate the accuracy of gradient estimation and parameter estimation, we begin with multi- variate Gaussian model that features fixed covariance. This model has flosed-form Fisher score, allowing us to directly compare the estimated gradients from FSM and KDE-SP against the ground truth. Further details and results of this experiment are presented in Appendix A.11
One key aspect of both the FSM and KDE-SP approach is the choice of the hyperparameters, specifically the perturbation constant in KDE- SP and the proposal variance in FSM. In Fig-
Perameter Error: |è 6
FSM Gradient
KDE-SP Gradlent
ure 3, we show the sensitivity of the gradient approximation quality different choices of this hyperparameter, as the simulation budget increases. Although the gradient approximation of both methods depends strongly on the choice of hyperparameters, we see that the FSM esti- mate is always able to match the accuracy of the KDE-SP estimate given sufficient simula- tion budget, even when the hyperparameters not favorably tuned. We provide the same ab- lation study in higher-dimensional settings in Appendix 1
100
10
10*
10a
Parameter Dimension
Figure 4: Parameter estimation accuracy of both the FSM and KDE-SP methods under increasing parameter dimensions, over 100 repeated optimiza- tion runs.
In Figure we show the quality of the result- ing parameter estimate for the same multivariate Gaussian model with increasing parameter di- mension while keeping fixed. While the FSM gradient is able to main- tain the quality of the parameter estimate, the KDE-SP struggles in higher dimensional param- eter spaces, likely due to the additional kernel density estimation required.
NLE
simulation budget
1.0
Parameter Error:
Prediction Error:
0.0
0.8
0.4
0.2
FSM
KDE-SP
FSM
KDE-SP
Figure 5: Parameter estimation and prediction ac- curacy of the NLE, FSM and KDE-SP methods.
7.2 LSST Weak Lensing Cosmology Model
In this example, we use the log-normal forward model proposed by Zeghal
2024 and Lanzieri
et al. 22025) which simulates the non-Gaussian structure in gravitational weak- weakk-lensing. Using the model in the full LSST-Y10 setting, this model is representative real-world weak-lensing data. Since the generated data are digh-dimensional tomographic convergence maps (5 x 256 x 256), we use trained ResNet-15 18 compressor Alsing al. [2018), producing a 6-dimensional summary statistic. As an additional benchmark beyond the KDE-SP method, we further implement a standard neural likelihood estimator (NLE) using the SBI package [Boelts et al. 2025), trained with the same total simulation budget given to both the KDE-SP and FSM gradient-based optimization methods. Evaluated the observations, NLE can directly optimized to obtain an approximate maximum likelihood estimator. Further details and results of this experiment are presented in Appendix(A.1] A.12

---

In Figure 5J we show both the parameter estimation and the accuracy of the prediction. Given the limited simulation budget available, we observe that sequential gradient-basee optimization methods outperform the more -inteniionnnteeivvee NLE approach and that the FSM approach is generally able achieve better performance with a smaller variance.
7.3 Generator Inversion Task
In this section, we tackle the canonical prob- lem of latent inversion of generator network
(Xia et al. 2022]. For a fixed generator G and query image Xp- the goal is to recover a latent vector z such that (z) s Xp- Al- though typically z is treated in this setting. we treat as a latent variable, z N (0, oz1) and focus 0 as the parame- ter of interest. Note the marginal likelihood
4
point estimate,
(z 0) dz Pio (x 0) ((x
Gm)
Figure 6: Images from different latent mean opti- mization procedure. Top row: FSM, Middle row: KDE-SP, Bottom row: Direct optimization
is intractable because the push-forward density under G. has no closed form. However, our FSM approach allows us to estimate the Fisher score Ve log (x 0) at Xo. enabling us to maximize the likelihood (0: xo) log Pw (Xo () without directly estimating the likelihood. Conceptually, this turns the generator inversion problem into likelihood-basee based inference. Alternatively, given a differentiable generator, we can directly optimize the reconstruction loss to obtain an estimated latent mean (referred to as direct optimization).
We train GAN model on 16 16 MNIST dataset and apply the generator inversion task, comparing the direct optimization approach,
Prodiction Error: IX- X61
Prodiction Error: [X X61
FSM, and KDE-SP method. From Figure 6 we
can see that while the FSM and direct optimiza- tion is able to recover the target observation, the KDE-SP struggles achieve the same pixel quality. This is also reflected in Figure 7 which shows the reconstruction loss for the different methods. More details and results for this exper- iment are provided in Appendix A.13
FSM
Opt
KDE-SP
FSM
Opt
KDE-SP
Figure 7: Prediction error for the FSM, KDE-SP and direct optimization method.
8 Conclusion
We introduced FSM-MLE, a novel likelihood-fiee maximum likelihood estimation technique based on local Fisher score matching. By directly estimating the Fisher score in simulation- based setting, our method circumvents the need to approximate the likelihood. This significantly reduces the com- plexity of existing approaches that either rely on kernel density estimation or train expensive neural density estimators. We further showed that under an isotropic Gaussian proposal, our local Fisher score matching estimator admits natural Gaussian smoothing interpretation, thereby inheriting robustness properties from well-studied Gaussian smoothing techniques in black-box optimization. Empirical results on synthetic examples, a cosmological weak-lensing model, and a generator inver- sion task highlight the simulation efficiency and robustness of our approach. Further work includes development of more principled selection of the proposal variance o2, richer parameterization of the Fisher score model beyond the linear model, and further investigation into better leveraging the smoothing behavior to tackle challenging likelihood optimization the SBI setting, as well as utilizing the approximate Fisher information matrix for uncertainty quantification. Since our method inherently sequential, a promising extension a "semiiammrtized" variant which would leverage a pretrained neural network encoder with training dataset, which would be coupled with our proposed linear FSM model during inference, thereby enabling a more expressive model while still preserving the benefits of fast, closed-form updates.
10

---

Acknowledgments and Disclosure of Funding
The authors thank the four anonymous reviewers for their constructive feedback. SK was supported by the EPSRC Center for Doctoral Training in Computational Statistics and Data Science, grant number EP/S02356//..
References
Justin Alsing, Benjamin Wandelt, and Stephen Feeney. Massive optimal data compression and density estimation for scalable, likelihood-free inference cosmology. Monthly Notices of the Royal Astronomical Society, 477 (3):2874-2885 2018.
Mark A Beaumont, Wenyang Zhang, and David J Balding. Approximate bayesian computation in population genetics. Geneticy, 162(4):2025-2033 2002.
Johanna Bertl, Gregory Ewing. Carolin Kosiol, and Andreas Futschik. Approximate maximum likelihood estimation for population genetic inference. Statistical Applications Genetics and Molecular Biology, 16 (5-6):291-3111 312, 2017.
Ayush Bharti, François-Xan Xavier Briol, and Troels Pedersen. A general method for calibrating stochastic radio channel models with kemels.. leee transactions antennas and propagation, 7066):3986-4001, 2021.
Jan Boelts, Michael Deistler, Manuel Gloeckler, Álvaro Tejero-Cantera Jan-Matthis Lueckmann, Guy Moss, Peter Steinbach, Thomas Moreau, Fabio Muratore, Julia Linhart, Conor Durkan, Julius Vetter, Benjamin Kurt Miller, Maternus Herold, Abolfazl Ziaeemehr, Matthijs Pals, Theo Gruner, Sebastian Bischoff, Nastya Krouglova, Richard Gao, Janne K. Lappalainen, Bálint Mucsányi, Felix Pei, Auguste Schulz, Zinovia Stefanidi, Pedro Rodrigues, Cornelius Schröder, Faried Abu Zaid, Jonas Beck, Jaivardhan Kapoor, David S. Greenberg, Pedro Gonçalves, and Jakob H. Macke, reloaded: toolkit simulation based inference workflows. Journal of Open Source Software, 10(108):7754, 2025.
doi: 10.21 105/josss.0777. 007554.
URI.
httpa: //doi. org/ 10. 21105/ 077544
George Casella and Roger Berger. Statistical inference. CRC press, 2024.
Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. Proceedings the National Academy of Sciences, 117748):30055-30002 May 2020. ISSN 1091-6490. doi: 10.1073/pnas. 1912789117. URL/http://oot doi org/10, 19tp:110011/////111111111110 1912789117
Katalin Csillery, Michael G. B. Blum, Oscar E. Gaggiotti, and Olivier Francois. Approximate bayesian computation (abe) in practice. Trends in Ecology and Evolution, 25(7):410 418, July 2010. ISSN 0169-5347.
doi: 10. tr016i.tee.22010.44..0
Peter J. Diggle and Richard J. Gratton. Monte carlo methods inference for implicit statistical models. Journal of the Royal Statistical Society, Series B (Methodological), 46(2):193-221 1984, ISSN 0035-9246.
John C Duchi, Peter L Bartlett, and Martin Wainwright. Randomized smoothing for stochastic optimization. SIAM Journal on Optimization, 22(2):67-771,, 2012.
Conor Durkan, George Papamakarios, and lain Murray. Sequential neural methods likelihood- inference. (arXiv:1811.087))) November 2018. URL htttp::/aaiii.tottt wwwpeiiii.rrr//aaa wwwgg//a/////ggg 08723 arXiv:181 NDV:1811.0872, [cs. stat].
Tomas Geffnner George Papamakarios, and Andriy Mnih. Compositional score modeling for simulation-bal based inference. International Conference on Machine Learning, pages 11098-11 111161 PMLR, 2023.
Aapo Hyvärinen and Peter Dayan. Estimation non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
E. L. lonides, C. Breto, J. Park, R. A. Smith, and A. A. King, Monte carlo profile confidence intervals for dynamic systems. Journal of The Royal Society Interface. 14(132):00170126. July 2017. doi: 1..09////...017..2222
Haoyu Jiang, Yuexi Wang, and Yun Yang. Simulation-basee inference langevin dynamics with score matching, 2025. URL
03883
Yanhao Jin, Tesi Xiao, and Krishnakumar Balasubramanian. Statistical inference for polyak averaged zeroth-order stochastic gradient algorithm, 2021. URL https://axxxx org/abs/2102 05198

---

Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optmizzation In Yoshua Bengio and Yann LeCun, editors, 3rd International Conference Learning Representations, ICLR 2015. San Diego, CA. USA, May 7-9, 2015, Conference Track Proceedings. 2015. URL htttt:://aaziv.org/bbbbas
6980
Denise Lanzieri, Justine Zeghal. T. Lucas Makinen, Alexandre Boucaud, Jean-Luc Starck, and François Lanusse. Optimal neural summarisation for full-field weak lensing cosmological implicit inference, 2025. URI.
/ttpp://arxiv...o
Shakir Mohamed and Balaji Lakshminarayanan. Learning in implicit generative models. 1ew.xiv:11100034555 February 2017.
doi: 10.48350Waxxi..6100.4444 URL http://arxiv http:://axxiw..rrabbbb////// http://bs///10.... 1610. 03483 arXiv:1610.03483 [es, stat].
Carl N. Morris. Parametric empirical bayes inference: theory and applications. Journal of the American statistical Association, 78(381):47-55, 1983.
Yurii Nesterov and Vladimir Spokoiny. Random gradient- free minimization of convex functions. Foundations of Computational Mathematics, 17(2):527-50 566,, 2017.
George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation.
Advances in neural information processing systems, 30, 2017.
George Papamakarios, David Sterratt, and lain Murray. Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows. In The 22nd international conference on artificial intelligence and statistics, pages 837 848. PMLR, 2019.
Joonha Park. On simulation- inference for implicitly defined models. (arXiv:2311.0m .09446), November 2023. doi: 10.48550/arXw..3111 1.09446. [math, stat].
htttp::/aarii..rg/abb////// 09446 1xXiv:331.00944
B. T. Polyak and A. B. Juditsky. Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization, (3888)888-855 July 1992. ISSN 0363-0129, 1095-7138. doi: 10.1 13703330333
Radhakrishna Rao. Information and the accuracy attainable in the estimation of statistical parameters. In Breakthroughs in Statistics: Foundations and basic theory, pages 235-247, Springer, 1992.
Herbert Robbins and Sutton Monro. A stochastic approximation method. The annals of mathematical statistics. pages 400 007, 1951.
FJ Rubio and Adam M Johansen. A simple approach to maximum intractable likelihood estimation. Electronic Journal Statistics, 1:1632-1654, 2013.
David Ruppert. Efficient estimations from slowly convergent robbins-monro process. Technical report, Cornell University Operations Research and Industrial Engineering, 1988.
Tim Salimans, Jonathan Ho, Chen, Szymon Sidor, and Ilya Sutskever. Evolution strategies as scalable alternative to reinforcement learning. arXiv preprint Xiv: /703.03868 2017.
Chad Schafer and Peter Freeman. Likelihood free inference in cosmology: Potential for the estimation of luminosity functions. Lecture Notes in Statistics, 209:3-19, 01 2012. doi: 10. 00.1007778----41144 1---444--520--11.
Louis Sharrock, Jack Simons, Song Liu, and Mark Beaumont. Scquential neural score estimation: likelihood-free inference with conditional score based diffusion models. In Proceedings of the 41st International Conference on Machine Learning, pages 44565-44602, 2024.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019.
Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon. Sliced score matching: A scalable approach to density and score estimation. Uncertainty in artificial intelligence, pages 574-584, PMLR, 2020.
J.C. Spall. Multivariate stochastic approximation using simultaneous perturbation gradient approximation. IEEE Transactions on Automatic Control, 37(3):332-341, 1992. doi: 10. 109/9.1 19632.
Andrew Starnes, Anton Dereventsov, and Clayton Webster. Gaussian smoothing gradient descent for minimizing functions (gsmoothgd). arXiv preprint arXiv:231 1. 1.0551,, 2023.
David Sterratt, Bruce Graham, Andrew Gillies, and David Willshaw. Principles of Computational Modelling in Neuroscience. 07 2011. ISBN 008-0-511-8779544 doi: 1110177/88097781119995888
12

---

Ce Sui, Shivam Pandey, and Benjamin D. Wandelt. Fisher score matching for simulation-base forecasting and inference, 2025. URL https://arxiv. org. hwwgg/abb/250..... 07833
Tijmen Tieleman. Lecture 6.5-rmsprop: Divide the gradient by running average of its recent magnitude. COURSERA: Neural networks for machine learning, 4(2):26, 2012.
Tina Toni, David Welch, Natalja Strelkowa, Andreas Ipsen, and Michael Stumpf. Toni welch d, strelkowa ipsen stumpf mpapproximate bayesian computation scheme for parameter inference and model selection in dynamical systems. J soe interface 6: 187-202. Journal of Royal Society, Interface the Royal Society, 187 202, 03 2009. doi: 10. 10.1098/rsff.2000 200..010007777770
Aad Van Vaart. Asymptotic statistics, volume Cambridge university press, 2000.
Pascal Vincent. A connection between matching and denoising autoencoders. Neural Computation, 23(7): 1661-1674, July 2011. ISSN 0899-7667, 1530-8853 8888.. doi: 10.1162/NECO 4000142.
Weihao Xia, Yulun Zhang, Yujuu Yang. Jing-Hao Xue, Bolei Zhou, and Ming- -suan Yang. Gan inversion: survey. IEEE transactions on pattern analysis and machine intelligence, 45(3):112133388, 2022.
Ling Yang, Zhilong Zhang, Yang Song, Shenda Hong, Runsheng Xu, Yue Zhao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang. Diflusion models: comprehensive survey methods and applications. ACM computing surveys, 56(4):1-39, 2023.
Justine Zeghal, Denise Lanzieri, François Lanusse, Alexandre Boucaud, Gilles Louppe, Eric Aubourg, Adrian Bayer, and The LSST Dark Energy Science Collaboration. Simulation-basee inference henchmark for Isst weak lensing cosmology, 2024. URL httL:httt://arrt httttp://axxd...or//a////////// orggabs//405 17975
13

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
See Section 3
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
The limitations are briefly discussed in Section and left for future work.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
14

---

Justification:
See Appendix A.8
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: See Section 7 and Algorithm T
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reprodueibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
15

---

Answer: [Yes]
Justification:
Only publicly accessible datasets are used and code is provided.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/t.ps..ccppullc//uuii /
www.SSubbisiooooiiiooooooooo for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
See Appendix A.11 A.12 A.13
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: All figures in Section 7 and Appendix A.1I
A.12 and A.13 include either
етror bars representing the 95% confidence intervals over 100 repeated runs, or boxplots that visualize the distribution of the results.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
16

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative еπхо rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
See Appendix A.11 A.12 and A.13
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform. in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
Justification:
The research presented in this paper fully complies with the NeurIPS Code of
Ethics
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- cration due to laws or regulations their jurisdiction..
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This paper is purely mathematical work and does involve direct societal applications.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
17

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This paper does not work on language models.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Lieenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: See Section
Guidelines:
The answer NA means that the paper does existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
18

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: This paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets,
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
This paper does not involve crowdsourcing and human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This paper does not involve crowdsourcing and human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
19

---

A Appendix Supplemental Material
A.1 Fisher score matching objective
Here, we provide a complete theorem and a proof for Theorem 3.1
Theorem A.1 (Local Fisher Score Matching), Let (WW be defined as Equation following assumptions:
Given the
p(x 0), q(O 0x) are differentiable with respect Sw(x is differentiable with respect to x
Vx E Rdck
lim p( 0) 0.) =
|101|
(W) can be rewritten (up to an additive constant w.r.t. W) as
;() - (x)122 2ww(x) Va log g(0) 0.)]
(6)
Proof. We denote the joint distribution over (x, 0) from the distributions p(x | 0) and q 6,) as p(x.8 0 0t). First, expand the square, and remove terms which are not dependent on the score model parameters W.
J(W)- 0p Vo logp(x ()) Sw (x))11
[ q(0 0t))
/ XER
p(x 0) V logp(x 0) (x()x) de
JOERS
9(0 (.)) /
p(x 0) logp(x 10)1 + [x]]1
XERRx
2 log O0T S) (x) dx de
q(B 0:) /
p(x 0) [[))(())) Voloop(x () Tw Sw (x)) dx do
DERC
XEE'
+ (constants w.r.t. W)
Next, by exchanging integrals and using the integration by parts tricks similar to Theorem 1 in Hyvärinen and Dayan 2005).
f q(8 q(110))
/
J(W)
(.))
p(x [x]]) d0
XXER
/ p(x 0) logp(x 0) T Sw (x)dx d0
DER q00 (0))))
XERd*
[ a(0 (0,) I
p(x 0)] (x)]1 do
JOER
XERRx
2
XERRx JOERR
q(0) 0.) Vep(x 0)T Sw (x)d0 dx
=f
9(0) (.)) /
p(x 0)] (x)]]
IXERAx
da
a
2
q(0) 0))
Σ
p(x ()) (x)d0 dx
XERd> JBE30
00
20

---

J(W)- 1
/
JOERS 9(0 ((0))0))
(+)
p(x
do
XEROx
dg
2 XEROx
s(t) W (x)
q(0 (e) 00,
p(x 0)d0 dx
AER
1=1
= JOER*
q(0
(8))
p(x ()] w (x) 2x dx
XER*X
de
s(t) w (x)
a u(( 0.)p(x 0)d0 dx
2
00,
XER
BER
Finally, by further simplification
J(W) = f JBER'S
/
q(0
(b))
p(x Sw (xx
XER1X
de
+
s(t) (x)
4(0 0c) 00,
log g(0 0,)p(x Oydd dx
XERe 1=1
DERd0
= -I JOERS
/
g(0 10a)
p(x
Sw] (x) 2x d0
XXERx
+ DERdo
9(0 (b))
/
( xx) W
log g(0 0.)p((x f))x d0
XER
00,
-f JOEN
9(0 (0)) /
p(x ((x10)
de
s (x) + 25 (x)
a
log q(0 0.)) dx do
00:
XEROX
do
E p(010))p(x00)
s(1) w (x)2 + +2 (x) 00,
log q( (.))
(x)]1 (x) Vo Volog (()) 0.)
A.2 Bayes-optimal solution to Fisher score matching objective
We present the complete theorem and proof for Theorem 3.2 here.
Theorem A.2. For a general differentiable function S Rdd + I2do
S = argmin S
E Va log 0) - 5(x)]
E p(8))
v log 0)
p(x
Proof. First, observe that since the function S only a function of x, we have
E S(x) E S(x) p(x(0.)
p(x,00,,)
We decompose the objective function by expanding the square, argmnin E V log 0) 5(x) || argmin E [[S(x)]|3-2233) T Vo logo 0))
s
s p(x00.)
Then, our objective can be equivalently expressed
E
[| 25(x)
E Vo logp(x 01] x())x0)))
p(x
Which has the optimal solution S*(x)
E V log 0)
p(0)x,0.)
21

---

A.3 Linear Fisher score model parameterization
Here, we provide details of the linear Fisher score model derivation.
Recall that the parameter and data space are IRd X
Rdx,, the linear score model weights are
(
x1
W E Rdxxdd, and we defined the data matrix as X,
Ruxdx and the corresponding
Gram matrix as G2 = x] Xy-
In practice, we include an intercept term in our regression by augmenting the data matrix with (
column of ones, i.e., x
E Rd and W as a (dx + 1) x do matrix. For simplicity, we omit this
intercept term in our derivation.
We start from the empirical version of the local Fisher score matching objective, Equation (replacing averages by sums for simplicity).
((W) = Σ ΣEUR Sw x2 + (x) Vo log 01)
Substituting our linear score model, S(x: 01) Tx
((W) - m
m
ΣΣ
W
(g) x 1| +
ΣΣ>
( x
WV a log 0:)la ocn)
j=
j=1k=
To obtain the first-order conditions, take derivative with respect to W, for each of the terms separately.
For the first term,
W
( x
Σ
(x((xW)T(x,WW)
j=1
Σ
tr(w N(WTx,x,JM
==
a Applying aw gives:
Σ J=1
2x X,W = Σ GjW
à For the second term, similarly apply aw to give:
701
m n
(J)
WV 1o0 0.)l0=001] = 2 Σ
V. log q(0 0.)
Oe
j=1
j=1 1=1mm Am
Combining the two terms, obtain
we (W) = Σ G,W +
( vo logo 0.
Setting this to 0 gives us the normal equations in Equation
If the sum of the Gram matrices, Tm
G, is invertible (otherwise, we may opt to use the ridge penalty).
can directly obtain the linear Fisher score matching estimator in Equation 5
22

---

Naturally, our linear score model setup can be extended to include Frobenius norm penalty 2 in the objective, leading a ridge-type solution:
W
Σ Gj + 1=1
Volog (0) 0+)1o=ar))
This stabilises the inverse and helps prevent overfitting in finite-sample regimes. In practice, implement the ridge-type linear Fisher score model.
A.4 Neural Network Fisher score model parameterization
An alternative to the linear Fisher score model provided in Appendix A.3] is neural network parameterization of the score model. In this setting, we denote S(x; 8, Sx(x) where S& is a neural network with parameters 6 for a fixed parameter iterate 8,. The parameters é can be obtained by optimizing the FSM objective J (ó; 0,) from Equation (2) (with its Monte Carlo estimate Equation (3) using standard neural network backpropagation. Thus, following Algorithm [ using neural network parameterization requires potentially costly inner optimization loop for each
parameter iterate 0,.
A.5 Fisher score matching proposal distribution
As discussed in Section 5.2 the theoretical optimal choice of the hyperparameter depends on the curvature of the log-likelihood function. In practice, the curvature is unknown, and thus selecting the optimal is challenging in general. We propose simple pilot calibration based grid search, before rumning the main FSM-MLE procedure, we execute a short pilot FSM-MLE procedure with
different candidate value 0% yielding corresponding candidate parameter estimate dz.. For each o we simulate data at Oa.. and select the candidate σı which minimizes the discrepancy between the observed data and the simulated data at the candidate parameter estimate. Thus, this procedure selects hyperparameters which can produce simulations most consistent with the observations.
simple annealing schedule which would reduce the hyperparameter over the course of the FSM-MLE optimization procedure could also be considered. This would ensure that the smoothing bias discussed Section vanishes asymptotically, however, such a procedure would be complicated by the increase in the variance of the FSM estimator and numerical instability when is too small. While we attempted to implement such an annealing scheme our experiments, we found that introduced additional complexity without meaningful performance gains over a fixed scheme.
When the likelihood function exhibits strong anisotropic curvature, an isotropic Gaussian proposal is suboptimal. Extending the calibration scheme to diagonal covariances, however, would make grid search scale exponentially with the parameter dimension, making the method computationally prohibitive for higher dimensional problems. Hence, it remains an open question on how to efficiently design scalable procedure to select more expressive proposal distributions.
A.6 Gaussian smoothing equivalence
We provide here a more detailed derivation of Theorem 5.1
Theorem A.3 (Equivalence as Gaussian Smoothing). Under isotropic Gaussian proposal, q(O =))=((0 Or, I² with the assumptions as Theorem A.1 the optimal score matching estimator is equivalent to the gradient of the smoothed likelihood
Vo, ((0,;;x)) -E0~p(0|x,00)
V logp(x 0)
where E(0,;x) = log r(x 0) ( 02) do and p(0 x 02) x p(x 0)q(0 02) is the induced posterior from the proposal distribution q( 0c)
Proof. First, define 7(0.)= = p(x 0))(0 0c)d0 such that ε(0,;x) = log Z(0.).
Now, observe that,
Vo, Z(0.)
Z(0.) / p(x (p(100VVo) (o) 0,)d0
Vo, ē(0,;x)
2(0,)
23

---

For an isotropic Gaussian proposal, q(0 0,) = =N(0) Bz. 12), we have that
Vo,4(0 0.) = Vog(( 0.)
Using the integration-by-part trick (similarly to the proof in Appendix A.6 we have,
f p(x () 0c)d0 =
p(x V) 0,)d9
/ Ver(x 0)₂(0 0,)d0
/ Vo logg(x 0)p(x 0)g(0 0.)d0
Substituting this expression into Va, x(0,;x), have,
Vo, [(f(:x)) = Z(6,) / p(x V) 02)d0
1
1
/ /Vo Vologr(x 0)p(x 0) 0,)d0
7(0,)
Vo log 0) p(x 0)g(O 0.)). d0
Z(e,)
Eg-p(0)x00)) Vo p(x ())
A.7 Bias of FSM
Theorem A.4 (Bias characterization the FSM estimator). Let be the true parameter, and denote
Xo Pg- random observations sampled from the model. Suppose there exists unique maximum likelihood estimator for this model, and that the log-likelihood is T-smooth. Recall that g(x;; 0.) Vo log p(xo o0 is the true Fisher score, S* (Xo; 0,) Be-p(((xx00) Vo log p(x 0) is the optimal FSM estimator. For a fixed parameter point or
The bias at Og is bounded by
Exol S* (*): 01) - g( 0:)1 s L Vdd Ex, (R(x))]
where R(x)
is likelihood ratio term and d is the dimension of the parameter space 02x
pix(0')
Proof. Recall that (xx; 0s) Vologp(x 0)le=a, is the true Fisher score and the optimal FSM estimator is defined as S (x; 02) = E E0-p(0)x00, V p(x ().
We first show the bias bound (x; 0.) g(x; 0.) at fixed data point x.
||5 (x: 0.) g(x: 0.)11 = Eampt0 Vo log () V logp(x 8 p)]
$ E0-p(0 0.)
V log ((x 0) Vo logp(x 0))==a
< / ||0 0ellp(8 x, 0.)d8
p(0 x, 0²)
q(0 0c) / 110 01ll9(0 0.)d0
< sup
< I p(0 0c)
q( 0e)
24

---

Denoting p 0,) f p 0) q( 0e) d0, note that sup, p(0xx..) p(x10xLLR(x)) ((11p.] p(x(02)
We now take expectation with respect to the true model, X ~ p(x 8²), and the only non-constant term is sup, p(0)x00,)
g(9)0.)
p(0 xo. 01)
p(xo (MLL(XXO) p(xo 0 0")xxo
Ex, sup 9
q(0 0+)
p(xo 0c)
p pxxo 0*)
p(xo dxxo
p X 0,)
p(xo| 8*)
Ex x 0c)
A.8 Convergence guarantees of FSM
In this section, we provide an asymptotic convergence analysis of the stochastic gradient method based on the local Fisher score matching gradient under Gaussian proposal distribution. Recall that, based on theoretical development Section [5.T and Appendix [A.6] we have shown that the FSM estimator, under isotropic Gaussian proposal distribution, targets smoothed log-likelihood,
Eo (0; x) = log / p(x (N 10 0,021)d8
Let the N independent and identically distributed observations be D = fx and the corresponding smoothed likelihood objective be
N
Eo (0;D)
(0:x.)
We define the smoothed maximum likelihood estimator for the dataset D as
0, argmax, (0:x,)
Equivalently, assuming the concavity of the smoothed likelihood function, characterize the smoothed maximum likelihood estimator with its first-order optimality condition.
N Σ
Vola(0;
St (xx
where S5x (x; & (f:x) is the Bayes-optimal FSM estimator.
In practice, however, we utilize the linear FSM estimator $(x; 0) diseussed in Section 3.2 and Appendix A.3] In order to reduce the variance of the resulting approximate maximum likelihood estimator, well provide stronger theoretical guarantees, use the averaged parameter estimate (Polyak and Juditsky 1992)
0
T
In the following, we state an asymptotic convergence result, which can be found in Proposition 2.1 of Jin et al. 2021], which is based on Polyak and Juditsky [1922]
25

---

Assumption A.1 (Smoothness and concavity of the true log- likelihood), The log- oog-likelihodd functiion
((0:D) = > p(x 0) is L smooth and μ strongly concave.
Note that this is a sufficient condition for the strong concavity of the smoothed log-likelihodd.
Assumption A.2 (Step Size Condition). The step-sizes Np > 0 satisfies for all t, ---1++1 o(ne)
and -t=1 17 (1+A/22-1/2 < 20
Assumption A.3 (Unbiasedness and Martingale Noise Control), Define the noise term & = ((0,-1, 0); - S" (0)-1;0), which is a martingale difference sequence with respect to Fi-1 a (Ul..... -1-1), where un {(Or,,, Xx,))==1 represents all the simulations used for the score model estimation at iteration t.
1. For all iterations t 2 1, the linear FSM estimator is unbiased:
B[$(61-1,4; 4; Fi-1l = (*(0--1;0)
2. Assume that there exists a constant K > 0 such that all t 2 1, almost surely:
1 Fe-1l+ + (0.-1;;)) S K(1 + [|0,-1 (0--1-00113))
Assumption A.4 (Hessian Bound). There is function (u) with bounded fourth moments, such that the operator norm of V (0, u) is bounded, VoS u)ll S (()) for all 0
Theorem A.5. Suppose Assumptions A.1 A.3 and A.4 hold and the sequence of step sizes fulfills A.2 Using the updates of the gradient descent 8(t1 < 0, + Si(: 8G, o), have that the averaged parameter iterates Or = + E Og satisfies as - oc:
1. OT ta. 00
2. VT -0,) to No,VV
where V = (,
D))
A.8.1 Relationship between smoothed MLE and the true MLE
Furthermore, we note here that we can establish an upper bound on the distance between the smoothed MLE 9, and the true MLE 9. For simplicity, assume that we only have a single observation in our dataset, x. Then, using the strong concavity in Assumption A.I L-smoothness of the log-likelihood as in Theorem A..4] and the result from A.7 for a fixed point x, and denoting the gradient of the true log-likelihood g(0; x) (g((0;x).
110, - 11| - 1 - g(0:x)) x)| = =1lo(00:0) x)]
|0((0xx)= - Sx; 0-)1l
Lo /d
p(0 x, 0o)
sup
q(0 0z) 8
Thus, we have shown that ||O,- 0| is approximately of the order O(a).
A.9 Uncertainty quantification of FSM
Here, we show the uncertainty quantification by leveraging the result from Appendix A.8 with classical MIEE theory. As before, we denote 0, + Et' O, as the averaged parameter iterate from the FSM-SGD procedure, and for clarity, we denote MMLE, as the MLE of the true likelihood based on N i..d. observations. Our goal will be to characterize the distribution of vN 0*))
26

---

First, we note that we can decompose 0r - 0* into both algorithmic and statistical uncertainty,
0r - (0T - BMLENN) (ŌMLE.N - 0*))
algorithmic uncortainty
sampling uncertainty
Multiplying by VN, we obtain the following.
V 0*) VN OMLE,N) + VN(ÓMLE, 0²)
Focusing on the algorithmic error, observe that
VN (Or MMLE,N) = T VT(0r 0MLENN)
Since by assumption we know that Appendix A.8. this implies that their product is
- and XNT = vT(Or OMLE,N) = Op(11 from
T
XNT VN OMLE,N) 7p N T
as N,T T + 20.
From classical MLE theory, under standard regularity conditions,
VN - 0²) -a N ((, I(0")=1)
where I(0") is the Fisher information matrix at the true parameter. Finally, to combine both results, using Slutsky's theorem,
VN 8*) - N(0,I(0)')
as N, T + 00
A.10
KDE-SP implementation
We implement the KDE-SP gradient estimator proposed in Bertl ct al. [2017]. combining a kernel
density estimate (KDH)-based likelihood approximation with a simultaneous perturbation stochastic approximation (SPSA). Specifically, at each iteration i, the approximate gradient of the log-likelihood is given by:
VE(0)= (00) - -((0)
2c,
where 0+ = + c₁ and 0- c, for random perturbation 5, This gradient estimate is then used in an SPSA update of the form
0, = 01-1 a VE (0,-1)
Following the specifications in Bertl et al. (017], we adopt the standard SPSA step size schedule:
04 (t + A)a
Ct tt
with a = 1, = 1/6, and = [0.17], where is the total number of iterations.
The constants and control the initial values of a and Ct. We tune both by performing grid search over pairs a, c ). For each candidate pair, we run short trial of the SPSA optimization, simulate data from the resulting parameter estimates, and measure prediction error relative to the observed dataset. We then select the pair ( (a, e that yields the lowest validation error. We also incorporate the KDE modifications proposed in Section 3.2 of Berti et al. (017], which refine the KDE-based likelihood approximation. These modifications help stabilize the KDE estimation for high-dimensional problems.
27

---

0.0001
0.001
0.01
4000
400
40
3000
300
30
Graciert ener Jj
2000
0- C- 0.5
0=c=1 1
0- 1.5
0=
1.5
0.8
0.4
areer
2.0
0.6
0.3
0.2
0.0
0.0
0.0
0.0
d=c=
0= C- 10
d Cm 50
Cm 100
0.0150
0.25
0.125
0.2
0.100
0.15
0.075
0.015
0.0075
0.050
0.010
0.0050
0.025
0.065
0.0025
0.00
0.000
0.000
0.0000
107
107
20
10a
10"
107
207
10a
20
10
10a
10"
207
10"
10
budiget
Simulatioon budget
budget
Simulation
Figure 8: Gradient accuracy of both the Fisher score matching (FSM) technique and the KDE- SP method for bivariate Gaussian likelihood for different choices of the proposal variance and perturbation constants. The error bars represent a 95% confidence interval for 100 repeated gradient approximations.
A.11
Additional details and results on numerical studies experiment
A.11.1 Additional results hyperparameter sensitivity
In Figures and we provide additional results on the ablation study showing the sensitivity of the gradient accuracy between the FSM method and the KDE-SP method for different choices of the proposal variance and perturbation constants, respectively. Figure 3 a subset of the results shown in Figure 8 Even in nigher-dimensiona. settings, we find that the FSM method can match the gradient accuracy of the KDE-SP method across a wide range of hyperparameter choices.
A.11.2 Additional results on parameter dimension scaling
Figure 10 includes an additional result with the neural network FSM method for the parameter dimension scaling experiment seen in Figure 4 Furthermore, Figure shows the scaling with parameter dimension, with increasing simulation budgets, complementing the results seen in Figure Generally, we find that the linear FSM method performs the best across different parameter dimensions and simulation budgets. The KDE-SP method performs worse in higher dimensions, likely due to the curse of dimensionality affecting the kernel density estimate. The neural network FSM method shows competitive performance in lower dimensions, but its performance quickly degrades in higher dimensions, possibly due to optimization and/or overfitting issues.
A.11.3 Additional results wall-clock time
We provide a comparison of the wall-clock time Figures T2 and 13 for repeated gradient estimation procedures for both the KDE-SP and FSM methods. As we can in Figure [2] the FSM scales
28

---

0.0001
0.001
0.01
(= 0,1
15000
1500
12.5
10.0
Gradint eroor
5000
500
5.0
2500
c= 0.5
0=
c= 1.5
0 cm 2
2.5
2.0
1.5
Gndder ererr
0.5
0.0
0.0
0m
0=c= 10
0= 50
0= 100
0.5
0.10
0.05
0.8
0.08
0.3
0.00
0.04
0.02
0.2
0.1
0.02
0.02
10"
10a
20"
100
10"
10'
10a
100
10
10a
10"
10"
10
10*
10F
10'
10'
100
100
Simulation budget
Simulation
miuusstion
Simulation budget
Figure 9: Gradient accuracy of both the Fisher score matching (FSM) technique and the KDE-SP method for 20 dimensional Gaussian likelihood for different choices of the proposal variance and perturbation constants. The error bars represent a 95% confidence interval for 100 repeated gradient approximations.
102
FSM Linear Gradient
Parameter Error: all
KDE-SP Gradient
101
FSM NN Gradient
100
10
101
102
Parameter Dimension
Figure 10: Parameter estimation accuracy of the Linear FSM. Neural Network FSM, and KDE-SP methods under increasing parameter dimensions, over 100 repeated optimization runs.
29

---

Parameter Error: [ 8*1
Parameter Error: j 0
d = 100
FSM Laue Coudiak ondire
Parame Error: [
15
FSN Jrsur Gradkeni d-st err
Paramener Error: [
60
roMuureu Studiat Gracker
100
FSR Unar Grudlsrt Gradert 12M
10
0.25
50
0.20
10a
10a
10a
10a
10"
10a
10a
10a
Number Simulations
Number Simulations
Number Simulations
Number of Simulations
Figure 11: Parameter estimation accuracy of the Linear FSM, Neural Network FSM, and KDE-SP
methods under increasing parameter dimensions and increasing simulation budgets, over 100 repeated
optimization runs.
SM Gradient
KDE-SP Gradient
10a
Wall clock time (ns)
10
Simulation Budget
Figure 12: Wall clock time comparison between FSM and KDE-SP estimation, over 1000 runs for increasing simulation budgets
favorably with respect to the increase in the number of simulation budgets. However, in Figure B the matrix inversion step of the linear FSM method grows cubically with the parameter dimension, and hence causes an increase the wall-clock time for the FSM method. We note that in practice reduce this cost considerably by employing faster linear solvers (e.g., conjugate gradient methods), which can greatly improve scalability higher dimensions.
A.11.4 Additional results on confidence interval construction
We also note that in Figure [4] we provide a simple validation test for the use of the FSM estimate for the Fisher information matrix estimation. This shows that can recover well-calibrated confidence interval even with the use a stochastic Fisher score estimate.
Further details of these experiments provided in the Appendix A.I1.
30

---

SM Gradient KDE-SP Gradient
10"
Wall clock time (ns)
107
10
10t
Dimension
Figure 13: Wall clock time comparison between FSM and KDE-SP estimation, over 1000 runs for increasing parameter dimension
0.9
Extimated Awaraped Coverage of
0.8
0.7
0.6
FSM Linear Gradient
0.3
FSM NN Gradient
Parametric Bootstrap (B-100)
0.2
95% Coverage
10"
10a
10a
Number of Observations
Figure 14: Estimated coverage of constructed confidence interval (averaged across all parameter dimensions) from the approximated Fisher information matrix estimation with FSM estimates, comparing repeating the FSM estimation procedure with nonparametric bootstrap
A.11.5 Experimental details
The gradient comparison experiment corresponding to the plots of Figures 3 and was carried out for a bivariate Gaussian mean model with 10 observations. Observations were generated from a true mean of (1.0, 1.0) (Figure 9 used 20-dimensional Gaussian with the same setting), and Fisher score estimates were taken at the observation means, which is also the maximum likelihood estimator. The uncertainty was obtained by repeating 100 runs of the score estimation for both methods.
The multivariate Gaussian parameter estimation accuracy Figure IT was performed with 100 obser- vations, and with parameter dimensions of d = 5. 20, 50, 100 for 100 optimization steps, using 100 re- peated runs as with the previous experiment. Figures and 10 were performed ing total simulation budget of 1000 and increasing parameter dimensions of d = 2, 5, 10, 20, 50, 100. The true parameters used generate the observations were similarly taken to be a vectors of ones with the previous experiment. The (a, c) hyperparameters for the KDE-SP gradient method was
similar way, by fix-
31

---

selected from a grid of [10-2.10-1, , ,,0,, x |10-2.10 1.0%.0' 10000'' .102, 10"). For the FSM-based estimation, the (a, 1) hyperparameters. corresponding to the proposal variance and step size, were tuned in the exact same as the KDE-SP gradient gradient hyperparameters (using the prediction b
error), but over a grid of [10-3; 10-2, 10--) x [10-2, 10-1, 10 instead. The Adam [Kingma and Ba 2015) optimizer was used for the FSM-based estimation, with averaging over the last 50 iterations of the parameter iterates.
For the wall-clock time comparisons in Figures and 13 each gradient estimation procedure was timed for 1000 runs on a bivariate Gaussian mean model with 10 observations. As both the FSM and KDE-SP gradient estimation was implemented in Python and the JAX package, best attempts were made equalize the comparison between the two methods. just-in-t (JIT) compilations for both methods were disabled for the wall-clock tests to remove compilation overhead.
For the confidence interval experiment of Figure 14 5 dimensional multivariate Gaussian mean model was used. A step size of 10-3 with a = 0.05 was used with the RMSProp (Tieleman 2012] optimizer. The final Fisher information matrix was estimated by simulating 100000 simulations from the resulting MLE estimate of the optimization run, which was used to construct the confidence interval. This was repeated for 100 runs to obtain an estimated coverage probability.
For the neural network- network-based FSM method, a standard feedforward neural network with two hidden layers of size 16 with ReLU activations used. Adam optimizer with a step size of 10 was used to train the neural network for 10 iterations, for each parameter iteration of the MLE optimization procedure.
All experiments in this section were performed on a standard consumer laptop, an Intel i7-1 1370H CPU with 64GB of RAM.
A.12 Additional details on LSST weak lensing experiment
For the weak lensing experiment in Figure 15 100 iterations of the gradient optimization method were used with both the KDE-SP and FSM estimators, with 100 simulations per iteration, giving a total simulation budget of 10000 simulations for the entire optimization process. The dimension of the parameter space is 6, and the dimension of the summary statistics used is 6 as well.
The same amount of simulations was provided to a neural likelihood estimater, which is a standard masked autoregressive flow model in the package SBI in Python [Boelts et al. 2025 To mimic general, uninformative prior. we used the priors for the parameters provided in Table 1 of Zeghal ct al. 2024]. which are all Gaussian priors, and converted them to a uniform prior by taking three standard deviations from the mean, Ui - represented as N N(H, o²). The NLE was trained with 10000 (parameters, data) pairs drawn from this prior, and 5000 iterations with a standard Adam optimizer were used to train the NLE. To optimize the NLE for a specific observational dataset, we evaluated the trained NLE the specific dataset, and directly differentiated through the NLE model, giving us deterministic gradient, which is used in a standard gradient-based optimization procedure. The likelihood is optimized until convergence, where there is no longer any change in the estimated likelihood with the NLE.
3 * a, f + 3 # or. where the original Gaussian priors are
The hyperparameters (a,c) e) for the KDE-SP gradient method were selected from -g of [10 10 10-3, 1,000-10--- 10-2] x [10-3,101 10-0, 10-1, 10"). For the FSM method, we set o = 10 and a step size of 10- 2, with parameter averaging over the final 50 iterations.
For the neural network-based FSM method, a standard feedforward neural network with two hidden layers of size 16 with ReLU activations used. Adam optimizer with a step size of 10 was used to train the neural network for 10 iterations, for each parameter iteration of the MLE optimization procedure. We find that the neural network-based FSM method did not perform well in this experiment, often suffering from high variance and instability during training.
An RTX 4090 GPU with 24GB of VRAM, 41GB of RAM was used in this experiment.
A.13 Additional details on generator inversion task
For the generator inversion task in Figure 16 we trained a standard GAN on down-scaled 16 x 16 MNIST dataset, giving a data dimension of 256 as no summary statisties were used. We used 500 iterations of the gradient optimization method with both the KDE-SP and the FSM gradient
32

---

100
250
500
0.8-
5.8-
Furamates
0.0
0.5
0.6
02
ISM-NN PSM-Lin KDE-SP
PSM-NN PSM-Lin KDE-SP
PSM PSM KDE-SP
PSM PSM-Lin KDESSP
14
1.8
122
1.0
5.8
0.6
0.4
PSM-NN PSM-Lin
KDE-SP
PSM-NN FSM-Lin KDE-SP
PSM PSM
KDE-SP
PSM-NN FSM Lim KDE
Figure 15: Parameter estimation and prediction accuracy of the NLE, FSM (linear and neural-network based) and KDE-SP methods for the LSST-Y weak lensing model, for increasing number of observations
estimation procedure, with 22500 simulations per parameter iteration used in the gradient estimation. The dimension of the parameter space is 50.
The ca hyperparameters for the KDE-SP gradient method was selected from a grid of 10-3, 5 x 10-3 10-2 5 10-2) x [10-1, 10-3 5 10-3 10-2, 5 10-2). For the FSM method, we set a = 0.2 and a step size of x 10-2, with parameter averaging over the last 300 iterations. The latent mean prior, o= set at 0.1.
The direct optimization approach was performed by directly minimizing a reconstruction loss (mean squared error in pixel space) between the generated images and the observations, and directly differentiating through generator network Gw. Specifically, we minimize the following loss function.
1 Σ (2, - xoll2
min (()
where ~ N0, 12). This is done with the Adam optimizer with step size of 5 10-1, and for 1000 iterations, with n = 100 simulations per iteration.
For the neural network-based FSM method, a standard feedforward neural network with two hidden layers of size 16 with ReLU activations was used. Adam optimizer with a step size 10 was used to train the neural network for 10 iterations, for each parameter iteration of the MLE optimization procedure. Compared to Appendix A.12 we found that the neural network-base FSM method performed better in this experiment, but still generally had subpar performance compared to the linear FSM method and with increased variance.
An RTX 4090 GPU with 24GB of VRAM, 41GB of RAM was used in this experiment.
33

---

Pred. Eiror: [X Xo1
Opt
escut raNiee
Opt
essil
Frod. Error: I Xall
Kit PSV Las PROC Ogr
K3o sP pax c06 AP P9A FS0 KM Cyl
198 So
Oyt
Kre op Piti PAEEN
Figure 16: Prediction error for the FSM (linear and neural-network based), KDE-SP and direct optimization method for the latent GAN inversion, each boxplot corresponds to different observation
34


Reviews:
[{"review_id": "VkFtkJdtfZ", "reviewer": "Reviewer_eobN", "summary": "The authors present a new likelyhood function gradient estimator by approximating the fisher score with a local linear approximation formed by local sampling of a parameterized likelyhood with a known distribution. This method is a variation of the previously discovered score matching technique, can apply to settings where users have access to samples from a model, but cannot write down a closed form for the model. THey connect this work to gaussian smoothing, as the sampled score approximates a smoothed version of the original score gradient. They prove results charaterizing bias and convergence of this score approximator, and evaluate it on toy data, gravitational lensing data, and generator inversion.", "questions": "- You use a local linear model, what changes if this local model is nonlinear? \n- Does this work when the parameter spaces are complex and the gaussian neighborhood sampling assumption doesent seem reasonable. How do the conclusions of the work change in this case.\n - How does performance scale with dimensionality of the parameter space", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- Clever mathematics to compute their score function gradient \n- Good applications to scientific domains like lensing\n- A nice theoretical contribution\n- Good supporting theorems like bias and convergence analysis\n- A simple idea that is reasonably quite general and of use to the community\n\nWeaknesses:\n- Hard to undersntad the magnitufe of improvement in fig 4. It seems pretty marginal from a quick glance. \n- Demonstrations of sample efficientcy seem to be on simulated data, a real data sample efficiency analysis would be more convincing\n- A little dense, a more direct introduction might be a big benefit for the general ML reader. Likewise an overview figure to help readers quickly grok the main idea", "was_revised": true, "final_justification": "Reviewers ran considerable additional experiments, andwering many of my questions. A good paper", "rebuttal": "Dear Reviewer eobN: \n\nWe thank the reviewer for their insightful feedback of our manuscript, and are encouraged by the reviewer finding our work to be a simple yet general approach, and useful to the broader community.\n\n> (Score Model Parameterization) “You use a local linear model, what changes if this local model is nonlinear?”\n\nThe choice of score model parameterization is indeed a central design choice, and all 4 reviewers astutely raised this point. We provide a unified response to address these concerns.\n\nOur initial research started with a neural network (NN) parameterization, but we ultimately opted for a linear model to align with our goal of developing a fast, robust, and lightweight MLE method. While we acknowledge that a NN model could potentially better capture the relationship of the Fisher score, this flexibility comes with significant computational costs in the form of an inner optimization loop and increased variance. In order to empirically validate this tradeoff, we implement the NN score model (2 layer MLP, using 10 steps of Adam per iteration) and extend the experiments in Section 7.1 and 7.3 of our manuscript. \n\nIn Table 1, we can see that as the dimension increases, the NN score model starts to become unstable, however the linear FSM is more robust and scales relatively well. In Table 2, on a more complex problem despite the additional complexity of the NN model, we are not able to outperform the linear model, and this comes with the cost of increased variance. \n\n\n| Param Dim   |               NN FSM|        Linear      FSM |             KDE-SP |\n|-------------|------------------|-----------------|-----------------|\n| 2           |  0.1235 (0.0065) | 0.1232 (0.0064) | 0.3013 (0.0138) |\n| 5           |  0.2107 (0.0071) | 0.2068 (0.0070) | 0.2324 (0.0086) |\n| 10          |  0.3051 (0.0079) | 0.2985 (0.0078) | 0.3455 (0.0085) |\n| 20          |  0.8285 (0.2324) | 0.4409 (0.0068) | 1.7109 (0.0232) |\n| 50          | 59.6913 (0.7021) | 0.7006 (0.0068) | 3.6558 (0.0433) |\n| 100         | 82.2349 (1.0051) | 0.9790 (0.0071) | 7.3584 (0.0457) |\n\nTable 1: Estimation Error vs. Parameter Dimension on multivariate Gaussian Model, std. err. in brackets\n\n\n|              NN FSM|       Linear       FSM |             KDE-SP |             Opt |\n|-----------------|-----------------|-----------------|-----------------|\n| 2.4963 (0.1635) | 2.0635 (0.0397) | 6.1476 (0.1594) | 1.6575 (0.0726) |\n| 2.2283 (0.1998) | 1.9612 (0.0533) | 6.2249 (0.0196) | 1.3495 (0.2286) |\n| 0.9394 (0.0634) | 0.7961 (0.0367) | 3.0341 (0.4033) | 0.6230 (0.0356) |\n| 0.8271 (0.0588) | 0.7863 (0.0427) | 1.5439 (0.1210) | 0.4170 (0.0142) |\n\n\nTable 2: Prediction Error on Generator Inversion Task, std. err. in brackets\n\nOur empirical findings show that a straightforward NN model is not able to match a linear model, and it remains an interesting research direction about how best to incorporate a NN model in our framework. One promising extension, inspired by the reviewer's comments, is a hybrid approach that incorporates both a NN and linear model. During an offline pretraining phase, we use our FSM method to train the NN model over a wide range of simulated data. Then, during inference time, we can leverage the NN model as an encoder and use the linear FSM model on the encoded observations, allowing us to maintain the computational efficiency of the linear model during inference. We will include a further discussion of this extension, as well as our discussion on the trade-offs between a linear and NN score model in a future revision of our manuscript.\n\n\n> \"Does this work when the parameter spaces are complex and the gaussian neighborhood sampling assumption doesent seem reasonable. How do the conclusions of the work change in this case.”\n\nWe would like to clarify that the “gaussian neighborhood sampling” is not an assumption imposed on the model, but a part of the Fisher score matching procedure that is used as a proposal distribution, dictating the local neighborhood of drawn parameter samples for which we perform the score matching. While our framework is general enough to accommodate proposal distributions beyond the Gaussian distribution (we require that the proposal distribution be fully supported, differentiable and with a known density function), the Gaussian proposal distribution is often the best choice in terms of its simplicity and convenience, as well as its deep connection to the Gaussian smoothing gradient estimator in black-box optimisation.\n\n\n> (Scalability of FSM) “How does performance scale with dimensionality of the parameter space”\n\nWe extend the numerical study in Section 7.1 of our manuscript to assess parameter estimation accuracy up to $d=100$, keeping the simulation budget fixed.\n\n| Param Dim   |              FSM |             KDE-SP |\n|-------------|-----------------|-----------------|\n| 2           | 0.1239 (0.0064) | 0.1257 (0.0066) |\n| 5           | 0.2108 (0.0072) | 0.2240 (0.0071) |\n| 10          | 0.3038 (0.0080) | 0.6999 (0.0325) |\n| 20          | 0.4438 (0.0069) | 1.1422 (0.0250) |\n| 50          | 0.7047 (0.0069) | 4.5179 (0.0281) |\n| 100         | 0.9898 (0.0074) | 7.6113 (0.0276) |\n\nTable 3: Estimation Error (MSE) vs. Parameter Dimension on multivariate Gaussian Model. Std. err. in brackets.\n\nAs Table 3 shows, the FSM estimator's error grows relatively modestly with the dimension compared to KDE-SP. The FSM estimator is more robust as it directly estimates the Fisher score, instead of incorporating a two-stage estimation process (KDE, and then finite-differences) for the KDE-SP baseline. This also is line with the theoretical analysis provided in Theorem 5.2 of our paper, where the bias of the FSM estimator scales as $O(\\sqrt{d})$. \n\n\n> (Clarity) \"A little dense, a more direct introduction might be a big benefit for the general ML reader. Likewise an overview figure to help readers quickly grok the main idea\"\n\nWe thank the reviewer for this valuable feedback, and will rewrite the introduction in a future revision of this manuscript to ensure that it is more accessible and clearer to a general ML reader."}, {"review_id": "qmLvPJLjPH", "reviewer": "Reviewer_e3CT", "summary": "The paper introduces a novel method for likelihood-free MLE by estimating the Fisher score via local score matching around each parameter iterate. A linear surrogate model yields a closed-form estimator. The method is connected to Gaussian smoothing and outperforms alternative approaches in simulation efficiency.", "questions": "- Have you tried more flexible (e.g., neural) score models?\n- How sensitive is the method to the choice of $\\sigma^2$ in high dimensions?\n- Your estimator is based on a local posterior $p(\\theta \\mid x, \\theta_t) \\propto p(x \\mid \\theta) q(\\theta \\mid \\theta_t)$ induced by the proposal — this resembles an empirical Bayes procedure where the prior adapts during optimization. Could you comment on this connection, and whether an EB perspective might offer additional insights or extensions?", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n\n- Novel, theoretically grounded estimator for likelihood-free Fisher score.\n- Closed-form linear model enables fast and stable optimization.\n- Strong empirical results across diverse SBI tasks.\n- Clear connection to Gaussian smoothing.\n\n**Weaknesses**\n- $\\theta_t$ and $\\theta_j$ notation may confuse readers due to overlapping roles.\n- Limited to linear score models; unclear if more expressive models would help.\n- Not amortized; optimization must be redone per dataset.", "was_revised": true, "final_justification": "**Paper:** *Direct Fisher Score Estimation for Likelihood Maximization*\n\n---\n\n1. **Flexible score models** – New experiments with a 2-layer MLP showed it matches linear FSM in low dimensions but becomes unstable and higher variance in high dimensions. Authors propose a hybrid NN + linear FSM “semi-amortized” extension.\n\n2. **Sensitivity to $\\sigma$** – Theory: bias–variance tradeoff worsens with higher $d$ (Theorem 5.2). Practice: no major tuning issues found, but calibration discussion will be expanded.\n\n3. **Empirical Bayes link** – Proposal $q(\\theta'|\\theta_t)$ acts as local prior; smoothed log-likelihood equals local model’s log marginal likelihood. EB view to be added in revision.\n\n4. **Notation** – Overlaps in $\\theta$ will be clarified.\n\n---", "rebuttal": "Dear Reviewer e3CT:\n\nWe thank the reviewer for their positive feedback. We are particularly encouraged that the reviewer found our work to be novel and theoretically grounded.\n\n> (Score Model Parameterization) “Have you tried more flexible (e.g., neural) score models?” and \"Limited to linear score models; unclear if more expressive models would help.\"\n\nThe choice of score model parameterization is indeed a central design choice, and all 4 reviewers astutely raised this point. We provide a unified response to address these concerns.\n\nOur initial research started with a neural network (NN) parameterization, but we ultimately opted for a linear model to align with our goal of developing a fast, robust, and lightweight MLE method. While we acknowledge that a NN model could potentially better capture the relationship of the Fisher score, this flexibility comes with significant computational costs in the form of an inner optimization loop and increased variance. In order to empirically validate this tradeoff, we implement the NN score model (2 layer MLP, using 10 steps of Adam per iteration) and extend the experiments in Section 7.1 and 7.3 of our manuscript. \n\nIn Table 3, we can see that as the dimension increases, the NN score model starts to become unstable, however the linear FSM is more robust and scales relatively well. In Table 4, on a more complex problem despite the additional complexity of the NN model, we are not able to outperform the linear model, and this comes with the cost of increased variance. \n\n\n| Param Dim   |               NN FSM|        Linear      FSM |             KDE-SP |\n|-------------|------------------|-----------------|-----------------|\n| 2           |  0.1235 (0.0065) | 0.1232 (0.0064) | 0.3013 (0.0138) |\n| 5           |  0.2107 (0.0071) | 0.2068 (0.0070) | 0.2324 (0.0086) |\n| 10          |  0.3051 (0.0079) | 0.2985 (0.0078) | 0.3455 (0.0085) |\n| 20          |  0.8285 (0.2324) | 0.4409 (0.0068) | 1.7109 (0.0232) |\n| 50          | 59.6913 (0.7021) | 0.7006 (0.0068) | 3.6558 (0.0433) |\n| 100         | 82.2349 (1.0051) | 0.9790 (0.0071) | 7.3584 (0.0457) |\n\nTable 3: Estimation Error vs. Parameter Dimension on multivariate Gaussian Model, std. err. in brackets\n\n\n|              NN FSM|       Linear       FSM |             KDE-SP |             Opt |\n|-----------------|-----------------|-----------------|-----------------|\n| 2.4963 (0.1635) | 2.0635 (0.0397) | 6.1476 (0.1594) | 1.6575 (0.0726) |\n| 2.2283 (0.1998) | 1.9612 (0.0533) | 6.2249 (0.0196) | 1.3495 (0.2286) |\n| 0.9394 (0.0634) | 0.7961 (0.0367) | 3.0341 (0.4033) | 0.6230 (0.0356) |\n| 0.8271 (0.0588) | 0.7863 (0.0427) | 1.5439 (0.1210) | 0.4170 (0.0142) |\n\n\nTable 4: Prediction Error on Generator Inversion Task, std. err. in brackets\n\nOur empirical findings show that a straightforward NN model is not able to match a linear model, and it remains an interesting research direction about how best to incorporate a NN model in our framework. One promising extension, inspired by the reviewer's comments, is a hybrid approach that incorporates both a NN and linear model. During an offline pretraining phase, we use our FSM method to train the NN model over a wide range of simulated data. Then, during inference time, we can leverage the NN model as an encoder and use the linear FSM model on the encoded observations, allowing us to maintain the computational efficiency of the linear model during inference. We will include a further discussion of this extension, as well as our discussion on the trade-offs between a linear and NN score model in a future revision of our manuscript.\n\n\n> “How sensitive is the method to the choice of $\\sigma^2$ in high dimensions?”\n\nFrom a theoretical perspective, we expect to find that the sensitivity of the smoothing hyperparameter $\\sigma$ increases as the dimension increases. That is because, as discussed in Section 5.2 of our manuscript, the $\\sigma$ governs a bias-variance tradeoff in the FSM estimator. Further, note that the dimension of the parameter space $d$ appears in the bias bound in Theorem 5.2. Hence, when the dimension increases, we expect the bias to increase, and hence $\\sigma$ will necessarily need to be decreased in order to control the bias. However, this inflates the variance of the FSM estimator, making the calibration of $\\sigma$ challenging in high-dimensional settings.\n\nDespite the theoretical challenges, we did not find the tuning of $\\sigma$ to be a practical issue, even for the relatively challenging, higher-dimensional Generator inversion task, with a parameter and data dimensions of $50$ and $256$ respectively. \n\nWe thank the reviewer for raising this important point. We will incorporate this into our discussion of the calibration scheme for $\\sigma$, and further investigate the sensitivity issues of $\\sigma$ in high dimensions in a future revision of the manuscript.\n\n> “Your estimator is based on a local posterior $p\\left(\\theta \\mid x, \\theta_t\\right) \\propto p(x \\mid \\theta) q\\left(\\theta \\mid \\theta_t\\right)$ induced by the proposal — this resembles an empirical Bayes procedure where the prior adapts during optimization. Could you comment on this connection, and whether an EB perspective might offer additional insights or extensions?”\n\nWe thank the reviewer for the insightful remark! Indeed, our FSM procedure can be interpreted as a sequential form of EB. At each iteration, the proposal distribution $q(\\theta \\mid \\theta_t)$ can be seen as a local prior, with the parameter iterate $\\theta_t$, the mean of the proposal distribution, acting as a hyperparameter. Furthermore, the smoothed log-likelihood objective $\\tilde{\\ell}(\\theta_t; x)$ as shown in Theorem 5.1 of our manuscript, is exactly the log marginal likelihood of the local model, that is used in EB\n\nWhile in our manuscript, we focused on the optimization viewpoint, focusing on interpreting FSM as Gaussian smoothing, the EB perspective provides the FSM with a complementary statistical interpretation. We will include a detailed discussion of this EB viewpoint in a future revised manuscript.\n\n> \"$\\theta_t$ and $\\theta_j$ notation may confuse readers due to overlapping roles.\"\n\nWe thank the reviewer for this feedback. We will revise the manuscript to ensure that there is no ambiguity in the notation for the thetas."}, {"review_id": "tbBN7zlrUr", "reviewer": "Reviewer_6bPt", "summary": "This paper tackles the problem of maximizing a log-likelihood $\\log p(x|\\theta)$, with data $x$ and parameters $\\theta$, in the simulation-based inference setting, i.e., when $p(x|\\theta)$ is intractable but model simulations can be obtained. The authors present an adaptation of Stochastic Gradient Descent (SGD) algorithm, where the gradient $\\nabla_\\theta \\log p(x|\\theta)$ is approximated via score matching at each iteration. To be more precise, given the current parameter value $\\theta_t$, they consider a Gaussian prior $p(\\theta|\\theta_t)$ centered in $\\theta_t$ (with arbitrary scalar variance $\\sigma$ defined as a hyperparameter) and propose to approximate $\\nabla_{\\theta=\\theta_t} \\log p(x|\\theta)$ by the conditional expectation of $\\nabla_\\theta \\log p(x|\\theta)$ over the posterior distribution $p(\\theta|x, \\theta_t)\\propto p(x|\\theta)p(\\theta|\\theta_t)$ (which I denote here as the \"Bayes estimator\"). Interestingly, the L2 loss induced by this conditional expectation can be written into a simulation-free loss similar to Implicit Score Matching (ISM), coined Local Fisher Score Matching, but with the advantage of bypassing the need of computing divergence of the approximation function. In practice, they minimize this loss over the class of linear functions wrt x, which results in an explicit minimizer and avoids intensive use of computational budget at each gradient step.\n\nOn the theoretical side, they show that the Bayes estimator of $\\nabla_\\theta \\log p(x|\\theta)$ can be easily interpreted as the score of a smoothed version of the original likelihood (namely, Gaussian smoothing which is widely used in zeroth-order optimization), and quantify the bias of this estimator with interpretable quantities. Interestingly, they demonstrate that the quality empirical estimator (i.e., the linear approximation of the Bayes estimator) is intimately induced by a trade-off on the scalar variance $\\sigma$. Finally, they present numerical experiments where they compare their method with a KDE-based SGD approach, and evaluate the robustness of their algorithm with respect to the hyperparameter $\\sigma$ and the dimension of the problem for various tasks where ground truth can be made available (cosmology, generator inversion..).", "questions": "- In practice (i.e., without access to ground truth), what is the strategy to tune $\\sigma$ ?\n- Have you considered (i) taking a diagonal covariance matrix to model the Gaussian prior $p(\\theta|\\theta_t)$ and (ii) optimizing the covariance on-the-fly (i.e., making it depend on $\\theta_t$) ?\n- Have you considered a neural network parameterization of the Bayes estimator combined with a joint SGD procedure (i.e, taking some local SGD steps wrt the FSM loss (2) at each global gradient step) ? I have the feeling that it won't be too much costly, while dramatically increasing the expressivity of the model.", "limitations": "The authors acknowledge two limitations of their work : (i) the potential lack of expressivity from the linear parameterization and (ii) the careful tuning of the hyperparameter $\\sigma$.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n- I find that the paper is very well written, I enjoyed very much the reading : the organization of the sections and the results is well chosen and easy to follow\n- although the presented approach may be simple (combining Gaussian prior with linear parameterization), the resulting algorithm is well justified in theory and very efficient in practice (notably due to the lack of divergence computation in the ISM-like objective) : the overall idea is very nice and has been thought to handle reasonable computational budget, which is key for practical applications\n- this paper present solid theoretical results (Theorems 5.1, 5.2, 5.3), which are interpretable and may help for applications\n- the proposed algorithm demonstrates better robustness to increasing dimension than previous work\n\n**Weaknesses** \n- in my opinion, the major weakness of this algorithm is the choice of the scalar parameter $\\sigma$ for practical applications. As explained above, it plays the central role in the quality of the Bayes estimator (a low $\\sigma$ would increase the variance, while a large $\\sigma$ would increase the bias) and has to be carefully tuned. As far as I understand, this hyperparameter has been tuned in the presented numerical experiments based on access to ground truth (hence, by minimizing the prediction error wrt the optimal parameter). However, in practice, such tuning would seem impossible since both ground truth and likelihood function are not available... Furthermore, this hyperparameter is fixed from the beginning of the optimization procedure, which may not be optimal in practice.\n- a second weakness (which is acknowledged by the authors) lies in the linear parameterization of the approximated gradient : the underlying assumption is that $x$ is distributed according to a Gaussian distribution given $\\theta$, i.e., $p(x|\\theta)\\approx N(x; f(\\theta), g(\\theta))$, which may be too simple for real-world applications.", "was_revised": true, "final_justification": "The authors addressed precisely my questions in the rebuttal, and plan to enrich the manuscript with the given explanations.", "rebuttal": "Dear Reviewer 6bPt:\n\nWe thank the reviewer for their thoughtful feedback.\n\n> (Tuning of $\\sigma$ hyperparameter) \"In practice (i.e., without access to ground truth), what is the strategy to tune $\\sigma$?\"\n\nWe would like to assure the reviewer that we *do not use* the ground truth parameters in the calibration of the smoothing hyperparameter $\\sigma$. Instead, we compare directly between simulated samples and the observation. Reviewer kgCt raised a similar point, and we provide a similar response here. \n\nWe use the following *$\\sigma$ calibration scheme*: In a pilot phase before the main optimization scheme, we perform grid search over a range of hyperparameter values. For each candidate $\\sigma$:\n\n1. We run the SGD algorithm to obtain a parameter estimate $\\hat{\\theta}_\\sigma$ corresponding to the candidate hyperparameter $\\sigma$\n2. Using $\\hat{\\theta}_\\sigma$, we draw simulated samples, and calculate the average Euclidean distance between the simulated samples and observations\n\nFinally, we select the hyperparameter $\\sigma^*$ that minimizes the average distance. Thus, our calibration procedure chooses the hyperparameter that can produce simulations consistent with the target observation. This principle of matching simulations to observations is a cornerstone of intractable likelihood methods.[1]\n\nWe would like to provide some insight into the *optimal selection of $\\sigma$*. The optimal choice of $\\sigma$ is intimately linked with the curvature of the log-likelihood function, and as discussed in Section 5.2 of our manuscript, dictates a bias-variance tradeoff for the FSM estimator. Using the definition of the Lipschitz constant $L$ in Theorem 5.2,\n* A *strong curvature* log-likelihood corresponds to a larger $L$ - in this regime, the bias term dominates, and a smaller $\\sigma$ is preferable to control the bias.\n* A *weak curvature* log-likelihood corresponds to a smaller $L$ - here, the variance term dominates, and a larger $\\sigma$ is preferred to reduce the variance.\n\nWe provide empirical evidence for this using a numerical study of a Gaussian mean likelihood model, $x \\sim N(\\theta, \\tau^2)$, where $\\tau^2$ governs the likelihood curvature. We provide the MSE between the FSM and the true Fisher score and the variance of the FSM estimator.\n\n|     σ |     MSE |   Var. |\n|-------|---------|-------------|\n| 1e-02 |  1.3981 |      1.3744 |\n| 5e-02 (*)  |  0.8202 |      0.2878 |\n| 1e-01 |  1.3539 |      0.1072 |\n| 5e-01 |  7.1327 |      0.0057 |\n| 1e+00 | 10.2868 |      0.0009 |\n| 5e+00 | 14.5088 |      0.0000 |\n| 1e+01 | 15.2296 |      0.0000 |\n\nTable 1: Strong curvature ($\\tau = 0.5$), optimal $\\sigma$ is small\n\n\n|     σ |    MSE |   Var. |\n|-------|--------|-------------|\n| 1e-02 | 0.9871 |      0.9839 |\n| 5e-02 | 0.1755 |      0.1739 |\n| 1e-01 | 0.1015 |      0.1003 |\n| 5e-01 | 0.0150 |      0.0146 |\n| 1e+00 (*)  | 0.0085 |      0.0055 |\n| 5e+00 | 0.0186 |      0.0004 |\n| 1e+01 | 0.0324 |      0.0001 |\n\nTable 2: Weak curvature ($\\tau = 2.0$), optimal $\\sigma$ is large\n\nFurthermore, we note that the variance of the FSM estimator decreases monotonically as $\\sigma$ increases, in line with our theoretical analysis. \n\nHowever, in practice, in the intractable likelihood setting, information about the curvature is unlikely to be available. Obtaining a reasonable estimate for the curvature, in the context of our goal of gradient estimation, would likely require further methodological development. While our proposed calibration scheme is more of a heuristic, even for black-box randomized optimization algorithms where the objective function is available, developing optimal, automatic tuning of the hyperparameters is a challenging and ongoing research area [2]. Our proposed calibration scheme is a pragmatic and data-driven heuristic, similar in spirit to hyperparameter tuning for many ML models that controls a bias-variance tradeoff.\n\nWe will ensure that the details of this calibration scheme and its rationale is fully clarified in our revised manuscript.\n\n> (Tuning of $\\sigma$ hyperparameter) \"Have you considered (i) taking a diagonal covariance matrix to model the Gaussian prior $p\\left(\\theta \\mid \\theta_t\\right)$\"\n\nUsing a diagonal covariance matrix for the proposal distribution would be advantageous, especially when the likelihood surface exhibits significantly different curvature across parameter dimensions. However, in practice, a calibration scheme using grid search scales exponentially with the number of hyperparameters, and using a diagonal covariance matrix introduces $d$ seperate hyperparameters, making the method computationally prohibitive for higher dimensional problems. We thus focus on an isotropic Gaussian proposal distribution in order to maintain the tractability and simplicity of our FSM procedure.\n\n> (Tuning of $\\sigma$ hyperparameter) \"(ii) optimizing the covariance on-the-fly (i.e., making it depend on $\\theta_t$ ) ?\"\n\nA typical choice of adaptation schedule would be to anneal the hyperparameter $\\sigma$ towards 0, which would ensure that the bias induced by the smoothing vanishes asymptotically, allowing recovery of the true MLE, and this is the standard approach in black-box randomized optimization [3]. However, in practical application, this approach is complicated by an increase in the variance of our FSM estimator and numerical instability when $\\sigma$ is decreased too aggresively. Furthermore, an adaptive schedule introduces new hyperparameters such as the decay rate of the schedule, increasing the complexity of the hyperparameter calibration. Our preliminary experiments with standard schedules did not produce a clear improvement over a fixed choice of $\\sigma$, and hence we opted towards a fixed $\\sigma$. This problem of finding a good schedule for $\\sigma$ is analogous to the tuning of the learning rate for SGD. A fixed learning rate is often a strong baseline, but many sophisticated learning rate schedules are available, and we are interested in exploring the incorporation of such schemes for an adaptive $\\sigma$. \n\nWe will include this discussion on adaptive and diagonal covariance matrices in our revised manuscript.\n\n> (Score Model Parameterization) \"Have you considered a neural network parameterization of the Bayes estimator...\"\n\nThe choice of score model parameterization is indeed a central design choice, and all 4 reviewers astutely raised this point. We provide a unified response to address these concerns.\n\nOur initial research started with a neural network (NN) parameterization, but we ultimately opted for a linear model to align with our goal of developing a fast, robust, and lightweight MLE method. While we acknowledge that a NN model could potentially better capture the relationship of the Fisher score, this flexibility comes with significant computational costs in the form of an inner optimization loop and increased variance. In order to empirically validate this tradeoff, we implement the NN score model (2 layer MLP, using 10 steps of Adam per iteration) and extend the experiments in Section 7.1 and 7.3 of our manuscript. \n\nIn Table 3, we can see that as the dimension increases, the NN score model starts to become unstable, however the linear FSM is more robust and scales relatively well. In Table 4, on a more complex problem despite the additional complexity of the NN model, we are not able to outperform the linear model, and this comes with the cost of increased variance. \n\n\n| Param Dim   |               NN FSM|        Linear      FSM |             KDE-SP |\n|-------------|------------------|-----------------|-----------------|\n| 2           |  0.1235 (0.0065) | 0.1232 (0.0064) | 0.3013 (0.0138) |\n| 5           |  0.2107 (0.0071) | 0.2068 (0.0070) | 0.2324 (0.0086) |\n| 10          |  0.3051 (0.0079) | 0.2985 (0.0078) | 0.3455 (0.0085) |\n| 20          |  0.8285 (0.2324) | 0.4409 (0.0068) | 1.7109 (0.0232) |\n| 50          | 59.6913 (0.7021) | 0.7006 (0.0068) | 3.6558 (0.0433) |\n| 100         | 82.2349 (1.0051) | 0.9790 (0.0071) | 7.3584 (0.0457) |\n\nTable 3: Estimation Error vs. Parameter Dimension on multivariate Gaussian Model, std. err. in brackets\n\n\n|              NN FSM|       Linear       FSM |             KDE-SP |             Opt |\n|-----------------|-----------------|-----------------|-----------------|\n| 2.4963 (0.1635) | 2.0635 (0.0397) | 6.1476 (0.1594) | 1.6575 (0.0726) |\n| 2.2283 (0.1998) | 1.9612 (0.0533) | 6.2249 (0.0196) | 1.3495 (0.2286) |\n| 0.9394 (0.0634) | 0.7961 (0.0367) | 3.0341 (0.4033) | 0.6230 (0.0356) |\n| 0.8271 (0.0588) | 0.7863 (0.0427) | 1.5439 (0.1210) | 0.4170 (0.0142) |\n\n\nTable 4: Prediction Error on Generator Inversion Task, std. err. in brackets\n\nOur empirical findings show that a straightforward NN model is not able to match a linear model, and it remains an interesting research direction about how best to incorporate a NN model in our framework. One promising extension, inspired by the reviewer's comments, is a hybrid approach that incorporates both a NN and linear model. During an offline pretraining phase, we use our FSM method to train the NN model over a wide range of simulated data. Then, during inference time, we can leverage the NN model as an encoder and use the linear FSM model on the encoded observations, allowing us to maintain the computational efficiency of the linear model during inference. We will include a further discussion of this extension, as well as our discussion on the trade-offs between a linear and NN score model in a future revision of our manuscript.\n\nReferences:\n\n[1] : Cranmer K, Brehmer J, Louppe G. Frontier of simulation-based inference. Proc Natl Acad Sci USA. 2020;117(48):30055-62.\n\n[2] : Ren K, Luo L. Parameter-free near-optimal zeroth-order algorithm for stochastic convex optimization. Proc 42nd ICML.\n\n[3] : Duchi JC, Jordan MI, Wainwright MJ, Wibisono A. Optimal rates for zero-order convex optimization: two function evaluations. IEEE Trans Inf Theory. 2015;61(5):2788-806."}, {"review_id": "2DC8z61nYT", "reviewer": "Reviewer_kgCt", "summary": "The paper presents an approach to estimate maximum likelihood (MLE) estimator is a simulation-based inference setting i.e. where we do not have access to the underlying likelihood distribution but can generate samples from it. The authors approach this by learning Fisher score with a local score matching technique. When the score is parameterized as a linear model, the authors derive a close-form solution for it. The authors present theoretical discussion on convergence and other properties of the estimator, and present various empirical results.", "questions": "Major questions--\n\n- While linear parameterization for score is convenient in that it gives a closed-form solution, I expect it can be severely constraining in some settings. Have the authors considered using a more flexible model for score and solve Eq. (3) numerically? Some discussion on this should be done in the paper. \n- Along above lines, it also seems that currently the authors fit a new score model (W) at every point, but can this be amortized? Similarly, will have a training dataset (as is generally done in SBI) help with more flexible score model? Again, some discussion on this will improve the paper.\n- The authors discuss the impact of $\\sigma_t$, but is there any guidance on how to choose it for a new problem? Along the same lines, do the authors consider varying/adapting $\\sigma$ over iterations (e.g. depending on where they are in parameter space or as they get closer to MLE).\n- How does the approach scale with dimensions, both parameter and data dimensions? If I understand correctly, the highest dimensions in current examples is 20. Again, some guidance along these lines is necessary for new applications. \n\nMinor points -- \n- I would recommend the authors include the algorithm-box in the main-text in the updated version (given the extra page). \n- Dimensionality of the parameters is not given in MNIST example\n- Convergence of FSM is difficult to see in Fig. 3, right panel.", "limitations": "While the paper does discuss some limitations, I believe it can be expanded.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The main text is well written and easy to follow. The approach is clearly motivated and presented, with discussion on convergence and other properties of the proposed estimator. \n\nRegarding weakness, see the questions section.", "was_revised": true, "final_justification": "The authors have satisfactorily responded to all the questions I raised, and the addition of new discussion and experiments further strengthens the paper.", "rebuttal": "Dear Reviewer kgCt:\n\nWe thank the reviewer for their thoughtful feedback.\n\n> Score Model Parameterization\n\nThe choice of score model parameterization is indeed a central design choice, and all 4 reviewers astutely raised this point. We provide a unified response to address these concerns.\n\nOur initial research started with a neural network (NN) parameterization, but we ultimately opted for a linear model to align with our goal of developing a fast, robust, and lightweight MLE method. While we acknowledge that a NN model could potentially better capture the relationship of the Fisher score, this flexibility comes with significant computational costs in the form of an inner optimization loop and increased variance of the gradient estimator. In order to empirically validate this tradeoff, we implement the NN score model (2 layer MLP, using 10 steps of Adam per iteration) and extend the experiments in Section 7.1 and 7.3 of our manuscript. \n\nIn Table 1, we can see that as the dimension increases, the NN score model starts to become unstable, however the linear FSM is more robust and scales relatively well. In Table 2, on a more complex problem despite the additional complexity of the NN model, we are not able to outperform the linear model, and this comes with the cost of increased variance. \n\n\n| Param Dim   |               NN FSM|        Linear      FSM |             KDE-SP |\n|-------------|------------------|-----------------|-----------------|\n| 2           |  0.1235 (0.0065) | 0.1232 (0.0064) | 0.3013 (0.0138) |\n| 5           |  0.2107 (0.0071) | 0.2068 (0.0070) | 0.2324 (0.0086) |\n| 10          |  0.3051 (0.0079) | 0.2985 (0.0078) | 0.3455 (0.0085) |\n| 20          |  0.8285 (0.2324) | 0.4409 (0.0068) | 1.7109 (0.0232) |\n| 50          | 59.6913 (0.7021) | 0.7006 (0.0068) | 3.6558 (0.0433) |\n| 100         | 82.2349 (1.0051) | 0.9790 (0.0071) | 7.3584 (0.0457) |\n\nTable 1: Estimation Error vs. Parameter Dimension on multivariate Gaussian Model, std. err. in brackets\n\n\n|              NN FSM|       Linear       FSM |             KDE-SP |             Opt |\n|-----------------|-----------------|-----------------|-----------------|\n| 2.4963 (0.1635) | 2.0635 (0.0397) | 6.1476 (0.1594) | 1.6575 (0.0726) |\n| 2.2283 (0.1998) | 1.9612 (0.0533) | 6.2249 (0.0196) | 1.3495 (0.2286) |\n| 0.9394 (0.0634) | 0.7961 (0.0367) | 3.0341 (0.4033) | 0.6230 (0.0356) |\n| 0.8271 (0.0588) | 0.7863 (0.0427) | 1.5439 (0.1210) | 0.4170 (0.0142) |\n\n\nTable 2: Prediction Error on Generator Inversion Task, std. err. in brackets\n\nOur empirical findings show that a straightforward NN model is not able to match a linear model, and it remains an interesting research direction about how best to incorporate a NN model in our framework. One promising extension, inspired by the reviewer's comments, is a hybrid approach that incorporates both a NN and linear model. During an offline pretraining phase, we use our FSM method to train the NN model over a wide range of simulated data. Then, during inference time, we can leverage the NN model as an encoder and use the linear FSM model on the encoded observations, allowing us to maintain the computational efficiency of the linear model during inference. We will include a further discussion of this extension, as well as our discussion on the trade-offs between a linear and NN score model in a future revision of our manuscript.\n\n\n> Amortization\n\nThe reviewer is correct that our method is sequential and not amortized. Our focus in this line of research was to develop a fast, simulation-efficient approach, and a non-amortized method allows us to maximize the simulation efficiency by only simulating where necessary. To benefit from a training dataset, the hybrid approach combining both a pretrained NN encoder model, and a linear FSM model for observations during inference time as we proposed in the previous paragraph would be a powerful \"semi-amortized\" hybrid scheme. A fully amortized score estimation approach was developed in [1] for full Bayesian inference. Due to length constraints here, we provide further elaboration in our revised manuscript.\n\n> Tuning of $\\sigma$\n\nWe use the following *$\\sigma$ calibration scheme*: In a pilot phase before the main optimization, we perform grid search over a range of hyperparameter values. For each candidate $\\sigma$:\n\n1. We run SGD to obtain a parameter estimate $\\hat{\\theta}_\\sigma$ corresponding to the candidate $\\sigma$\n2. Using $\\hat{\\theta}_\\sigma$, we draw simulated samples, and calculate the average Euclidean distance between the simulated samples and observations\n\nFinally, we select the hyperparameter $\\sigma^*$ that minimizes the average distance. This principle of matching simulations to observations underpins many intractable likelihood methods.[2]\n\nWe would like to provide some insight into the *optimal selection of $\\sigma$*. The optimal choice of $\\sigma$ is intimately linked with the curvature of the log-likelihood function, and dictates a bias-variance tradeoff for the FSM estimator (Section 5.2). Using the definition of the Lipschitz constant $L$ in Theorem 5.2,\n* For a *strong curvature* log-likelihood (Larger $L$), the bias term dominates, and a smaller $\\sigma$ is needed to control the bias.\n* For a *weak curvature* log-likelihood (smaller $L$),the variance term dominates, and a larger $\\sigma$ is needed to reduce the variance.\n\nWe provide empirical evidence for this using a numerical study of a Gaussian model, $x \\sim N(\\theta, \\tau^2)$, where $\\tau^2$ governs the curvature. We provide the MSE between the FSM and the true Fisher score and the variance of the FSM estimator.\n\n|     σ |     MSE |   Var. |\n|-------|---------|-------------|\n| 1e-02 |  1.3981 |      1.3744 |\n| 5e-02 (*)  |  0.8202 |      0.2878 |\n| 1e-01 |  1.3539 |      0.1072 |\n| 5e-01 |  7.1327 |      0.0057 |\n| 1e+00 | 10.2868 |      0.0009 |\n| 5e+00 | 14.5088 |      0.0000 |\n| 1e+01 | 15.2296 |      0.0000 |\n\nTable 3: Strong curvature ($\\tau = 0.5$), optimal $\\sigma$ is small\n\n\n|     σ |    MSE |   Var. |\n|-------|--------|-------------|\n| 1e-02 | 0.9871 |      0.9839 |\n| 5e-02 | 0.1755 |      0.1739 |\n| 1e-01 | 0.1015 |      0.1003 |\n| 5e-01 | 0.0150 |      0.0146 |\n| 1e+00 (*)  | 0.0085 |      0.0055 |\n| 5e+00 | 0.0186 |      0.0004 |\n| 1e+01 | 0.0324 |      0.0001 |\n\nTable 4: Weak curvature ($\\tau = 2.0$), optimal $\\sigma$ is large\n\nFurthermore, we note that the variance of the FSM estimator decreases monotonically as $\\sigma$ increases, in line with our theoretical analysis. \n\nHowever, in practice, in the intractable likelihood setting, information about the curvature is unlikely to be available. Obtaining a reasonable estimate for the curvature, in the context of our goal of gradient estimation, would likely require further methodological development. Our proposed calibration scheme is a pragmatic and data-driven heuristic, similar in spirit to hyperparameter tuning for many ML models that controls a bias-variance tradeoff. We will ensure that the details of this calibration scheme and its rationale is fully clarified in our revised manuscript.\n\n> Adaptation of $\\sigma$ \n\nA typical choice of adaptation schedule would be to anneal the hyperparameter $\\sigma$ towards 0, which would ensure that the bias induced by the smoothing vanishes asymptotically, allowing recovery of the true MLE, and this is the standard approach in black-box randomized optimization [3]. In practice, however, this approach is complicated by an increase in the variance of our FSM estimator and numerical instability when $\\sigma$ is decreased too aggresively. Furthermore, an adaptive schedule introduces new hyperparameter such as the decay rate of the schedule, increasing the complexity of the hyperparameter calibration. Our preliminary experiments showed that standard annealing schedules did not outperform a fixed $\\sigma$, so we opted for the simpler approach. This problem of finding a good schedule for $\\sigma$ is analogous to the tuning of the learning rate for SGD. A fixed learning rate is often a strong baseline, but many sophisticated learning rate schedules are available, and we are interested in exploring the incorporation of such schemes for an adaptive $\\sigma$. We will include this discussion on adaptive $\\sigma$ in our revised manuscript.\n\n> Scalability of FSM\n\nWe extend the numerical study in Section 7.1 of our manuscript to assess parameter estimation accuracy up to $d=100$, keeping the simulation budget fixed.\n\n| Param Dim   |              FSM |             KDE-SP |\n|-------------|-----------------|-----------------|\n| 2           | 0.1239 (0.0064) | 0.1257 (0.0066) |\n| 5           | 0.2108 (0.0072) | 0.2240 (0.0071) |\n| 10          | 0.3038 (0.0080) | 0.6999 (0.0325) |\n| 20          | 0.4438 (0.0069) | 1.1422 (0.0250) |\n| 50          | 0.7047 (0.0069) | 4.5179 (0.0281) |\n| 100         | 0.9898 (0.0074) | 7.6113 (0.0276) |\n\nTable 5: Estimation Error (MSE) vs. Parameter Dimension on multivariate Gaussian Model, std. err. in brackets\n\nAs Table 5 shows, the FSM estimator's error grows relatively modestly with the dimension compared to KDE-SP. The FSM estimator is more robust as it directly estimates the Fisher score, instead of incorporating a two-stage estimation process (KDE, and then finite-differences) for the KDE-SP baseline. This is also in line with the theoretical analysis provided in Theorem 5.2 of our paper, where the bias of the FSM estimator scales as $O(\\sqrt{d})$. \n\nReferences:\n[1] Sharrock L. et al. “Sequential neural score estimation: likelihood-free inference with conditional score-based diffusion models,” ICML 2024.\n\n[2] Cranmer K. et al. “Frontier of simulation-based inference,” Proc. Natl. Acad. Sci. USA, vol. 117, no. 48, pp. 30055-30062, 2020.\n\n[3] Duchi J.C. et al. “Optimal rates for zero-order convex optimization: two function evaluations,” IEEE Trans. Inf. Theory, vol. 61, no. 5, pp. 2788-2806, 2015."}]

Author rebuttal:

Meta-review: The authors introduce a general method for maximizing likelihood when it is intractable to compute, in the simulation-based inference setting (when samples from any p_theta are obtainable), based on Fisher score matching (with linear parameterization). The paper is supported with both theory and experiments. After discussion, all authors voted for acceptance. The main strength is that the paper gives a new, very generally applicable method for a generally important problem and demonstrates its effectiveness. The discussion centered around answering reviewer requests for experiments (in the high-dimension setting and in comparison to neural network models) and hyperparameter tuning, both of which reviewers found were answered satisfactorily.
Decision: Accept (spotlight)