Dimension reduction via score ratio matching
Authors: Ricardo Baptista
Michael Brennan
Youssef Marzouk
Venue: Accepted by TMLR ()

Abstract: Gradient-based dimension reduction decreases the cost of Bayesian inference and probabilistic modeling by identifying maximally informative (and informed) low-dimensional projections of the data and parameters, allowing high-dimensional problems to be reformulated as cheaper low-dimensional problems. A broad family of such techniques identify these projections and provide error bounds on the resulting posterior approximations, via eigendecompositions of certain diagnostic matrices. Yet these matrices require gradients or even Hessians of the log-likelihood, excluding the purely data-driven setting and many problems of simulation-based inference. We propose a framework, derived from score-matching, to extend gradient-based dimension reduction to problems where gradients are unavailable. Specifically, we formulate an objective function to directly learn the score ratio function needed to compute the diagnostic matrices, propose a tailored parameterization for the score ratio network, and introduce regularization methods that capitalize on the hypothesized low-dimensional structure.  We also introduce a novel algorithm to iteratively identify the low-dimensional reduced basis vectors more accurately with limited data based on eigenvalue deflation methods. We show that our approach outperforms standard score-matching for problems with low-dimensional structure, and demonstrate its effectiveness for PDE-constrained Bayesian inverse problems and conditional generative modeling.

Published in Transactions on Machine Learning Research (04/2025)
Dimension reduction via score ratio matching
Ricardo Baptista"
rsb@callech.co edu
California Institute of Technology Pasadena, CA 91125 USA
Michael C. Brennan'
mebrenn edu
Massachusettss Institute of Technology Cambridge, MA 02139 USA
Youssef Marzouk
ymars@mit edu
Massachusetts Institute of Technology Cambridge, MA 02139 USA
Reviewed OpenReview: https:// openrreveee, net/ forum? id= mub mXXZBa
Abstract
Gradient-based dimension reduction decreases the cost of Bayesian inference and probabilis- tic modeling by identifying maximally informative (and informed) low-dimensional projec- tions of the data and parameters, allowing high-dimensional problems to be reformulated cheaper low-dimensional problems. A broad family of such techniques identify these projections and provide error bounds on the resulting posterior approximations, via eigen- decompositions of certain diagnostic matrices. Yet these matrices require gradients or even Hessians of the log-likelihood, excluding the purely data-driven setting and many problems of simulation-based inference. We propose a framework, derived from score-matching, to extend gradient-based dimension reduction to problems where gradients are unavailable. Specifically, we formulate an objective function to directly learn the score ratio function needed to compute the diagnostic matrices, propose tailored parameterization for the score iat network, and introduce regularization met that capitalize on the hypothe- sized low-dimensional structure. We also introduce a novel algorithm to iteratively identify the low-dimensional reduced basis vectors more accurately with limited data based on eigen- value deflation methods. We show that our approach outperforms standard score-matching for problems with low-dimensional structure, and demonstrate effectiveness for PDE- constrained Bayesian inverse problems and conditional generative modeling.
1 Introduction
A central aim of Bayesian computation is to develop efficient algorithms for characterizing conditional dis- tributions, e.g., the posterior =X\===" of a parameter X E R", given an observation y" E IRm and their joint distribution xx,y The cost such procedures can become prohibitive with growing n and m, making inference difficult for high (or possibly computational burden is dimension reduction.
dimensional problems. One approach for mitigating this
In this paper, we interested in two types of low-dimensional structure that appear in many Bayesian inference problems, and more generally in probabilistic modeling. First is the notion that the target dis- tribution can be well approximated as a low-dimensional update of a dominating reference distribution. Second is the notion that. the conditioning variables can be replaced with low-dimenssooaal projections or summaries, A recent line of work (Brennan et al., 2020; Zahm et al., 2022; Cui & Zahm, 2021; Baptista et al., 2022; Cui Tong, 2022) has developed gradient- buused methods for ideutifying and exploiting both
* *Atthors contrihuted equally this work.

---

Published in Transactions on Machine Learning Research (04/2025)
types of low-dimensional structure in non non-Gaussian settings. These approaches construct specific diagnostic malrices containing averaged gradient or Hessian information of the posterior or joint log-density. Eigen- decompositions of these diagnostic matrices reveal informative projections of the conditioning variables and informed projections of the parameters, and also yield bounds on the crror of the resulting approximations to the posterior in terms of the rank of these projections. Specifically, the spectra of these matrices indicate to what extent given problem has the posited low-dimensional structure.
Here we propose framework derived from score-matching to extend these ideas to problems where gradients are unavailable. In particular, we introduce learning problem to approximate the gradient of the log-ratio of two densities, which we term a score ratio function. We show that score ratios are central objects within the diagnostic matrices described above, and that learning the score ratio can itself take advantage of low- dimensional structure in the problems at hand. Doing so yields novel architectures and algorithms that outperform standard score-matching approaches, and ultimately enables more efficient approximations of the targeted conditional distributions. Our main contributions are as follows:
1. We propose algorithms for uncovering two types of low-dimensional structure in probabilistic models (low-dimensional updates of a reference measure and low-dimensional conditioning) based on score ratio matching ($3).
2. We introduce a novel training objective, network parameterizations, and regularization methods lailored to our dimension reduction goals (54).
3. We develop a new algorithm that iteratively identifies a basis for the desired reduced subspaces more accurately with limited data ($5).
4. We demonstrate that our score ratio matching method better reveals low-dimensional structure com- pared to standard score matching, and that it enables more accurate and efficient high-dimensional approximate inference ($6).
Related work: Many previous papers have highlighted the benefits of gradient-based based dimension reduction in Bayesian computation, for methods ranging from MCMC (Cui & Tong, 2022; Cui et al., 2014; Constantine et al., 2016) to SVGD (Chen & Ghattas, 2020) to normalizing flows (Brennan et al., 2020; Cui et al., 2023; Radev al., 2020) ensemble filtering (Le Provost et al., 2022). The present. work is concerned with realizing these benefits in gradient-free settings, which include Bayesian inverse problems with complex forward models (Kaipio & Somersalo, 2006; Biegler et al., 2010), goal-oriented inference (Lieberman & Willcox, 2013; Berger et al., 1999), simulattoon-aase inference (SBI) more broadly (Cranmer et al., 2020). and the purely data-driven setting where only fixed set of samples from Tx2Y are given.
The second form of dimension reduction we seek (low-dimensional conditioning) is related to SBI-focused efforts of creating summary (ideally sufficient) statistics of the observations to improve the quality and efficiency of inference algorithms (Fearnhead & Prangle, 2012; Joyce & Marjoram, 2008; Nunes & Balding, 2010). So far, most of these methods have not been gradient- based, and do not explicitly take advantage of low-dimensional structure when learning statistics. For instance, Radev et al. (2020) learns an (arbitrary) summary neural network to compress the observation variable before it enters inference network. Brehmer et al. (2020); Alsing et al. (2018) define data summarics based on locally sufficient statistics, but only around one reference parameter value. More generally, hand-crafted summary statistics are often the norm in SBI (Sisson et al., 2018). These statistics often require expert knowledge to form or use specific selection criteria that are not tied to posterior quality, and hence may lead to issues with "poor quality fcatures," " as seen in Lueckmann et al. (2017). In contrast, the approach proposed here seeks optimal summaries that provide error guarantees on the resulting posterior approximation.
2

---

