The Quotient Bayesian Learning Rule
Authors: Mykola Lukashchuk
Raphaël Trésor
Wouter W. L. Nuijten
Ismail Senoz
Bert de Vries
Venue: NeurIPS 2025 poster ()
Keywords: Bayesian methods
Fisher information
Heavy-tailed distributions
Information geometry
Natural gradient
Riemannian optimization
Variational inference
Primary area: probabilistic_methods
TL;DR: The Quotient Bayesian Learning Rule extends Bayesian updates beyond exponential families via quotient geometry.

Abstract: This paper introduces the Quotient Bayesian Learning Rule, an extension of natural-gradient Bayesian updates to probability models that fall outside the exponential family. Building on the observation that many heavy-tailed and otherwise non-exponential distributions arise as marginals of minimal exponential families, we prove that such marginals inherit a unique Fisher–Rao information geometry via the quotient-manifold construction. Exploiting this geometry, we derive the Quotient Natural Gradient algorithm, which takes steepest-descent steps in the well-structured covering space, thereby guaranteeing parameterization-invariant optimization in the target space. Empirical results on the Student-$t$ distribution confirm that our method converges more rapidly and attains higher-quality solutions than previous variants of the Bayesian Learning Rule. These findings position quotient geometry as a unifying tool for efficient and principled inference across a broad class of latent-variable models.