Published in Transactions on Machine Learning Research (04/2025)
2 Background
2.1
Gradient-based dimension reduction
We provide a bricf description of two gradient-based dimension reduction methods for which we will develop gradient-free versions. These methods are dala-free certified dimension reduction (CDR)* (Gui Zahm, 2021) and conditional-muttue information based dimension reduction (CMIDR) (Baptista et al., 2022). Both methods decompose the parameter and observation into two subspaces,
U, X U,X,, + UI where X, U x, X_ U[x
(1)
= V,Y, + 1 where =, V, Y, Y1 = V
(2)
where U= = [U, Ui] E Rax" and V = [V, V, E Rmxx, s S m.
are unitary matrices, and U, E R r S n, and
Intuitively each method seeks transformations U and V such that the projected variables X, and Y, capture where the parameter is most informed and where the observation is most informative, respectively. Given decomposition in (2), Definition 2.1 below introduces low-dimensional approximation of the posterior that departs from known reference distribution using low-dimensional subspaces for the parameter and for the observation.
Definition 2.1. Let p be a chosen reference density on R Given unitary matrices U E Ruxa and E R and integers r S n, s SI m, let (U.V) denote set of distributions with densities of the form (xlr(=ly) f(U, V.UTx.VVTvvlll
for some f: IR.+ + R2y, where U, Raxr contains the first r columns of and V, E mxx contains the first columns of V. The class of distributions where only parameter dimension reduction is considered is denoted Dr, (U), and the class where only observation reduction is considered is denoted (V).
Propositions 2.2 and 2.3 below are key results for selecting the optimal transformations U and V, respectively. The transformations arise from minimizing an error bound for the resulting posterior approximation of the form in Definition 2.1. To derive the error bound, Proposition 2.3 assumes that the joint distribution *x,Y satisfies a subspace log-Sobolev inequality; we provide a definition in Appendix A.1. We note that multivariate Gaussian distributions, Gaussian mixtures, and uniform distributions on compact and convex domains all satisfy subspace log-Sobolev inequality; discussion. Conditions on the forward model of an inverse problem that are sufficient for πx, to satisfy a subspace log-Soboles inequality can be found in Baptista et al. (2022, Example 2).
Zahm et al. (2022, Assumption 2.5) for more
Proposition 2.2 (Modified from Section 3.3 of Cui & Zahm (2021)). Let p be the standard Gaussian density, and define the porameter diagnostic mafria
xx]x(x)u) p(x)
(x|r(211)) p(x)
HO Exx,Y
Va log (
V log
Raxn
(3)
Let (A ux) Rz0 IR" be the k-th eigenpair of H H with Ax 2 ... 2 Ax and take U [ui...., unl Then there erists x E (U) such that
Exy [ S Tr (II- ).,UJ)HHDD))
Σ
(4)
Proposition 2.3 (Modified from Theorem of Baptista et al. (2022)²). Define the observation diagnostic matric
w)) p(x)
HY Exx,y V log
V V.V,
V.Vlogg
log
mm
p(x)
We drop 'data-free' in the acronym will consider the dimension reduction problem initially presented in Zahm et al. (2022) for the posterior corresponding single realization of the observation. WWiil Baptista al. (2022) present joint, parameter and observation reduction results, focus here observation dimension reduction for of presentation.
3

---

Published in Transactions on Machine Learning Research (04/2025)
Assume that Xx, satisfies the subspace logarithmic Sobolev inequality with constant C(xx,)). Let (A Ub) E R2 Th be the k-th eigenpaör of H with X 2 ... 2 Am and lake V = ( -um]- Then, there exisls xxIr E (V) such that
Bzy [xx((xxy|[xxx)) & C(Tx,y)) Tr ((T V.V. V.V,") Hom)) C(EXY)) Σ
(6)
Note that there is no explicit assumption of subspace log-Sobolev inequality in Proposition 2.2 because, in the setting of parameter dimension reduction, we need this assumption only on the reference measure p. and since p here chosen to be standard Gaussian, the inequality is satisfied with log- Sobolev constant of one.
Propositions 2.2 and 2.3 have several practical implications. Given the diagnostic matrices H and H we should choose U, and V, as the leading cigendirections of these matrices. (As shown in Zahm et al. (2022), solving these eigenvalue problems minimizes a more general upper bound.) With this choice, we also have upper bounds on the approximation error (in KL divergence) incurred by parameter and observation roduction, which can be used select the reduced dimensions and s.
2.2 Approximating score functions
Score matching has recently appeared powerful unsupervised learning framework with applications to generative modeling (Song & Ermon, 2019; Song et al., 2021; Ho et al., 2020; De Bortoli et al., 2021) and Bayesian inference (Zhang et al., 2018; Pacchiardi & Dutta, 2022). The core task is to approximate the score function, the gradient of a log-density function, for various rownstream tasks. We direct readers to Song & Ermon (2019); Song et al. (2021) for an overview of score matching, especially the derivation of the objective function for learning the score, network training strategics, and application to generative modcling using Langevin sampling. In this section, we focus on conditional score-matching, where one approximates the score function of the conditional distribution Vo (xly)
Let wa: 1 R" be the neural network approximation of v log ( were 0 denotes the learnable parameters of the neural network. Score matching seeks to minimize a weighted L error between the true and approximate scores
1 J'wo)) = 2 BTx, V log x (aly) wo v)1
While minimizing this objective is intractable as requires access to the score function, approaches known as implicit score matching (Hyvärinen Dayan, 2005) and denoising score matching (Vincent, 2011) introduce reformulated objective functions for J (0) to learn the score based only on samples from given distribution, such xx,.. We will focus on implicit score matching in this work, although this is not limitation. Under mild regularity assumptions on the true and approximate score functions (see Assumption 3.1), Hyvärinen & Dayan (2005) showed that the objective J can be written as
| [1loo(, y)1 +Trll++Tr(r.... Tr(V (VVV.ol( v)) +C,
(wwo Exx.Y
where C' is a constant with respoct to the score network. This objective only depends on æ via the expoctation, and in practice one can learn wa(x, y) by minimizing an objective that estimates the expectation using joint parameter-observateeee samples {2(1), y b Tx,.. That is,
J(w)) = =Σ -olooooo (10)) Tr(V (((()))(()) (()).
As discussed in Song et al. (2020), directly evaluating the term Tr(V, ". (2(1), prohibitively expensive for even moderate dimensions n. practice, this is alleviated using the Hutchinson's 's trace estimator (Hutchinson, 1989).
(,:ti,,'' وزر) in the objective function is
4

---

Published in Transactions on Machine Learning Research (04/2025)
3 Score ratio matching
Here we propose a tailored score matching approach to construct the diagnostic matrices in Section 2.1. To do so, we consider approximations of the score ratio function
w(x, v) v. log (xxI (dlu)/oo))
where p is chosen to be tractable reference density (e.g., standard normal) and present several results related to the dimension reduction approaches of Section 2.1. First we note that both diagnostic matrices HeDn and НmMI can be expressed in terms of the score ratio function. Indeed, HEDR is defined explicitly in terms of w(x,)) y), while HE depends on the mixed gradient function
V log (Tx,y y (z, v)/p(a)) = vyw(z,y).
Therefore, approximating the score ratio allows us to construct the two diagnostic matrices and perform dimension reduction of both parameters and observations.
A naive strategy would be approximate V log (xl() and use to compute the score ratio as the difference,
V ( = V, log 1xg Va log p(z).
Instead, we take a different approach that leverages the (possible) low-dimensional structure of xx-y- Under the hypothesis that T is well approximated within VU, for some choices of U and V, we expect the score ratio, rather than the score itself, to be well approximated by a ridge function (Pinkus, 2015), i.e., a function that is constant for z, E Im(U1) x Im(V)).
Section 4, we describe a parameterization of ua((z; y) and a regularization method that are tailored to learning this low-dimensional structure.
Theorem 3.2 provides objective function that allows for direct approximation of the ratio using score ratio network We: Ru+m R". This and subsequent results rely on the following mild assumptions on the posterior and the score ratio approximation.
Assumption 3.1. Let wy denote the score (or score ratio) approzimation. We assume (i) wi is differentiable with respect x; (ti) (xly is differentiable with respect x; (ti) E ||we||? ca for all 0; (tv) E Exx, |l Vr (Fl))//))) < oo; and (v) (xly)=)
||e|| + for all 0 and y-
Theorem 3.2. Let wa: Ra + R" be the score ratio approzimation. Under the conditions of Assump- tion 3.1, we have the following equivalence of objectives:
1
[o - Va Vzlog (x]r(
|1²
Bxx,
E
2
p(x)
1
= =Enx,y
[u)( w0( + Tr(V 4) + V. log p(x) wa(2,)) v)
+ +C,
(7)
2
where C is a constant that only depends on the densities xxp and p.
We give a proof of the theorem in Appendix A.2. As in implicit score-matching, Theorem 3.2 shows that the approximation can be learned with an objective that does not explicitly depend on the true score ratio.
In practice, we replace the expectation on the right-hand side of (7) with a Monte Carlo estimate using joint samples (a((0, y(o) and thus define the optimization objective,
J(wo)= =
Σ wo(x60) ,()) wo(x (0) + Tr(V (1 (()) + V. Pog (ool0))TToooo wo( y(0
(8)
Once we identify the score ratio approximation wy that minimizes this objective, the diagnostic matricos can be approximated as
HO Rs AO = Exxx [wo(s, y)wo []
HY ~ M = Exx. Vywa (z, y) [V) [ywo(x,))]

---

Published in Transactions on Machine Learning Research (04/2025)
Estimators for the parameter and observation transformations U and V then given by the leading eigenvectors of H and H respectively. Given a reduced posterior that is constructed using the estimated transformations, a question that naturally arises is how the accuracy of the reduced posterior is affected by error in our approximation of the score ratio. The follow theorem provides an answer specifically for CDR.
Theorem 3.3. Let wo be an approsimation to the score ratio satisfying
|1
2
Bxx,y wo(x,y) () - V log
Sc.
p(z)
Define the approrimate parameter diagnostic matrix H Bax. [wo(2,, y)) y) Let ui) E
Rz0 R" be the i-th eigenpair of H&DR with A1 2 ... Xn and take U = .u.,...... url Then there crists
Exp E (U) such that
Exy [ Y [ + Σ An-
See Appendix A.3 for the proof. We note that the approximation error of the score ratio network e is generally unknown in practice. Theorem 3.3, however, establishes stability of the error in the approximate posterior distribution Txiy, which depends on the subspace U, with respect error in the score ratio used to identify that subspace.
Choice of the KL divergence the error metric We note that our methods specifically seek projec- tions of the parameter and observation that minimize bounds on the posterior approximation error mcasured by the data-averaged KL divergence E,y || directly following the results of Propositions 2.2 and 2.3. While the data-averaged KT. divergence is a natural and widely used choice, alternative error metrics such as a worst-case KL divergence (over Y), Wasserstein distances, or total variation distances may provide additional insights or be better suited for specific inference tasks. Developing methods that minimize these alternative metrics and studying the differences in the posterior approximations could be an interesting direction for future work, particularly in cases where robustness to tail behavior accuracy of specific posterior expectations are important.
Comparison with other common dimension reduction methods Here comment on the differ- ences between the present work and several common dimension reduction methods. The classical method of principal component analysis (PCA) seeks directions that maximize the variance of the data. Nonlinear methods such as t-SNE (Van der Maaten Hinton, 2008) and UMAP (MeInnes et al., 2018) can effec- tively preserve local neighborhood structures, making them powerful tools for visualization and marginal dimension reduction. However, these methods do not explicitly consider the joint structure between param- eters and observations- whichhhhhhhhh is crucial for inference problems, where preserving conditional relationships is essential. We reiterate that our framework seeks projections of the parameter and observation that are informed and informative, respectively, making it particularly suited to inference. In comparison, methods that act marginally (e.g., performing PCA, t-SNE, UMAP etc. on X and Y individually) do not capture the dependence structure between the parameter and observation and thus do not explicitly consider the underlying inference problem.
We note that in the case of linear-Gausian likelihood models, the projections proposed by Propositions 2.2 and 2.3 are equivalent to those found by canonical correlation analysis (CCA) based covariance information between X and Y; see Proposition 5 in Baptista et al. (2022). Insofar as the bradeett-based methods proposed in Baptista et al. (2022) can be viewed as a generalization of CCA to non-Gaussian likelihood models (c.g., by using gradient information of a possibly nonlinear forward model), our methods further generalize to problems where gradients of these likelihood functions are not available.
Choice of reference distribution We emphasize that the choice of the reference distribution our method is a degree of freedom, provided that its score function casily evaluated. Our method identifies directions in parameter space where the posterior differs from the reference distribution on average, and is

---

Published in Transactions on Machine Learning Research (04/2025)
most effective when these directions are contained within low-dimensional subspace. The marginal distri- bution of the parameter (i.e., the prior in Bayesian settings) is thus the natural first choice, when tractable. In cases where the prior can be easily mapped to a standard Gaussian, we recommend reparameterizing the problem accordingly for computational convenience (as done in Sections 6.1 and 6.2). However, when such a "whitening" is not easily performed, as in Sections 6.3 and 6.4, we suggest using a standard Gaussian reference distribution after centering and scaling the parameter. Future work will explore other methods for picking the reference distribution when the parameter marginal not available.
4 Structure-exploiting networks and regularization
We now describe a parameterization for (ng v) and regularization method that uncover possible low- dimensional structure in the target distribution. Recall that for x E (U, V), we have
(x f(U, x, V v)((z),
for function f: Rh+ Rso, where U, E Raxx contains the first columns of U, and V, Rmxrr contains the first columns of V. Then note that the score ratio and its observation gradient take the specific form
T Exir((z/l) U.g(U, v), Vy V. log p(x)
Txir( p(x)
V x log
where g Un) = Va. log (Ir,, and h(z,, U.) V V
log f Ua). Observe that the range of the score ratio lies within the subspace spanned by Ur, and the range of V. log(Fxx((elv)/o(z)) lies within the subspace spanned by V.. We encode these structural ansatzes into the parameterization of the score ratio network
we y) Wo WW x, W. y),
where Wa E TRax"' W, E R.xxs tha IR'''' - R²' a typical conditional score network for some r' S n and s' S m, and θ (0, W,, Wg) denotes all trainable model parameters. If W, and Wy converge toward matrices that have low (effective) ranks (much smaller than n and m) during optimization, the ranges of we(x,9) 3) and Vywe(z, y) are restricted accordingly. We note that this parameterization also allows to constrain the highest possible ranks of the estimated diagnostic matrices. That is, the rank of the parameter diagnostic matrices is bounded by r' and the rank of the observation diagnostic matrix is bounded by s'.
To promote We and Wy be low-rank when Txjy is expected to have low-dimensional structure, we pe- nalize the nuclear norms of W. and W, during optimization, as is commonly used for low-rank matrix estimation (Fazel, 2002). This leads to our final objective function
J(O) = J[we) A WW/]+ + Ae W,
where ||- is the nuclear norm, '1, A2 are regularization parameters, and J(we) is the objective function defined in (8). Algorithm 1 presents the complete score ratio matching procedure to estimate the parameter and observation Iransformations as well as the dimensions of the reduced variables based on a specified tolerance for the posterior approximation error.
5 Deflating score-ratio matching
Given that dimension reduction acts as pro-processing step before performing inference, we wish to reduce the cost of constructing the diagnostic matrices as much as possible. To this end, we would like to use relatively small sets of training samples and small networks that can be optimized easily. For problems with large parameter and observation dimensions, we find that a given score ratio network often estimates diagnostic matrices that accurately capture the leading eigenvectors of the true diagnostic matrix, but that higher-indexed eigenvectors inaccurate.

---

Published in Transactions on Machine Learning Research (04/2025)
Algorithm Single network score ratio dimension reduction
1: Input: Target data (00),y0)3 ~ xX,,, and user tolerances Ex,Ey > Solve
min J(0, Wy) 0 o,W.WWy
to obtain the score-ratio approximation w2( y)-
3: Estimate the diagnostic matricos
N wo(x(1) , (01))o0(x(1) 100)
N
j=1
N
HOMI
Σ
v wow y(0) T vewo(zl) yo0)
N
j=l
4: Compute the eigenpairs: (A ū.) E Ryo x R of HODR and (A ū.) E Rgo x R of HOMI
5: Pick r and 8 so that
xx <sx, Σ <
2 k>r
and set U, [ū ... 7,1, Vg 10 T.] output: U,, V.
To improve the estimation of more eigenvectors of the two diagnostic matrices, we propose an itcrative method inspired by eigenvalue deflation methods (see Saad (2011, Chapter 4) for comprehensive review). Rather than finding all of the columns of U, and V, using a single score ratio network, we construct a sequence of score ratio networks designed to capture progressively more vectors that form the transformations.
The following lemma defines deflated matrix that can be used to reveal higher-index index eigenvectors. Then, Proposition 5.2 shows how to construct modified score ratio to compute deflated diagnostic matricos.
Lemma 5.1. Let (Ai, Φi), for i = 1....,d, ,d, be the eigenpairs of symmetric matrir A E Rdxd with A A 2 Aa and ||oill2 = 1. Define the subunitary matric D, = = Φr] and projector P = 1 - 2, Then the deflated malriz Ā = PAP has eigenpuirs (1, for i= =l,,..,, and (2, for i = r
d.
See Appendix A.4 for the proof.
Proposition 5.2. Let pX E TRax and PY E Rmxm be orthogonal projectors. Define the deflated score ratio
wr(x,y) = PXV. log (xpr 4)
(9)
p(x)
Then the diagnostic matrices computed using lake the form
HCDR == Exx.y
v) = Px Hpx
N Eaxy [V,up(z,v) vywp(z,y)) v)] px H'px
where H. E Raxx and H' E pmxm are symmetric matrices. Nolably, this implies ker(PX) c ker( (HeDR) and ker (PY) ker (HY
See Appendix A.5 for the proof. Let PX and PY be orthogonal projectors onto the span of a few previously computed eigenvectiors, e.g,, those computed from the diagnostic matrices estimated in Algorithm 1. Then, Proposition 5.2 allows us to learn a new score network that produces deflated diagnostic matrices whose lcading eigenvectors are orthogonal
previously computed eigenvectors. This process can be then repeated by defining new projectors onto the span of larger collection of eigenvectors. Algorithm describes the complete numerical procedure using multiple deflation steps.
8

---

Published in Transactions on Machine Learning Research (04/2025)
Discussion computational costs Here, discuss the computational costs and scalability of the methods described in Algorithms 1 and 2. First, we emphasize that our framework functions a pre- inference step: we only require joint prior samples, {z(i), y ~ TX,.. Also, rather than needing sufficient sample size and network expressivity to preciscly recover the score (as one might want for postcrior sampling), empirically we find that much less needed to capture the low-dimensional structure underlying the problem, i.c., the dominant eigenspaces of the diagnostic matrices. This allows us to leverage that structure during inference, reducing the overall computational burden of the problem.
In general, approximating the score ratio encounters challenges similar to those of traditional score matching techniques. Assuming the score of the reference distribution can be evaluated efficiently, the computational cost of evaluating the traditional score matching loss and the loss function derived in Theorem 3.2 are the same up to a constant. That said, our network parameterization and regularization scheme are specifically designed exploit low-dimensional structure in the target problem. As demonstrated in Section 6.1, our approach exhibits better sample efficiency than traditional score matching. The iterative deflating method (Algorithm 2) can further improve sample efficiency, as each score ralio function is restricted to learning very low dimensional subspace.
Algorithm 2 Iterative-deflated score ratio dimension reduction
1: Input: Target data (I(0), 1 ~ xx,y, number of deflation steps T, number of eigenvectors to keep at each step & S <r',8' Initialize deflating orthogonal projectors
Px =x= P Im
3: for = 1........ T do
4: Parameterize the projected score ratio network the form of eq. (9)
5: Obtain the leading reduction vectors U. and Ve via Algorithm
6:
Update the reduction basis vectors Ú -U Uzl, V + [V Vil
7: Update the orthogonal projectors
px - px Ü,Ü,, p t px V.V.
8: end for
9: output: Ü, V
6
Numerical examples
We now present several numerical experiments to show the utility of our methods. Table I in the Supple- mentary Material summarizes all network and training hyperparameter choices for each numerical example. For the problems presented in Sections 6.1 and 6.2, the score ratio is tractable and thus we can construct the true diagnostic matrices H&D and H We thus evaluate the accuracy of our mothod by comparing the posterior approximation errors for the optimal basis transformations (4) and (6) with the following error bounds achieved using the learned bases Ü and V3;
EC FDR = 1 TT
EC = V.V. HH
We emphasize that this analysis is meant to validate our method and is only possible when the true diagnostic matrices are computable. For the problems of Sections 6.3 and 6.4, the true diagnostic matrix is not computable. In these cases we validate our method by showing we achieve better inference fidelity for the reduced problems as compared to the non-reduced problems.
We use tilde to distinguish learned bases, l.e. eigenvectors of the diagnostic matrices computed using score (ratio) network, from the eigenvectors the diagnostic matrices.

---

Published in Transactions on Machine Learning Research (04/2025)
6.1 Distribution with planted low-dimensional structure
First, we consider the "embedded banana" distribution where the data-generating process is given by
Xi 1),, X312 Modf.)), X3 N),
(10)
and X = RX', where R E 10 a random rotation matrix that sampled by computing the QR factorization of random matrix with standard Gaussian entries. In this case, the distribution лх is not function of any observation, and so we only consider parameter dimension reduction. Given that xx only departs from standard Gaussian along the coordinates (x), we expect our algorithm to find the subspace spanned by the two leading columns of R. In this example, we are able to compute the score ratio analytically and directly estimate the true diagnostic matrix HCDn-
In Figure Ia we plot the error bound E,DR for three different bases: (1) the eigenbasis of the true diagnostic
matrix; (2) the eigenbasis of the diagnostic malrix computed with our score ratio approximation; and (3) the eigenbasis of the diagnostic matrix computed with standard score approximation (i.e., a score network approximating V. log xx (x)). Both the score ratio and standard score networks were trained with N 1000 samples. For our method, we see that the error bound sharply drops at 2 to less than 10-3, We also see that our method yields considerably lower errors at each r compared to standard score matching. For a visual representation of the results, Figures 1b and lc show a scatter plot of additional held-out samples from лх (which were not used during training) and these samples rotated into our discovered basis Ú when taking d = In the learned basis, non-Gaussia in the problem has been concentrated to the first two directions, and the third direction is now essentially independent of the first two.
10'
Tue score
SR necwork lour mathod) standard network
Ertor bourdd decay
10-14
-5566 251001550
10-1
10-4
Reduced parameter dimension
-2560 25 53 75 25312150
-0
-2 0
1322 25 50 1000
(a)
(b) Before rotation
(c) After rotation
Figure 1: Embedded banana problem: (a) idcal parameter reduction error bound, score matching error bound (standard net work), and score ratio matching error bound (our method). Our method beiter captures the subspace where the target distribution म deviates from the reference distribution. Middle and right: scatter plots of held-out samples from the embedded banana distribution with d 3 (b) before (in the original random basis) and (c) after rotation by the learned basis U for X. We observe that non non-Gaussianity has been concentrated in the first two directions, and the third direction is now essentially independent of the first two.
6.2 PDE-constrained inverse problem
Next, we consider an inverse problem where the forward model involves an elliptic partial differential equation (PDE). The inference parameter describes the permeability field on a two-dimensionall domain D and the observations are pointwise measures of the corresponding pressure field. This so-called "Darcy flow" problem is widely-used test case in the literature on nonlinear Bayesian inverse problems (Stuart, 2010; Iglesias al., 2014; Cui et al., 2014). The permeability field el and the pressure field u are related by the following Poisson equation (Neuman & Yakowitz, 1979; McLaughlin & Townley, 1996; Carrera et al., 2005; Sun, 2013):
10

---

Published in Transactions on Machine Learning Research (04/2025)
. (eeVu) = =0, in D = [0, 1]2
0) 0 u(Er, 1)
(11)
Su = for & E (0.1}
8m
Neumann boundary conditions on the left and right boundaries impose a zero-flux condition modeling im- permeable layers either side of aquifer, and Dirichlet. boundary conditions on the top and bottom boundaries correspond to fixed pressures.
We endow the pogperrmaayyltyy & with a Gaussian prior N(0,E) where E is the Matérn covariance defined by the differential operator
C (ōl - ~
where A is the Laplacian operator, and 8 and control the variance and correlation of the prior realizations (Lindgren et al., 2011).
The inference parameters are the leading n = 100 coefficients of the log-permcability in the Karhunen- Loève expansion of the Gaussian prior. Given that the solution operator mapping the permeability field the pressure field, e ~> u, nonlinear, the inverse problem is nonlinear even without accounting for the exponential in the parameterization of the permeability, Let F denote the forward model mapping the parameter to 100 values of collected
uniform 10 x 10 grid [0.1, 0.9]?. Observations for the parameter follow the model Y = F(x) + e, where ~ N(0, 10-31,) and define the likelihood function Ty|x (-Iz). The resulting distribution of interest in this problem is the posterior xxp. We take the prior covariance parameters to be = 0.5 and 7 0.1. Figure 2 shows three realizations of the log- log-permeability field drawn from the prior, the corresponding pressure fields, and observations.
2.0
1.5
0.5
0.5
0.0
05
0.5
-15
-2.0
0.0
0.0
0.5
1.0
0
0.5
1.0
(a) Log- nog-permeability fields drawn from the prior
1.0
1.0
0.8
0.8
0.6
0.6
0.6
0.5
0.5
0.4
0.0
(b) Corresponding pressure fields and observations
Figure 2: Darey flow problem: example prior realizations of the log- perreeabiiitty fields and corresponding pressure fields and observations
Here we fix the sample size to N = 90000 and apply both the single network (Algorithm 1) and iterative- deflated (Algorithm 2) versions of score ratio dimension reduction, taking the number of eigenvectors ex- tracted at each deflation step to be = 1,2, or 3. Figures 3a and 3b compare the error bounds achieved by the optimal parameter and obsorvation bases identified by our methods. We see that parameter reduc- tion error achieved by the single score-ratio network separates from the optimal error earlier than the error achieved by the iterative deflated methods, and that overall e = performs the best by a small margin. Fach
11

---

Published in Transactions on Machine Learning Research (04/2025)
of the networks performs similarly well for observation reduction. Figure compares the four leading true
and approximated parameter and observation basis vectors for the & = iterative-deflated method, which visually match well (up to an immaterial change in sign).
true score
single SR network =
f=1
true score
single SR network 1
==2
f=2
t=
10a
Error bound decay
Error bound decay
100
10
6xx00
10
12
14
Reduced parameter dimension
Reduced observation dimension
(a)
(b)
Darcy flow problem: error bounds for parameter (a) and observation (b) subspaces identified by various methods. l is the number of basis vectors extracted at each deflation step.
Figure 3:
6.3 Flux, a quantity-offinterest inference example
In this example, we consider a simulation-based inference problem related to the Darcy flow problem of Sec- tion 6.2. The observation model is the same, and though the permeability remains uncertain, the parameter of interest is now the log of the flux across the lower boundary of the domain scalar-valued function of the permeability and the corresponding pressure. The log-flux defined
log / e ds R.
(12)
In
where IB- 00, 1] {0} denotes the bottom boundary. For this problem, we introduce source term on the right hand side of the Darey flow PDE,
V V. Ve f, in D= [0, 112 0) =0 u(f). 1) =1
(13)
Ou 0 for E {0,1},
where the selected source term f(E) = exx 20 cll?), with (0.2, 0.2), models groundwater recharge.
The posterior distribution interest in this problem is TQIY, which follows from the joint distribution TQ.Y. We note this case that evaluating the unnormalized posterior density (and thus evaluating the score function) is not tractable; more details are given in Appendix B.3. We seek reduce the dimension of the observations Y via Algorithm 1 using N 10000 samples. Figure 5a shows bounds posterior approxima- tion error that are estimated using learned score-ratio network. The fast decay in this bound suggests that the posterior should be well approximated using relatively low-dimensional obscrvation projections, c.g. = 4. To validate the utility of the corresponding learned basis V., we perform inference using a conditional normalizing flow (specifically an unconstrained monotonic neural network (UMNN) as described Wehenkel & Louppe (2019)) trained from samples of aQ,:: see Appendix B.2 for the implementation details. We build two such flows: one depending on the call-dimensional observations and the other depending on the reduced observations of dimension = 4 defined by the basis V. We compare the performance of the UMNN flows
12

---

Published in Transactions on Machine Learning Research (04/2025)
True 1
True 2
True
True 4
6.75
0.5
5.75
075
0.0
0.0
0.5
600
6.5
0.00
0.
3.00
0.5
0.00
-2.35
0.50
4.50
4.50
-.50
0075
-0.5
-5.7
-279
000
Estimated
Estimated 2
Estimated
Estimated
1.0
2.0
10
1.00
10
1.00
1.
1.0
0.50
0.50
642
025
0.25
0.5
-0.2
-1.25
-.50
0,3
6.5
1:
120
%00
100
10
1.00
0.3
05
100
(a)
True and estimated parameter reduction vectors
True
True
True
True
3.2
200
1.0
1.
8.50
1.0
100
6.50
0.5
5.50
6.25
025
0.25
025
-0.5
-1.25
-22
-0.0
-050
-1.50
-15
0.5
2.0
-100
05
-1,00
0.0
10
-1.0
Estimated 1
Estimated 2
Estimated 3
Estimated 4
0.5
9.75
0.5
6.66
05
0.00
35
5.50
0.5
4.00
-035
0.77
-6.35
-17
0.2
00
10,0
05
0.0
(b) True and estimated observation reduction vectors
Figure 4: Darcy flow problem: the four leading parameter and observation reduction vectors from the true diagnostic matrix and our score ratio estimates of them. We note that reduction vectors are equivalent up to change in sign.
learned with an increasing number of map training samples (independent from the samples used to learn the score ratio approximation) to true posterior predictive samples computed via MCMC,
Figure shows posterior samples generated using MCMC and approximate posterior samples generated by the conditional normalizing flows with and without dimension reduction. Samples from the reduced- dimensional inference procedure more closely match the true posterior samples, across the full range of training set sizes. For quantitative comparisons, generated samples, respectively, and let F and F denote their respective empirical cumulative distribution functions. The Kolmogorov-Smmnn Smiinov (KS) statistic is defined as ((2, Q) = sup, FF(g) F(g)). Figure 5b shows the decay of the KS statistios as a function the number of flow training samples, averaged across 10 different realizations of the observation variable Y. We see that reducing the observation dimension prior performing inference yields lower KS statistics on average for any number of training samples.
let. Q and Ô denote sets of MCMC samples and flow-
We perform MCMC high-dimensional latent. parameters X descrihed in Section 6.2, and then compute log-flux samples Q via (12).
13

---

Published in Transactions on Machine Learning Research (04/2025)
Full-rank, mean
Reduced,
Error bound decay
100.
0.8
std
std
KS-statistic
0.6
10-4
0.4
0.2
0.0
10-8
0
2
6 8
10
102
103
104
Reduced observation dimension
Number of training samples
(a) Flux: bound on the posterior error
(b) Flux: KS statistic vs. map training set size
Figure 5: Flux problem: (a) decay of the estimated error bounds that may used to select the dimension of the reduction observation; (b) average KS statistic versus the number of training samples. Reducing the dimension of the observation enables higher-fidelity inference.
Reference solution (MCMC)
Reduced map True realization
Full-rank map
5
5
4-
4-
3-
2
1-
1
0
0
1.0
0.5
1.0
0.5
0.0
0.5
1.0
1.0
0.5
0.0
0.5
1.0
(a) Full dimensional map,
(b) Full dimensional map, 1000 training samples
(c) dimensional map,
100 training samples
10000 training samples
5
4
4
3-
3
2-
2
1
0
1.0
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
(d) Reduced map, 4,
(e) Reduced map, 4,
(T) Reduced map, 4,
100 training samples
1000 training samples
10000 training samples
Figure Flux problem: comparison of samples from the full dimensional maps (without dimension re- duction) and reduced maps (with dimension reduction) to reference solution computed using MCMC for different amounts of map training data (left to right).
14

---

Published in Transactions on Machine Learning Research (04/2025)
6.4 Energy price modeling
We now consider conditional generalive modeling problem related to energy market modeling. In the United States, there are seven independent system operators (ISOs) that operate competitive wholesale eloctricity markets where gencrators and resellers can buy and sell power. For renewable generators and microgrid operators, these markets provide important revente opportunities that can improve the economic viability of new and existing projects. A generator's ability to captare these financial opportunities, however, is dependent on their ability to forecast wholesale electricity prices days ahead of time. Wholesale electricity prices vary by location and market and can be docomposed into the sum of three components: the energy price, the congestion price, and the cost of transmission losses. Energy prices are related to the cost of generation across an entire ISQ and are constant over the entire region. Congestion costs vary by location and are nonzero when the physical constraints of electrical components, like transmission lines or transformers, are reached. As a step toward probabilistic forecasting, we investigate the conditional relationship between energy prices of the PJM ISO, which services much of the north-mid Atlantic United States, and temporal and weather- covariates.
Our training dataset contains N 20000 energy price and observation realizations from Jamuary 2020 through February 2023. Energy prices regularly take values across several orders of magnitude. As a pre-processing step, we take the natural log of the energy prices, then shift and scale the samples by its mean and standard deviation. Each observation realization comprises temperature and cloud coverage data from 43 weather stations around the serviced region (see Figure Sa), as well as a temporal encoding T = (sin (ein(22Do///a0 /365), cos(2xDoY /365). sin(2:HoD/224 (os(22ooD/245 /24)), where DoY and HoD denote the day of the year and hour of the day of the prediction. This time encoding allows and daily patterns of the energy market. The parameter X and observation Y (i.e., temperature, cloud coverage and r) dimensions are 1 and
model the seasonal
90, respectively.
We use this dataset to estimate the score ratio for πxÿy and thus the diagnostic matrix HÊ via Algorithm 1. In Figure 7a, we observe a sharp drop in the estimated posterior error from observation dimension reduction, suggesting that the leading
modes should capture most of the information carried by the full observation about the parameter. Figures 8b and 8c plot the two leading basis vectors for the observation, and labels the observation components with the highest magnitudes. We see that the first vector primarily captures the temperature predictions at few key citics. The second vector gives the most weight to the temporal encoding. Figure 9 plots the temperature and cloud coverage contribution of the leading four basis vectors on a map of the United States. We see that many vectors have interpretable structures. For example, the temperature component of 03 soems to be a weighted average of the temperature prodictions across the weather stations.
As before, we validate the learned reduction vectors by comparing the performance of approximate inference using full- and reducee--dimensoonal conditional normalizing flows as function of the number of training samples from XX,.. While the ground-truth posterior distributions are not available in this problem, we do have the "true" " realized value of the energy price. Thus we report the continuous ranked probability score (CRPS),, (Gneiting & Raftery, 2007) for the approximate posterior samples, a predictive metric used to evaluate the performance of probabilistic forecasts, on average over 100 held-out energy-price samples. Figure 7b shows that the reduced- dimenssooll flows are more predictive on average than the full-dimensional flow. We also note that the CRPS is stable with as few as 100 training samples for the map. Figure 10 shows six histograms of posterior samples generated from full- and reduucedddinensional flows trained with 10000 samples for different realizations of the observation variable. In two of the examples, we see that the reduced flows may capture underlying multi-modality in the posterior that the full-dimensional flow (i.c., without observation reduction) does not.
dimeesioall
5Te CRPS defined as CRPS(x", F)= ( )dz where x" the true parameter value and F the cumulative distribution function of the approximate one-dimensional posterior distrihution.
15

---

Published in Transactions on Machine Learning Research (04/2025)
1.3
100
Full-rank Reduced
Error bound decay
1.2
Prior
1.1
CRPS
10-4
1.0
0.9
10-8
0.8
0.7
0
2
6 8
10
102
103
10
Reduced observation dimension
Number of training samples
(a) Energy price: bound on the posterior error
(b) Energy price: CRPS vs. map training set size
Figure 7: Energy price problems: (a) decay of the estimated error bounds that may used to select the dimension of the reduction observation; (b) average CRPS versus the number training samples. Reducing the dimension of the observation enables higher-fdeeityy inference.
1.0
0.5
0.0
pmmam
JEMP-KRDU
45
0.5
JEMPATOL
TEMP-KLM
-1.0
JEMP-KEWR
60 Data variable index
(b) V1
1.0
CCACRRy
cos tay sin tod
0.5
0.0
-0.5
CC-KAC
-100
-05
-90
-85
-80
-75
-70
-65
-1.0
20
40
60
80
(a)
Data variable index
(c)
Figure 8: Energy price problem: (a) the location of weather stations, (b/c) the two leading reduced vectors plotted as line graphs with the five highest magnitude components labeled in red.
16

---

Published in Transactions on Machine Learning Research (04/2025)
2.50
9.5
45
150
3.50
525
3.55
200
2.00
6.25
-0.00
-000
-35
-10
105
-90
-75
(a) Temperature contribution, 0i
(b) Temperalure contribution,
2.0
4.00
375
3.75
350
40
-1.00
-3.00
(c) Temperature contribution, 03
(d) Temperature contribution, Va
50
5.50
-0.25
-0.25
-5.00
-5.00
-0.5
-2.75
1.00
3.00
(e) Cloud coverage contribution, V/
(f) Cloud coverage contribution, 0'2
2.00
5.5
45
2.50
335
3.25
40
2.00
0.00
0.25
-0.00
1000
-35
-00
700
2000
-40
-75
70
(g) Cloud coverage contribution,
(h) Cloud coverage contribution,
Figure 9: Energy price problem: the temperature and cloud coverage contributions of the leading four basis vectors for the observation space
17

---

Published in Transactions on Machine Learning Research (04/2025)
Full-rank map
Reduced map
True realization
0.6
1.0
0.6
0.5
0.5-
0.8
0.4
0.4
0.3
0.3
0.2
0.4
0.2
0.1
0.2
0.1
2
2
1
2
(a) Example posterior
(b) Example posterior
(c) Example posterior 3
1.0
1.2
0.5
1.0
0.4
0.6
0.3
0.6
0.4
0.2
0.4
0.2
0.0
0.0
0.0
(d) Example posterior
(e) Example posterior 5
(f) Example posterior 6
Figure 10: Energy price problem: six estimated posteriors (corrosponding to different realizations of the observation) using both the full and reduced dimensional flows, trained with 10000 samples. The reduced flow samples seem more predictive of the true realization for each realized energy price. In (c.d), we observe that the reduced flow possibly captures bi-modality in the posterior.
7 Conclusions and future work
We introduce a framework for gradient-based dimension reduction of conditional distributions based score ratio matching. In the Bayesian setting, for example, our methodology identifies low-dimensional subspaces of the parameter and observation spaces that best capture how the posterior departs from a chosen reference distribution and how the observations can be reduced. These subspaces are accurately identified from a score ratio network using a paramelerization that exploits the score ratio's gradient structure, logether with low- rank matrix regularizers and an eigenvector deflation technique. The discovered ooy-dimensional structure is shown to improve the accuracy of amortized inference and prediction in gradient-free settings for diverse applications, including PDE-constrained inverse problems and simulation-based inference of energy prices.
We outline several promising directions for future work. While the posterior approximation error guarantees in this work are currently limited to reduced subspaces of the parameters and obscrvations, future work will investigate extensions to nonlinear dimension reduction; see Bigoni et al. (2022) for related work in a function approximation setting where gradients are available. Tn addition, will be valuable understand how the intrinsic dimension of the inference problem, along with the smoothness of the underlying densities, affect the number of samples N needed to estimate the score ratio and to uncover this low-dimensional structure. For example, if one can characlerize how the error in an estimated score ratio wi-
xx y (xll)) (alv)) p(x)
)²
Exx, wi (x. v) V. log
in expectation over the sample, scales with an appropriate intrinsic dimension and sample size, then The- orem 3.3 could be reformulated using user-specified parameters. We also note that Theorem 3.3 currently addresses parameter dimension reduction but not observation dimension reduction. Future work will aim establish an analogous result for observation reduction, thereby mirroring the duality between Theorem 2.2 and Theorem 2.3.
18

---

Published in Transactions on Machine Learning Research (04/2025)
Acknowledgments
The authors extend their gratitude to Olivier Zahm for his insights on the methodology and theoretical results this work, and Matt. Parno for his help constructing the energy price problem and interpreting Caltech,
the results of our method. RB acknowledges support from the von Kármán instructorship the US Department of Energy AEOLUS center (award DE-SC0019303), the Air Force Office of Scientific Research MURI on "Machine Learning and Physics- Based Modeling and Simulation" (award FA9550-200 0358) and the Department of Defense (DoD) Vannevar Bush Faculty Fellowship (award N00014-22-1-2777 held by Andrew M. Stuart. MB and YM acknowledge support from the US Department of Energy, Office of Advanced Scientific Computing Research, M2dt MMICC center under award DE-SC0023187.
References
Justin Alsing, Benjamin Wandelt, and Stephen Feeney. Massive optimal data compression and density estimation for scalable, likelihood-free inference in cosmology. Monthly Notices of the Royal Astronomical Sociely, 477(3))2874 2885, 2018.
Ricardo Baptista, Youssef Marzouk, and Olivier Zahm. Gradient-based data and parameter dimension reduction for Bayesian models: information theoretic perspective. arXi::2207.08670, 2022.
James O Berger, Brunero Liseo, and Robert L Wolpert. Integrated likelihood methods for eliminating muisance parameters. Statistical science, 14(1):1-28 1999.
Lorenz Biegler, George Biros, Omar Ghattas, Matthias Heinkenschloss, David Keyes, Bani Mallick, Youssef Marzouk, Luis Tenorio, Bart van Bloemen Waanders, and Karen Willeox (eds.). Large-Scale Inverse Problems and Quantification of Uncertainty. Wiley Online Library, 2010. doi: https:/ htps::/doo. 10..//0000000 9780470685853.
Daniele Bigoni, Youssef Marzouk, Clémentine Prieur, and Olivier Zahm Nonlinear dimension reduction for surrogate modeling using gradient information. Information and Inference: Journal of the IMA, 11(4): 1597-1639, 2022.
Johann Brehmer, Gilles Louppe, Juan Pavez, and Kyle Cranmer. Mining gold from implicit models to improve likelihood-free inference. Proceedings of the National Academy of Sciences, 117(10))5242-5544)) 2020.
Michael C. Brennan, Daniele Bigoni, Olivier Zahm, Alessio Spantini, and Youssef Marzouk. Greedy inference with struehure-exploiting lazy maps. Advances in Neural Information Processing Systems, 33:8330 8342, 2020.
Jesús Carrera, Andres Alcolea, Agustín Medina, Juan Hidalgo, and Luit J Slooten. Inverse problem in hydrogeology. Hydrogeology journal, 13:206-222 222, 2005.
Peng Chen and Omar Ghattas. Projected Stein variational gradient descent. Advances in Neural Information Processing Systems, 33:1947 1958, 2020.
Paul G Constantine, Carson Kent, and Tan Bui-Thanh. Accelerating Markov chain Monte Carlo with active subspaces. SIAM Journal on Scientific Compuling, 38(5):A2779 A2805, 2016.
Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-basee inference. Proceedings of the National Academy of Sciences, 117(48):30055 -0062, 2020. doi: htw.:///paa..117789111. URL attps://vwe pnas. org/ /t///do///////////// http///...10////a... pnas. 1912789117.
Tiangang Cui and Xin T Tong. A unified performance analysis of likelibhood-nffrmmed subspace methods.
Bernoulli, 28(4):2788-2 2815, 2022.
Tiangang Cui and Olivier Zahm. Data-free likelihood-infommt dimension reduction of Bayesian inverse problems. Inverse Problems, 3744004450,, 2021.
19

---

Published in Transactions on Machine Learning Research (04/2025)
Tiangang Cui, James Martin, Youssef M Marzouk, Antti Solonen, and Alessio Spantini. Likelihhoodinffommmll dimension reduction for nonlinear inverse problems. Inverse Problems, 30(11):114015, 2014.
Tiangang Cui, Sergey Dolgov, and Olivier Zahm. Scalable conditional deep inverse Rosenblatt transports using tensor trains and gradient-based dimension reduction. Journal of Computational Physics, 485: 112103, 2023.
Valentin De Bortoli, James Thornton, Jeremy Heng, and Arnaud Doucet. Diffusion Schrödinger bridge with applications to scoro-based generative modcling. Advances in Neural Information Processing Systems, 34: 17695 17709, 2021.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. nflows: normalizing flows in Py- Torch, November 2020. URL https://doi.org/// 5ttp///o..21//eoooo.. 4296288.
Maryam Fazel. Matrix rank minimization with applications. PhD thesis, PhD thesis, Stanford University, 2002.
Paul Fearnhead and Dennis Prangle. Constructing summary statistics for approximate bayesian computation:
semi-automatic approximate bayesian computation. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74(3):419 474, 2012.
Tilmann Gneiting and Adrian E. Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477):359 378, 2007.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusson probabilistic models. Advances in Neural Information Processing Systems, 33:6840-6851, 2020.
Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in critistics-Simlaatioo and Computation, 18(3):1059-1076, 1989.
Aapo Hyvärinen and Peter Dayan. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
Marco A Iglesias, Kui Lin, and Andrew M Stuart. Well-posed bayesian geometric inverse problems arising in subsurface flow. Inverse Problems, 30(11):114001, 2014.
Paul Joyce and Paul Marjoram. Approximately sufficient statistics and Bayesian computation. Statistical applications in geneties and molecular biology, 7(1), 2008.
Jari Kaipio and Erkki Somersalo. Statistical and computational inverse problems, volume 160. Springer Science Business Media, 2006.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR, 2015. URL http://arre .rrg waaa///2..... //////112..9.99/
Mathieu Le Provost, Ricardo Baptista, Youssef Marzouk, and Jeff Eldredge. Regularization of the ensemble Kalman filter in elliptic inverse problems. In progress, 2022.
Chad Lieberman and Karen Willcox. Goal-oriented inference: Approach, linear theory, and application to advection diffusion. STAM Review, 55(3):493-519, 2013.
Finn Lindgren, Håvard Rue, and Johan Lindström. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society; Series B (Statistical Methodology), 73(4):423 498, 2011.
Jan-Matthis Lueckmann, Pedro J Goncalves, Giacomo Bassetto, Kaan Öcal, Marcel Nonnenmacher, and Jakob H Macke. Flexible statistical inference for mechanistic models of neural dynamics. Advances in neural information processing systems, 30, 2017.
Y. Marzouk, T. Moselly, M. Parno, and A. Spantini. Sampling via measure transport: An introduction. In Handbook of Uncerlainty Quantification, R. Ghanem, D. Higdon, and H. Owhadi, editors. Springer, 2016.
20

---

Published in Transactions on Machine Learning Research (04/2025)
Leland Mclnnes, John Healy, and James Melville. Umap: Uniform manifold approximation and projection for dimension reduction. arXiv preprint or arXiv: NiN::1802.... 03426, 2018.
Dennis McLaughlin and Lloyd Townley. A reassessment of the groundwater inverse problem. Water Resowrces Research, 32(5):113111161, 1996.
Shlomo P Neuman and Sidney Yakowitz. statistical approach to the inverse problem of aquifer hydrology: 1. theory. Water Resources Research, 15(4):845 860, 1979.
Matthew A Nunes and David J Balding. On optimal selection of summary statistios for approximate Bayesian computation. Statistical applications in genetics and molecular biology, 9(1), 2010.
Lorenzo Pacchiardi and Ritabrata Dutta. Score matched neural exponential families for likelihood-free inference. Journal of Machine Learning Research, 23(38):1-71, 2022.
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshmi- narayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1 64, 2021.
Allan Pinkus. Ridge functions, volume 205. Cambridge University Press, 2015.
Stefan T Radev, Ulf K Mertens, Andreas Voss, Lynton Ardizzone, and Ullrich Köthe. Bayesflow: Learning complex stochastic models with invertible neural networks. IEEE transactions newral networks and learning systems, 33(4):1452-1466, 2020.
Yousef Saad. Numerical methods for large eigenvalue problems: revised edition. SIAM, 2011.
Scott A Sisson, Yanan Fan, and Mark Beaumont. Handbook of approzimate Bayesion computation. CRC Press, 2018.
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 density and score estimation. In Uncerlainly in Artificial Intelligence, pp. 574 584. PMLR, 2020.
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations, 2021.
A. M. Stuart. Inverse problems: Bayesian perspective. Acta Numerica, 19:451-555 2010.
Ne-Zheng Sun. Inverse problems in groundwaler modeling, volume 6. Springer Science & Business Media, 2013.
Laurens Van der Maaten and Geoffrey Hinton. research, (11),, 2008.
Visualizing data using t-sne. Journal of machine learning
Pascal Vincent. A connection between score matching and denoising autoencoders. Neural computation, 23 (7):1661 1674, 2011.
Antoine Wehenkel and Gilles Louppe. Unconstrained monotonic neural networks. Advances in neural infor- mation processing syslems, 32, 2019.
Olivier Zahm, Tiangang Cui, Kody Law, Alessio Spantini, and Youssef Marzouk. Certified dimension reduc- tion in nonlinear Bayesian inverse problems. Mathematics of Computation, 91(336):1789-118) 1885- 2022.
Cheng Zhang, Babak Shahbaba, and Hongkai Zhao. Variational Hamiltonian Monte Carlo via score matching. Bayesian Analysis, 13(2):485550) 2018.
21


Reviews:
[{"review_id": "DRRfqRdbZi", "reviewer": "Reviewer_n4op", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n* The paper proposes an approach to generalize score-based dimension reduction methods to cases where the gradient is not available in an analytical, or at least straightforward, form. To my knowledge, this is new.\n* The paper offers an enlightening theoretical analysis. In particular, Theorem 3.3 provides an upper-bound on the expected KL-divergence, both in terms of the error introduced by the score network and the low-dimensional representation.\n* The proposed algorithm is clear and not overly complicated. As far as I can tell, there are not too many tuning parameters to deal with. The main choices are: $\\varepsilon_X$, $\\varepsilon_Y$, architecture of score network. I find it very compelling that the ranks of $U$ and $V$ are not set manually, but computed to match the tolerances $\\varepsilon_X, \\varepsilon_Y$.\n* The numerical experiments are diverse and analyzed in a good level of details. I appreciated that the authors examin how the performance of their algorithms changes with the size of the training set. Overall, I found these extended \"case studies\" quite insightful.\n* Another strength of the method is that it provides an estimable error bound, which can be computed (as in Figures 6(a) and 6(b)).\n* Identifying projections of $X$ and $Y$ is insightful, as showcased in Section 6.4.\n* The paper is very well written.\n\nWeaknesses:\n* The error bounds are in terms of $\\mathbb E_{\\pi_Y} [D_{KL}]$, which is a reasonable choice, though other measures of disagreements between the approximation and the target may be consider. For example total variation and Wasserstein distances, which provide bounds on errors in the estimation of expectation values and quantiles. I recommend that the authors justify their choice of error measure (in terms of theoretical guarantees and mathematical convenience).\n* Theorem 3.3 seems incomplete: can the authors add an inequality which accounts for the low-rank of $V$ (thereby mirroring the duality between Propositions 2.2 and 2.3), and another inequality which accounts for $\\varepsilon$, $V$, and $U$ simultaneously? If this is not straightforward, I recommend adding a sentence or two as to why, and an acknowledgment that this could be future work.\n* The steps in Sections 4 and 5 are interesting, and make sense but, as far as I can tell, the paper does not evaluate the benefits in terms of computation cost and the error potentially introduced in Section 4.\n* In Figure 3, there does not seem to be a clear a difference between the trajectories (except for the \"true score\" benchmark). The results seem noisy. Can the authors add error bars or error shades? From Algorithm 1, it seems the main source of variation would be the data drawn from $\\pi_{X,Y}$, so maybe the experiment could be repeated with different realizations of the data.\n* In Section 6.4, the authors remark that the reduced map finds multiple modes, compared to the full-rank map. Why is this better? How do the authors know whether the modes are spurious or not, and is there a way to interpret these modes?\n\n* The examples in Sections 6.2 and 6.3 admit a tractable score. For the example in Section 6.4, the authors do not provide enough details for me to assess whether it is straightforward to evaluate the score. If this falls within an SBI framework, this should be more clearly emphasized.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies the application of gradient-based dimension reduction, in settings where the gradient is not (directly) available. The proposed method approximates a posterior $\\pi_{X|Y}$, by reducing the dimension of $(X, Y)$ using the transformations $(U, V)$. The proposed method proceeds in three steps: (i) approximate a score ratio, (ii) approximate diagnostic matrices, and (iii) estimate $U$ and $V$ using the leading eigenvectors of the estimated diagnostic matrices. \n\nThe authors bound the error of the approximate posterior, based on the error in the score ratio approximation and the rank of the transformation $U$. As far as I can see, there is no analogous result for the transformation $V$.\n\nNext, the authors refine their method and propose learning a neural network for the score ratio, by taking advantage of the assumed low-dimensional structure of the posterior. Finally, the authors propose using an eigenvalue deflation method to iteratively learn a score-ratio network and corresponding columns for $U$ and $V$. \n\nThe paper examines the proposed algorithms on three examples: a multivariate distribution with low-dimensional non-Gaussianity, an inverse PDE problem, and energy pricing problem.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See the above discussion on weaknesses.\n\nI also have a few minor comments and questions:\n* In Section 3, why is it \"naive\" to approximate $\\nabla_x \\log \\pi_{X|Y}(x|y)$ and then compute the score ratio?\n* In Theorem 3.3, the R.H.S on the inequality, $\\lambda_k$ can be evaluated but, as far as I can tell, $\\varepsilon$ cannot. This prevents us from computing the upper-bound in practice and this limitation should be acknowledged.\n* How does the procedure in Section 4 change the error in the score-function network?\n* In Section 6.2, the model is $Y = \\mathcal F(x) + \\varepsilon$ and $\\varepsilon \\sim \\mathcal (0, 10^{-3} I)$. I'm curious: how is the $10^{-3}$ scaling determined?\n* In Section 6.4, why are the authors believe there are $s = 4$ leading modes (say as opposed to $s = 3$ or $s = 5$)?\n* I recommend splitting Figure 6 into two figures, and trying to order the figures based on when they are mentioned in the text.\n* Can the authors report the formula for the continuous ranked probability score?", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "gqxhgb6Jni", "reviewer": "Reviewer_Sz6n", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths:**  \nThe paper is generally well-written and clear, with only a few minor typographical errors (listed below). The proposed method is interesting, and its application to both simulated and real datasets effectively showcases its utility.  \n\n**Weaknesses:**  \n- It is unclear how sensitive the method is to certain factors, such as:  \n    - The robustness of results to changes in the neural network training dataset. How robust are the projections when multiple training datasets are considered?  \n    - The choice of the reference distribution. What happens if a different reference is considered?\n    - The number of training samples used. While Figure 5 addresses sample size to some extent, the authors could expand on its implications. Under what conditions is a large sample size essential for the method to work well?  \n- If I understand correctly, the method performs linear projections through the eigendecomposition of diagnostic matrices; comparing its results with those from standard dimensionality reduction techniques like PCA, T-SNE, or UMAP would be helpful. Highlighting strengths and weaknesses relative to these methods would clarify the method's applicability.  \n- Could the author elaborate more on the computational cost implied by their method?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper explores two types of low-dimensional structures in conditional distributions, such as posterior distributions in Bayesian settings. The first type involves approximating the target distribution as a low-dimensional update of a reference distribution, while the second type replaces conditioning variables with low-dimensional projections or summaries. Detecting these structures is achieved by constructing specific _diagnostic matrices_ that capture information from the gradient or Hessian.  \n\nThe paper extends these methods to scenarios where gradients of distributions are unavailable, employing a score-matching procedure to approximate the gradient of a log-density. Specifically, using a neural network, the authors focus on approximating the gradient of the logarithm of the ratio between the target conditional distribution and a reference standard Gaussian distribution. This shift in focus to the score-ratio function is motivated by the ability to express diagnostic matrices in its terms. An iterative deflation method is also introduced to reduce the computational cost of constructing these diagnostic matrices.  \n\nThe proposed methodology is evaluated on various datasets, including both real and simulated examples, to demonstrate its applicability and effectiveness.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. Address the weaknesses outlined above.  \n2. Explicitly refer to the Supplementary Material in the text where appropriate (e.g., Section 6: \"Table 1 _in the Supplementary Material_ summarizes all network and training hyperparameter choices for each numerical example\").  \n3. Use a different letter for variables in Lemma 5.1 to avoid confusion with the network \\( w(x, y) \\).\n  \n4. I list some minor inconsistencies and typos:\n  - Section 2.2: The sentence \"contribution to the latter sections\" seems extraneous and should be removed.  \n  - Figure 1: Ensure consistency between figure labels and text. In the figure, the axes are labeled \\( X_0, X_1, X_2 \\), while the text uses \\( X_1, X_2, X_3 \\).  \n  - Double-check the bibliography: Address duplications (e.g., Cui et al., \"Likelihood-informed dimension reduction for nonlinear inverse problems\" is listed twice).  \n  - Equation (8): I assume the 1/N multiplies the entire expression, as MC is applied to approximate Equation (7) - use parenthesis to avoid confusion.", "broader_impact_concerns": "No one detected.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Yo0zw6mtJ9", "reviewer": "Reviewer_XpJ3", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### **Strengths**\nThe framework extends dimension reduction to gradient-free settings, addressing a critical gap in existing methodologies.\nIt also has applicability to diverse areas, including PDE-constrained Bayesian inference and conditional generative modeling.\n\nRegularization techniques effectively exploit the presumed low-dimensional structure of data.\nThe eigenvalue deflation algorithm enhances the identification of low-dimensional projections with limited data.\n\nExperiments demonstrate the method’s superiority over standard score-matching techniques in relevant scenarios.\n\n### **Weaknesses**\nThe score ratio learning and iterative deflation process may be computationally intensive for large-scale problems.\n\n While the experiments are thorough, additional real-world applications could better demonstrate practical impact.\n\nThe choice of regularization parameters might heavily influence performance, requiring careful tuning.\n\nThe method’s scalability to very high-dimensional or extremely large datasets remains unclear.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces a novel framework for dimension reduction in scenarios without gradient information, addressing limitations of traditional gradient-based methods. It derives an objective from score-matching to learn score ratio functions, enabling computation of diagnostic matrices without gradients.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Comparison with other dimension reduction techniques, such as t-SNE\n\nDiscussion about scalability and computation overheads\n\nMore real world demonstrations", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers agreed that the paper's claims were supported by evidence and the content would be of interest to the TMLR community, and I agree with their assessment.

Please remember to amend the text to conform to the TMLR style (in particular, change red text used for highlighting diffs back to black text) for the publication-ready version.
Decision: Accept (poster)