The Quotient Bayesian Learning Rule
Mykola Lukashchuk
Raphaël Trésor'
Wouter W. L. Nuijten'
ismail Şenöz² Bert de Vries
Department of Electrical Engineering, Technical University of Eindhoven, the Netherlands
Lazy Dynamies, Utrecht, the Netherlands
GGN Hearing, Eindhoven, the Netherlands
{m. lukashchuk, r. wre.rr.rrrr... w. w. 1. nuijten, bert. de. vries) nl 1senoz@lazydynami com
Abstract
This paper introduces the Quotient Bayesian Learning Rule, an extension natural-gradient Bayesian updates probability models that fall outside the expo- nential family. Building on the observation that many heavy-tailed and otherwise non-exponential distributions arise as marginals of minimal exponential families, we prove that such marginals inherit a unique Fisher-Be Rao information geometry via the muniient-manifold construction. Exploiting this geometry, derive the Quotient Natural Gradient algorithm, which takes steepest-descent steps in the well-structured covering space, thereby guaranteeing parameterizatiionnnnaaiia optimization in the target space. Empirical results on the Student-l distribution confirm that our method converges more rapidly and attains higher-quality solu- tions than previous variants of the Bayesian Learning Rule. These findings sition quotient geometry as unifying tool for efficient and principled inference across broad class of latinl-variable models.
Introduction
Statistical models with heavy-tailed likelihoods are indispensable when data contain outliers or ex- treme values that violate Gaussian assumptions. A prime example the Student-t distribution: its gngrees-of-freedone parameter lets the tails stretch or contract, providing the robustness practitioners require.
Fitting such models is considerably harder than specifying them. The latent-scale representation that makes the Student-t analytically convenient also renders Angectation-Maximizatii painfully slow in high dimensions, while naïve gradient methods stumble on the strong curvature induced by heavy tails. We therefore seek an algorithm that preserves the full tail flexibility of the Student-f and (ii) exploits the well-behaved geometry enjoyed by exponential-fame family (EF) distributions.
The Bayesian Learning Rule (BLR) of Khan and Rue (023) offers natural starting point: it frames inference as gradient ascent distribution space and, when the candidate posterior is an EF member, replaces ill-conditioned Euclidean steps with natural-gradient updates that follow Fisher geodesics. In its manifold formulation [Lin et al. (2020b), the EF's natural parameters form a Riemannian mani- fold equipped with the Fisher information metric, yielding both elegant theory and fast convergence. Unfortunately, Student-t distributions lie outside the exponential family, so standard BLR cannot be applied directly.
The novel extension of the BLR, The cases by using group actions, but the Lie-group BLR framework has yet to be extended to multi- variate settings a sigmificant limitation that our work specifically overcomes while maintaining the desirable information-geometri properties of the original BLR formulation.
BLR" Kiral et al. 2023 addresses some non-EF
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

The central insight motivating our work comes from fundamental property of the Student-t distri- bution: can be represented as the marginal of Normal- Wishart distribution, the so-called scale- mixture structure, first studied by [ndreww and Mallows 1974), where posterior candidates are parametrized through a latent "scale variable" that transforms an arbitrary base distribution. Normal- Wishart is an Exponential Family distribution. Moreover, Normal-Wishe minimal exponential family parametrization. This representation has been leveraged various con- texts, from mixture modeling (Peel and Mclachlan 2000) to robust regression [Lange et al., 1989]. primarily to facilitate EM-style algorithms through data augmentation.
distribution possesses the
We take this insight in new direction by exploring its implications for the geometric structure of the parameter space. Specifically, we observe that this marginalization relationship naturally induces a quotient manifold structure, where the Student-t manifold can be viewed as a quotient of the Normal-Wishan manifold under an equivalence relation defined by identical marginalized distributions.
Our key theoretical contribution lies in showing that the Fisher-Rao metric, which defines a natural Riemannian structure on statistical manifolds, can be extended from the Normal-Wishal manifold the Student-t manifold through this quotient relationship. Furthermore, by carefully choosing a base measure and family of scaling distributions in the scale-mixture, a wide range of non-EF models can be captured in this unified framework BBandorff--Nelsenn et al. 1982], enabling robust Bayesian updates generalizing our approach beyond the Student-t.
More precisely, prove that if a distribution is marginal of minimal exponential family, then its parameter space inherits a unique Fisher information metric structure as a quotient Riemannian manifold.
Building on this theoretical foundation, we propose an extension of the BLR that leverages the scale mixture representation and the quotient manifold structure. This insight leads us to develop the "Quotient Natural Gradient" algorithm, which efficiently optimizes on the Student-t manifold using horizontal lifts between manifolds. Our approach computes steps in the well-structured Normal- Wishart space and maps them appropriately to the Student-t parameter space through the established quotient relationship. In the remainder of this paper, we formalize these concepts, develop the nec- essary mathematical framework, and evaluate our approach empirically. We compare the Quotient Natural Gradient against both standard EM and naïve manifold optimization, demonstrating its ad- vantages in terms of convergence speed and solution quality. Our results highlight the practical value of this geometric perspective and suggest broader applications to other statistical models with similar latent variable structures.
2 Background and problem setup
2.1 Bayesian learning rule
Given a model parameter space Z and a loss ((z), the Bayesian Learning Rule (BLR; Khan and Rue 2023 ) optimizes over distributions rather than point estimates
q"' arg min E.y [1] [g arg min CEQ
-[(g), .Ee
where Q = { & E = c RRd parametrizes candidate posteriors, Hi[g] denotes the Shannon
entropy, and r > 0 a temperature. In other words, BLR minimizes negative ELBO, or maximizes EL.BO (we stick to maximization convention); just for the re-use in the future, we will define
Clq] rH[[g] E, [4]
The key component of the BLR the use of the natural gradient Amari [1998] in place of the naive Euclidean updates. Euclidean gradients ignore the underlying geometry of the set Q. Nutural-gradien descent Amari [1998] instead preconditions the gradient by the inverse Fisher information F (E), yielding steps of constant KI length and trajectories that are invariant to re- parameterization. More formally, the natural gradient update given by
F(C) Eg [V log ge (z) Ve log qe (z)
(3a)
7el = F(E) Ve((E),
(3b)
f+1 = & ar V (c((E))
(3c)

---

The obstacle in (3b) is computing and inverting F(E), an O(n²) operation that quickly becomes prohibitive (for full d-dimensional Gaussian, n ( and the cost (dd)))
For the BLR objective [ this cost disappears when QA member of the minimal, regular ex- ponential family, which means that A is open subset of the Fuclidean space and the sufficient statistics maintain independence (Jordan and Sejnowski 2001 Chapter 3J. The distribution be- longs the exponential family if
a (z) h(z) exp (A T(2) - A(A))
(4)
holds, where h is the base measure; 7 is the sufficient statistic; are the natural parameters, and A the log- log-partition function
A() = log / h(z) exp (AT)())))))))))) dz,
(5)
ensuring that is probability distribution. The dual (expectation) coordinate
Vx)(A)
(6)
yields the following gradient identity
VoC.
where (0) ((x(0) is the objective expressed in expectation parameters [Khan and Nielsen 2018 Thm. 1J. No matrix inversion required: one simply computes an ordinary gradient in 0. Unless explicitly stated otherwise, we assume that all exponential families considered this paper are minimal and regular.
2.2
Reparamcterization through marginalization
To make the general concepts of our construction more visible to the reader, we will refer to run- ning example, the so-called Normal-Gamma distribution, which serves as univariate preparation for the multivariate Normal-Wishal that is the main example underlying our experiments. Read- ers seeking an even simpler introduction may first consult Appendix [A] where the two-dimensional Negative Binomial example illustrates the quotient geometry with transparent visualizations.
The Normal-Gamma distribution is defined
T N (0-
(8a)
G(a,
(8b)
The Normal-Gamma distribution is a four-paran four-parameter distribution that defines a joint over two vari- ables and т. This four-parant parameter distribution defining joint over variables and exemplifies the broader class of scale-mixture distributions [Andrews and Mallows) 1974) that forms the foundation of our approach. More importantly, in the current context, is that the Normal-Gamma distribution (8) is a minimal exponential family distribution and its marginal over is a Student-t distribution, which lies outside of the exponential family. The mapping from the standard parametrization to the natural parametrization given by
12
0
2
2
The complete exponential family representation of the Normal-Gamma distribution is provided in Appendix B.5] Equation (54).
More generally, consider a joint exponential family density on (z, zv
qx (Zex)) h((zax) exp A T(Zex) A())),
A
where z Rd, zv Rdv d + dy. Marginalizing over zy defines
qe = (2,22 dzv.
(10)
See Amari 2016
Chap. 6] for thorough treatment of the dual affine structure.

---

and hence a surjection
π A - A A
(11)
For the running the natural-parameet vector À (A1, Ag, λ3, A4) [ef. Eq. (9), the projection maps the four- dimensional Carmal-Gamma space onto the three Student-t parameters (μ, o² 1) via
example (8) we have a correspondence z = 2 and zy = T. Writing
{ = ( =
d 2A
-Ae+A 2A (Aa+1/2) 2 +1)
(12)
Because many distinct A's yield the same triplet (1e, o2, (2, v), this exemplifies that in general the is not a bijection. But we obtained a minimal exponential family reparameterization our marginal that lies outside of the exponential family.
3 Marginal quotient structure
Natural-gradine steps are cheap in the joint exponential-famin space A but expensive in the marginal coordinates E, because the Fisher inverse can be avoided due to relation (7). Our plan is therefore to run the natural gradient scheme completely in A and afterwards marginalize the result of our procedure X by sending it back to x (A) E.
However, our aim is to minimize the BLR objective ①1) in the marginal parameter space E, rather than in the full natural-paramen
space A. This raises two questions:
(i) Is the outcome of the gradient scheme independent of the choice of representative A E π (()?
(ii) Does running natural-gradient descent in A and marginalizing each X (where & is the interate of the gradient scheme) actually minimize -L(E) in the marginal coordinates E?
The resolution hinges on quotient topology. Specifically, the marginal parameter space E can be viewed as the quotient set A/ ~' defined by the following equivalence relation:
A
T(>i) क (A2), A1, E A.
(13)
The equivalence classes (elements) of A/ ~' usually called fibres. We will align with this convention.
The quotient manifold theory ensures that (ii) is resolved if A/ ~' Riemannian quotient mani- fold and we project the gradient on the horizontal space with respect to F(A) Bouuma ] 9.9 and Def. 9.24]. A small background quotient manifold theory is provided in Appendix B
is an open subset of Euclidean space, and, by the moment parametrization assumption on (see Definition [ 3\t is an embedded submanifold of Euclidean space. Under this assumption, π smooth map between two embedded submanifolds, so the horizontal subspace can be simply expressed as
Hx* = (ker D (A))
where (T) is the differential of the smooth map between two Fuclidean spaces (see Boumal| 2023 Proposition 3.35J), and the orthogonal operator taken according the Riemannian metric
(Fisher metric) of the manifold A [Boumal 2023 Def. 3.10].
Definition 1 (Moment-per parametrized family). Let Q {ge E E} be k-dimensional family of probability densities on measurable space Z. We call moment--paameetio if there exist measurable moment functions m....., my 2 - R such that
(i) E is embedded k-dimensional submanifold of TR
(ii) For every & E the expectations e' (E) Ege [m,(z)] exist and are finite;
(iii) The mapping e E R* + (e) (E..... (E)) is a smooth bijection whose Jacobian has full rank k everywhere on E.
We refer to E (or m(E)) as the moment coordinates of qç-
Definition can be understood as labeling of each distribution by the values of finitely many expectations (e-g. the mean, the variance, the skewness.... where those expectations vary smoothly and uniquely according to =.

---

Many common families families-including including all Student-t's with degrees of freedom > 1-fit the pattern of Definition [ but some heavy-tailed laws such as the Cauchy (v 1) do not because their first moments are undefined.
Theorem resolves (i) from our problem statement because the theorem states that E is the quotient manifold of A. A proof of Theorem is given in Subsection B.3 of Appendix B
Theorem 1 (Marginalization yields a smooth quotient manifold). Let 9A be a minimal, regular exponential famity with parameter space A C. Rd. Suppose partition Zes Zv) is chosen so that the marginal family {qe heez obtained via π A - E is momentt-araaetrizeld (Definition 4
Then E is the quatient manifold of A induced by
Point (ii) is settled by Theorem 2 Theorem 2 shows that E is the Riemannian quotient manifold of A under the Fisher-Ras Rao metric and that the induced quotient metric coincides with the Fisher information metric of the marginal family itself. Putting the pieces together, Theorem shows that running the natural gradient in A projected on the horizontal space H. is equivalent running the natural gradient in E. The full proof is given in Subsection B.4 of Appendix B
Theorem 2 (Induced Fisher-Rao metric). Assume the setting of Theorem and equip the natural- parameter space A with its Fisher information metric Fx, Then:
(i) The map n, that project the Riemannian manifold (A, Fx) on E, induces a Riemannian quotient manifold structure E:
(ii) The Riemannian quotient metric on E is then the Fisher metric of
Summarizing, our approach replaces minimization of L by the natural gradient descent in E by natural gradient in A where at each step the gradient vector is projected onto the horizontal space as follows:
P (ker DR(A)) LE Va
(14)
where P is the orthogonal projection for metric (:; 1)x-- Theorems and matical equivalence between the two approaches.
ensure us of mathe-
The quotient Bayesian learning rule
We now translate the quotient-manifold theory developed Theorems 2 into concrete opti- mization procedure for evidence-lower-bond
(ELBO) maximization. Throughout this section let (z(z), & E denote the marginal variational family which we ultimately seek an optimum, and pick Eo E such that the prior factor of the model can be written p(z) (z).
Assume that the marginal distribution qg (z) arises by marginalizing minimal, regular exponential family [z, 2v), parameterized by natural parameters X A, over the extended latent variable Zexi zz zv) (see partition (TOp). The map π A E. A + E. induced by this marginalization
precisely the marginalization map defined earlier.
Choose representative Ào π (Eoo of the prior and define
N
((A) = E,y [log Dso (z, zv) log qs (z. 2v) + Σ E [ollooppz 22)]
(15)
Because C(X) is constant along every fibre π (E), moving within a fibre-thar is, in a "vertical" di- rection belonging to ker ( only the parameterisation, not the marginal distribution. By projecting each gradient step onto the Fisher-orthogonal complement HA via the operator (16). ensure that every update alters À solely through its image { イ (A). Hence the optimization trajectory produced in the joint space coincides exactly with the one obtained by running natural- gradient ascent on C(E) in
Let (AA) denote the expectation parameters of qu- For minimal exponential families the or- dinary gradient V. ((X) coincides with the natural gradient in the joint space; see Khan and Nielsen 2018). We therefore
compute VaL(A),
5

---

(ii) identify - (x)(X) V via the duality between 0 and X,
project V onto the horizontal subspace HA (ker D (A))
(iv) take a step of size Be in that horizontal direction.
Because the horizontal space is orthogonal to the fibers π (E), each update stays within a single equivalence class A, thereby realizing the quotient-natural-graade
flow guaranteed by Theo- rem 2 The procedure terminates when the horizontal component of the gradient falls below tol- crance e. At convergence, the optimizer X is mapped back to the marginal space via 5" = π (\*), yielding the desired posterior approximation q (z). The complete routine is summarized in Algo- rithm Il which we call the quotient Bayesian learning rule (QBLR).
An immediate question is how to make the step (IIII in the above scheme efficient. Let V23 ker (A) be the vertical sub-space and the differential of the marginalization map Jn (A), then its right null-space is the vertical subspace of our quotient. Pick some matrix (()) that forms basis of the Vx. Then the projection on the horizontal space (in the Fisher-Rao geometry) can be formed by
(A = I-X(A)(K(A)
F(N) F(A) K[A] K(A) F(X).
(16)
A short algebraic derivation of the identity (16) is provided in Appendix B Subsection
B.2
Crucially, the inversion involves only the dim VA x dim V1 matrix KJ FK; for the
case dim V, = 1 (Appendix C Subsection C.2) so (16) collapses to single scalar divide, and no
full Fisher inversion ever required. The general computational analysis of the expression (T6)
given in Appendix E.1
Algorithm 1 The Quotient Bayesian Learning Rule
Input: lifted prior parameters Ao canonical projection π : A E, data set = fz; ) ELBO defined in the lifted space (()) (15), step-size schedule (8.)120. tolerance
1: A + Ao
D initialize in the lifted (joint) space
2: repeat
3:
9e + VoL(A)
b compute natural gradient through the dual coordinates Eq. (6) b project onto the horizontal space, defined in Eq. (16)
4:
g + Projker + (90)
5:
X X + B.
natural-gradien ascent step
6:
until 11g I|2 <E
7:
E"
π (A))
8:
return marginal variational posterior qe- (-)
Student-t via Normal-Wisharr representation
In this section, we present an alternative approach heavy-tailed posterior approximation using the Normal-Wishat scale mixture representation. While Lin et al. [2020a] developed updates for Student-/ distributions through a curved exponential family formulation using the Normal-Inverse Gamma scale mixture, our approach leverages the quotient manifold structure induced by the marginalization map from the Normal-Wishall to the Student-l manifold. We first introduce Mishar parameterization and derive the explicit marginalization mapping the Student-t distribution. Then, we develop natural gradient updates that exploit the geometric structure this mapping, avoiding the need for reparameterization tricks. We demonstrate how method retains the computational efficiency exponential family updates while capturing the heavy-tailed nature of the Student-t distribution, comparing our approach with Lin's Normal-Inverse Gamma formula- tion both theoretically and empirically.
5.1 Comparing parameterization approaches
The fundamental difference between our approach and that of Lin et al. [2020a] lies how we represent the Student-f distribution. Lin's approach reparameterizes the Student-t as a curved ex- ponential family, ensuring a one-to-one correspondence between the scale mixture parameter space and the distribution space. Their key insight was finding specific parameterization that maintains this one-to-one correspondence, but the cost of working with a curved (non- monnmiii exponential family.

---

In contrast, our approach begins with the Hisha--Wsharrt distribution, which is minimal exponen- tial family distribution (see Appendix Cp. When marginalized, this yields the multivariate Student-t distribution through a many-to-one mapping, creating a quotient manifold structure. We can work directly the unconstrained minimal exponential family space, leveraging its well-understood ge- ometric properties. The quotient structure allows us to handle the redundancy in parameterization through the horizontal space projection.
The critical trade-off' between these approaches can be summarized as follows:
Lin et al. 2020a Curved Exponential Family + Student-f
(17a)
Minimal Exponential Family quotient
Ours
Student-t
(17b)
Mathematically, these approaches are represented as
[p(x/w) N(,,, wx) (p(w) ~ InvGamma(v, v)
Lin (NIG):
one to one
Talp E,
(18)
Ours (NW):
[((])S) ~ NN,, 1) (P(S) Wishart(r/
quotient
T ((/111
1 - d
-d
(19)
Instantiating the generic QNG-VI template (Algorithm D with the Normal-Wishan lift yields Al- gorithm 2 the algorithm is provided in Appendix C Following the scalar-NIG construction of Lin ctt al. [2020a], we apply the Bonnet- and Price-theorem analogues developed in Appendix C.4 to the parameters μ4, K, and V obtaining unbiased stochastic natural gradient the corresponding quotient manifold. For the shape parameter v, we construct an unbiased gradient estimator with the Implicit Reparameterization Trick of Figurnov et al. (2018).
For each sample 2n we draw an auxiliary scale matrix A, Wa(u; V), couple with the latent vector Zn, accumulate the data-fit gradients in the natural parameters (A1:4), add the analytic prior terms, and convert the result to expectation-space via the chain-rule identities Eqs. 77a) (70d).
The stochastic natural gradient is then projected onto the horizontal subspace (Alg. 2 Step 3) before a single ascent step in (A1:4) is back-transformed to (μ, V, K, v/).
Because every intermediate quantity depends on V and к only through the quotient-invariants
V-1 +
pp
7= == S=
K
the update is invrssentatioo-imvriinnn any smooth reparameterization that preserves the marginal Student-t-c.g. the joint rescaling (V, x) (V/e, ex) with e 0-produces the identical step on the Student-t manifold. The resulting trade-offs vis-à-vis the curved-NIG scheme of Lin et al. (0200a are summarised Table
6 Experimental validation
The full, Prrsion-pinned codebase that recreates every number in Table is archived at https: // anonymous. 4open. wcienccer//MRRB--CCT C735, A line-by-line description of the training pipeline, hardware, and hyper-parameters is given in Appendix [ all information needed for exact re- execution therefore lives in one place and does not clutter the main text. All experiments were conducted on a MacBook Pro (2021) equipped with an Apple MI Pro chip and 32 GB of memory.
We benchmark three variational-inferon (VI) optimisers that operate the same Bayesian logistic-regression model:
1. BBVI-NS the correfunction-fii free black-box VI variant of Roeder et al. 2017];
2. NG-LIN the natural -radiett approach of Lin al. [2020a];
3. NG-Ours the quotient intural-gradien optimizer introduced in this work, using a Nor- mal-Wiihaar marginal representation.
We run the methods for four different datasets that are taken from the UCI/OpenMI repostoory:
Breast Cancer Wisconsin (Diagnostic) 569 samples, 30 features (Wolberg al. 1993]
7

---

Aspect
Lin al. (scalar NTG)
Ours (quotient-NG. NW)
Scale-mixture lift
[( μ, (G((,,
(( ( ((5)
Minimality of joint EF
curved, rank-3
minimal, rank-4
Parameter-invariant
only to lincar re-labelling of same coords
any smooth paramctrisation (log-scale, NG, etc.)
Tail expressiveness
one scalar w = isotropic kurto- sis
per-direction (matrix) kurtosis
Need explicit F-1
no (mean- grad trick)
no (mean-grad trick + 2-scalar projection)
Extra work vs.
one outer-product (O(d))
Limit 00000
behavior unknown
smoothly becomes Gaussian
NG
Table I: Compact comparison of Lin's Student-t update and representatioon-ivaaiann quotient natural-
gradient step.
Pima Indians Diabetes 442 samples, features [Smith al. 1988).
Sonar (Mines vs. Rocks) - 208 samples, 60 features Gorman and Sejnowski 1988).
Spambase 4 6001 samples, features [Hopkins et al. 1999..
Each dataset split 80:20 (stratified) and Seature-standardizee using training statistics only.
For every (dataset, method) pair report test-set accuracy of the posterior mean together with the empirical standard deviation estimated from ten posterior samples; see Table
NG-Ours matches or surpasses BBVI-NS on three of the four benchmarks while requiring roughly one-tenth as many optimization iterations. The advantage is most striking on Sonar, where the richer Normal-Wishen marginal representation lifts accuracy significantly higher over BB and NG- LIN, confirming the benefit of a geometry- awwre update coupled with a more expressive variational family. Moreover, BBVI-NS marginals collapsed, so we do not benefit from the Bayesian procedure; we did obtain collapsed estimate.
Method
Metric
Breast cancer
Diabetes
Sonar
Spambase
Mean
0.9314 + 0.0210
0.7494 + 0.0473
0.7951 + 0.1760
0.8894 + 0.0078
BBVI-NS
Sample
0.9314 0.0210
0.7494 = 0.0473
0.7951 0.1760
0.8894 0.0078
Entropy
0.0000 + 0.0000
0.0000 = 0.0000
0.0000 = 0.0000
0.0000 0.0000
Mean
0.8919 + 0.0391
0.7022 = 0.0356
0.7476 = 0.0502
0.8904 0.0089
NG-LIN
Sample
0.9214 + 0.0209
0.7526 = 0.0260
0.8150 = 0.0116
0.8906 = 0.0090
Entropy
0.0696 + 0.0021
0.0026 + 0.0040
0.0905 + 0.0010
0.0112 0.0019
Mean
0.9711 + 0.0194
= 0.7292
0.9095 = 0.0417
0.8891 0.0124
NG-Ours
Sample
0.9599 0.0110
0.7791 + 0.0232
0.9142 0.0178
0.9057 0.0076
Entropy
0.1751 + 0.0153
0.1490 = 0.0100
0.1863 = 0.0011
0.1046 + 0.0060
Table 2:
Comprehensive evaluation of Bayesian logistic regression performance
four UCI/OpenML.
datasets. Each entry shows mean + standard error across 10 train-test splits with adaptive learning rates. Mean: test accuracy using posterior-maan weights (MAP estimation); Sample: test accuracy averaged over 100 posterior weight samples (capturing parameter uncertainty): Entropy: predictive entropy over test outputs nats (higher values indicate greater prediction uncertainty). BBVI-NS is the score-function-free black-hox VI of Roeder et al. [2017]; ; NG-LIN is the natural-gradien method of Lin et al. [2020a]: ; NG-Ours is the quo- tient natural-gradient optimizer introduced in this work. Note that BB VI-NS collapses near-point posteriors (entropy 0), while NG-Ours maintains the highest uncertainty quantification and achieves superior sample- based accuracy.
7 Discussion
Why horizontal-space projection matters. Properly removing the vertical component of the stochastic natural gradient stabilizes training: with projection, the ELBO converges to higher ELBOO values, whereas without it the optimization drifts and eventually blows up (Fig. TT. This empirical

---

lifted ELBO
-2
0
500
1000
1500
2000
step
Figure 1: Comparison of lifted ELBO convergence with and without Fisher-orthogonal projection in the Poisson-Gamma of Srrative-Binomil target (detailed Appendix A).. initialize five different repre- sentatives on the same fiber and optimize for 2000 iterations, estimating the lifted ELOO (15) with 5000 Monte Carlo samples each step. Curves display the reoossrepreseetttive mean with 11 standard deviation ribbon; the y-axis clipped the 2-98% quantile range to suppress outiers. With horizontal projection (blue), optimization remains stable and attains higher ELBO values; without projection (red), the flow drifts along fiber and eventually becomes unstable. The step-size schedule follows the Riemannian distance-over-gnane grrdggnadnnnss optimizer Dodd ct al. (2044 with initial distance estimate 0.005.
result matches the theoretical analysis of $ staying in the horizontal subspace keeps every iterate inside single marginal equivalence class, preventing spurious motion along the gauge orbit.
Position within the BLR landscape. Conditional-EF methods Lin et al. [000aa] rely on non- minimal embeddings and bespoke per-family updates, whereas quotient Bayesian learning rule (QBLR: see 4 uses a minimal embedding and a single closed-fone form natural gradient for all Nor- mal-Wishart scale mixtures. Lie-group BLR [Kiral et al.) 2023) enforees manifold constraints through group actions while keeping the Fisher geometry implicit; the published instantiation han- dles diagonal covariances, although a full-covariance extension is, in principle, conceivable but has not yet been demonstrated.
Toward mixture models. Student-t mixtures (GST-MMs) handle multimodal or heterogeneous data [Meitz al.) 2018 Revillon al. (0017] A drop-in combination of our horizontal projection with the variational mixture update of Minh et al. [2025) would yield fully natural-gradient GST- MM: per-component Normal- formal-Wishart factors follow our update, while the mixing weights use Minh et al.'s 's rule. Derivations and large-scale experiments are deferred to future work.
Breadth of applicability and open challenge. Scale mixtures, first systematised by Andrews and Mallows [1974] and greatly expanded by Barnndorf- et [1982). include the Laplace, exponential-povee power, and many other heavy-tailed families [West] 1987] Whenever the scale kernel admits a regular, minimal exponential-fame family lift, the quotient structure Eq. [13 emerges and
QBLR applies unchanged. The principal remaining challenge is to construct such lifts for exotic priors-ea e.g. skewed or asymmetric heavy-tailed laws-so so that our template can used out of the box.

---

Concluding remarks. A single geometric ingredient the Fisher-or orthoonnll projection onto the horizontal space- turns natural-graan -gradeeet BLR into stable, representation-free optimiser for broad class of heavy-tailed Bayesian models. Respecting the quotient structure is therefore not a pedantic luxury but a practical necessity for reliable optimisation.
8 Conclusions
We introduced the Quotient Bayesian Learning Rule (QBLR), which extends natural- -rradieet vari- ational updates to distributions that fall outside the exponential family yet arise as marginals of minimal exponential families. By casting the marginal parameter space Riemannian quotient, we showed that inherits a unique Fisher-Rao metric and derived the associated quotient nat- ural gradient (QNG). The algorithm performs steepest descent the well-conditioned covering space, projects the update horizontally, and thereby preserves parameterization invariance. A closed- form Dermmal-Gamm/Studenn
example makes the construction concrete, and empirical results on Bayesian logistic regression demonstrate faster convergence and superior predictive calibration compared with earlier BLR variants. The same geometric template is readily transferrable to wide class of scale-mixture priors and their mixture extensions, opening a path toward robust, heavy-tailed Bayesian learning at scale. While our method demonstrates strong geometric properties, its main limitation is computational complexity in high dimensions, which we suggest addressing through structured covariance proposals in Appendix E.2 see valuable future work.
Acknowledgements
We gratefully acknowledge financial support by the Dutch Ministry of Economic Affairs (PPS fund- ing), by the Dutch Research Council (NWO) and by hearing aid manufacturer GN Hearing, under contracts SLL-TTSM/21.0611/21111 (project: Auto-AR) and KICH3.LTP.20.06 (Project: RO-
BUST).
References
NIST digital library of mathematical functions. https:/tt/////mmf mist.gov.
P.-A. Absil, R. Mahony, and R. Sepulchre. Optimization algorithms on matrix manifolds. Prince- ton University Press, Princeton, N.J. Woodstockk 2008. ISBN 978-0-691-1329833 OCLC: ocn174129993.
Shun-ichi Amari. Natural Gradient Works Efficiently in Learning. Neural Computation, 10(2):
251-276, January 1998. ISSN 0899-7667. doi: 101116/09996669330111446 URL https: //doi.org/1n 10. 11220089769883007776
Shun-ichi Amari. Information Geometry and Its Applications, volume 194 of Applied Math-
ematical Sciences. Springer Japan, Tokyo, 2016. ISBN 978-4433-555977-1 978-4-431- 55978-8. doi: 10. 01779978--431-557788 URL htttt::////itt www.://.ssiiigg..oom//// wtt.tttg///aaa....pgggr....////a
778-44331-5997--
D. F. Andrews and C. L. Mallows. Scale Mixtures of Normal Distributions. Journal of the Royal
Statistical Society Series B: Statistical Methodology, 36(1 36(1):99-102, 909-102, 199--002,Sepeember 1974. ISSN 1369-
7412, 1467-9868. doi: 10.1111/j.2517-666 1974.tb00989.x URI. com hwtjjjssbaariiccee wtw.tiieee/////////00000000
O. Bamdorff-Nielsen, J. Kent, M. Sprensen, and M. Sorensen. Normal Variance-Mean Mixtures and z Distributions. International Statistical Review Revue Internationale de Statistique, 50(2): 145. August 1982. ISSN 03067734. doi: 10.2307/1402598 URL https:://ww. http:///www....rrrroooo org/
#tabiee/1025887orrginccceess
R. H. Bartels and G. W. Stewart. Algorithm 432 [C2]: Solution of the matrix equation AX XB C [F4]. Communications of the ACM, 15(9):220-82) 826, September 1972. ISSN 0001- 0782, 1557-7317. doi: (0.11455331153336611 URL https://dd.. acm. www./oogg//oo//oeoo 1145/
acm.
011573.361580
Georges Bonnet. Transformations des signaux aléatoires travers les systèmes non linéaires sans mémoire. Annales des Télécommunications, 19:203-220, 1964. doi: 10. 10.1007/BF03014773)
10

---

Nicolas Boumal. An Introduction to Optimization on Smooth Manifolds. Cambridge University Press. 1 edition, March 2023. ISBN 978-1-00-9166164 978-1-00-91661777 978- 170-11-00--
7. doi: 111017/9781009966111
URI. https://www.cmmrrddge..oo
(bentiffier/9711001166666 bopee//boooeeteeoobboo/booo/
Lawrence D. Brown. Fundamentals of statistical exponential families with
ap-
plications in statistical decision theory.
SPIE, January 1986.
doi: 10.
1214/lnms/12 216667577.
URL
https::/projecteucllie htw.t//eediid.rrr//eoooo//////aa wewellissswwhhhhhhe wwwh-aappliaaiinss--oooop ietennniiinntttssttiiieeee
wtt.etiisssiicceeeeeeeeetttttait eISBN-/10 2100011224//mmm//1546677555
Fundamentals-
Carlos A. Coelho and Ding-Geng Chen, editors. Statistical Modeling and Applications: Multivari-
ate, Heavy-Tailed Skewed Distributions and Mixture Modeling. Volume 2. Emerging Topics in
Statistics and Biostatistics. Springer Nature Switzerland, Cham, 2024. ISBN 978-3-0331-6621-- 0ttp://00000//778--01-------3330 --1-
978-3-031-66622--0 CRI.https: h/ttps://iink.... httt:::://iinnssprrnnggg. spwinger.n..... cwmm 1007 90000//97833311
Daniel Dodd, Louis Sharrock, and Christopher Nemeth. Learning-rate- stochastic optimization over Riemannian manifolds. In Proceedings of the 41st international conference on machine learning, ICML'24, Vienna, Austria, 2024. JMLR.org.
Morris L. Eaton. tics, volume 53, pages 302-334.
Chapter 8: The Wishart Distribution. In Multivariate Statis-
Institute of Mathematical Statistics, January 2007.
URL https::/ppoojectuucla. / institute-of wttp://eoeooactcaaaaaaaaaaaaatta catpatisttte-lltccuuu-ooooeeeeoo mon.uurrrr..oeeee-ooooooo...... mww.graph- series/ Mult hhlliiaariatte-staacssiic/c/eeea 1214/1mms/1196285114
doi: 10. 10.1214//mm//1112285111
Chapter 8- The Wishart- 1805 [cs, stat], May 2018. URL Disi---------iiiaaaaiiiiiittat Dit.ssarrrt-istttruuttoooo/o// 10. 1805.08498.
Michael Figurnov, Shakir Mohamed, and Andriy Mnih. Implicit Reparameterizztion Gradients.
08498 08498 arXiv:
Paul Gorman and Terrence Sejnowski. Analysis of hidden units in a layered network trained to classify sonar targets. Neural networks, 1(1):75-89, 1988.
Magnus R Hestenes, Eduard Stiefel, and others. Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards, 49(6)4009-436, 1952.
Mark Hopkins, Erik Reeber, George Forman, and Jaap Suermondt. Spambase. UCI Machine Learn- ing Repository, 1999. DOI: http:///oo..rr110.2443433//
Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in MkatistissSSmmlaatii and Computation, 1059-1005- 1989.
Maor Ivgi, Oliver Hinder, and Yair Carmon. DoG is SGD's best friend: A parameter- free dynamic step size schedule. In International conference machine learning, pages 14465-14499, 144... 2023.
Michael Irwin Jordan and Terrence Sejnowski, editors. Graphical models: foundations of neural computation. Computational neuroscience. MIT Press, Cambridge, Mass, 2001. ISBN 978-0- 262-60042-2.
Mohammad Emtiyaz Khan and Didrik Nielsen. Fast yet Simple Natural-Gradient Descent for Vari- ational Inference in Complex Models. arXiv: 1807.04489 [cs, math, stat], July 2018. URL
httt:://axivv http://arri..oo//abb////aaa 1807 04489 arXiv: 1807.04489
Mohammad Emtiyaz Khan and Hävard Rue. The Bayesian learning rule. Journal of Machine Learning Research, 24(281):1-466 2023.
Eren Mehmet Kiral, Thomas Moellenhoff, and Mohammad Emtiyaz Khan. The Lie- Bayesian Learning Rule. In Francisco Ruiz, Jennifer Dy, and Jan-Willem van de Meent, editors, Proceed- ings of The 26th International Conference on Artificial Intelligence and Statistics, volume 206 of Proceedings of Machine Learning Research, pages 3331-3352, PMLR, April 2023. URI. ccctps:/pprceeddingg mlr. wwwppp/vessp.... wtwp://///w.a....... html

---

Kenneth L. Lange, Roderick JA Little, and Jeremy MG Taylor. Robust statistical modeling using the t distribution. Journal of the American Statistical Association, 84(408):881-89 1989.
John M. Lee. Introduction Smooth Manifolds, volume 218 of Graduate Texts in Mathe- matics. Springer New York, New York, NY, 2012. ISBN 978-1-4419-998111 978-1-4419- 9982-5. doi: 10. 00.10007978-11440 1-9998---4419-992-55 URL https: http//iiik.... http://Lmm..prrrmeerr..o//o.o...
978- 11-11-4155 9982
Wu Lin, Mohammad Emtiyaz Khan, and Mark Schmidt. Fast and Simple Natural-Gradient Vari- ational Inference with Mixture of Exponential-family Approximations, November 2020a. URL
attt:///rxiv..oom www..rgg/bbbba. abs/ 1906. 02914 arXiv: NtVV:1906.0999 [stat].
Wu Lin, Mark Schmidt, and Mohammad Emtiyaz. Khan. Handling the Positive-Definite Constraint in the Bayesian Learning Rule. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine mlr.
Learning Research, pages 6116-6126. PMLR, July 2020b. htttttt::///poooo
gress/v119/1in221 html
Wu Lin, Mohammad Emtiyaz Khan, and Mark Schmidt. Stein's Lemma for the Reparameterization Trick with Exponential Family Mixtures, February 2025. URL http://aaivvvoor////////// 13398, arXiv: 1910.113999 [stat].
Mika Meitz, Daniel Preve, and Pentti Saikkonen, mixtrre autorerresseeeeeeeeeeee model Sis-distribution, May 2018. URL [econ].
model based on Student's
Tåm Le Minh, Julyan Arbel, Thomas Möllenhoff, Mohammad Emtiyaz Khan, and Florence Forbes. Natural variational annealing for multimodal optimization. arXiv preprint arXiv:2501.0460 2025.
Victor H Moll. Special integrals of gradshteyn and ryzhik: the proofs-volume Il, volume 2. CRC Press, 2015.
Pierre Del Moral and Angele Niclas. A Taylor expansion of the square root matrix functional, January 2018. URL http:://axxiii.i wtt.//g///////.ee.. hts///bs//1/ 1705. 08561 1eXiv:1705.0855 [math].
Kevin P Murphy. Conjugate bayesian analysis of the gaussian distribution. def. 1(2o2))16, 2007.
Atsumi Ohara, Nobuhide Suda, and Shun-ichi Amari. Dualistic differential geometry of positive def- inite matrices and its applications to related problems. Linear Algebra and its Applications, 247: doi: 10. 001101600224-3755594)0003833
31 31-53, November 1996. ISSN 0024-3795.
URL |https:
/ //www.sc www.ncedireett..... wtt.//ccccccce.eeeeeeeecee.eeeee
a Date: 2025505--22709::72333666
Read_Statuss: To
Read Read_Status Date:
D Peel and G Mclachlan, Robust mixture modelling using the t distribution. Statistics and Com- puting, 10(4):339-344 2000.
William D Penny. Kullback-lieble divergences of normal, gamma, dirichlet and wishart densities. Wellcome Department of Cognitive Neurology, 2001.
Guillaume Revillon, Ali Mohammad- -ijafari, and Cyrille Enderli. A generalized multivariate Student-t mixture model for Bayesian classification and clustering of radar waveforms, July 2017.
URL http::/araxxii..a////// org/ 1707 09548
Geoffrey Roeder, Yuhuai Wu, and David K Duvenaud. Sticking the landing:
Simple, lower-varianee gradient estimators for variational inference. Advances in Neural Information Processing Systems, 30, 2017.
Jack W Smith, James E Everhart, William C Dickson, William C Knowler, and Robert Scott Jo- hannes. Using the adap learning algorithm to forecast the onset of diabetes mellitus. In Proceed- ings of the annual symposium on computer application in medical care, page 261, 1988.
Michael K Tippett, Stephen E Cohn, Ricardo Todling, and Dan Marchesin. Conditioning of the stable, discrete-time Lyapunov operator. SIAM Journal on Matrix Analysis and Applications, 22 (1):56-65, 2000.
12

---

Mike West. On Scale Mixtures of Normal Distributions. Biometrika, 74(3):646-645 648, 1987.
William Wolberg, Olvi Mangasarian, Nick Street, and W. Street. Breast Cancer Wisconsin (Diag-
nostic), UCI Machine Learning Repository, 1993, DOI: https//dbb.oog//022444///5v
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: The abstract asserts that (i) heavy-tailed marginals of minimal exponen- tial families inherit Fisher-Rao geometry via a quotient manifold, (ii) this yields the parameterisationnnnaariaa QBLR update, and (iii) QBLR outperforms prior BLR variants on Student-t tasks. Claim (i) is proved in |3} (ii) is implemented in $ and exemplified in [5 (iii) is confirmed experimentally in $ Assumptions and limits are stated with Def- inition [ and revisited in 7 Hence the introductory claims accurately match the paper's
scope and results.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as it 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 scope and limitations of our method follow directly from Definition and are examined in greater detail in the Discussion (Section
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are 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 au- thors should reflect on 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 a few datasets or with few runs. general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the ap- proach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-ter system might not be used reliably 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 to ad- dress problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't it 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 a complete (and correct) proof? Answer: [es]
14

---

Justification:
This paper presents two main theorems (Theorems || and 2 with proofs provided in their respective subsections of Appendix [ Additionally, we introduce several auxiliary theorems that establish the gradients forms of the Algorithm ① as applied to the Norhal-Wishart distribution. These auxiliary results are proven in Appendix 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 theo- rems.
The proofs can either appear in the main paper the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be comple- mented by formal proofs provided in appendix 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 experimental results of the paper to the extent that it affects the main claims and/or conclu- sions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: The complete experimental setup is detailed in Section 6 the accompanying
code link is provided there, and implementation specifies appear in Appendix [
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 a dataset and/or model, the authors should describe the steps taken make their results reproducible verifiable.
Depending on the contribution, reproducibility can accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or the contribution is a specific model and empirical evaluation, it may be necessary either make it possible for others replicate the model with the same dataset, or provide access to the model. In general. releasing code and data also be provided via
often one good way accomplish this, but reproducibility detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of model checkpoint, other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all sub- missions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or a way to re- produce the model (c.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which thors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it 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 reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
15

---

Answer: [Yes]
Justification: The paper provides the link the anonymous repository at the beginning of Section 6
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" is acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g,, for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //nips.ou cc/ wwwppcpppicccccccccccccccccccccc public/ gwtp//uiiess//eg CodeSubmi wwwsssssssiiiioooeooo onPolicy) for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc. 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 to the paper) recommended, but including URLs to data and code 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: All required information to reproduce experimental results reported in Section is provided in Appendix D
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropri- ate information about the statistical significance of the experiments? Answer: [Yes]
Justification: Table 2 validates our approach and reports error bars, which are explained in its caption.
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 are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall
run 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).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should prefer- ably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
16

---

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 of range (e.g. negative
error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text. 8. Experiments compute resources
Question: For each experiment, does the paper provide sufficient information the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Answer: [Yes]
Justification: Section 6 specifies that we conducted all experiments MacBook Pro (2021) equipped with an Apple M1 Pro chip and 32 GB of memory.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as 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 't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics Answer: [Yes]
Justification: The authors declare to have done their best to adhere to the NeurlPS Code of Ethics. The research does not include human subjects, sensitive data, or societal dangers.
Guidelines:
The answer NA means that the authors have not 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- eration 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: There is societal impact of the work performed. 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 or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact spe- cific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a 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 intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitiga- tion strategies (e.g., gated release of models, providing defenses in addition attacks,
17

---

mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
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)? Answer: [NA]
Question:
Justification: The paper does not pose such a risk. Guidelines:
The answer NA means that the paper poses no 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 re- quiring that
adhere usage guidelines restrictions access the model 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. Licenses 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: [Yss]
Justification: All datasets used in our paper are properly cited in Section 6
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include a URL.
The name the license (c.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets released, the license, copyright information, and terms of use in the pack- age should be provided. For popular datasets, papprssitthooeeeeo com/ datasets has curated licenses for
datasets. Their licensing guide can help determine the li-
cense of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is 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 documenta- tion provided alongside the assets?
Answer: [NA]
Justification: The 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/modell 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 used.
At submission time, remember to anonymize your assets (if applicable). You can either create 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 pa- per include the full text instructions given to participants and sereenshots, applicable, as well as details about compensation (if any)?
18

---

Answer: [NA]
Justification: The paper does not involve crowdsoureing nor research with 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 contri- bution of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethies, workers involved in data collection, cura- tion, or other labor should be paid at least the minimum wage 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: The paper does not involve crowdsourcing nor research with human subjects. Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research
with human subjects.
Depending on the country in which research conducted, IRB approval (or equiva- lent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurlPS Code of Ethics and 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.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The concept of the reasearch does not involve LLMs as the core methods. Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components. Please refer to our LLM policy (https: cc/ co/Conf eeenneeeeces//e (tcpe/ceee///2///// 2025/LLM) for what should or should not be described.
19

---

in their standard form. Whether or not a lift can be A Illustrative Example
We illustrate the QBLR algorithm with low-dimensional example: negative binomial distribution family lifted to a (3)-dimensional scale-mixture distribution family, the so-called Poisson-Gamme distribution. The exponential hierarchical lifts for non-exponential family distributions like the Neg- ative Binomial distribution are plentiful in the literature (see Section 7. but they are often curved (non-minimal) in their form: the classical Poisson- Gamma parameterization in (r,p) ties together two natural coordinates and thereby lives on 2D submanifold of a 3D joint. In this section, we present an ad-hoc-yy -et fully constructive way to uncurve that representation by introducing one free scale, yielding minimal 3- 3-parameter exponential- family lift in natural coordinates. This is exactly the setting required by our quotient-manifole theory.
We first show how to build a minimal marginal exponential representation lift for the Negative Binomial distribution from the Poisson-Gamma distribution in Subsection [A.] and then we show how to use the established lift to instantiate Algorithm| for this specific scenario in Subsection A.2 Finally, we discuss the limitations of the same uncurving trick for the Laplace case in Appendix A.3
We thank the anonymous reviewer (bBHv) who proposed this section.
Note (erratum). In our rebuttal we incorrectly stated that the Laplace distribution could also be re- covered using this construction technique. However, the resulting marginal family is actually richer than the Laplace family alone. Whether the Laplace distribution can be obtained as a minimal lift through a different non-minimal representation, or requires a fundamentally different construction, remains an open question. We apologize for this oversight; see Appendix A.3 for details.
A.1 Building the minimal marginal lift
The framework of application of QBLR is quite general as many non-exponential family distribu- tions have some joint exponential family representation. However, our work is limited by an even stronger assumption: the existence of a minimal parameterization of these joint exponential family representations. While assume that at least some lift to joint exponential family is given, such representations are often curved inon-minine found is then crucial question for us.
By investigating this point, this section enables us to understand where and why our QBLR Algo- rithm should be applied.
The heavy-tailed distributions families that have an exponential family scale mixture representation are well documented in the literature. Regarding the particular case of being a scale mixture of normal distributions, (Adrews and Mallows (1974) provides necessary and sufficient conditions in their paper which is applied to Student-t, Laplace. and Logistic distributions. More examples can be found in Coelho and Chen (2024).
Many scale-mixture joints found in the literature come in a curved form-that is, their sufficient statisties are linearly dependent, so the family is not minimal. The textbook parameterisations of both the Normal-Wishart and the Normal- Normal-Exponential joints fall into this category. (By contrast, the Normal-Wishart lift we use-set see Appendix C-is is explicitly minimal; the distinction is made concrete in the Laplace example that follows.) Discrete over-dispersed families admit analogous lifts; in particular, the Negative-Binomial (NB) arises as Poisson-Gamma mixture.
Because a curved exponential family violates the minimal-regular assumption, it cannot serve as lift for QBLR unless one first "uncurves" it by adding extra, independent natural parameters. We now explain why this step is necessary and how those additional degrees of freedom restore minimality.
Curved vs. minimal lift for the Negative-Binomial distribution. The textbook mixture
L.a Poisson[k|] A)
NB r.p)
Gamma r, B
T-p Yp
(NB-curved)
Poisson
Gamma (rate)
curved when viewed as a joint EF (k; A): the joint density has three sufficient statisties, Ti log X, T2- = X, T3- = k, but only two free parameters (r, p) (equivalently, T3's 's's natural parameter is fixed at zero). Hence the Jacobian of the ceffccinnt-ttaisttinn map has rank 2 and Brown's minimality
criterion fails [Brown] 1986, Prop. 1.5].
20

---

Uncurving with one extra degree of freedom. the Poisson mean:
Introduce an independent positive scale >> on
J. Poisson cxa)
NB(k r.p)
Gamma A r. B 0) dx.
(NB-minimal)
scaled Poisson
Gamma (rate)
Now the joint admits a minimal EF representation with three independent natural parameters
77 (71 77, 1r-1, 1, -(3 c)2 log c),
sufficient statistics T (log A, X, k). base measure h(k;; 1-111 [x=0] > /k!, and fog-partition
A(n) log (1) + 1) (71 + 1) log( 1/2 - c D (m > 11 72 + ns 0}.
Crucially, integrating out À gives, for every D, the Negative-Biit Binomial marginal with parameters
nB
gn = ( k 1 (1-
r= 7/1 > 0, Equivalently, the marginalisation map written in natural coordinates is the smooth surjection x (1,; 72, 77) (r = m + p emh /(-
E (0, 1),
72
π R x (0, 1),
so = = ~ forms a quotient manifold with a one-dimensional fibre and the rank-one projector used in Algorithm T
~
(see Appendix A.2)
Open question. We do not know whether every curved exponential family can be "uncurved" by judiciously adding degrees of freedom; establishing necessary and sufficient conditions remains, to knowledge, open problem in exponential-famity theory. In particular, applying the same uncurving strategy to the Laplace family via a Normal- pormal-Exponential lift restores minimality in the joint but yields marginal family that is strictly richer than Laplace and therefore does not reproduce Laplace globally across the natural domain. We thus present our "add-a-fiee frree--hpeerprrameee trick as empirical recipe, a theorem; see Subsection A.3 for details.
A.2 Instantiation of the QBLR
Let (k, A) E 1,2,,,.... Rzo and define the minimal, regular exponential family
A
dn (k, A) l(k:; A) expf m log 122 + nak A(r)), h(kk A)
Ke! 10s00
(20)
Minimality is immediate: a; log À A + agk const on {(k, A)}. then varying & forces ds and varying forces 01 = a2 = O. Summing over &s and integrating in A gives the log part-titiiion
expfm m log A + ded + nokk dd = ( X expl (12 + e² ) ds
7(n)
k!
= =[(m + (- + e )²
(21)
which converges on the open domain
A, {πER? R3: mn > -1, 1/2 Th < <0}.
Hence
A(n) log I'mm + +1) (m + 1) logl 7/2 - en²).
(22)
Marginalization map.
For each fixed k; integrate out A:
(k)) 1/0 Oy (k, A) e A(r) elisk
k
ena d.
kl
l'(k + m 1)
(-112 to
11
e
(23)
T(m 11)
7/2
1/2
Writing
== +1 +1 > 0, =
E (0, 1),
1/2
21

---

(where p E (0, 1) follows from e < 0), we obtain
k + -K+++1 ( p)"
(k)
k: = 0,1,
(24)
Lee. the Binomial NB(r,P) p) for every oj E Anp Thus the marginalization map in natural coordinates is the smooth surjection
Ag = = R= x (0, 1),
T(/1 12, 18)) (r (r= 1/1 + p e -=-"3\(-72))
(25)
is visibly surjective: given any (r, p) E E, take mn = r - and, for arbitrary c > 0, set 7/2 = -c
and 17s = log(pe) then E A, and (n) (r. p).
Fibres and rank-one projector in natural coordinates.
The Jacobian of 25 is
0
0
DT(1)) -
e
7/2
-12
so ker Dr(n) = spann (n) with the vertical vector
k (n) (0, 12 1)² DT(1)) (n))))))))))))))))))))))))))
(26)
Hence each fibre is the smooth ID curve
Fn = { (m, e' 12, + + t tER
which leaves (r, p) invariant because r = 7n + and p = t (-1/)) Equipping the lift with its
Fisher metric F(T) A2)(m), the horizontal projector is rank-one:
Kip (1)) F(n)
P=(1) = (1) F(n) (T) (n).)
(27)
matching the quotient geometry used throughout (Theorems and
Figure illustrates the practical benefit of using the QBLR algorithm. Panel (a) shows the Fuclidean
gradient field -V (r.pp
- KT (gr.pl ((mme) computed via finite differences in the marginal coordinates
(r,p); these directions ignore the Fisher-Rao geometry and can yield poorly conditioned trajectories. Panel (b) displays the quotient natural gradient: at each (r,p) lift an arbitrary representative E π - (r. p), compute the natural gradient V (n))(() Tuse) in the 3-parameter Poimson-Gamma space, project horizontally, and push it forward through DTr(1). The resulting arrows respect the quotient manifold structure and are invariant to the choice of lift within each fibre, thereby guaranteeing parameterization-T optimization on the Degative-Binomine manifold itself.
A.3
When does the uncurving trick fails?
The Negative-Binomial example showed that introducing free scale parameter can uncurve Poisson-Gamma mixture and enable QBLR. Does the same strategy work for the Laplace distri- bution's Normal-Exponential representation? As we demonstrate below, adding variance-scaling parameter к > 0 does restore minimality, but the resulting marginal family is strictly richer than the standard ro-parameter Laplace: the uncurved lift introduces degrees of freedom that survive marginalization. This cautionary example illustrates that uncurving is an empirical recipe whose validity must be verified case by case, not universal construction.
Le
A= (21, Az) A Ae) As <0, d2 22 parameterize the (minimal) Normal-Exppoenntal variance-mixture lift with latent variance T E Reo:
exp
2r
X1 + 22= + AaT A(2)).
(x(7,) = /2
exp
>> 0.
(28)
Erratum:: In the rebuttal response we mistakenly the joint without the Gaussian base factor
(2xT) 1/2 exp -2² =2r()))) which makes the z-integral non non-oormalizable The corrected minimal EF is given
by Eqs. 28 29
22

---

(a)
10
(b)
14
54
Figure 2: Comparison of gradient flows on the Negative Binomial parameter space (r.p) toward the true distribution (red dot at 5.0,p 0.7). (a) Euclidean gradient field: arrows show the naive gradi- ent (rap) KL(gr.ellgau
computed via finite differences. This approach ignores the underlying statistical geometry and can lead to inefficient trajectories. (b) Quotient natural gradient field: arrows represent the horizontally projjccce natural gradient obtained by lifting to the Poisson Gamma representation, computing the natural gradient in the 3-parameter exponential family, and projecting onto the horizontal space. The quo- tient natural gradient respects the Fisher-Rao geometry of the marginal Negative Binomial manifold, yielding parameterizationinnarii
descent directions that follow geodesies fields converge to the same optimum, but the quotient approach provides more stable and geometry-awan up- dates.
information- -egaainnn---- sense. Both
The log partition function is
A(A) log (
8: -3 0, 7: - 53) 00
(29)
with the continuous boundary extension A(X) log 3 at y 0. Here K1 denotes the modified Bessel function the second kind
noa (10.32.10) defined by the integral representation
2
Ki
4 exp
>> 0.
(30)
41
Note that the form of A(X) trivially follows from the integral identity Moll 2015 7.2.12]. Equation
minimal, regular exponential family: the sufficient statistics (z/T, 1/T, r) are linearly independent, and A is an open subset of R33
(8
Why the marginal is not Laplace. To see why the uncurved lift fails, we compute the z-marginal by integrating from (28 Completing the square in the exponent gives 2² 2λ1
2
(z >11² 27, where (A2 + ( > O. The Bessel integral then yields
(=) x
exp - -23 v (=- \1 21)
(31)
21 27
This is not Laplace distribution. The standard two-parameter Laplace has the form
Lap(> b)
26 exp -
26 exp (x))/\)=m))/
involving only v z By contrast, (31) includes an additional constant 2y under the square root. This extra degree of freedom survives marginalization, producing three-parameter family
strictly richer than Laplace.
The uncurved Normal- PermallExpooential lift does define valid quotient manifold via Theorems 2 but the marginal family overshoots the target: we obtain a generalized symmetric distribution rather than Laplace. The standard Laplace embeds a constrained wwo-dimensional submanifold (requiring const) within this richer structure. Consequently, while the geometry correct, the probabilistic target wrong- iilit---ttttatnngggg that uncurving is case- aase-by-case recipe. not universal construction.
B Marginal quotient manifold theory
This appendix gathers the piiecrentiil-geommete material needed for Sections 2 of the main text.
We first recall just enough quotient-manifole theory to set notation (Subsection B. 1. derive the hor- izontal projector used in Algorithm] I and then give the proofs of Theorems 1-2 (Sections B.4
23

---

The exposition assumes familiarity with embedded submanifolds and the basic vo- cabulary of Riemannian geometry. Readers new this topic may find Boumal (2023 Chapter 3] concise primer before diving in.
Prerequisite.
B.1 Quotient manifold theory
A quotient manifold arises when we identify points in a manifold M according to equivalence relation ~. However, not every equivalence relation on M defines a quotient manifold. The con- ditions under which equivalence relation yields quotient manifold structure have been studied extensively in differential geometry [Absil et al. 2008 [Section 3,4 Quotient manifolds].
Formally, the quotient space / ~ consists of equivalence classes [a:] {u E M : ~ v). The canonical projection : .M + M/ ~ sends each point its class, π (x) [æ).. the pre-image of that class is defined as
The fibre through
F2 = IT ([=)) = {u E M :()) = π())} = {y E y
(32)
Throughout work with embedded submanifold M C R and projection M + N whose image C IR' is itself embedded. In this context, the general quotient manifold criterion reduces simple test:
is smooth and rank (T dim N Ye E M.
If condition (*) holds, then called a smooth submersion; every fibre - submanifold, and the quotient inherits a unique d-dimensional smooth structure. Consequently, we can define
1 ([x]]) is an embedded
dim(M/~) dim
secure in the knowledge that this integer is well-defined by the constant rank condition.
Under condition (+), each fibre π ([x]) is an embedded submanifold and there is a unique smooth structure on M/ that makes smooth submersion into M/ ~ [Absil et al., 2008 Prop. 3.4.2]. Moreover, the quotient M/ ~ is automatically Hausdorff and second second-counnnnn coountable.
Vertical space.
The tangent space of F3 the kernel of the differential
T.F, = ker Dx(I)) S T.M.
(33)
We call this subspace the vertical space and write Ve := ker Dx(x).
Horizontal space. Let (s, -)x be Riemannian metric on M. The orthogonal complement of V, is the horizontal spuce
H₂ = {ue TT : w). oY Vz).
(34)
key property of quotient manifolds is that Riemannian metric on M induces unique metric on M/~ if it is invariant along fibers. Specifically, if for any ~ and any horizontal vectors Hr and w E H, with [x) = Dw(y) we have (u, u). = (v, v) then we can define
well-posed metric on the quotient
(5. [p) (E. 2)x
(35)
where & and are the horizontal lifts of tangent vectors E,C E (M/~>). This makes M/ ~ a Riemannian quotient manifold.
Readers seeking a more concrete treatment of these abstract concepts may refer to Appendix B.5
where we examine them in the context of the Normal-Gamma distribution.
Throughhut, we treat the equivalence class [x: as a point the quotient manifold M/ ~e; the fibre is the full pre-image of that point. We write Ti / ~²) (with parentheses) only to emphasise that the tangent taken after the quotient, not the quotient of the tangent space at point æ. The shorter Tix or I M/ //~ can be used whenever confusion arises.
24

---

B.2 Orthogonal projection onto the horizontal space
At every point À A the tangent space splits as TxA H₂ Va, where V\ == ker (T is the vertical subspace (directions that leave the marginal unchanged) and H2 is its F(A)-orthogonal complement (horizontal directions that do change the marginal). For gradient based optimisation we need fast way to remove the vertical component of an arbitrary vector E TxA.
To do so, let K (A) RRdr be any matrix whose columns span VA (where dim VA = =). Write the desired horizontal part 9X = Ka for some coefficient vector E R.. Imposing F())- orthogonality to every vertical vector Ku gives the normal equations
KT FTA)() Ka) = = a [K] F())K] KT F())g. FFF))))
Substituting this a yields the explicit projector
Px(A)
K(A) (K(NTPAA)KON ((A) F(A),
(36)
so that 91 = g.
The matrix to be inverted is only with dim V\- In Normal-Wishart example 1; 36 then reduces to single scalar division, completely sidestepping the (dd) cost of inverting the full Fisher matrix.
B.3 Proof of Theorem 1
This section is devoted to the proof of Theorem T Before proceeding with the proof, we recall the notation established in the main text.
Theorem (Induced Fisher-Rao metric), Assume the setting of Theorem and equip the natural- parameter space A with its Fisher information metric Fx. Then:
(i) The map n, that project the Riemannian manifold (A, Fx) on E, induces Riemannian quotient manifold structure on E:
(ii) The Riemannian quotient metric on then the Fisher metric of
Setting and notation. In this paragraph, we restate the symbols used in the main text, all in one place. We work on a measurable product space Zext = Zu x Zv and write žext = (zu, zv) for generic element. The block æy collects the coordinates whose distribution we ultimately care about, whereas zv will be integrated out. Therefore, to shorten our notation, we refer to Zy and zy as Z and 2, respectively.
The marginal family that defines distribution over is parametrized by E. Its ambient "parent" is minimal, regular exponential family with open natural-paramette parameter space A that defines a distribution over Zext-
The key connection between E and a marginal relation
qe (z) = f gx q2(2, 2v) dzy.
(37)
This relation naturally defines function A E.. And the function naturally defines the corresponding equivalence relation ~ep on in the following way:
A A2 x x(A)) = TT22).
(38)
That is, two points in are equivalent they yield the same marginal distribution.
Note that, for minimal regular exponential family, the log- aog-parttiion function A(X) is infinitely dif- ferentiable, and its derivatives correspond the moments of the sufficient statistics. The marginal- ization can be expressed in terms of these moments, which inherit the smoothness properties of A()).
We remind the reader that in our setting to prove Theorem it suffices to show that is smooth submersion (see the condition (+) For convenience, we quote the theorem we are about to prove in the notation fixed above.
25

---

Theorem (Marginalization yields a smooth quotient manifold), Let be a minimal, regular exponential family with parameter space A C. R Suppose a partition 2 = (Z, Zv) is chosen so that the marginal family {ge)ee tee obtained via π A A is moment-parametrried (with dim E = k) (Definition)) 77. Then E is the quotient manifold of induced by n.
Proof. Recall the marginal relation
/ Zv
25 (z)
as (a(x)2y)day ay Zu--
Step O. Set-up and notation. Write T (Ti,..., .Ta) and (X², (A-,...,AA 14). For any bounded
measurable ş Zu + R set
(0) gE) = / Zv
(lzz) (ZU) dzu.
Step 1. A finie-ddimensinnal probe of the marginals. Because the marginal family moment- parameterized with dim E = k, it is attached to k integrable functions m....... m, E 1 (2y):
A - R²
ea()) == (m; (a())) (i = 1.....,)) kk).
Writing the marginal single integral gives the equivalent form
{ mz(ze (gx(Zex)) (Zea)) dzeat-
(A)
Zen
The goal is to show that is a smooth submersion constant rank r = rank (x(A) (the same r for every X). Once established, each fibre -'y) is automatically embedded submanifold of A.
Step 2. Computing the derivative of e. Fix À E and a tangent vector v = (c',...,,)) E TyA. Differentiate under the integral (dominated convergence allows this):
D,et(A)=
pilzu (Zeex) dz
ON
Z
j=1
d
vd
pi(zu)
Oxi Qu (Zea) dzext.
#=1
Zoa
Because 8X3
9x3 (x)) (Tz
OA(A)
a2 (z), we obtain the exact Jacobian entry
80
det
Jiy(A) = Ox (A)
(39)
(z((exs) - 0A(A)
mi (z(
9x (Zexa) dzex
(40)
Zea
B.g (Ty)
Note, that from (39 smoothness trivially follows from the fact that A(A) E (A)
A convenient way to rewrite equation 40 is with covariances (up to additive constant which does not change the rank):
Jig (A) = Cov,, ( (Zea))
(41)
So each column j of J(A) stores the covariances between the function my and the statistic T under the joint distribution
Step 3. Why the rank is constant. Minimality of our exponential family guarantees that the d d covariance matrix F(A)= == Covg, (T.(Zea)) Ty (Zext) is positive definite for every À [Brown 1986
Theorem 4. Denote by
H()) = J(2) J(X) E R
26

---

Because Fx > O, H()) is positive semidefinite for every X and satisfies
det (A)
rank J() k.
Suppose that det H (Ap) > at some point Ao Then the matrix H()a) = F(
is positive definite, and in particular, the Jacobian matrix J(Au) has full rank k. This means the map À + e(X) has full rank at Xo- Since both F()) and J()) are real-analytic functions of X, the composition (A)
((A)) is also real-analytic. Consequently, det (A) is real-analytic scalar function on the parameter space. By basic property of real real-analytic functions, if det H (A) not identically zero, then its zero set has empty interior. Since det H (Ao) > O, the function cannot be identically zero, and hence there exists an open neighborhood of X where det H (2) > O. Thus, rank (J()) = in a neighborhood of X Therefore, the rank of J()) cannot drop in any open neighborhood where det (A) is positive. If rank were to drop at some point, this would force det (A)
at that point, contradicting the real- callannayyt and strict positivity nearby. Hence, the rank remains full wherever it is full once.
Step 4. Submersion embedded fibres. We now consider the smooth map e A + R which we have shown to have constant rank k. By the imite-dimensional constant-rank theorem [Lee] [Theorem 5.12]. it follows that each fibre (y) is embedded submanifold of A of codi- mension & (i.e.. of dimension d - k) and these fibres vary smoothly with y. forming a regular foliation of A. Moreover, since e(A) depends on À only through the marginal distribution qg. we have:
(e))) = (X' E A (m(x)) = (m())) = π = ((m()),
where ग A Q= denotes the map sending À to its marginal distribution (=(A)- Therefore, each marginal pre-image an embedded submanifold of A, and the space of parameters decomposes smoothly according to level sets of the marginal.
B.4 Proof of Theorem 2
This section is devoted the proof of Theorem 2 statement of the theorem from the main text. We use the notation established in the previous section.
Before proceeding with the proof, recall the
Theorem (Induced Fisher-Rao metric). Assume the setting of Theorem parameter space A with its Fisher information metric Fx. Then:
1 and equip the natural-
(i) The map n, that project the Riemannian manifold Fx) on E, induces Riemannian quotient manifold structure on E:
(ii) The Riemannian quotient metric on then the Fisher metric of
Proof. Consider Ç and À ग. (E), then let f(z, denote the log-density of the distribution (z))). and let f(zext A) be the log-density of the distribution q (Zex).
Because π A + E is a smooth submersion of constant rank (proved in Theorem D the Local Section Theorem [Lee 2012) Theorem 4.26) guarantees that for every E E E and every À E म -1() there exists an open neighbourhood U C E of and a smooth map a U -> A such that mo idg and (E) X. Patching these local sections with a smooth partition of unity yields smooth global idg.
section = + A satisfying
The log-density of the distribution ge (z) can be related to f X) in the following way:
f Zy
fz. { log
exp(F(zen, dzv.
(42)
Consider the following helpful function
((z, A) : /
exp fz, dzv.
(43)
Zr
27

---

Then we can express f(z, {) in the following way:
f(z, E) = log (()).
(44)
Now, we differentiate both sides of the identity 44 with respect to & and we get the following:
( / exp (i)
0g.5(a,5)
(Leibniz + dominated convergence)
(<(2) 2v
(5(z) / exp (1)
Σ
O2 f(z, zv, (E)
(chain rule)
Bv
==1
Ezvl=
Σ Fx zv, A)
(recognize conditional density).
The last identity can be re-written in vector form in the following way:
vef(z,E) = Do()) Bz.j2 Vxf(z, (
(46)
Introduce the conditional joint score
s(z, Zv.
A) ZV.
2v,
g(2, = Egul z).
With this notation (46 reads
V vef(z,E) =Do(EE ozz
Taking the outer product and integrating over z ~ 9E
FE(() = E. = Do(()) Eggg Do(E), A= oE))
(47)
MM((A)
Write F(A) E[ss"] for the Fisher matrix in A and C(A) E& [Varls =]] for the average conditional covariance. Then using the total variance decomposition, we obtain the following:
F(A) C(A) + M(A) = M(A) C(A).
(48)
Let
R(A) [AAT J(A) = DT(
(49)
Because every residual score 81 the vertical space ker J: consequently
E(s z] satisfies JsL = 0. the matrix C(A) acts entirely in
RJ C(A)R = 0.
(50)
Inserting the facts 48
50
into the pullback formula (7). obtain the following:
FE(() = R(X) ((A)R(M) J(A) = Fs(E).
(51)
Hence the Fisher information of the marginal family {Eg] } is obtained from the full Fisher on A simply by pushing it forward- eeuiialenlli pulling it back tthrouhh the Jacobian J Dz
This metric compatibility (5I) fulfills exactly the hypotheses of the Riemamian-quotient theorem, so
all conditions of [Boumal 2023 Theorem 9.35] satisfied: A/~ inherits the unique Riemannian metric that turns π into a Riemannian submersion, that is compatible with relation (51).
28

---

B.5 Univariate Student-t as a quotient manifold of Normal-Gamma
This section is dedicated to a concrete example Riemannian quotient manifold theory applied to marginalization. Even if the formal derivation of the mathematical objects introduced in Sec- tions [.I and B.2 is not required to implement our main result: Algorithm I it offers intuition on how quotient manifold and QBLR work.
Narmal-Gamma distribution a joint distribution Iiu,o-1,00 .0,0) 3) (x, 7) over a random variable (2; r) defined by the following relationship:
x/T < NH, (o ()
(52a)
Gamma (3)-
(52b)
It is straightforward to rewrite Q11ar11 .) (I,, 7) into the minimal exponential family representation
qx(2, T) 22x
exp (T X A(2)),
(53)
where the natural parameters, sufficient statistics, and the logpartition are:
o 1'1p²
1
x = 0
2
2 -
EA-RxR.
x
R4
2
x R_,
(54a)
T(x,T) = (IT, 2T T,log T,
(54b)
1 A(A) = logI' Aa + 2)
1
log( (2x2)
(28 + log - +
A 412
(54c)
The marginalization over - defines mapping from the Normal-Gamma parameter space À (A1, A2, λ3, X4) to the Student-t parameter space & = (u,," via the following marginalization (quotient) map:
2222 2As + 1/2) 2Ax 27s+1) +1)
A1 - + \2222)2)
T(A)
(55)
Our construction starts with V> ker D (A) the vertical space; to obtain it, we need to compute the differential of our quotient map. Then we are left compute Hy = (V
Differential the quotient map.
The Jacobian J(A) = DT(A) E 1x344 is
A
0
0
22a
2
A(A)
A
442d (A3(2+1122)
A--422A4 2Ax
1
(56)
113(0++1/2)
2>a(As+1/2)
0
2
0
Vertical space. With the Jacobian of the quotient map in natural coordinates (Eq. (56)) ), the vertical space À is simply its kernel:
V> ker DT(A) span k(A)) k(A): (4), 4x2, 0, X 4
A direct row-by-row multiplication shows DT(A) k(A) 0. k(A) lies in the kernel. Since (T(A) has full row rank 3 for every À E A, the kernel is one dimensional and dim Vx = At this stage, no Riemannian metric is needed the vertical space is determined purely by the quotient map π..
Horizontal space. To define the horizontal space, we must specify Riemannian metric, as differ- metries generally yield different horizontal spaces. In our case, we employ the Fisher-Rao metric, which for regular minimal exponential families equals the Hessian of the logpartition function. For Normal-Gamma distribution specifically, the Fisher information matrix takes the following form
the Normal- Gamma
ZAga
b A7h
4
à 2
+
A 8A
2Aya
A
DAga
A1b
b,b
X/b 16A
A7
X²
xab
223a
2X
2X
42z0
F(A)
(57)
A2
As)
22za
4
2
A
X
1
A+
4
29

---

where a Xi - A3. and e is the trigamma function. Equipped with the Fisher in- IAp
formation matrix in natural coordinates (57). the horizontal space is defined as the F(A)-orthogonal (:)): a tangent vector (V1, U2, V3, v2)" is hori-
complement of the vertical line Vi- span (
zontal iff
(v. k())rr R(A) F(A) v =
The oothogonayyy condition (A))F(() v = 0 is equivalent to requiring v to be orthogonal (in the Euclidean sense) the single vector
n(A): == F(A)
short calculation with the entries of (()) in 57 gives
2
11 4)2d4
Th a²
Abbda
4bd224
n2 = 2
n4
Aza?
a²
Provided 0a 0 (true the admissible domain A), we have T4 the linear constraint n u = 0 can be solved explicitly:
711
112
13 n
14
21
1/2
23
72
n4
Choosing v1, 12, vy successively as the standard basis vectors produces an F-orthogonal basis the horizontal space:
(A) (1, 0, -m//na,
(A) = (0, 1, 0, -22/
(A) = (0, 0, 1, -na/ -ns/na). na).
Any natural gradient g now be decomposed as 9 911 + 91 with
k(A)
9n
(g.n)
k(A) x()) and 91 9 9/
so that Hy is the direction used in Algorithm
C
C.1 Definition and properties
random variable Sz follows a multivariate Normal- Wishart distribution with parameters
(j, V k. if
S ~ N(u, (RS) S W(v, V)
(58)
(59)
where μ IRd the location parameter, scaling parameter, and v > d - is the degree of freedom parameter.
dxdd is positive definite scale matrix, 0 is
The joint probability density function of the Normal-Wishart distribution is given by
2lz, Slp, V =
(60)
C.2 Marginalization and the Multivariate Student-t
Wishart distribution S Walp, 1 S NM, (x:S)-1) marginalizes multivari- ate Student-t (see MMurphy| 2007 Section 9])
p(z μ. V, K,
IA,, To(z (,
(61)
where Σ K(P- I+ v/ = v d + 1. Hence the mapping
(n, K, - Tdl, V - /[x[v d + 1)) -dd
many-to-one: a fixed (j,, v) and fixed product T determines a unique Student-l. So by changing (x, V) while keeping their product fixed, yield the same Student-t.
30

---

C.2.1
Canonical exponential family form
The Wishar- distribution can be written in exponential family form in the following way:
p(2, exp (ATT(2,S) -A(A)).
(62)
where the sufficient statistics are
Sz S
S)z
(63)
2 Sz log det S
and the natural parameters are defined trough the standard parameters (H, V, (2, 1/) in the following way:
A1
1 1 -1 + upT))
A2
A
X
=
(64)
y 2
A
Then the log-pau pog-partiioon function is
d
d d
A(A) = log( (222222)))
log det
(65)
2
+ d(d + 224)
+ 2\4
d
log(2) + logl'a
+ log(2x).
2
2
2
As established equation
7 the mean parameters are given by the gradient of the bog-partitinn
function
A(A)).
(66)
For the Wishalt distribution, these parameters are
01 = dA(x)
E [Sz) = MW
dA(A)
02 dA2
E [S] =VV,
(67)
dA(A)
03 dAg
[== Sz] un V K
04 dA(X)
[log det S] log det vp log 2 t/4 (
dA
where via the multivariate digamma function.
Proof. 02 is computed using conditional expectation,
Ep.S) [Sv] = Es [ST [u]] = Es [Sn] Nap.
(68)
The value of O2 is directly the moments of the Wishart distribution Faton 2007 Proposition 8.3] and 0g can be derived from them as follows:
31

---

E
S2
2
E ( [
(69a)
1.j
>
Es (Sh. Euis [=]=]]]
(69b)
Es [S.. (Cov((C))Zz)) + Ras [2s] Buis [z]]]]
(69c)
Σ
Es [5 (KS)) + [Wmas]]
(69d)
Es [s,s()) + Σ Es [S
(69c)
Es Lrr (S(S) + Σ Lللlill as SeS
(69f)
d + vp wp..
(69g)
n
04 is direcly the log-expectation of Wishart distribution given by Penny (2001].
C.3 Derivation of the NGD update
Let's consider (SS) Wa(S)u, v) and q(2)8) = Nizla, (x:S)-
We denote the log-likelihood for the n' ''th data point by fn fn(z) == log p(Dn with a Normal- Wishart prior with parameters μ = 0, v I, 1, and the degree of freedom parameter Vo. We use the lower bound defined in the joint distribution, p(D, z. S)
((A) = 1y((,99 [log pDD, 2, S) log ((D, S.]
N
log 2)) + log N (Z](0a, S-1) N (2/u, (x:))-1) + log
Wa(S Po, La) Wa(SIu, VV)
Σ
(2)
Our goal is to compute the gradient of this ELBO with respect to the expectation parameters O (de- fined in (67). Because is an invertible re- parameterization of the standard parameters (J4, V Ki 2), their gradients are related by the chain rule as follows:
2 vo,C= d
VnC,
(70a)
Vo,L -v -1-1 V F 24 CosC,
(70b)
VoL VV v L) - vVLL
(70c)
bu (///(/22)
VorL - 1 - V v v< D)pT v(Vo, C)p -Ve,)-(V
(70d)
note that by 3i2 we denote the multivariate trigamma function.
32

---

The ELBO gradients for the standard parameterization can be obtained as follows
Val(A) - Σ
V. g(x [fn(z)] - vd
(7la)
N
V.C(A) = Σ
V g( [f((z]]
dR
(71b)
2
1
Vel(A) Σ
g( [fn(2]]
2
μ
+
2
(7lc)
i=1
d
1
V-
V (( [fn(z]] V&C(A) Σ
2
Up
v/ 2
(71d)
2
4
where di the multivariate trigamma function.
The last thing to instantiate Algorithm [ for the Normal-Wishart is to implement the projection onto the horizontal space (see Appendix B.2). For the Wishart lift, the vertical space is one-dimensional, so the vertical subspace at any X is VA span(k(A)) with
k(A) = vec(A)),) 0)
the last natural coordinate \4 always effect the marginal. Given the natural gradient g V,L, its orshorool projection is obtained by removing the component along k(A)
b(A)TF(A) k(2). k()) F(A)k(A)
9X = 9
(72)
Because k(A) is single vector, the denominator is scalar; evaluating (72) therefore requires only one call to the Fisher-matrix-vetto product and one scalar division; no inversion of the full Fisher matrix is ever needed.
Using the derivations this section, we can now summarize our algorithm. Specializing the generic quotient-natura--nnnde loop (Algorithm 1 to the Norhal-Wishart lift (ji, V, K, v) gives a fully explicit routine:
1. computes the stochastic data-fit gradients in the standard parameter space (gdati 99th "o .g .g""");
2. adds the analytic prior terms 7la 7ld):
3. converts the result to the expectation coordinates (go, go.) via the chain rule
(70a) 70d);
4. removes the vertical component with the rank-one projector (72
5. performs natural-gradient ascent step of size 3r in the horizontal direction and back- transforms to (p, V, K. v)).
The whole procedure, including the projection 72 is collected in Algorithm below.
33

---

Algorithm One step of the quotient natural-gradient update for Normal-Wishart parameters
Input: current standard parameters (1,, Mr Ki v), minibatch Bi, dataset size N, step size B, Data-fit contribution
N
data
Ix
data
- Iv
data
data
V (mn,. (o)
nEB.
Add prior terms
(Eqs. 71a (7d
94 -- 4p
data
UB
d
9x. 9k.
2
K
14l1 )
gv 9
2
d
1
-1
Iv Iv
data
2
Up
Vid v/22
2
Chain rule
(Eqs. 70al (70d)
(go,, 902: 90 go.)
ChainRule(gp In. 9v
Horizontal projection (rank-ond one)
(Eq. 72
A(A)T F( k(A)T F(A)
9e t ge a((A
b Natural-gradden update in X Sppcce
d + A Bc Be
i= 1:4
b Back-transform to standard parameters
(Eq. (64)
K 22x,, f /IRR, V v-t -222 -xppp vt 2X + d.
C.4 Path-gradients for Normal-Wishart
In Section D we implement Algorithm for Student-t distribution through the Normal-Wishart marginal representation. For a concise implementation of the algorithm, refer to Appendix C.3 This implementation requires unbiased gradient estimators v μC, ÖxL; & C. For the Normal -isshar variational family, these estimators can be obtained from the general gradient form provided in Theorem for function f : Rd + IR.
In the following statements, we will use the so-called Lyapunov operator
TY : S + gd TY AY + YA.
(73)
We denote with - the inverse Lyapunov operator, defined by:
T [B]
(74a)
with AY + Y B, A A BES
(74b)
According to Bartels and Stewart 1972 A 0 is sufficient condition for T-1 to be correctly defined. We will also refer to the operator Sym that associates a matrix to the sum of its transpose and itself follow:
Sym A E Tec
A+AT, V'.
(75)
Theorem 4 (Gradient Identities for the Normal-Wishart Distribution). For a dimension d 2 1 and gd ++ v > d + 1, consider the joint density of the Normal Kormal-Wishart
parameters P E Rd x v
distribution
34

---

(u.,v,(2 =M (x.))) Wi(S v, v
pp,s(*)
War
Let f: RRd R. be a twice-differentiable function that is integrable with respect to qu. vzz and whose first and second derivatives also integrable. The ensuing gradient identities are valid:
1. Gradient with respect μ (Bonnet identity):
Vp Eg [f()) = Egg [V,((z)
2. Gradient with respect to k (Price identity):
Eng [[(=)) = 2x² Ey [1((s-
ŌR
3. Gradient with respect vp (Price identity):
V4 [x)] = 2x
T-1 (S
where B W(v,T) TT and z N (1s, (N (xb-1/236-1122)-)
The gradient of any real function in the mean parametrization (including ELBO) can be straightfor- wardly deduced from the equations of Theorem Detailed proofs for each identity are provided in Lemmas and 4 respectively.
Lemma (Bonnet identity for the Normal -Wihart lift). Under the conditions of the theorem 4 the following identity holds
V [x[((=)] = [=[(((())]
(76)
Proof.
V By [f(2)] Eur,
V T. s(=) [f()]]] (dominated convergence + Fubini)
(77)
= ( (B),,(=) V [=)]] (by Lin et al. 2025 Theorem IJ)
(78)
= [o))].
(79)
The proof above employs the vanishing surface term, mirroring the classical Bonnet proof (in
French) [Bonnet 1964). more contemporary explanation of the same finding is provided Lin
et al. 2025 Theorem I]).
Lemma (R- Price identity for the Normal-Wishee lift). Under the conditions the theorem 4 the
following identity holds
[,]] = 2x² Rg (irS--V2s(0))
(80)
Proof. Let (2))(2) = N2 n, Σ) with E = (x:)) then
OfE.
Eur.
O Bo [f]
(Fubini's theorem)
Eur..
(Or Vs... [/])]
(chain rule)
Eus.
x (-x Een.s!
(by [Lin ct al. 2025 Theorem 4])
2x Eyg [(((S-)\(()))
The first line exchanges the differentiation operator and integration operator, which is possible be- cause the derivative of Up. can be bounded from above by integrable function. The second applies the chain rule. The third uses two facts: (1) Σ x -5- implies O, -5-1 and
(2) the classical Price formula [Lin et al. 2025 Theorem 4] Ez = [ The final line simplifies using the trace inner product, the linearity of the trace, and the expectation.
D
35

---

h($) Φ 1-1/2
84
8
+X
[
84
1//2 1- 1-1211
1/22x
1/21
Figure Commutativity diagram Lemma The following commutative diagram illustrates the Fréchet differentiability of the inverse square-root map h + + 4 --12 on the space of symmetric positive-definte matrices 44... The vertical arrows represent perturbations in the input space and the corresponding linearized response the output space symmetric perturbation X E 3d to the input & corresponds, under the linearization of h. to symmetric output given by the Lyapunov operator. The boltom the diagram means that effect of first perturbing and then applying h, versus applying then diffcrentiaaing, yields the same result to first order in X.
derivative Dh( This diagram expresses the fact that applying small
represents this linear transformation. Commutativity of
Lemma 3 (Fréchet differential of the inverse square-root). Let & $4 The map + -- is Fréchet differentiable with:
Dh( X - f
(81)
22
Figure illustrates the commutative structure of the differential relationship provided in the Lemma
Proof. The Fréchet differentiability of the square root in 8 is a direct implication of Moral and
Niclas 2018 Theorem 1.1]. In
1 Id,
(82)
we substitute the functions with their respective Taylor expansion at point & in direction X E $d. With that substitution, we get the following:
- + Dh( + + 4 + X) ( + D + = Id
(83a)
+ क - 1 + - -
Db( + x(I ld
(83b)
+ - 4 D + D) x[Ixx] = 0.
(83c)
Given that > 0, equation 83ch implies that we can define [X] as T--1- t X $ which is precisely the statement of the lemma.
$2
Lemma 4 (Φ Prririce identity, Lyapunov version). Under the conditions of the theorem 4 the follow- ing identity holds
1
Va Eg[s(=]
E.
T #b
(84)
2x
where B ~ W(v,)), z/B ~ ~Nu,, (nd
1/2 112 B 1 and Sym
1 := AT).
36

---

Proof. We compute the gradient of + + 1xtth by treating it as the composition of two functions: first, + (i0b-1/223 1/26)1/1BB11/3)- and second, a E H 1/
(x))
Va = (V
(85)
where (D)( represents the adjoint operator of ()
Under the conditions on f from Theorem 1 we can apply Lin et al. 2025 Theorem 4] to obtain the following: Ve(D) = ===((ll [=)]
We express & as the composition of three functions:
1 E S+ + $-
with differential Dé (+): : X e ES-4 - T Φ x #]
(86a)
02 S++ + ABA with differential Dø A) X S ABX XBA,
(86b)
4 S E S++ + (x:5)-² with differential (S) X E 8 -
(86c)
We recall that any Riemannian metric on the manifold S can be expressed as (X Y) ur(# ( X$-1y) where S is isomorphic to the tangent space of S + and X, Y E E (see Ohara et al. 1996]) ) Based on the form of the Riemannian metric on the tangent space of + we can state that for any A.A E S++ the differentials of (A2 and Dég (A) self-adjoint. According to Tippettt et al. 2000]. T-1 is also self-adjoint, making DΦj (4) self-adjoint for any & E S++. The differential DØ(
S4+ then self-adjoint as the composition of self-adjoint operators. The gradient of a 6 can be expressed using the formula (85) as follows:
Vx(ao = (Dés o o2 o
(87a)
(Døy(T= B b) ( ( D (+) 'IV [ooaaa
(87b)
o Diz2 D61(0)
))* Døs( BΦ V
(87c)
= -
-T-
Sym( - BA V )
Vo -s #) +)
(87d)
-
1T-1
BBp
(87c)
45
x-1
- [+]( ot B-1)
T
(87f)
We can apply our formula (77) directly under the expectation over B ~ W(Id, v) and under the linear operator T-1 to obtain final gradient follows:
VgE.D [f()]] = Eg
V E [=)]]
(88a)
= EEg
V + a (a o Φ))]
(88b)
EB
K
T
Sym( - - -
(88c)
2
42
E
T-11
(88d)
2
47
D Experimental setup and reproducibility protocol
Benchmarks.
Pre-processing and splits. (1) 80/20 stratified train-test split with random state=42; (2)
feature- wise standardisation using training means/variances only.
37

---

Models and inference schemes. All tasks use Bayesian logistic regression (BLR), We compare
three -infreenceeeeeeeeeeeee schemes:
Abbrev.
Variational family
Optimiser
BBVI"
Student-r
Black-box VI Roeder al. 2011
NG-LIN
Student-r
natural-gradient Lin al. 2020a
NG-Ours
Normal- Wishart (lift)
Quotient Natural Gradient (Alg.
Optimisation schedules (parameter-free). To eliminate hand-tuned learning rates, we use the Bestance-ove-Graulle (DoG) rule lfgi al. (2023] and its Riemannian generalisation (RDoG) [Dodd et al. 2024]. Both schedules set the step size adaptively from quantities the algorithm can
measure on-the-fly.
Euclidean DoG (for BBVI*). Let wt be the parameters and gr. the Euclidean gradient. Maintain
Ft max e, max
xoll22 soll
G.
>
1
a
and set
Ft
T1+1
Th. 9t-
GG
We use = 10-3
Riemannian DoG (for NG-LIN and NG-Ours). Replace Euclidean norms by natural-gradient norms and the Euclidean distance by the geodesic distance d!(,.) associated with the Rasher-Rao metrie:
Fz max 6, max d(I,, 20)), G. = Σ ist
Ft
est
G.
2++1 (-
We set = 10-3 and, unless noted, use the non-positive curvature correction Sx (i.c. A 0). For NG-LIN, d(.,,) -) is approximated by the symmetric KL between two Student-t distributions, estimated with a fixed set of 64 Monte-Carlo draws anchored at the start point to reduce variance. For NG-Ours, d(::,) .) is the exact KL in the lifted minimal exponential family (Normal-Wishart) available in closed form.
Safety of the lift-based distance. Let A - E be the marginalisation map from the lift to the marginal parameters. The quotient-metric result (Theorem 2 implies that, for X with =A.).
KI(gs, 4x, |l 2xa) 2 KI(9s, (l.)..
Thus the lifted KL we plug into RDoG an upper bound on the (unknown) marginal KL. Be- cause DoG/RDoG chooses ?/, Ft/ S/. (Fr) G,, larger distance yields (mildly optimistic) larger step. A tighter, future alternative is the abre-minimised lift distance, infx,e=-1(E) KL l14x,) KL (l i.c. the true quotient metric.
(shared).
Epochs = 8,000; mini-batch size = 32.
Step sizes: hand-tuned learning rale. DoG/RDoG schedules determine we with 10-3; curvature correction disabled by default (x 0).
Monte-Carlo samples: 10 per update for BBVI*; for NG variants (gradients), plus 64
fixed draws for the NG-LIN symmetric-KI. distance used by RDoG.
Software environment. Python 3.11 (CPU-only): jax 0.6.0. numpy 2.2.4. scikit-learn 1.6.1, torch 2.6.0, pandas 2.2.3. A version-pinned pyproject.. toml is included in the repository.
Hardware and runtime. All experiments run on a single CPU-only machine (no GPU/TPU). End-to-end wall-clock time to regenerate Table is 2 hours.
38

---

Reporting,
For every (method, dataset) pair we report: (1) posterior-men mean accuracy (acc,),, (2) its
standard error of the mean (SEM). Results are produced by python run_ comparissn py.
Reproducibility assets.
Code https: wtt.ccoo//iiaa http:///tttub.cmmm//aaaaaaa////a One command reproduces all num-
bers and figures.
Determinism. NumPy, JAX and scikit-leam PRNGs fixed to 42; JAX in deterministic
mode.
Environment capture. pyproject. toml and a generated requirements-loch txt
freeze packages; a Markdown "compute card" records CPU model, cores, OS, and energy draw.
E Complexity Analysis
In Algorithm [ two operations have a non-trivial computational complexity: the natural gradient computation and its projection onto the horizontal space. Based the current literature, we analyze these complexities. In Appendix E..I we explain why projecting the natural gradient is negligible compared its estimation. In Appendix E.2 in the context of the Normal-Wishen example, propose a methodology to reduce the complexity of the natural gradient estimation.
We thank anonymous reviewers for raising this question.
E.1 Complexity Analysis
Projection computational cost. The projection operator (>) defined in Equation (6) is used to compute the horizontal component of the natural gradient. For clarity, we denote the projection of the natural gradient go IddimA follows:
9и = PH (A)
naive approach to compute gu includes the inversion of the symmetric positive definite matrix
RV(A) == K(>) F(A) ((A) RRd.xd.
where K(A) is the matrix representation of basis of the vertical space ker Dr A
33
(with
- is the marginalization map (0)). ), F(X) the Fisher information matrix (3 d
dim Va
(with Vx the vertical space). This approach would require computing the matrix inverse, with cost
Odd)) However, we avoid materializing the inverse. Instead, we solve the linear system
RV(A)u K(A) F(X) go
(89)
using a small dense linear solver: due to the fact that RV(A) is symmetric positive definite, the conjugate gradient solver by Hestenes al. 1952 is applicable. This reduces the computational cost to (d
In practice, well-constructed lifting ensures that d, < dim A, making the cost of projection negligible relative to the Fisher-vector product F(X) go. which has cost Oddim A²).
In our main Normal-Wishare case, the vertical space is one- me-dimensional (dg 1),, so the projection reduces a single scalar division.
E.2 Future Improvement
Quadratic time is acceptable up to a few hundred dimensions on commodity hardware, but larger problems call for additional structure. In the following, we propose plan reduce computational complexity.
Let us denote by d the dimension of sample. We propose two different factorization methods to reduce the O(d)) computational cost of the natural gradient while preserving its geometric property.
Structured covariances
Restrict the scale B blocks diag( II...., V B) with sizes de-
Sampling: E. O(d)); purely diagonal v needs only O(d) Gamma draws. Fisher products Hessian trace both factor block-wise, yielding the same E, Odf)) and O(d) in the diagonal case.
39

---

Low-rank with diagonal factorization Decompose the scale matrix as
V = LL" + diag(u),
L E Rd k < d,
so that only kd + free parameters are stored instead of jd(d + 1).
Sampling. Each column of L is drawn from matrix- -ormall and the diagonal entries of from independent Gammas. The two draws require kd and d random numbers, respectively, hence (kd) time and memory.
Fisher-vector products. In the horizontal projector need (diag(u) +LL ) x. Compute as
y - diag(u) + x).
O(d)
(kkd)
Both multiplies with L cost O(kd), so the total is O(kd) per Fisher product-linearr in for any fixed rank k.
Hessian trace (two Price identities). Both the x-Price and the 4-Price terms require tr(S (2f). Rather than forming the dense Hessian, we use the Hutchinson estimator Hutchinson 1989)
tr(S-1v2f) = Σ " (((S-1 ~ [+1]
Each term needs Hessian product (HVP) and one multiplication with S-1. The HVP is model-specific; the 5-1- multiply uses the Woodbury identity:
S-ix Dー D- 4(4 + L -'
D diag(v),
which is again (kd)) Choosing probes keeps the overall trace cost bounded by (l22)).
Both variants leave the vertical space one-dimensional, the horizontal projection stays single scalar divide.
Summary
Dense QBLR is Odd)); with a diagonal or block-diagonal scale drops to (dd), and with rank-k plus diagonal it is Okk²d) (linear in d for fixed k). These paths scale QBLR to far larger latent spaces without sacrificing its geometry or requiring matrix inversions.
40


Reviews:
[{"review_id": "PwywAXJOrg", "reviewer": "Reviewer_Jttp", "summary": "The paper proposes  the Quotient Bayesian Learning Rule, especially on the heavy-tailed distribution models. Introducing some definitions, the paper proves that using some projections, the marginals inherit a unique Fisher–Rao information geometry via the quotient-manifold construction, and hence the optimization methods, for example, SGD can be used.", "questions": "Please see Strengths And Weaknesses.", "limitations": "N/A.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "1. More discussions and examples on Def 1 and Thm 1 would be much helpful.\n2. What is your benefit to use your new method over the original Bayesian learning rule in Khan and Rue 2023.\n3. What is measurable moment functions in Def 1? Are they the same as $e$ in your supplement material? More discussion on the existence and uniqueness of $e_i$ will be helpful. \n4. Still for the Def 1. Is (iii) the bijection too strong? But it seems that the target of using manifold and Def 1 is only for the local smoothness of the parameter space. That is strong (although I admit that is often ignored by many papers). So essentially, to get a better SGD, we need a better (ie regular and smooth) parameter space. I may miss something important. So you just compute the gradient, find the dual parameter and do horizontal projection, and find the optimizer and corresponding optimizer for $\\xi^{*}$. What is the benefit? See quite normal ways to do variational inference. \n5. I don't understand why you mention Lie-group in the paper. \n6. The paper only uses logistic regression models which have quite analytical expressions. This is a little bit simple. Have you tried other models with more complicated structures and multi-dim outputs?  Have you scale the questions up?\n7. For the applications, the computational complexity seems high. In real problem, finding projection and inverse of matrix leading to higher computation cost and less accuracy. \n8. Have you tried to find the reason for the divergence of ELBO in Figure 1 without horizontal space projection? Have you tried other experimental setting?\n\nMinor: some typos like in Algorithm 1.", "was_revised": true, "final_justification": "Thank the authors for the response. I have read the rebuttal and will keep my score.", "rebuttal": "We thank **Reviewer Jttp** for the thoughtful and constructive feedback. Below we address every point in turn.\n\n## Question 1\n\n\"More discussion and examples for Def. 1 and Thm. 1.\"\n\n## Answer 1\n\nWe will insert a self-contained example in the final text—the **Laplace** distribution viewed as a Normal–Exponential scale mixture—in Section 3 and a new figure for it. The fibres of the quotient map appear as 1-D curves, and the horizontal natural gradient is visualised as moving orthogonally (under the Fisher–Rao metric) to those fibres. This concrete picture should make Definition 1 (moment-parametrised family) and Theorem 1 (quotient manifold) intuitive.\n\n## Question 2\n\n\"What is the benefit over the Bayesian learning rule (BLR) of Khan & Rue 2023?\"\n\n## Answer 2\n\nBLR assumes that the *variational* family is itself minimal and regular exponential family, which rules out many heavy-tailed or over-dispersed targets (e.g. Student-$t$, Negative-Binomial, Laplace, ...). Our QBLR only requires an exponential *lift*; the marginal can be non-exponential. Hence QBLR:\n\n1. extends natural-gradient methods beyond exponential families;\n2. preserves cheap Fisher pre-conditioning (no explicit inverse) for the gradient.\n\n## Question 3\n\n\"What are the measurable moment functions in Def. 1? Are they the same as in the supplement? Please discuss existence and uniqueness.\"\n\nWe greatly thank the reviewer for this question, which contributes a lot to the clarity of our work. We will include it in Appendix B and if space would permit in the main text.\n\n### A3.1 Notation inconsistency\n\nThanks for pointing this out, we will clean this notational inconsistency in the final text. Yes, $m_i$ from the main text are identical to the $e^i$ from the Appendix B.3.\n\n### A3.2 Uniqueness\n\nAny collection of functions $m = (m_i)_{0 \\leq i \\leq \\dim \\Xi}$ such that each expectation $\\mathbb{E}[m_i]$ is finite and the Jacobian of the collection $Dm(\\xi)$ has full rank at all $\\xi \\in \\Xi$ suffices; we do *not* require uniqueness. The Laplace example will illustrate our point, where we will show two different and valid expectation parameterizations.\n\n### A3.3 Existence\n\nExistence is inherently *model-specific*. For *scale-mixture* families such as the Laplace and Student–$t$, the latent–scale decomposition supplies moment functions whose Jacobian is full rank. A closely related construction applies to discrete heavy-tailed laws: the Negative-Binomial distribution admits a Poisson–Gamma representation (Greenwood and Yule, 1920), and the Gamma moments furnish the required $m_i$ (see our rebuttal to reviewer Y5N7 Question 1 for more details on this specific example). We do *not* claim that our algorithm is applicable to any distribution family; providing at least one full-rank moment map is a prerequisite for our method.\n\n## Question 4\n\n\"Is the global bijection in Def. 1(iii) too strong? What is the practical benefit?\"\n\n## Answer 4\n\nThe existence of the moment functions is a limitation of our approach, but this limitation is significantly weaker than the limitation of the BLR of Khan & Rue 2023 which requires the family to be a member of the exponential family. Indeed, every member of the exponential family has an expectation parametrisation (the so-called Amari's dual parametrization, see equation (6)). Due to our weaker constraints, the QBLR algorithm uses the advantage of the BLR algorithm of Khan & Rue, but on a much wider range of distribution family.\n\n## Question 5\n\n\"Why mention Lie groups?\"\n\n## Answer 5\n\nWe cite as Kiral et al. (2023) an *approximate* alternative that uses Lie-group geodesics when the exact Fisher geometry is unavailable. Our quotient construction gives the *exact* Fisher–Rao metric for the marginal family, subsuming their intuition; we will clarify this in the final paper. If the reference is really distracting, we are ready to discuss possible removal of it in the camera-ready but we would rather save it.\n\n## Question 6\n\n\"Only logistic regression—any harder models or larger $d$?\"\n\n## Answer 6\n\nThe computational story is already laid out in our reply to **Reviewer Y5N7**, see the paragraph entitled *\"Routes to larger $d$\"*. There we show that two standard factorizations—block-diagonal and low-rank $+$ diagonal—turn every $\\mathcal{O}(d^2)$ step of QBLR into $\\mathcal{O}(d)$ or $\\mathcal{O}(kd)$, so the method remains tractable even when $d \\gg 100$.\n\n**Harder example (planned).** The rebuttal policy forbids updating the experimental link, but we have already derived QBLR updates for a more challenging *Bayesian neural network* with roughly 1,861 weights (60 → 30 → 1 architecture) that replaces the logistic regressor.\n\nWe will:\n1. include the full derivation in the camera-ready appendix;\n2. benchmark this network before final submission; and\n3. release the associated JAX code so that reviewers and the community can easily reproduce and extend the experiment.\n\n## Question 7\n\n\"Projection and inverse seem expensive.\"\n\n## Answer 7\n\nThe projection is not particularly expensive because the inverse can be avoided in the algorithm. Please excuse us; in order to provide a deep and extensive answer to your question while respecting the amount of characters allowed by rebuttal, we invite you to read the answer to question 1 of reviewer SfwL. This also prevents us from a copy paste.\n\n## Question 8\n\n\"Why does the ELBO diverge in Fig. 1 without projection? Have you tried other settings?\"\n\n## Answer 8\n\nThank you for raising this point; the behaviour of the lifted ELBO was indeed expected and observed over all our experiments.\n\nEquation (15) is constant on every fibre $\\pi^{-1}(\\xi)$. If the update is not projected, it acquires a vertical component that moves $\\lambda$ inside the same fibre while leaving $q_\\xi(z)$ unchanged. Because no maximum-entropy distribution exists within most fibres, the first KL term $\\mathbb{E}_{q_\\lambda}[\\log q_{\\lambda_0} - \\log q_\\lambda]$ can be driven to $\\pm\\infty$; the lifted ELBO therefore explodes, whereas the marginal ELBO stays flat. This is the runaway trace in Figure 1.\n\n## References\n\nGreenwood, Major, and G. Udny Yule. \"An inquiry into the nature of frequency distributions representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents.\" *Journal of the Royal statistical society* 83, no. 2 (1920): 255–79.\n\nKhan, Mohammad Emtiyaz, and Håvard Rue. \"The Bayesian Learning Rule.\" *Journal of Machine Learning Research* 24, no. 281 (2023): 1-46.\n\nKiral, Eren Mehmet, Thomas Moellenhoff, and Mohammad Emtiyaz Khan. \"The Lie-Group Bayesian Learning Rule.\" *Proceedings of The 26th International Conference on Artificial Intelligence and Statistics* (2023): 3331–52."}, {"review_id": "MexT8hfFrT", "reviewer": "Reviewer_bBHv", "summary": "The paper introduces an efficient algorithm for performing natural gradient updates to minimise the Bayesian learning rule (BLR) objective with respect to non-exponential distributions that can be represented as marginals of minimal, regular exponential distributions. The proposed method works by performing cheap natural gradient updates on the lifted parameter space, followed by marginalisation to the original parameter space. The approach is justified by quotient manifold theory arguments, giving it the name \"Quotient BLR\" (QBLR). QBLR is evaluated on a variational inference problem for a Bayesian logistic regression model with a Student-t variational posterior, performing favourably against two baselines on four UCI/OpenML datasets.", "questions": "- Could you expand on how big the challenge of finding the minimal, regular exponential lift (Section 7) is? Are there any examples in the literature for other scale mixtures?\n- Could you provide the missing ELBO trace plots? Do you have trace plots with ELBO (y-axis) against runtime (x-axis)? In Section 6, you mention that, even for the datasets where QBLR and BBVI-NS have similar performance, QBLR takes roughly one-tenth as many optimisation steps than BBVI-NS. Given that the single BBVI-NS updates are probably cheaper, ELBO vs runtime would be a better indication of the actual speed-up.\n- I found a bit surprising that BBVI-NS outperforms the NG-LIN baseline across all datasets, given that [1] showed the opposite for vanilla BBVI; is this mismatch due to the BBVI-NS being the score-free variant? Do you expect that standard BBVI would underperform NG-LIN in the same experiments?\n- Appendix E seems to indicate that all methods compared in Section 6 uses the fixed step size $\\beta_0 = 10^{-6}$. How was this step size chosen? Given the different geometric derivation of the three methods, would one potentially expect the optimal learning rate for each method to be different? Do you expect that learning rate scheduling could further improve the performance of your method?\n\n[1] Lin, W., Khan, M. E., & Schmidt, M. (2019). Fast and simple natural-gradient variational inference with mixture of exponential-family approximations. In International Conference on Machine Learning (pp. 3992-4002). PMLR.", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "**Strengths**:\n- Fast and scalable Bayesian learning with non-exponential distributions (such as several heavy-tailed ones) is a valuable problem\n- The presentation is generally clear, with the main ideas and relevant background properly introduced\n- The method proposed for natural gradient optimisation over distributions that admit a minimal, regular exponential lift is sensible, and its geometric justification seems solid\n\n**Weaknesses**:\n- Mismatch between very general presentation and very specific evaluation:\n\t- QBLR applies to distributions that admit a minimal, regular exponential joint representation; however, the only example given of such a construction is the Normal-Wishart parameterisation of the multivariate Student-t distribution (Section 5); in Section 7, it is noted that scale mixtures (like the Student-t) include many other interesting heavy-tailed distributions, for which QBLR could be applied once a minimal, regular exponential lift (like the Normal-Wishart for the Student-t) is found; while this is mentioned as a challenge, it is crucial to the out-of-the-box applicability of QBLR, and it is not clear from the text how big of a challenge it actually is.\n\t- QBLR is initially motivated for general BLR objectives (Equation 2), which is then changed for the more specific ELBO objective (Equation 15) in Section 4 without notice.\n\t- Among ELBO objectives, only Bayesian logistic regression VI training is considered in the experiments (Section 6).\n- ELBO trace plots are missing: in the text, it is mentioned that they are in Appendix E, but I couldn't find them there; ELBO trace plots are an important visual tool to complement Table 2, especially in the cases where black-box VI matches QBLR.\n- Minor clarity:\n\t- Section 2.2: while the scalar Student-t/Normal-Gamma example is useful to make things more concrete, it is mentioned as a *running* example, but it is not actually used in any other section of the main text. Possibly, I think it could be recalled in Section 5, to discuss the changes in the associated joint exponential parameterisation when going from scalar to multivariate Student-t distribution.", "was_revised": true, "final_justification": "My main concerns were addressed by the authors' rebuttal. After reviewing the answer to my and other reviews, I have increased my score to a 4. I have decided not to increase it further because, in its current form, I still consider the experimental evaluation a bit limited for the generality of the presentation of the proposed method.", "rebuttal": "## Question 1\n\n\"Could you expand on how big the challenge of finding the minimal, regular exponential lift (Section 7) is? $\\dots$\"\n\n## Answer 1\n\nOur work is indeed limited by the existence of a minimal, regular exponential lift between a distribution family and its extended exponential family representation. Whether or not a lift can be found is then a crucial question for us, and we thank the reviewer for highlighting it.\n\nThe heavy-tailed distributions families that have an exponential family scale mixture representation are well documented in the literature. Regarding the particular case of being a scale mixture of normal distributions, Andrews and Mallows (1974) provides necessary and sufficient conditions in their paper which is applied to Student-$t$, Laplace, and Logistic distributions. If you are interested in more examples, we suggest a book by Coelho and Chen (2024).\n\nMany scale‑mixture joints found in the literature come in a *curved* form—that is, their sufficient statistics are linearly dependent, so the family is *not* minimal. The textbook parameterisations of both the Normal–Wishart and the Normal–Exponential joints fall into this category. (By contrast, the Normal–Wishart lift we use—see Appendix A—is explicitly minimal; the distinction is made concrete in the Laplace example that follows.)\n\nBecause a curved exponential family violates the minimal‑regular assumption, it cannot serve as a lift for QBLR *unless* one first \"uncurves\" it by adding extra, independent natural parameters. We now explain why this step is necessary and how those additional degrees of freedom restore minimality.\n\n**Curved vs. minimal lift for the Laplace distribution** The textbook Normal–Exponential mixture\n\n$$\n\\text{Lap}(z\\mid\\mu,b ) = p(z\\mid\\mu,b)\n= \\int_{0}^{\\infty} \\mathcal N\\bigl(z\\mid\\mu,\\tau\\bigr)\n  \\operatorname{Exp}\\Bigl(\\tau\\Bigm|\\lambda=\\tfrac{1}{2b^{2}}\\Bigr)\n  \\mathrm d\\tau\n\\tag{L‑curved}\n$$\n\n*is curved*: the joint density has **three** sufficient statistics,\n$T_1=\\tau(z-\\mu)^2$, $T_2=\\tau$, $T_3=\\log\\tau$,\nbut only **two** parameters $(\\mu,\\lambda)$.\nHence the Jacobian of $T$ has rank 2 and Brown's minimality criterion fails (Brown, 1986, Proposition 1.5) minimal representation is not possible and the family is internally curved.\n\n**Open question.** We do *not* know whether *every* curved exponential family can be uncurved by judiciously adding degrees of freedom. Establishing necessary and sufficient conditions is, to our knowledge, an open problem in exponential‑family theory. We therefore present our \"add‑a‑free‑hyperparameter\" trick as an *empirical recipe*, not a theorem.\n\n### Uncurving with one extra degree of freedom\n\nWe introduce an independent positive scalar $\\kappa>0$ that scales the latent variance:\n\n$$\n\\text{Lap}(z\\mid\\mu,\\kappa b ) = p(z\\mid\\mu,\\kappa,b) = \\int_{0}^{\\infty} \\mathcal N\\bigl(z\\mid\\mu,\\kappa\\tau\\bigr)\n  \\operatorname{Exp}\\Bigl(\\tau\\Bigm|\\lambda=\\tfrac{1}{2b^{2}}\\Bigr)\n  \\mathrm d\\tau .\n\\tag{L‑minimal}\n$$\n\nNow the joint has three *independent* parameters $\\kappa$, $\\mu$, $\\lambda$, all varying over an *open* domain $\\mathbb R\\times\\mathbb R_{>0}^{2}$. The rank of the sufficient-statistic map matches the parameter dimension and natural parameter representation is possible, and it is only a question of the algebraic manipulation to find it, which for the matters of space we avoid here, so (L-minimal) is a **minimal, regular exponential lift** suitable for QBLR.\n\n**Practical Consequence.** Our \"uncurving\" trick adds parameters to the family representation, but compared to the $O(d^{2})$ number of parameters, this $O(d)$ addition of parameters is negligible — for example, for Normal-Wishart.\n\n## Question 2\n\n\"Could you provide the missing ELBO trace plots? $\\dots$\"\n\n## Answer 2\n\nThank you for catching this oversight. The ELBO–iteration traces were generated but were inadvertently left out of the PDF; in addition, we fully agree that ELBO–*runtime* plots give a more honest picture of QBLR's efficiency.\n\nUnfortunately the rebuttal policy for this year does not permit us to insert new figures during the author‑response period, so we cannot attach the plots here. We will therefore:\n\n1. include both ELBO–iteration and ELBO–runtime panels (one per dataset) in the camera‑ready appendix;\n2. release the plotting script alongside our code so reviewers and readers can reproduce or extend the traces immediately.\n\nWe apologize for the omission and appreciate the suggestion to add the runtime comparison.\n\n## Question 3 \"I found a bit surprising that BBVI-NS $\\dots$\"\n\n## Answer 3\n\nThank you for flagging this inconsistency. In hindsight, we believe the reversal is an artifact of the *single, fixed* step size ($10^{-6}$) that we applied to *all* methods. The very next answer (Q4) describes our new, parameter‑free learning rate schedule; with that adaptive choice, NG‑LIN gains its expected advantage over BBVI‑NS. Under identical tuning, we also expect vanilla BBVI (with score term) to lag behind NG‑LIN.\n\nIrrespective of the final acceptance decision, we see this clarification as a valuable contribution to the community.\n\n## Question 4\n\nAppendix E seems to indicate that all methods compared in Section 6 uses the fixed step size $\\beta_{0} = 10^{-6}$. How was this step size chosen? Given the different geometric derivation of the three methods, would one potentially expect the optimal learning rate for each method to be different? Do you expect that learning rate scheduling could further improve the performance of your method?\n\n## Answer 4\n\n**Static choice (in the submission).** We conducted a grid search on the Diabetes dataset. For $\\beta_{0}>10^{-6}$ BBVI‑NS diverged; for $\\beta_{0}<10^{-6}$ NG‑LIN and QBLR converged too slowly. To keep the comparison fair, we fixed $\\beta_{0}=10^{-6}$ for *all* methods and datasets.\n\n### Adaptive schedule study\n\nMotivated by your questions, we tried to explore options to remove the single hand-tuned hyper-parameter still present in our experiments—the learning rate—we switched to the *Distance‑over‑Gradient* (DoG) rule of Ivgi et al. (2023), together with its Riemannian generalisation RDoG (Dodd et al., 2024). Both schedules are *parameter‑free*: at every iteration they set the step size to the largest value that has been \"empirically safe\" so far, using only distances that the algorithm can measure on the fly.\n\n**Euclidean DoG.** For parameters $x_{t}$ with Euclidean gradient $g_{t}$ we keep\n\n$$\n\\bar r_{t}=\\max\\Bigl(\\epsilon,\\max_{s\\le t}\\lVert x_{s}-x_{0}\\rVert_{2}\\Bigr), \n\\quad\nG_{t}= \\sum_{i\\le t}\\lVert g_{i}\\rVert_{2}^{2},\n$$\n\nand set\n\n$$\n\\eta_{t}= \\frac{\\bar r_{t}}{\\sqrt{G_{t}}},\n\\quad\nx_{t+1}=x_{t}-\\eta_{t} g_{t}.\n$$\n\nNo learning‑rate constant needs to be supplied. We applied this directly to **BBVI‑NS**.\n\n**Riemannian DoG (RDoG).** Replace Euclidean norms by natural-gradient norms and Euclidean distance by the geodesic distance $d(\\cdot,\\cdot)$ generated by the metric $g$:\n\n$$\n\\bar r_{t}= \\max\\Bigl(\\epsilon,\\max_{s\\le t} d(x_{s},x_{0})\\Bigr),\\qquad\nG_{t}= \\sum_{i\\le t}\\lVert g_{i}\\rVert_{g,x_{i}}^{2},\n$$\n\n$$\n\\eta_{t}= \\frac{\\bar r_{t}}{\\sqrt{\\zeta_{\\kappa}(\\bar r_{t})G_{t}}},\n\\quad\nx_{t+1}= \\exp_{x_{t}}\\bigl(-\\eta_{t} g_{t}\\bigr),\n$$\n\nwhere $\\zeta_{\\kappa}$ is the standard curvature correction ($\\zeta_{\\kappa}=1$ in non-positive curvature). We used RDoG for\n\n- **NG‑LIN** — symmetric KL distance between two Student‑$t$ distributions, approximated by $64$ Monte‑Carlo samples drawn once from the starting point and re‑used at every iteration;\n- **QBLR** — exact symmetric KL distance *in the lift*, available in closed form for any minimal, regular exponential family (here: Normal–Wishart).\n\n**Safety of the lift‑based distance.** Theorem 2 implies that for any two lifted points $\\lambda_{1},\\lambda_{2}$ in different fibres,\n\n$$\n\\mathrm{KL}\\bigl(q_{\\lambda_1}\\|q_{\\lambda_2}\\bigr)\n\\geqslant\n\\mathrm{KL}\\bigl(q_{\\xi_1}\\|q_{\\xi_2}\\bigr),\n\\qquad\n\\xi_i=\\pi(\\lambda_i).\n$$\n\nThus the distance we plug into RDoG is an *upper bound* on the unknown marginal KL. Because DoG chooses the step size $\\eta_{t}=\\bar r_{t}/\\sqrt{\\zeta_{\\kappa}(\\bar r_{t})G_{t}}$, a larger distance $\\bar r_{t}$ translates into a *larger*, i.e. mildly optimistic, step. A tighter alternative—still open for future work—would minimize the lift KL over each fibre (this is the true metric on the quotient):\n\n$$\\inf_{\\lambda_i\\in\\pi^{-1}(\\xi_i)}\n   \\mathrm{KL}\\bigl(q_{\\lambda_1}\\|q_{\\lambda_2}\\bigr)\n=\\mathrm{KL}\\bigl(q_{\\xi_1}\\|q_{\\xi_2}\\bigr),\n$$\n\nwhich would shrink $\\bar r_{t}$ and yield an even safer schedule.\n\n**Outcome.** With these parameter-free schedules the ranking becomes $\\text{QBLR} > \\text{NG‑LIN} \\gtrsim \\text{BBVI‑NS}$ on three of the four datasets, and the QBLR–NG-LIN gap widens on Sonar ($d=60$). We will report the full RDoG results in the camera-ready and release the scheduler code for reproducibility.\n\n## References\n\nAndrews, D. F., and C. L. Mallows. ‘Scale Mixtures of Normal Distributions’. Journal of the Royal Statistical Society Series B: Statistical Methodology 36, no. 1 (September 1974): 99–102. https://doi.org/10.1111/j.2517-6161.1974.tb00989.x.\n\nBrown, Lawrence D. Fundamentals of Statistical Exponential Families with Applications in Statistical Decision Theory. SPIE, 1986. https://doi.org/10.1214/lnms/1215466757.\n\nCoelho, Carlos A., and Ding-Geng Chen, eds. Statistical Modeling and Applications: Multivariate, Heavy-Tailed, Skewed Distributions and Mixture Modeling, Volume 2. Emerging Topics in Statistics and Biostatistics. Cham: Springer Nature Switzerland, 2024. https://doi.org/10.1007/978-3-031-69622-0.\n\nDodd, Daniel, Louis Sharrock, and Christopher Nemeth. ‘Learning-Rate-Free Stochastic Optimization over Riemannian Manifolds’. arXiv Preprint arXiv:2406. 02296, 2024.\n\nMaor Ivgi, Oliver Hinder, and Yair Carmon. ‘DoG Is SGD’s Best Friend: A Parameter-Free Dynamic Step Size Schedule’. In International Conference on Machine Learning, 14465–99. PMLR, 2023."}, {"review_id": "VfEQEuv3m7", "reviewer": "Reviewer_MVLt", "summary": "This paper addresses the problem of learning the parameters of a distribution that is the marginal of a joint distribution, where the latter (but generally not the former) lies in an exponential family.  It develops a detailed analysis of how one should properly extend the so-called Bayesian Learning Rule methodology for exponential families to this scenario, relates the result to existing alternative approaches, and provides experimental validation and comparisons.  Further discussion of the results is also included.", "questions": "- The paper develops how to preserve Fisher-Rao metric for the scenario of this paper, and Figure 1 provides evidence for its value.   But this is just one example.  Is it representative?   More generally how sensitive are the results to preserving this either exactly or approximately?   In what settings/regimes is such preservation important in practice?    And when it is not, are there simpler alternatives?  There are some comments in the discussion section, but they are quite brief.   \n\n- In the experiments, sometimes the new method is significantly better than existing ones, but in other cases it is comparable?   What explains why the methodology matters more in some cases than others?   More generally, how might practitioners recognize in advance scenarios where they should anticipate a significant improvement from the proposed approach?\n\n- It would be helpful to understand the scaling implications of the results (when d is large).  In particular, the empirical results show statistical performance, but is the computational complexity of the proposed methodology similar to those being compared?  More importantly, if the methodology is applied to much larger problems, how should one expect the computational and sample complexity to scale (relative to other methods or baselines).", "limitations": "Yes.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n- The paper is very well written:  clear exposition and well organized.  \n\n- The analysis yielding the resulting methodology is carefully developed and quite thorough, together with accompanying intuition.  A nice contribution to the statistics literature.\n\n- The numerical evaluations and comparisons (for the Bayesian logistic regression case) are meaningful.\n\nWeaknesses/Issues:\n\n- The examples and comparisons are a bit limited.\n\n- There is limited discussion of computational complexity and scaling behaviors in high-dimensional settings.\n\n- The demonstration and discussion of the importance of horizontal space projection step is bit limited", "was_revised": true, "final_justification": "Based on the response to my reviews, and on the other reviews and their responses, I will keep my score unchanged.", "rebuttal": "We sincerely thank the reviewer for a careful reading and for the three insightful questions. Below we address each point in turn.\n\n## Q1\n\n\"Only one example for the Fisher–Rao projection—how representative is it, and when does exact preservation matter?\"\n\n## A1\n\nFigure 1 is merely the most compact illustration. Whenever the lifted family contains fibres on which no maximum-entropy member exists—a common situation for heavy-tailed posteriors—*any* naïve optimisation in the joint space can drift vertically and blow up the lifted ELBO. We observed the same runaway trace on all four datasets whenever the horizontal projection was disabled.\n\n## Q2\n\n\"Why is QBLR sometimes much better and sometimes only comparable? How can a practitioner tell in advance?\"\n\n## A2\n\nThe gain depends on two factors:\n\n(a) **Tail heaviness of the posterior.** The more mass the posterior puts in the tails, the larger the vertical drift in the lift and the bigger QBLR's advantage. This explains the pronounced gap on Sonar (heavier tails, $d=60$) and the smaller gap on Diabetes (lighter tails, $d=8$).\n\n(b) **Condition number of the log‑likelihood Hessian.** When the Hessian is ill‑conditioned, the natural‑gradient pre‑conditioner matters, hence QBLR $\\gtrsim$ NG‑LIN $\\gg$ BBVI‑NS; when the Hessian is well‑behaved, all three converge similarly.\n\n*Rule of thumb.* Start with the cheap runs of BBVI‑NS. If with a large gradient norm, only a little ELBO progress is observed, QBLR is likely to yield a substantial improvement.\n\n## Q3\n\n\"How does the method scale with large $d$?\"\n\n## A3\n\nThe raw cost of QBLR is $\\mathcal{O}(d^{2})$ per iteration—the same as NG‑LIN—and is dominated by two Fisher–vector products in the lifted space. The \"Routes to larger $d$\" paragraph in our reply to Reviewer Y5N7 details two standard factorizations that retain the quotient geometry:\n\n- *Block‑diagonal scale:* $\\mathcal{O}\\bigl(\\sum_{b} d_{b}^{2}\\bigr)$; purely diagonal: $\\mathcal{O}(d)$.\n- *Low‑rank + diagonal:* rank $k$ yields $\\mathcal{O}(kd)$ memory and $\\mathcal{O}(k^{2}d)$ time (linear in $d$ for fixed $k$).\n\nThese structures reduce wall-clock time to that of BBVI-NS while retaining the faster convergence of a natural-gradient method. A concise guide and timing table will be added to the camera-ready version.\n\n**Once again, thank you for the thoughtful feedback and encouraging assessment. Your comments helped us sharpen both the theory and its practical guidance.**"}, {"review_id": "9QDDVq7vZK", "reviewer": "Reviewer_SfwL", "summary": "Natural gradient descent (NGD) is more suitable for Bayesian machine learning, and prior works have shown that for exponential family distributions, the natural gradient over natural parameter is the same as the gradient over mean parameter, which allows efficient computation of NGD. This paper aims to extend this result into probability distributions that fall outside the exponential family. Denote the random variable of interest as z, this paper proposes to find a joint distribution (z, zV) which lies in the exponential family, compute the natural gradient on the joint distribution, and then map the update back into the variable of interest z. The authors show theoretical guarantee of the proposed method and evaluate the effectiveness on Student-t distribution.", "questions": "- How expensive is the projection onto the horizontal space (Eq. (16))? How to choose K(lambda)? What is F(lambda) in this equation? \n\n- Line 272, by making posterior mean italics, does it mean the results are computed by making prediction with the posterior mean value? This seems very strange to me. Shouldn't the results be reported by (1) sample weights from the posterior (2) making prediction with the sampled weights (3) report the mean and std?\n\n- How fast is the proposed method compared with previous methods?", "limitations": "See weakness and questions above.", "rating": 5, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength:\n- This proposed method allows efficient and stable learning for Bayesian models with distributions lie outside the exponential family. The proposed method seems principled and elegant (I lack background in geometry and optimisation to verify the solidity, but from reading the paper, it appears to be like this).\n- The paper is very well written and easy to follow.\n\nWeakness:\n- The proposed method is very distribution-specific, as in we need derive everything for any new distribution. In the paper, the authors showed results for student-t distribution, I don't know how easily applicable the proposed method is for other non-exponential family distributions.\n- The evaluation is very limited, and the improvement is also not very obvious. The experiment results are reported on a single train/test split of the UCI datasets. Given the relatively small size of these datasets, the results would be more convincing if averaged over multiple random splits. I would be interested to see the result over say, 5 splits. Also, 10 seems to be a way too small number for sampling.", "was_revised": true, "final_justification": "The rebuttal has addressed my concerns and hence I increased my score to 5.", "rebuttal": "### Question 1\n\nHow expensive is the projection onto the horizontal space Eq. (16)? How to choose $K(\\boldsymbol{\\lambda})$? What is $F(\\boldsymbol{\\lambda})$ in this equation?\n\n### Answer 1\n\n\n\nWe thank the reviewer for raising this insightful question. We plan to incorporate the answer into the final version of the paper, as it significantly enhances the clarity and precision of our work. We are grateful for the reviewer’s valuable contribution to improving the manuscript.\n\nTo aid clarity, we first restate the projection operator, (Equation (16) in the paper)\n$$\n    P_{\\mathcal H}(\\boldsymbol{\\lambda}) = \\operatorname{I} - K(\\boldsymbol{\\lambda})\\left[K(\\boldsymbol{\\lambda})^{\\top} F(\\boldsymbol{\\lambda}\\,K(\\boldsymbol{\\lambda})\\right]^{-1} K(\\boldsymbol{\\lambda})^{\\top} F(\\boldsymbol{\\lambda}),\n$$\n which projects a vector from the tangent space $T_{\\lambda}\\Lambda$ onto the horizontal space $\\mathcal{H}_{\\lambda} := \\text{ker} D \\pi(\\boldsymbol{\\lambda})$ in the Fisher–Rao geometry.\n\n\n\n**What is $F(\\boldsymbol{\\lambda})$ in this equation?** \n We start by answering the last question to make our exposition clearer. $F(\\boldsymbol{\\lambda})$ is the Fisher information matrix defined in Equation (3a). For further exposition, it suffices to note that it is a symmetric positive-definite matrix of dimensionality $\\text{dim} \\Lambda \\times \\text{dim} \\Lambda$, where $\\text{dim} \\Lambda$ is the dimensionality of the extended distribution parameter space. The extended distribution is defined in equation (10).\n\n**How to choose $K(\\boldsymbol{\\lambda})$?** \n\nWe remind the reader that $K(\\boldsymbol{\\lambda})$ is the matrix representation of a basis of the vertical space $\\ker D\\pi(\\\\boldsymbol{\\lambda})$ (36) where $\\pi:\\Lambda \\rightarrow \\Xi$ is the marginalization map (10).\nSo, $K(\\boldsymbol{\\lambda})$ is a matrix of dimensionality $\\text{dim} \\Lambda \\times d_{v}$ where $d_v$ is the dimensionality of the vertical space. \n\nA generic method for computing $K(\\boldsymbol{\\lambda})$ is based on the singular value decomposition  [Klema and1302\nLaub, 1980]  of the Jacobian representation  $J_\\pi(\\boldsymbol{\\lambda})$ of $D\\pi(\\boldsymbol{\\lambda})$. Note that the Jacobian representation $J_\\pi(\\\\boldsymbol{\\lambda})$ is justified in the paragraph starting at _line 135_ in the Section 3. $K(\\boldsymbol{\\lambda})$ can then be extracted from the singular value decomposition of $J_\\pi(\\boldsymbol{\\lambda})$:  the group of vectors that correspond to the singular value $0$ is the basis of the vertical space $\\ker D\\pi(\\boldsymbol{\\lambda})$. To compute the singular value decomposition, the strategy divide-and-conquer can be applied. A notable example is the randomized Krylov method by Musco and Musco [2015]. For a comprehensive understanding of these algorithms and their analysis, we recommend the tutorial by  Li et al. [2019]. The time complexity of the best algorithms is $\\mathcal{O}(\\dim\\Lambda\\dim\\Xi^2)$.\n\n\nUsually, in practice, an analytical expression for $K(\\boldsymbol{\\lambda})$ can be derived, and the singular value decomposition can be bypassed. In our main Normal-Wishart distribution example, we derived an analytical expression for $K(\\boldsymbol{\\lambda})$, which can be found in Appendix A.3 at _line 829_.\n\n\n**Projection onto the Horizontal space complexity?** \n\nThe projection operator defined at the beginning of the answer  is used to compute the horizontal component of the natural gradient. For clarity, we restate the projection applied to the natural gradient $ g_{\\theta} \\in \\mathbb{R}^{\\dim \\Lambda} $ as follows:\n$$\ng_{H} = P_{\\mathcal{H}}(\\boldsymbol{\\lambda})[g_{\\theta}].\n$$\n\nA naïve approach to compute $g_H$ includes the inversion of the symmetric positive definite matrix \n\n$$\nRV(\\boldsymbol{\\lambda}) := K(\\boldsymbol{\\lambda})^{\\top} F(\\boldsymbol{\\lambda})K(\\boldsymbol{\\lambda})\n\\in \\mathbb{R}^{d_v \\times d_v}\\,,\n$$\n\nwhere $ d_v = \\dim V_\\lambda $. This approach would require computing the matrix inverse, with cost $ \\mathcal{O}(d_v^3) $. However, we avoid materializing the inverse. Instead, we solve the linear system\n\n$$\nRV(\\boldsymbol{\\lambda}) v = K(\\boldsymbol{\\lambda})^{\\top} F(\\boldsymbol{\\lambda})g_\\theta,\n$$\n\nusing a small dense linear solver: due to the fact that $RV(\\boldsymbol{\\lambda})$ is symmetric positive definite, the conjugate gradient solver by  Hestenes et al. [1952] is applicable. This reduces the computational cost to $ \\mathcal{O}(d_v^2) $.\n\nIn practice, a well-constructed lifting ensures that $ d_v \\ll \\dim \\Lambda $, making the cost of projection negligible relative to the Fisher--vector product $ F(\\boldsymbol{\\lambda}) g_\\theta $, which has cost $ \\mathcal{O}(\\dim \\Lambda^2) $.\n\nIn our main Normal--Wishart case, the vertical space is one-dimensional ($ d_v = 1 $), so the projection reduces to a single scalar division.\n\n\n### Question 2 \n\n_Line 272_, by making posterior mean italics, does it mean the results are computed by making prediction with the posterior mean value? This seems very strange to me. Shouldn't the results be reported by (1) sample weights from the posterior (2) making prediction with the sampled weights (3) report the mean and std?\n\n### Answer 2\nThank you for catching this point and for reminding us of the proper Bayesian evaluation protocol.\n\n\n**Policy note.**\nConference rules this year prevent us from adding new tables or updating our experiments during rebuttal, so we summarize the outcome verbally here.\n\n\n**Why we originally used the posterior mean.** \n The experiment is conducted on the Generalized Student-$t$ distribution family. For Generalized Student-$t$ distribution, the mean parameter is equal to the maximum a posteriori estimator. Reporting test accuracy through a maximum a posteriori estimator provides a noise-free benchmark. This approach is commonly used in logistic regression studies that prioritize decision boundaries over calibrated probability estimates.\n\n\n\n**Posterior‑predictive check (100 Monte‑Carlo draws).** \nFollowing your suggestion, we re-evaluated each inference method by drawing 100 weight samples from the variational posterior, predicting with each sample, and averaging the test accuracy. Across all four datasets, results \\emph{confirm} the previous ranking: QBLR remains best, widening its margin on the Sonar task. We will include this metric in the final table of the results.\n\n\nWe hope this clarifies our choice of metric and demonstrates that QBLR’s advantage extends from the mode to the full predictive distribution.\n\n### Question 3\n\nHow fast is the proposed method compared with previous methods?\n\n### Answer 3\n\nNG--Ours and NG--LIN each cost $\\mathcal{O}(d^{2})$ per iteration, whereas BBVI--NS is $\\mathcal{O}(d)$.  \nThe difference stems from algorithmic order: BBVI--NS is a first--order Euclidean method that evaluates only the score--function gradient, while the two natural--gradient variants pre--condition that gradient with the Fisher information; applying the (dense) Fisher block entails two Fisher--vector products and is therefore quadratic. Empirically, the higher iteration cost is offset by faster convergence: NG--Ours reaches its final ELBO in roughly one--tenth the iterations required by BBVI--NS and about half those of NG--LIN, so overall wall--clock time remains comparable while the accuracy gains reported in Table~2 are preserved.\n\n\n\n### References\n\nKlema, V., and A. Laub. ‘The Singular Value Decomposition: Its Computation and Some Applications’. IEEE Transactions on Automatic Control 25, no. 2, 1980: 164–76. https://doi.org/10.1109/TAC.1980.1102314.\n\nMusco, Cameron, and Christopher Musco. ‘Randomized Block Krylov Methods for Stronger and Faster Approximate Singular Value Decomposition’. In Advances in Neural Information Processing Systems, Vol. 28. Curran Associates, Inc., 2015. https://proceedings.neurips.cc/paper_files/paper/2015/file/1efa39bcaec6f3900149160693694536-Paper.pdf.\n\nLi, Xiaocan, Shuo Wang, and Yinghao Cai. ‘Tutorial: Complexity Analysis of Singular Value Decomposition and Its Variants’. arXiv, October 2019. https://doi.org/10.48550/arXiv.1906.12085.\n\nHestenes, Magnus R., and Eduard Stiefel. ‘Methods of Conjugate Gradients for Solving Linear Systems’. Journal of Research of the National Bureau of Standards 49, no. 6 (1952): 409–36."}, {"review_id": "MvRfQqjgRP", "reviewer": "Reviewer_Y5N7", "summary": "The paper introduces the Quotient Bayesian Learning Rule (QBLR), an extension of natural-gradient variational inference used for probabilistic models whose variational families lie outside the exponential family. Drawing on the observation that many heavy-tailed distributions (e.g. Student-t) arise as marginals of minimal exponential families (e.g. Normal–Wishart), the authors show that the latent family’s Fisher–Rao geometry induces a well-defined Riemannian structure on the marginal parameter space via a quotient-manifold construction. Building on this geometry, they give Algorithm 1, which (i) lifts the variational objective into the exponential-family space, (ii) computes ordinary (Euclidean) gradients in the joint natural parameters, (iii) projects those gradients onto the horizontal subspace orthogonal to the fibres of the marginalisation map, and (iv) takes natural-gradient steps that are parameterization-invariant. Specialising to the multivariate Student-t via its Normal–Wishart scale-mixture representation, they derive a quotient natural-gradient update that avoids the Fisher inverse. Empirically, they benchmark QBLR against (a) score-function–free black-box VI (BBVI-NS) and (b) the curved-exponential-family method of Lin et al. (2020) on Bayesian logistic regression across four UCI datasets.", "questions": "* Are there other examples that you could have considered where the quotient Bayesian rule would be expected to display a significant improvement over existing approaches? Related to this, how does QBLR perform on higher dimensional examples?\n* Would it be possible to include a section, or perhaps a figure, in the paper to help readers gain some intuition for how the quotient manifold works?\n* The framework seems to be quite general, but the authors do not seem to propose a wide range of models which QBLR would be appropriate for. Could more information be given on other model classes beyond the heavy tail setting?", "limitations": "The authors do not adequately address the limitations of their work. In the case of societal issues, there are no obvious issues that the authors would need to address.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Quality\n\nThe paper is presented with a high degree of theoretical rigour, with both the manifold-structure theorem and the induced Fisher–Rao metric theorem stated clearly and proved in full detail. The algorithmic exposition is similarly precise: Algorithm 1 and its specialisation to the multivariate Student-t family are given in closed-form, and practical implementation details, such as constructing bases for vertical subspaces and avoiding large Fisher inversions, are addressed in the paper. Moreover, by exploiting the Normal–Wishart lift, each update reduces to small-scale linear solves, ensuring computational efficiency which means the proposed approach is potentially viable for moderate-dimensional models. On the downside, the empirical scope is relatively narrow: all benchmarks are confined to Bayesian logistic regression on UCI datasets, leaving open questions about performance in higher-dimensional settings or even for other models in general. In addition, the framework’s reliance on identifying a minimal, regular exponential-family representation means that extending QBLR to other heavy-tailed or skewed distributions may be nontrivial and is not yet demonstrated.\n\nClarity\n\nThe paper is well-organised in a logical, (mostly) reader-friendly sequence. The paper begins with motivation and background and then moving on to geometric constructions and the algorithm itseflf. Plus, the experiments are well-presented, as well is the running Normal–Gamma example. The results are also clearly presented. However, the paper is quite dense and the information-geometric notation (quotient manifolds, horizontal lifts, fibre bundles, etc.) would be a steep learning curve for those without a differential-geometry background. Some additional (and quite useful) technical definitions are relegated to the appendices. The heavy reliance on the supplementary material does disrupt the narrative flow for readers who would find it easier to read a self-contained version of the paper within the main text.\n\nSignificance\n\nBy framing heavy-tailed variational inference as an optimisation problem on a quotient manifold, the authors provide a new geometric template that subsumes classical Bayesian learning rules for exponential families and earlier non-EF extensions. This general framework offers a more robust, parameterisation-invariant inference scheme in the presence of outliers or multimodal posteriors, potentially broadening the applicability of variational methods to real-world problems where heavy tails are common. Nonetheless, the paper stops short of validating this promise on more challenging models—such as high-dimensional latent spaces or deep hierarchies or really anything other than logistic regression—and does not compare against state-of-the-art non-natural-gradient VI techniques (for example, advanced black-box methods with control variates). The results that are stated are also not significantly better than competing methods, perhaps with the exception of the Sonar data example. As a result, the full practical significance of QBLR remains to be established.\n\nOriginality\n\nThe introduction of the quotient-manifold perspective for variational inference represents an interesting geometric insight, and the horizontal-space projection mechanism ensures parameterisation invariance across arbitrary smooth reparameterisations. These contributions stand out as original within the information-geometry literature. At the same time, the paper could engage more deeply with related work—particularly Lie-group BLR approaches that also aim to preserve geometric invariances under group actions—and clarify the precise distinctions. Moreover, although the Student-t instantiation is useful, it follows straightforwardly from the general theory and may strike some readers as an incremental extension rather than a wholly new algorithmic advance.", "was_revised": true, "final_justification": "The authors have provided some excellent additions to the paper, reflecting the suggestions of multiple reviewers, which I believe will enhance the quality of this work.", "rebuttal": "We appreciate the reviewer's thoughtful questions and constructive feedback.\n\n> We apologize if some equations do not render perfectly; we have done our best to ensure compatibility with the OpenReview markdown compiler.\n\n## Question 1\n\n\"Are there other examples that you could have considered where the quotient Bayesian rule would be expected to display a significant improvement over existing approaches? Related to this, how does QBLR perform on higher dimensional examples?\"\n\n## Answer 1\n\n### Are there other examples that you could have considered where the quotient Bayesian rule would be expected to display a significant improvement over existing approaches?\n\nHowever, we do not have a complete answer to this question; we can illustrate our intuition with the following example. Consider the negative binomial distribution introduced by Greenwood and Yule (1920); it's a distribution over $\\{0, 1, 2, \\dots \\}$ with the following probability density:\n\n$$\nf(k| r, p) = \\binom{k+r-1}{k} (1-p)^k p^{r},~r>0, p \\in (0, 1).\n$$\n\nThis distribution has the following integral representation, the so-called Poisson-Gamma mixture:\n\n$$\nf(k| r, p) = \\int_{0}^{\\infty} \\mathcal{P}(k|\\lambda)\\mathcal{G}\\left(\\lambda \\big{|}r, \\frac{p}{1-p}\\right) \\mathrm{d} \\lambda,\n$$\n\nwhere $\\mathcal{P}$ and $\\mathcal{G}$ represent the Poisson and the Gamma densities, respectively. The negative binomial distribution is not within the exponential family. However, the joint distribution (the under integral expression) over $(k, \\lambda)$ is an exponential family distribution, which can be represented as a minimal and regular exponential family distribution. The negative binomial distribution is known to be within the exponential family once $r$ is fixed.\n\nThis example shows a contrast between the initial definition of the negative binomial distribution and its restricted definition ($r$ fixed). Natural gradient descent cannot be performed on the negative binomial distribution family, but it can be performed on the exponential family structure of the Gamma-Poisson mixture (the integral representation above). The QBLR method enables one to transform the natural gradient descent on the Gamma-Poisson mixture into an efficient gradient descent into the negative binomial distribution family at low cost.\n\nIn the restricted case, the family is an exponential family distribution, so natural gradient descent optimisation is already possible. A natural gradient on the Gamma-Poisson representation (the integral representation above) will then be equivalent and will not produce any benefit but will introduce a computational burden: we will run optimisation in the higher-dimensional space.\n\n### Related to this, how does QBLR perform on higher‑dimensional examples?\n\nOur current benchmarks already span almost an order of magnitude in feature dimensionality—from $d=8$ (Diabetes) to $d=60$ (Sonar)—and therefore in parameter dimensionality vary from 72 to 3660. Across this range, QBLR delivers consistent gains over both BBVI‑NS and NG‑LIN; however, we acknowledge that, by modern standards, these are still small-scale problems.\n\n**Computational cost.** With a *dense* Normal–Wishart lift, one QBLR iteration is driven by three quadratic terms (see Appendix A.3):\n\n1. *Wishart draw.* Sampling a full $d \\times d$ scale matrix costs $\\mathcal{O}(d^{2})$.\n2. *Fisher products.* A Fisher–vector products $F(\\boldsymbol{\\lambda})g_{\\theta}$ costs $\\mathcal{O}(d^{2})$ but *never* require a matrix inverse.\n3. *Hessian trace.* The Price identities (see Appendix D) that involve $\\mathcal{O}(d^{2})$ when the Hessian is formed explicitly.[^1]\n\nQuadratic time is acceptable up to a few hundred dimensions on commodity hardware, but larger problems call for additional structure. In the following, we propose our plan to reduce computational complexity.\n\n**Routes to larger $d$.** Two common factorizations reduce *all three* $\\mathcal{O}(d^{2})$ components while preserving the quotient geometry:\n\n1. **Structured covariances.** Restrict the scale to $B$ blocks $\\Psi=\\operatorname{diag}(\\Psi_{1},\\ldots,\\Psi_{B})$ with sizes $d_{b}$.\n   - **Sampling**: $\\sum_{b} \\mathcal{O}(d_{b}^{2})$; purely diagonal $\\Psi$ needs only $\\mathcal{O}(d)$ Gamma draws.\n   - **Fisher products & Hessian trace** both factor block‑wise, yielding the same $\\sum_{b} \\mathcal{O}(d_{b}^{2})$ and $\\mathcal{O}(d)$ in the diagonal case.\n\n2. **Low‑rank + diagonal factorisation.** Decompose the scale matrix as\n\n   $$\n   \\Psi = L L^{\\top} + \\operatorname{diag}(v), \\qquad L\\in\\mathbb{R}^{d\\times k}, k\\ll d,\n   $$\n\n   so that only $k d + d$ free parameters are stored instead of $\\frac{1}{2} d(d+1)$.\n\n   - **Sampling.** Each column of $L$ is drawn from a matrix‑normal and the diagonal entries of $v$ from independent Gammas. The two draws require $k d$ and $d$ random numbers, respectively, hence $\\mathcal{O}(k d)$ time and memory.\n\n   - **Fisher–vector products.** In the horizontal projector we need $y = \\bigl(\\operatorname{diag}(v) + L L^{\\top}\\bigr) x$. Compute it as\n\n    $$y = \\operatorname{diag}(v)x + L \\bigl(L^{\\top}x\\bigr).$$\n\n     Both multiplies with $L$ cost $\\mathcal{O}(k d)$, so the total is $\\mathcal{O}(k d)$ per Fisher product—linear in $d$ for any fixed rank $k$.\n\n   - **Hessian trace (two Price identities).** Both the $\\kappa$-Price and the $\\Phi$-Price terms require $\\operatorname{tr}(S^{-1}\\nabla_{z}^{2}f)$. Rather than forming the dense Hessian, we use the Hutchinson estimator (Hutchinson, 1989):\n\n     $$\n     \\operatorname{tr}(S^{-1}\\nabla^{2}f) = \\frac{1}{R} \\sum_{r=1}^{R} u_{r}^{\\top} \\bigl(S^{-1}\\nabla^{2}f\\bigr)u_{r}, \\quad u_{r}\\sim\\{\\pm1\\}^{d}.\n     $$\n\n     Each term needs one Hessian–vector product (HVP) and one multiplication with $S^{-1}$. The HVP is model‑specific; the $S^{-1}$–vector multiply uses the Woodbury identity:\n\n     $$\n     S^{-1}x = D^{-1}x - D^{-1}L \\bigl(I_k + L^{\\top}D^{-1}L\\bigr)^{-1} L^{\\top}D^{-1}x, \\quad D=\\operatorname{diag}(v),\n     $$\n\n     which is again $\\mathcal{O}(kd)$. Choosing $R\\le k$ probes keeps the overall trace cost bounded by $\\mathcal{O}(k^{2}d)$.\n\nBoth variants leave the vertical space one‑dimensional, so the horizontal projection stays a single scalar divide.\n\n**Summary.** Dense QBLR is $\\mathcal{O}(d^{2})$; with a diagonal or block‑diagonal scale it drops to $\\mathcal{O}(d)$, and with rank‑$k$ plus diagonal it is $\\mathcal{O}(k^{2}d)$ (linear in $d$ for fixed $k$). These paths scale QBLR to far larger latent spaces without sacrificing its geometry or requiring matrix inversions.\n\n## Question 2\n\n\"Would it be possible to include a section, or perhaps a figure, in the paper to help readers gain some intuition for how the quotient manifold works?\"\n\n## Answer 2\n\nThank you for the question. We will definitely include one in the final text. We found an interesting example for the quotient manifold theory: the Laplace distribution. This distribution was not initially described precisely in the paper. The Laplace distribution is a two-parameter, non-exponential family distribution that we can represent as the Normal-Exponential scale mixture, which is a three-parametric minimal and regular exponential family. This example is visually pleasing, as we will represent fibers as one-dimensional curves. To strengthen the intuition of the QBLR algorithm, we will even show how the projection of the natural gradient (the quotient natural gradient) moves perpendicular (in the Fisher-Rao metric) to the fibers in the Normal-Exponential scale mixture space.\n\n## Question 3\n\n\"The framework seems to be quite general, but the authors do not seem to propose a wide range of models which QBLR would be appropriate for. Could more information be given on other model classes beyond the heavy tail setting?\"\n\n## Answer 3\n\nOur initial goal was to extend natural gradient descent to the heavy-tail distribution family. However, as you understood, the solution is more general and could be applied to the non heavy-tail distributions family as far as we can define the target distribution family as a marginal of a minimal and regular exponential family distribution. The example of the negative binomial distribution, developed above, shows such a case. We also believe that our result could be applied to the so-called Gaussian Variance-mean Mixture family introduced by Lin et al. (2019), which is a generalization of the scale-mixture family.\n\n## References\n\nGreenwood, Major, and G. Udny Yule. ‘An Inquiry into the Nature of Frequency Distributions Representative of Multiple Happenings with Particular Reference to the Occurrence of Multiple Attacks of Disease or of Repeated Accidents’. Journal of the Royal Statistical Society 83, no. 2 (1920): 255–79.\n\nHutchinson, Michael F. ‘A Stochastic Estimator of the Trace of the Influence Matrix for Laplacian Smoothing Splines’. Communications in Statistics-Simulation and Computation 18, no. 3 (1989): 1059–76.\n\nLin, Wu, Mohammad Emtiyaz Khan, and Mark Schmidt. ‘Stein’s Lemma for the Reparameterization Trick with Exponential Family Mixtures’. arXiv, February 2025. https://doi.org/10.48550/arXiv.1910.13398. \n\n[^1]: For GLM models the Hessian is rank 1, so this step is already $\\mathcal{O}(d)$."}]

Author rebuttal:

Meta-review: The paper extends natural-gradient variational inference to models outside the exponential family that can be obtained by marginalising a minimal exponential-family model. The main example considered is the Student-t distribution. This marginalisation induces a quotient-manifold structure. The key contribution is to show that the Fisher-Rao metric, a natural Riemannian structure on statistical manifolds, extends from the Normal-Wishart to the Student-t manifold, and to use this insight to develop a new algorithm for variational inference with heavy tail variational families.

Reviewers agreed that the theoretical part provides interesting geometric insights and is rigorously and clearly presented. The experiments on Bayesian logistic regression required some clarification, which the authors provided in their rebuttal. It was noted that, while the experiments effectively illustrate the method, the empirical evaluation could be strengthened by going beyond logistic regression. Overall, the paper makes nonetheless a valuable and technically solid contribution.
Decision: Accept (poster)