Least squares variational inference
Authors: Yvann Le Fay
Nicolas Chopin
Simon Barthelmé
Venue: NeurIPS 2025 poster ()
Keywords: variational inference
computational statistics
gradient-free methods
kullback-leibler divergence
bayesian statistics
stochastic optimization
Primary area: probabilistic_methods
TL;DR: Gradient-free variational inference within exponential families using least squares regression with applications on discrete distributions and likelihood-free models.

Abstract: Variational inference seeks the best approximation of a target distribution within a chosen family, where "best" means minimizing Kullback-Leibler divergence. 
When the approximation family is exponential, the optimal approximation satisfies a fixed-point equation.
We introduce LSVI (Least Squares Variational Inference), a gradient-free, Monte Carlo-based scheme for the fixed-point recursion, where each iteration boils down to performing ordinary least squares regression on tempered log-target evaluations under the variational approximation.
We show that LSVI is equivalent to biased stochastic natural gradient descent and use this to derive convergence rates with respect to the numbers of samples and iterations.
When the approximation family is Gaussian, LSVI involves inverting the Fisher information matrix, whose size grows quadratically with dimension $d$.
We exploit the regression formulation to eliminate the need for this inversion, yielding $O(d^3)$ complexity in the full-covariance case and $O(d)$ in the mean-field case.
Finally, we numerically demonstrate LSVI’s performance on various tasks, including logistic regression, discrete variable selection, and Bayesian synthetic likelihood, showing competitive results with state-of-the-art methods, even when gradients are unavailable.

Least squares variational inference
Yvann Le Fay', Nicolas Chopin' Simon Barthelme
ENSAE, CREST, IP Paris
GIPSA-Lab, CNRS {yvann. www.lleaany,nlllll micolas chopin)@ensae.ff simon. www..aaeeeaeeeeeeaaaaaaaaaaae
Abstract
Variational inference seeks the best approximation of a target distribution within chosen family, where "best" means minimising Kullback-Leibler divergence. When the approximation family is exponential, the optimal approximation satis- fies fixed-point equation. We introduce LSVI (Least Squares Variational Infer- ence), a gradient-free, Monte Carlo-based scheme for the fixed-point recursion, where each iteration boils down performing ordinary least squares regression on tempered log-target evaluations under the variational approximation. We show that LSVI is equivalent biased stochastic natural gradient descent and use this to derive convergence rates with respect to the numbers of samples and iterations. When the approximation family is Gaussian, LSVI involves inverting the Fisher information matrix, whose size grows quadratically with dimension d. We exploit the regression formulation to eliminate the need for this inversion, yielding (dd) complexity in the full-covariance case and (d) in the mean-field case. Finally, we numerically demonstrate LSVI's performance on various tasks, including lo- gistic regression, discrete variable selection, and Bayesian synthetic likelihood, showing results competitive with state-of-the the-rr methods, even when gradients are unavailable.
Introduction
This paper focuses on parametric variational inference (VI, [1-3]). Given an (unnormalised) target density n, we aim at finding the distribution that minimises the (reverse) Kullback-Lebler diver- gence:
arg min KT. (q क) = f log (q/#) qEQ
(1)
where Q is user-chosen parametric family (e.g., Gaussians), and ग = JT. This approach has become de facto standard in probabilistic machine learning recent years and implemented in various software packages, such as STAN, NumPyro, PyMC3, and Blackjax [4-7]. The min- imisation typically carried out through gradient-based procedures using automatic differentiation, either stochastic gradient descent (SGD, (8te- applied after reparameterising the target distri- bution [9, 10]-or or its faster alternative natural gradient descent (NGD, [11-14]). This convenient for users, they only have provide the function == log x to the software.
These procedures use different gradient estimators; some require log be amenable automatic differentiation, which is the case when using reparameterisation, while others only require gradient estimators of expectations under the variational distribution via the log- deriviaie trick [15]. The gradient estimator for expectations usually suffers from high variance, and practical implementations
"Corresponingg author
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

rely on the reparameterisation trick, which is not possible in several important cases, for instance when discrete distribution, or when ı is intractable or non-differentiable (as in likelihood- free inference). Additionally, convergence of SGD is sometimes slow and/or tedious to assess, and requires careful step sizes tuning [16] while naive implementation of NGD) requires costly matrix inversions.
1.1 Outline and contributions
We introduce practical algorithms for VI within exponential families when gradients of log are unavailable. These algorithms involve taking biased stochastic gradient descent steps. but we show both theoretical convergence and good performance in non-toy problems. In Section 2, we derive an exact, but intractable, iteration we call LSVI, that boils down to performing successive least squares (OLS) regression. We highlight connections to NGD and discuss its convergence properties. Sec- tion 3, we introduce stochastic variant that updates the OLS estimate using multiple draws from the current approximation. Importanll, under standard smoothness and relative convexity assump- tions the objective, and bounded-momest assumptions convergence guarantees and rates with respect to the numbers of draws and iterations, conditioned on high-probabilltyy events. In addition, we provide an adaptive method to calibrate step sizes by controlling the linear regression residuals. Section 4 focuses on the Gaussian variational family; we propose reparametrisstion of the linear regression such that the OLS procedure requires inver- sion of the Fisher information matrix (FIM). These schemes tailored to Gaussian distributions are cost-efficient: methods scale linearly with d in the mean-field case, and in the full-covariance case, the cost matches the cost of computing d x matrix products, i.e., (dd)). J. In Section 5, we extensively illustrate the performance of our methods compared to other inference procedures, in- cluding gradient-based VI and exact Bayesian inference procedures. Limitations are discussed in Section 6. We provide a Python package supporting GPU parallelisation via JAX to replicate the experiments: https:/ / github. cwwcccmm/lloo.oo cowmylefayLLss//
the variational family, we establish
2
Exact LSVI
Let : X' + F some unnormalised target density, with X C Rd. It will be convenient work with an exponential family Q of unnormalised densities:
(x)x)) == exp(n s(x)),
NEV= = {n: Z(an) &o}
where 77 E V is the natural parameter associated to 9, E Q 2(g) = Sx q denotes the partition function, and X > R. the extended statistic function defined
s sxx) = )
In words, we include an intercept in s to make the family closed under multiplication by a positive scalar. For 7 (1)(0) .77 where ()) denotes the first component of 17, let yn be the normalised version of 9n, (which therefore depends only 1j): q; = 4n/2,, using the short-hand Z, == Z(g,)- Notation [-1[1] means a properly normalised expectation, i.c. [,, Sx h/ Likewise, we replace the standard Kullback-Leibe objective with divergence for unnormalised densities [17], which is defined by
uKL(q m) m):=/ f
f gloo 2 + Z(T) Z(g),
for any density q absolutely continuous with respect to म. addition, we assume the variational family Q is minimal and regular, which is a standard assumption in VI [12, 18 18-20]. and is met by any standard exponential families (e.g.. Gaussian, Beta, Poisson, Bernoulli, etc., [21, Table 3.1.). These assumptions ensure 7) V 9, is injective and the log-partition function is differentiable everywhere [21, Prop.
Assumption 2.1 (Minimality and regularity of Q). The components of s are linearly independent (minimality), and the set of natural parameters V is open (regularity).
The next proposition shows that the critical points of the uKI. divergence are also critical points of the KL divergence, and vice versa. In words, nothing is lost by considering the uKL instead of KL.

---

Proposition 2.2. Let = (1) . " E V. Vf g the reciprocal holds: V, (L, π) and 0,100 u
T) = 0 then V. K म) = 0, and
x then V, UKL π
The first-order condition of the ukL minimisation problem is given by the following proposition. Proposition 2.3. Let f = log n be the (unnormalised) log target density. Let 1 = (1710) 57 ) V.
V, (g, m) and only if {Ensss)) [}} 1 = E,,]Ss). Furthermore, if V. uKL T) = 0,
then n70) = ( #)) + log (Z(T)/ Ix exp (7)
2.1 The exact LSVI scheme
The first-order optimality condition is equivalent to the fixed point equation: ((1) with
<(n) = FT Enp Fn == Essss J] on = [f,ss,
(3)
and Fú is the Fisher information matrix (FIM) associated to qn- Salimans and Knowles [22] remark that is the ordinary least squares regressor (OLS; [23]) of f(X) with respect to s(X) when X 9n;
(n) = argmingeRT En +x) - ((x))]
(4)
A nice property of Φ when is in the variational family with π = 9n-, is that for any 17 E V, ((n) '7" i.e., exactly recovers л. However, in general ((n) may not be in V. and naively performing a fixed-point scheme can lead to unstable variational approximations or, worse, result in non-oormali densities (i.e., ((n)) & V). To address this, we consider a relaxation of the fixed-point scheme obtained via a momentum fixed-point iteration [24]:
M++1 = e26(n)) + (1 62),,,
t 0
(5)
where >> 0 is such that 7 is in V. Such e, necessarily exists because V is open (Assump- tion 2.).. Since iteration (5) assumes that one has access to expectations under the variational family (which in general is not the case), we refer to (5) as the exact Least Squares Variational Inference
(LSVI) iteration. This relaxation has a natural interpretation in this specific context: 77-1 in (5) is the solution of the least squares objective (4) when I exp f is replaced by the tempered (annealed) 1--c, 61 Tt..
density qny
2.2 LSVI as natural gradient descent (NGD) and mirror descent (MD)
This subsection summarises a well-established connection between NGD and MD in the variational inference literature [12, 18, 20] but generalised to the unnormalised KL divergence.
Let us define the (unnormalised) moment parameter mapping w : 17 E V H -,, sxx) and let W ((V) be the set of moment parameters. We denote by 7) : W - the inverse mapping of w : W + V. whose existence is guaranteed under Assumption 2. [21, Ch. 31. Define as the unnormalised KI. divergence (2). When expressed in natural parameters, we write 1 n E V
uKL(g, T), when expressed in moment parameters, we write | : w E W + uKL(gw T), and similarly for expectations: E.,- == B,(w)-
The following proposition states that LSVI iteration (5) is a NGD iteration on the uKL divergence in the natural space of parameters, and equivalently MD in the moment space [25, 26, Ch. 3]. Proposition 2.4 (LSVV is NGD which is equivalent to MD, [20, Lemma 1J). Under Assumption 2. and provided the sequence (1r) defined by (5) is in V. (n)) satisfies the dynamic (NGD).
M++1 = 17- - F ((n())/Zn.
(6)
or equivalently.
Me41 17- Ex v.l(((n))).
Furthermore, let Wo E W and define for 2 0 (MD).
W+41 == argminwew { V. + - (w,ww))
(8)
where Dz- is the Bregman divergence [27] with respect to Z the Legendre transform of Z: Z (w(w) = argmin,,er( w- Zn Then, the sequence (Tr) defined by (5) with yoo = 7(ww)) satisfies for all L 2 0, 1). = n(wx). In words, LSVI performs a natural gradient step in the space of natural parameters, which corresponds to a mirror descent step in the dual (moment) space.
3

---

Algorithm Generic LSVI (any family Q)
Require: E V, N 2
f1o Yoo while not converged do
X1...., ... qi
Fa. - N Ei
s ((X))
2/, - N D (Xx))((X))
16+1 t ti 24,
> ordinary least squares estimator (OLS)
Et +- stepsizel Fr 2 f/-11 X)
îh+1 < E 1 + f/
end while
Proposition 2.4 allows us to leverage known convergence results for MD under standard smoothness and convexity assumptions the ukL objective [in the VI literature, see, 12, 28, 29).
Assumption 2.5. The uKL objective w W + uKL(yw T) is L-smooth, justtrongly convex relative to Dz-.
Under Assumption 2.5, MD is known to converge with rate O(1/k) for sufficiently small and linearly decreasing step sizes: E, = (L + at) for S a < μs, [see, e.g., 30, 31, Theorem 4.5 and Lemma 4.8, Theorem 4]. The non-strongly convex case (μ+ = 0) exhibits O(1 /V//)) convergence rate for specific choice of step sizes [see, e.g., 30, Corollary 4.6]. In practice is trivial to set the to obtain a O(1/k) rate, as the relative strong convexity parameter if it exists, might be unknown and eventually very small.
Remark 2.6. The wrong-convexity/conexiig assumptions are rarely verified practice, however, such assumptions are standard for analysing convergence of optimisation algorithms (including NGD and MD), to ensure a unique minimiser and tractable rates [see, e.g., 32, Ch. 5]. While non-conjugate VI objectives may not be globally convex [20] (but it holds when the variational family contains the target), local convexity near optima often suffices for local convergence to hold.
See [33, 34] for provable smoothness guarantees on the KI. objective.
Practical algorithms and their analysis
The exact LSVI mapping assumed that one has access to expectations under the variational family. In practice, exact computation those expectations is intractable for general target log-density f. In this section, introduce practical algorithm in which these expectations estimated via Monte Carlo, and we study the impact of the Monte Carlo error on the convergence guarantees.
3.1 Generic LSVI
Our first algorithm comes down to replacing the two expectations in (3) with Monte Carlo estimates:
F, (X.)ss(X)T,
in N f(X))sX)),
where Xy...., XN 9n- The counterpart the exact iteration (5) is then
Lid.
i+ := fo
1-1 2n. + =
(10)
with to 7/n- At any iteration 1> 1, the step size can, in all generality, depend on the current state the algorithm via a function stepsize. We discuss one possible choice Section 3.2. This leads naturally to generic LSVI Algorithm 1, whose one-iteration cost is O(m3 + m²N).
Iteration (10) replaces the exact computation of F, with a Monte Carlo estimate --1 This approx- imation introduces bias in the estimation of the inverse FIM, and consequently, in the estimation of the natural gradient involved in (6). Further analysis of the statistical properties of the sequence

---

(ft). in particular, its convergence toward a neighbourhood of the optimum, requires careful con- trol of the bias. When admits uniformly bounded fourth-order moment, and the spectrum of Fy is bounded away from zero, the bias conditioned to high-probability event can be controlled.
Assumption 3.1. The sufficient statistic admits uniformly bounded fourth-ordes moments:
J44 == sup max (En [|s(x):11)) < 00, == sup wEW IS
(Ew [lu s(X)1))
sup
<0.
WEW I| =1,u3
(11)
Assumption 3.2. The smallest spectral value wnf |1 is strictly positive.
Both assumptions are verified if i) W is compact set, and ii) x(X) admits fourth-order moments, for ~ qu and for all w E W. While W is generally not compact set, it should not be considered as limiting assumption in practice, and can be lifted, see [35-37]. We further assume that f admits uniformly bounded second-ordes moment as this is required to control the norm of &w.
Assumption 3.3. m2 sup E [
We derive the convergence in expectation to the minimum of the KL loss conditioned on the event that the estimated FIMs are well-conditioned.
Theorem 3.4 (Explicit convergence rates for LSVI), Assume 2.1, 2.5, 3.1, 3.2, and 3.3. Let k 2 O, and let io. f1 ............ fis be given by (10), with بe (fj:) for 0 S t S k. Let A₁ nf A())))(a)))))))))) with A(w)= [|| Full < || ||=""'11'') Further assume that at each iteration
the quantities Fa. and 8n, are computed using two independent sets of samples. Let o. = -c-187-1 11
for t 2 Co 1, Cx Σ Ct-1- Let Wik c. EE=rct-ru ct-1 be the weighted average of the iterates, and let w be the minimiser of
Fix 8 E (0. 1), provided VN 2 2 + (Vlo(mm))) + 22 log(m)) for some constant Co 0, Ap happens with probability at least 1 6.
2. Conditioned Ak.
E(((a)) Ax] - 1(w") (5
11) uKL(qw- 1ww))
1 N
= 0
+
+0
Ck
N
C
(12)
where the big-O terms are independent of k.
3. Let Let L+ for some > 0. The RHS of (12) has asymptotic convergence rates that depend on a compared to the strong-convers parameter When a > p, the sequence (cx) is strictly decreasing, and the rate is O (k-p/aa + O (N-1). When sequence (c+) is constant, and the rate
p, the
(k)k- 1-1 +0 (N-1). When a < μ, the sequence (cx) is strictly increasing, and the rate is O (k-#/a) +0 (N-1N-1)+ O ON--)..
Remark 3.5. Our proof follows similar strategy to that of Hanzely and Richtárik [30], extending their mirror descent lemma to biased estimates. We control both the bias and the variance of the FIM estimate, conditionally this occurs with high probability when N is sufficiently large, using concentration inequalities for dositive-definite matrices [38].
the event that the estimated FIM is well-conditioned (A&)). We show
Remark 3.6. The convergence guarantees can be decomposed in three terms. The first term is due to initialisation and vanishes & + oo, the third term is the Monte Carlo bias and vanishes as N + 0o, and the second is a cross term and vanishes whenever k + oo or + oo.
Remark 3.7. The OLS estimate to the regression problem uses a single set of samples to compute both Fr. and 2q,, contrary to the estimate introduced in the previous theorem. Additionally, for many exponential families, closed-form expressions for F, are known. Since the OLS optimal with respect to the variance, it exhibits lower variance compared to others estimates. Importantly, it inefficient to use two distinct set of samples or to replace the estimated FIM with the exact FIM [22].
5

---

3.2 The choice of the Ey's
Setting & to a small enough and linearly decreasing sequence of step sizes ensures convergence of the sequence (10) to a neighbourhood of local minimizer 17" [30, 31J, see Theorem 3.4. However, the smoothness and strong-convexity parameters (4) of the KL objective, they exist, are rarely known in practice [33, 34]. For these reasons, choosing the can be tedious task, as in any stochastic optimisation scheme [16]: step sizes that are too large lead unstable behaviours while too small step sizes lead to slow convergence.
Let 77 E V and 17 = ((1)) be the OLS, consider the following linear regression objective,
f(X,) TT s(X,) + Vis X XN
Lid. On
(13)
where is is the residual of the regression. Then (13) implies that for any E (0.
ef,),) + (1-e)n) =XX)) (so(n) + (1 -))) (X.) + EUS.
(14)
The previous equation (14) shows that descending toward the direction of the OLS with step size multiplies the variance of the residuals v1,.... UN. x² by g² Let u² be some upper bound on the variance of the residuals, and let & S u/v, then the residuals have variance less than u². This remark, combined with a backtracking procedure parameters, yields an adaptive schedule for choosing the step sizes (Algorithm 4 in Appendix B), which we have found to be robust against noisy iterates and slow descents.
ensure that the iterates remain in the of natural
4 Gaussian families
The two most commonly-used families Q in variational inference are the full-covariance Gaussian family (Na(p, Σ) with arbitrary μ and E > 0), and the mean-field Gaussian family (E diagonal). A single iteration of LSVI requires inverting the Fisher information matrix (FIM) F, which is too expensive to be practical in high-dimension; i.e., Omm)), with m = ((d (resp. m = (dd)) in the mean-field (resp. full-covariance) case.
Attempts lessen the computational complexity of inference procedures over Gaussian distribu- tions either rely on access to cheap gradient estimates in the space of moments [13, 18, 20, 39], on single draw updates making the FIM estimate cheap invert but noisy [13, 40], or on restrictive assumptions on the target density [41, 42]. We derive closed-form formulae for the natural gra- dient descent iteration whose cost, in the full-covariance case, essentially amounts to the cost of computing d x d matrix products, that is (dd). In the mean-fied feen-field case, the cost is (d).
Full-covariance Gaussian family Let be the family of (unnormalised) Gaussian densities of dimension d. The sufficient statistic is s(x) = (1, x (vee (xax" )) T E IRm with m d(d-1) 1) 1, where vee (xx*) denotes the vector obtained by vertically stacking the columns of xx and we denote by unvec the inverse operation. Consider a natural parameter (r) (1),, (,,(2)T) e V
with (0) E IR, r7 I& and r7(2) E Idd", then it defines unique Gaussian distribution with mean and covariance matrix given by
1 (2) 21 - (- (1) (1) unv
(n, Σ) =
(15)
We reparameterise the linear regression of f(X) with respect to s(X). where X N(, E), into a regression of f(u + CZ) with respect to t(2), where Z ~ (0, and C = Chol(E)) is the Cholesky of Σ, and
t(2) = 1, T 2 -1
2/ -
2/ -
21
Zdd
22
(16)
V2
V2
V2
and
:= argmin. E [[>T((2)
s( CZ))
(17)
In brief, is one-to-one transformation such that the output vector has un-correlated components: ([)( (Z)) = T. That makes possible the estimation of 7 without inverting the FIM. The explicit mapping from ? 77 depending on (n, Σ) is given by the next theorem.

---

Theorem 4.1 (LSVI mapping ó for coll-covariance Gaussian distributions). Let E V defines Gaussian distribution X (H, E), and let Chol(E) be the Cholesky of L. Then, B o(n) is defined recursively from bottom to top by
3
(0)
Σ d
i=1 Tú, [(1).T 14 ((2).T
vec uu
3 =
B01
x-1,(11 (1)1-21 -011-111 2p ((2)
(18)
1833
vec (c-IrC-T)
where = [L((Z))(( + CZ) has subcomponents 7 (,(0), ,,), ,,, (T)TT) 7(0)) E R. 7() E RR,, E 12(d+1)/2 and where I is the symmetric malrix given componentwise by [
(2) - /2, Ti,++ = 7 -1-1122242--0( f S S d and S S d - In
addition, if has second-orde derivatives such that [|Ex Vfll < oc and -Ex V fl then
6(1)) defines a Gaussian distribution with mean and covariance given by
(( " (VXx)())) X ~~N(,E,
(19)
Theorem 4.1 gives the regressor with respect to s(X) of f(X), as a function of (p, C*) and y. Furthermore, all the involved operations have cost dominated by the computation of which the same as computing products of d x d matrices, (dd)).
Mean-field Gaussian family The family of mean-field Gaussian distributions treated similarly to the previous one by removing the cross-terms ž;zj in the sufficient statistic. The total cost of the OLS computation is (dd) See Appendix D.3 for the explicit regression procedure.
4.1 Stochastic schemes tailored to Gaussian distributions
We now take advantage of the reparametrisation tricks previously introduced to derive tailored im- plementations of LSVI for Gaussian variational families, with optimal one-iteration cost in d. An unbiased estimate of the OLS (17) is given by
N
Ltd. 3 = N (Z)f((+ CZJ,, Zi...., Zn ~ (O,
(20)
1=1
We define î as the estimate obtained by plugging & into (18) Theorem 4.1. The mean-field case treated in a similar manner. See Algorithms 2 and 3.
Algorithm 2 LSVI-MF (mean-field Gaussian family)
Algorithm 3 LSVI-FC (full-covariance Gaus- sian family)
Require: (No, og) on,i > 0, é E [1,4], N
Require: Mo, Eo 0, 2 (Ao, Eo) + (no, 2o)
(fo, 0ū) + (HO, ot)
no 6- (-x0, 1-x0--- - μ1 - vec 5-11
No + (-200 -p/00- 276
while not converged do
while not converged do 2i...., 10,,
Ć + Cholesky(2,, (Σ,)
Zi.... ZN (O.
941 + + 1 + G. & 2,)
(Z) Compute 1/411 using (43)
241 + 1 ΣN f(Z,f + C.Zi)
Compute 144 using (18)
Et f- stepsizeefpee...
f4 it, Z1.N)
E, + stessize (eriiizeeff11,1f R4+1 ir. Zi.N)
fi+ < e1f1+111 (2).-11 (1)
(1
Í + E + ft
f41 2//441 1/11
P t- - - 16+1 1/11
(2), (1)
0b_1 <- 21/111
(2).-1
Σ f- (unvec(fi+1]))
end while
end while
7

---

5 Numerical experiments
We consider three examples: one where SGD may be used minimise the KL objective, and two where it may not, because the reparameterisation trick is not possible: distributions in Q are discrete, log æ not differentiable, or because the log-derivative trick yields noisy estimates [15].
In the first example (logistic regression)), we compare all three LSVI' instances with other gradient- based KL minimisation procedures, including ADVI, NGD, and a gradient-free procedure for Gaus- sian mixtures. In the second and third examples (variable selection and Bayesian synthetic likeli- hood, BSL), since SGD is not available, assess the approximation error of LSVT relative the true posterior using exact Bayesian inference.
5.1 Logistic regression
Given data (I,, Us) E (-1,1, model is: r(8) x IT= F(yer) (3) where F(x) 1/(1 + e e=*) and p(3) is a (typically Gaussian) prior over the parameter B. This type of posterior is often close a Gaussian, and is popular benchmark in Bayesian computation [43]. See Appendix C.2 for summary of the considered datasets and the priors.
the posterior distribution of logistic regression
Whenever applicable, we compare LSVI (Algorithms 1, 2. 3) with NGD and ADVI. For NGD, the gradients are obtained via JAX autodifferentiation [44] and the FIM is estimated via Monte Carlo. For ADVI, we use the standard implementations given by PyMC3 [6] and Blackjax [7] with default step size schedules (that is, modification of Adam and RMSProp for pyMC, and comparable fixed step sizes for Blackjax). In addition, we provide a comparison of LSVI (Algorithm 1) in low dimension with the gradient-free iteration for Gaussian mixtures (GMMVI, [45]) which is a fair comparison since GMMVI and LSVI Algorithm 1 have the same complexity in this case. In addition, we illustrate the compatibility of our proposed methods with subsampling procedures for large datasets [46, 47] to reduce the cost of the log-likelihood evaluations.
Figure summarises this comparison for the Pima dataset (full-covariance case). One sees that LSVI (Algorithm I) converges essentially in one step, LSVI-FC (Algorithm 3) converges in less than 100 steps for linearly decreasing step sizes. For such a low-dimensional dataset (d 9). LSVI remains competitive with LSVI-FC since it converges faster, and the matrices it needs to invert are small. LSVI performs comparably to NGD and GMMVI, but is less noisy (with or without an adaptive schedule 4). We consider larger and more challenging datasets as recommended by [43]. Figure 2 (left) does the same comparison for the MNIST dataset (mean-field covariance), In Appendix C, Fig- ure 4 for the Sonar dataset (full-covariance) and Figure 7 for the Census-Income Income dataset (mean-field covariance with subsampling). This time, inverting the FIM too costly (e.g., 2015 X 2015 for Sonar), so we only use the tailored schemes LSVI-MF and LSVI-FC. Section C.2 contains extra details and results for all datasets in Table 2. including runtimes and memory usage (Table 1), aver- age cost time per iteration with respect to N (Figures 4 and 5), loss elapsed time and classification performance (Figure 6). details on the considered schedules for the step sizes (Table 3).
5.2 Variable selection
Given dataset D = (Iz; V1 Us):=i....... I, Rd Us R, the variable selection task in Bayesian linear regression may be modelled as y/ x. diag(y)8 + oE,, NO, 1). where E {0, 13dd is vector of inclusion variables, which is assigned prior distribution that is a product of Bernoullipp); e.g., P- = 1/2. If (3,0²) is assigned a conjugate prior, the marginal posterior distribution (>) (with 3, o² integrated out) admits a closed-form expression, the support of which is {0, 134. It is
therefore natural to set Q to the family of Bernoulli products, i.e. q(>) = II 9, (1 q:) --m with
q. E [0, 1] for i = 1...., d. This family is discrete, which precludes a reparametrisation trick, and the application of ADVI.
Figure 2 (right) compares the posterior inclusion probabilities, i.e. ক( = 1|D) approximated either through LSVI (Algorithm 1), or the Sequential Monte Carlo (SMC) sampler of Schäfer and Chopin [48], for the concrete dataset (d 92). This dataset challenging as it generates strong posterior correlations between the 7. Despite this, LSVI gives reasonable approximation of the
Python package: https: /ttps://ee hewgthhb..o/ylleftteeee
8

---

745
11
GMMVI (default)
720
NGD
375
ADVI (pyMC, defaul))
500
Zoom 5
374
475
373
GMMVI (default)
LSVIFC. sch. Alg.
368
LS, sch.
367
NGD
350
325
357 20
200
40
100
Iteration
Iteration
Figure Logistic regression, Pima data, full-covariance approximation. KL divergence (up to an unknown constant) between the variational approximation and the posterior, function of the number of iterations. Left: truncated from iteration 20 for better readability. Right: fo- cus on GMMVI, I.VI and NGD. Mean over 100 repetitions and one standard deviation interval (jax www.nummy..ttt.
(blackjax,
2000
NGD, scht
LSVI-MF, sch. 1 (with Alg.41
0.8
divergence
1500
LSVI-MF, sch. (with Alg.41
LSVI
2000
0.4
0.0
Elapsed time seconds
SMC
Figure 2: Left: Logistic regression, MNIST, mean-field approximation. KL divergence between the variational approximation and the posterior, as function of the elapsed time. Truncated from itera- tion 10. Mean over 100 repetitions and one standard deviation interval. Right: Variable selection example, posterior marginal probabilities (<; 1[D): LSVI vs SMC. (LSVI: 100 repetitions, the min-max intervals are reported with arrows, SMC: 3 repetitions).
true posterior. To the best of our knowledge, this is the first time variational inference is implemented for variable selection using the Bernoulii product family. See Section C.3 for extra numerical results and more details on the prior, the data, and the implementation.
5.3 Bayesian synthetic likelihood
BSL is a popular way perform likelihood-free inference, that is, inference on parametric model which is described only through simulator: one is able to sample Y ~ Po. but not to compute the likelihood p(y(0), see Frazier al. [49] for review.
BSL requires to specify syy), a low-dimensional summary of the data and assumes that syy) ~ N (b(O), E())), leading to posterior density =0) p(0) (s(y); b(0), I(O)), where p(0) is the prior. Since functions b and are unknown, they are replaced by empirical moments b(0). ε()), computed from simulated data. This makes BSL, and in particular its Markov Chain Monte Carlo (MCMC) implementations, particularly CPU-intensive, as the data simulator must be run many times. Furthermore, each evaluation of x is corrupted with noise, making it impossible to differ- entiate log म. Note that, in general, the data simulator is too complex implement some form of reparametrisation trick, the application of automatic differentiation procedures.

---

We consider the toad's displacement example from [50], which has been considered in various BSL papers [49, 51]. The model is parameterised by (a. 7.Po) E Rt x R x [o, 1]. See Section C.4 for more details on the model. We implement both LSVI-MF and LSVI-FC. For the former, we use family of truncated Gaussian distributions, while for the latter, we re-parametrise the model in terms of Ç f(), where is one-to-one transform between Q and IR The panel of Figure 3 shows that both LSVI algorithms converge quickly. The bottom panel shows that the full-covariance LSVI approximation matches the posterior obtained via MCMC, at fraction of the CPU cost, see Table 1. Again, we refer to Section C.4 for more details on the implementation of either LSVI or MCMC.
1.8
0.8
0.7
Scaled
0.4
.22
0.3
4/100
0.6
0
10
20
30
40
50
1.6
1.8
35
40
0.6
0.8
Iteration
Figure 3: Left: Variational approximations of each coordinate of 0 with one standard deviation interval, normalised. Truncated Gaussian: solid line. Full covariance Gaussian: dashed line. Right: Full-covariance Gaussian variational approximation (blue), MCMC approximation (orange).
6
Limitations
The current approach is limited to exponential families; mixture exponential families may be tackled by adapting the inrectationmmaxmmiaatiio approach of Arenz al. [45], or by building on existing applications of NGD VI methods to mixtures of exponential families [19, 22]. For Gaussian approximations, if the posterior contains directions that strongly non-Gaussian, then conditional- Gaussian strategies like integrated nested Laplace approximations may be applied [52]. In discrete exponential families, independence can be lifted by considering tree-structured dependencies, which are quite flexible, see, e.g., Wainwright and Jordan [21].
Acknowledgments and Disclosure of Funding
The first author gratefully acknowledges partial support from the Magnus Ehrnrooth foundation. The authors thank Sam Power, Mohammad Emtiyaz Khan and anonymous reviewers for insightful remarks on preliminary version.
References
[1] Michael I. Jordan, Zoubin Ghahramani, Tommi S. Jaakkola, and Lawrence K. Saul. An intro- duction to variational methods for graphical models. Machine Learning, 37(2):183-233, 1999. ISSN 1573-0565. doi: 10. 10.1023//4100766507777 URL https://doi. /00 1rg/10.1023/AA. 1023/A: 1007665907178.
[2] David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational inference: review for statisticians. Journal of the American Statistical Association, 112(518):59-877, 2017. ISSN
0992-149,,1537-244 doi: 10.1080/0162) 1099000061199.2217711211811 10. 10.1008/01621459 2017. 1285773.
URL
[3] Cheng Zhang, Judith Bütepage, Hedvig Kjellström, and Stephan Mandt. Advances in vari- ational inference. IEEE Transactions Patiern Analysis and Machine Intelligence, 41(8): 2008-2026, 2019. doi: 10.1109/PAAM1.20888889777 URL htt.hhttp://eee :/doi. mecii.rg//10 otg/10..1000 1109/ TP AMI. 2018. 2889774.
10

---

[4] Bob Carpenter, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Mareus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A prob- abilistic programming language. Journal of Stafistical Software, 76(1):1-32, 2017. doi: (0.18637//ss.v077..it URL https://wwwt jww.tssfttt..ttt.....oo.. org/ ptrg//ndex.ce pttp:/ppppsssscplplppspo view/v076101.
[5] Du Phan, Neeraj Pradhan, and Martin Jankowiak. Composable effects for flexible and acceler- ated probabilistic programming numpyro, 2019. URL https://arxiv http:///axvvvv....aaa.// org.oog///ao.. htt//bbs/19111/ 11554.
[6] John Salvatier, Thomas V. Wiecki, and Christopher Fonnesbeck. Probabilistic programming in python using PyMC3. Peerl Computer Science, 2:e55, 2016. doi: 10.7717/peerj-cs.55 URL https://doi org/10. 7717 htw.//perr-cs..000
[7] Alberto Cabezas, Adrien Corenflos, Junpeng Lao, Rémi Louf, Antoine Carnec, Kaustubh Chaudhari, Reuben Cohn-Gordon, Jeremie Coullon, Wei Deng, Sam Duffield, Gerardo Durán- Martín, Marcin Elantkowski, Dan Foreman-Mackey, Michele Gregori, Carlos Iguaran, Ravin Kumar, Martin Lysy, Kevin Murphy, Juan Camilo Orduz, Karm Patel, Xi Wang, and Rob Zinkov. Blackjax: Composable bayesian inference in jax, 2024, URL https: /ttp:://axvvvvio org/ abs/2402. 10797.
[8] Rajesh Ranganath, Sean Gerrish, and David Blei. Black Box Variational Inference. In Samuel Kaski and Jukka Corander, editors, Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, volume 33 of Proceedings of Machine Learning Research, pages 814-822, Reykjavik, Iceland, 22-25 Apr 2014. PMLR. URL html.
https://prrs cttps://proceedinggs mlr
[9] Alp Kucukelbir, Dustin Tran, Rajesh Ranganath, Andrew Gelman, and David M. Blei. Auto- matic differentiation variational inference. J. Mach. Learn. Res., 18(1):430-774, 2017. ISSN 1532-4435. URL http: //jmmr/mm////////// org.oogg//....ooo paapersv18//ii 107 html.
[10] Michalis Titsias and Miguel Lázaro-Gredine Gredilla. Doubly stochastic variational bayes for non- conjugate inference. Eric P. Xing and Tony Jebara, editors, Proceedings of the 31st In- ternational Conference on Machine Learning, volume 32 of Proceedings of Machine Learn- ing Research, pages 1971-1979, Bejing, China, 22-24 Jun 2014. PMLR. URL https: l/proceedings.mll mlr press/ https//////tttitiiieeeee html.
[11] Shun-ichi Amari. Natural gradient works efficiently in learning. Neural Computation, 10(2): 251-276, 02 1998. ISSN 0899-7667. doi: 10.1 N0.116228997699301177666 URL https: //doi. 10. (1112//00887766690000777777
[12] Mohammad Khan and Wu Lin.
Variational Inference : Converting
Variational Inference in Non- Non-Conjugatt Models to Inferences Conjugate Models. Aarti Singh and Jerry Zhu, editors, Proceedings the 20th International Conference Artificial In- telligence and Statistics, volume 54 of Proceedings of Machine Learning Research, pages 878- 887. PMLR, 20-22 Apr 2017. URL https: http://poceeen...cccom mlr haaass/vv4//aann html.
[13] Mohammad Emtiyaz Khan and Didrik Niclsen. Fast yet simple natural-grading descent for variational inference in complex models. In 2018 International Symposium Information Theory and Its Applications (ISITA), page 31-35. IEEE Press, 2018. doi: 10.2391 /9//TT.. 2018. 80866664.... URL https: org/00 10919/IIITA. 2018 8664326.
[14] Mohammad Emtiyaz Khan and Håvard Rue. The bayesian learning rule. Journal of Machine Learning Research, 34(281)::14 2023. URL http://jmmr... .or...org//e wwt.ppapers/v24/2220000 html.
[15] Ronald J. Williams. Simple statistical gradient-following algorithms for connectionist rein- forcement learning. Machine Learning, 8(3):229-256, May 1992. ISSN 1573-0565. doi: 10. 10.1007/BF0099266%. URL https://doi: org/ 1ttp:///////.1.///////00000000/0

---

[16] Manushi Welandawe, Michael Riis Andersen, Aki Vehtari, and Jonathan H. Huggins. A frame- work for improving the reliability of black-box variational inference. Journal of Machine Learning Research, 25(219):1-71, 2024. URL http://jjmml jtrr///mrooooooo .org wtrgpappers/25/2222 html.
[17] Thomas P. Minka. Divergence measures and message passing. Divergence measures and message passing, 2005. URL https://miat: inrae.... At/AATot//bbbloss/T-22055 pdf.
[18] Mohammad Khan, Didrik Nielsen, Voot Tangkaraut, Wu Lin, Yarin Gal, and Akash Srivas- tava. Fast and scalable Bayesian deep learning by nright-perturbation in Adam. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Ma- chine Learning. volume 80 of Proceedings of Machine Learning Research, pages 261 1-2620. PMIR, 10-15 Jul 2018. URL aattps://pppl https/////pooeees...... .ml........... press, /tte//80//aaaaaa 8a. html.
[19] Wu Lin, Mohammad Emtiyaz Khan, and Mark Schmidt. Fast and simple natural-gradient vari- ational inference with mixture of exponential-fan and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Ma- chine Learning. volume 97 of Proceedings of Machine Learning Research, pages 3992 4002. PMLR, 09-15 Jun 2019. URL https:://ppeee https://proceedings mlr prwpppressssspeeep press/v97/1in1900 html.
approximations. In Kamalika Chaudhuri
[20] Kaiwen Wu and Jacob R. Gardner. Understanding stochastic natural gradient variational infer- ence. Journal of Machine Learning Research, 2024. URL https://dl. acm. org/doi/100 doi/10. 5555/3692070. 3694258.
[21] Martin Wainwright and Michael I. Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning, 1(1-2):1-305, jan 2008. ISSN 1935-8237.
[22] Tim Salimans and David A. Knowles. Fixed-Fomm Variational Posterior Approximation through Stochastic Linear Regression. Bayesian Analysis, 8(4):837 882, 2013. doi:
10.1214/13-BAB BA858.. URL http:: stttpa://doi.. org/ 10. 1214/12111111 BA858.
[23] R. Penrose. On best approximate solutions of linear matrix equations. Mathematical Proceedings of the Cambridge Philosophical Sociery, 52(1):17-19, 1956. doi: 10.1017/ $03055044000039999, URL https://ddaa:: http:///d....r//.00/00000/660000
[24] Heinz H. Bauschke and Patrick L. Combettes. Fejér Monotonicity and Fixed Point Iterations, pages 91-109. Springer International Publishing, Cham, 2017. ISBN 978-3-319-48311-5, doi: 5.
UtL.httpss:/bbi.orrg 1tt.:///r.....00////------------
[25] Arkadij Semenovič Nemirovskij and David Borisovich Yudin. Problem Complexity and Method Efficiency in Optimization. Wiley, New York, 1983. Originally published in Russian in 1979.
[26) Haihao Lu, Robert M. Freund, and Yurii Nesterov. Relatively smooth convex optimization by first-order methods, and applications. SIAM Journal on Optimizationn 28(1):333-354 2018. doi: 10.1137/16M10995446 URL bttpps///del http:///ooooogg.../// htt////i...r....../111111/111111 100.0111//161119956000
[27] L.M. Bregman. The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming. USSR Computational Math- ematics and Mathematical Physics, 233)2200-117, 1967. ISSN 0041-5553. doi: 10.1016/ 0041-5553(67)90000--7 URL https:://ool oof//ooooor///00
[28] Mohammad Emtiyaz Khan, Reza Babanezhad, Wu Lin. Mark Schmidt, and Masashi Sugiyama. Faster stochastic variational inference using proximal-gradient methods with gen- eral divergence functions. In Proceedings of the Thirty-Second Conference on Uncertuinty in Artificial Intelligence, UAI'16, page 319-328, Arlington, Virginia, USA, 2016. AUAI Press. ISBN 9780996643115, 15. URL http::tt/:/s/e //ds acm. wtgg/doi/10..om 10..0..555553300 /0255/3 3020982.
12

---

[29] Mohammad Emtiyaz Khan, Pierre Baque, François Fleuret, and Pascal Fua. Kullback-leibler proximal variational inference. In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Gar- netr, editors, Advances in Neural Information Processing Systems, volume 28. Curran As- sociates, Inc., 2015. URL. https::/p/// http://aapppe.eeeeee .nips. cc/paper__aa wttpeeefflleepppppeeeaeee / 321 224466B444c666577777166666622222 www.:ra5552224-----a-aaeee Abstract. html.
[30] Filip Hanzely and Peter Richtárik, Fastest rates for stochastic mirror descent methods. Compu- tational Optimization and Applications. An International Journal, 79(3):717-766, 2021. ISSN 1926-6000, doi: 10. 1011007/11058---01-00224444 URL 00007/s10899021--1 021-00284-5.
[31] Pierre-Cyril Aubin-Frankowski, Anna Korba, and Flavien Léger. Mirror descent with rel- ative smoothness in measure spaces, with application to sinkhorn and EM. In Proceed- ings of the 36th International Conference on Neural Information Processing Systems, NIPS '22. Red Hook, NY, USA, 2024. Curran Associates Inc. ISBN 978-1-713-87108-8 URL /ttts://dl... acm. acm.org/ www.oogggo//ooooooooooooo /10//10..00 50..55555/36555 00000277036005555
[32] Francis Bach. Learning Theory from First Principles. Adaptive Computation and Ma- chine Learning. The MIT Press, Cambridge, Massachusetts London, England, 2024. ISBN 9780262049443. "The aim of this book is provide the simplest formulations that can be derived from first principles' with simple arguments".
[33] Justin Domke. Provable smoothness guarantees for black-box variational inference. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 2587- 2596. PMLR, 13-18 Jul 2020. URL https: http:////coccceee//ee mlr. war..ppes//vie /vt///v11//////// prress/v19/domkeeeeee /d0kee00... html.
[34] Justin Domke, Robert Gower, and Guillaume Garrigos. Provable convergence guar- antees for black-box variational inference. In Oh, T. Naumann, son, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Informa- tion Processing Systems, volume 36, pages 66289-66327. Curran Associates, Inc., 2023. URL httpps://proceeding..ttt mew.eeerddnnnnenenneeeeeee https//se..eeeree..pppppeeeeeeee sateerrfille/pppppeee httperrrilleeeappee/2/////////// /pppp//p////fflle/// (0lcf666Bbb555003355577776666 Paper- pdf.
Glober-
[35] Kevin Scaman and Cedric Malherbe. Robustness analysis of non-convex stochastic gradient descent using biased expectations. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 16377-16387. Curran Associates, Inc., 2020. URL paper_files
wewerise...pp.... cc/ Paper pdff.
[36] Kevin Scaman, Cedric Malherbe, and Ludovic Dos Santos. Convergence rates of non-convex stochastic gradient descent under generic lojasiewicz condition and local smoothness. In Ka- malika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 19310-19327. PMLR, 17-23 Jul 2022. URL https: httpa/////ccceeesss..om mlr. prw.s/ressseee vvtess/vv12////// scaman22a..... html.
[37] Bastien Batardière, Julien Chiquet, Joon Kwon, and Julien Stoehr. Importance sampling-based gradient method for dimension reduction in poisson log-normal model, 2024. URL https: //arxiv. org/ abs abs/2410. 2410.00476.. 00476.
[38] Richard Y. Chen, Alex Gittens, and Joel A. Tropp. The masked sample covariance estimator: an analysis using matrix concentration inequalities. Information and Inference: A Journal of the IMA, 1(1):2-20, 20, 05 2012. ISSN 2049-8764. doi: 10. 10.109//maiaa/aas001... URL https: //doi. 100 1093/ http://mmmaaiiaaceeo 1as001.
[39] Linda S L Tan. Analytic natural gradient updates for cholesky factor gaussian variational approximation. Journal the Royal Statistical Society Series B: Statistical Methodology, page gkaf001, 01 2025. ISSN 1369-7412, doi: 00.1093/jsssh//kaa0000. URL httph: //do:://dooooooo 10. 1093/ wtt.//////ssseee jrssssssssgeeee qkaf0011.
13

---

[40] D. Nguyen A. Godichon Godichon-Baggion Baggioni and M.-N. Tran. Natural gradient variational bayes with- out fisher matrix analytic calculation and its inversion. Journal of the American Staristi- cal Association, 000):1-12, 2024. doi: 10. 10.10800012214592222233333
URL. https:
//doi. www.:////...10010222221
2024. 2392904.
[41] David Barber and Christopher M. Bishop. Ensemble learning for multi-layer networks. In Michael I. Jordan, Michael J. Kearns, and Sara A. Solla, editors, Advances in Neural Infor- mation Processing Systems 10, [NIPS Conference, Denver, Colorado, USA, 1997], pages 395- 401. The MIT Press, 1997. URL https:: /papp////pppeaapee nips ppapper_tlleepppprr//SSa
[42] Manfred Opper and Cédric Archambeau. The variational Gaussian approximation revisited. Neural Computation, 21(3):786-792, 2009. ISSN 0899-7667,,5330-88888 doi: 10.1 16.1162/neco. URL https://doi. ort/ 11e1//1..111//1....e 2008 08-07-592.
[43] Nicolas Chopin and James Ridgway. Leave Pima Indians Alone: Binary Regression as Benchmark for Bayesian Computation. Statistical Science, 32(1):64 87, 2017. doi: 10.1214/
16-STS581. URL https:://oi.orreeeeeee 11111.1211/1--SSSS0
[44] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Minn and Qiao Zhang. JAX: composable transformations of Python+NummPy programs, 2018. URL http: // //github... com/jax jax-m/jjeeeeeeeee ml/jax.
[45] Oleg Arenz, Gerhard Neumann, and Mingjun Zhong. Efficient gradient-free variational in- ference using policy search. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 234-243. PMLR, 10-15 Jul 2018. html.
URL https: htttttt//ppoceeeeeeee /ttp//eeecceee....
mlr. prress/v00/aree
[46] Matthew Hoffman, David M. Blei, and Francis Bach. Online learning for latent dirich- let allocation. In Proceedings of the 24th International Conference on Neural Information Processing Systems- Volume 1, NIPS' 10, page 856-864, Red Hook, NY, USA, 2010. Curran Associates Inc. URL https: /dl. acm. org/doi org//oi/100 /10. 0655/2997189 2997285.
[47] Matthew D. Hoffman, David M. Blei, Chong Wang, and John Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14(40):1303-134) 2013. URL http: //jmlr.org hwtppp///////aaaaeeaaeappppaaeee
html.
[48] Christian Schäfer and Nicolas Chopin. spaces. Stat. Comput.. 23(2):163-884, 2013. ISSN 0960-3174,1573-13777 doi: 10.1007/ s1 s11222-011-9299-2 URL https://tti..ooo
Sequential Monte Carlo large binary sampling
1007/ http:/111122221 011 9299-2.
[49] David T. Frazier, David J. Nott, Christopher Drovandi, and Robert Kohn. Bayesian inference using synthetic likelihood: Asymptotics and adjustments. Journal of the American Statisti- cal Association, 1186544):2221-28822 2023. doi: 10. 00018800161115522222006611222 URI.
btttss://oi.orgg// 10t0//0001180//121155555 2022. 2086132..
[50] Philippe Marchand, Morgan Boenke, and David M. Green. A stochastic movement model reproduces patterns site fidelity and long-distance dispersal in a population of fowler's ss toads (anaxyrus fowleri). Ecological Modelling, 360:63 69, 2017. ISSN 0304-3800. doi: 10.101 Ecolmodel.2017.066.00 URL https: //d::////oo// org/10. 1016/ (0016//eecoomoeel... www.model. 2017 025.
[51] Ziwen An, David Nott, and Christopher Drovandi. Robust Bayesian synthetic likelihood via a semi-parametric approach. Statistics and Computing, 30(3):543-557, 2020. ISSN 0960- 3174,1573-1375. doi: 10.1007/s11222 12207//|11222--19999944 011-222-11--9904-2
URL
[52] Håvard Rue, Sara Martino, and Nicolas Chopin. Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations. Journal of the Royal Statistical Society: Series (Sattistical Methodology). 71(2):319-392, 2009, doi: 10.1111/ 1467-9868. 2008 008.00000.x.
14

---

[53] Nicolas Chopin and Omiros Papaspilliopollos SMC Samplers. Springer International Publish- ing. Cham, 2020. ISBN 978-3-030-47845-2.
[54] I.-C. Yeh. Modeling of strength of pigh-erooooee concrete using artificial neural networks. Cement and Concrete Research, 28( 28(12):1797-1808, 1998. ISSN 0008-8844 doi: 10. 10.1016/ S0008- 00008--884699800065-33 3. URL http::://bb..org//...... 1001110001080808 (844((88) 00165-3.
[55] Edward George and Robert E. McCulloch. Approaches for Bayesian variable selection. Statistica Sinica, 7(2):339-3777 1997. ISSN 1017-0405.
[56] David I. Warton. Penalized normal likelihood and ridge regularization of correlation and co- variance matrices. Journal of the American Statistical Association, 103(481))340-39,, 2008. ISSN 01621459. URL ptttp://wwer wwwp:wwww..oroooooo.oo.ooo.....w jstor. . org/ atg/staable/2764044.
[57] Tropp Joel Aaron. An introduction matrix concentration inequalities. Foundations and Trendy in Machine Learning. 8(1-2):1-230, 2015. ISSN 1935-8237. URL https://doi.
org/10. /561/220000048..
[58] Léon Bottou, Frank E. Curtis, and Jorge Nocedal. Optimization methods for large-scale machine learning. SIAM Review, 60(2):23-311, 2018. doi: 10.1137/16M1080173. URL https: http::////oooor///.... 11011111137/1
15

---

Contents
For any definite positive matrix Σ, denote by C Chol(E) the unique lower triangular matrix such that CCT Σ. A Notations
16
B Adaptive schedule algorithm
17
C Extra details on numerical experiments
17
C.1 Runtime analysis
17
C.2 Logistic regression
17
C,3 Variable selection
19
C.4 BSL and toads displacement model
20
D Proofs
22
D.I First order condition and critical points of the uKL objective
22
D.2 The exact LSVI is natural gradient descent
23
D.3 The exact LSVI mapping for mean-field Gaussian distributions
23
D.4 The exact LSVI mapping for Gaussian distributions
25
D.5 Concentration bounds for the Fisher matrix in the compact case
27
D.6 Convergence analysis of the stochastic L.VI algorithm
28
A
Notations
For any vector u RP, we denote by u-1 E R² the componentwise inverse of u, and u & u' the componentwise product
with u'. For any matrix U E RPXq, we denote by vec(U) the X 9 vector obtained by vertically stacking the columns of U, and by unvec the inverse operation satisfying unvec( vec (U)) U, while U 6) U" denotes the Kronecker product between matrices. For any square matrix U E TPPx, let diag(U) (U) be the p vector composed of the diagonal components of U and let |U]]| be the spectral norm of U.
For any set A, U(A) denotes the uniform distribution over A. N Σ) denotes the Gaussian distri- bution with mean z and covariance matrix Σ, and (N,, with 0² = (of...., o.) denotes the Gaussian distribution with mean j and diagonal covariance matrix diag(
The O is the usual big-O notation, i.c., An = O(B,) for some sequences A., B,, let it be reals, vectors or matrices, if there exists a constant C' such that for N large enough and all n > N. ||A.]]l < C' ||. If A.]] > C. we write A,, I(B.). We write A. Op(1) for sequence of random variables (A,) such that, for any A( > B) < & for n large enough.
there exists constant B such that
16

---

B Adaptive schedule algorithm
Algorithm 4 Variance control and backtracking strategy
Require: 00 v/ E Rd, 1, XI...., XN iid
On u
1: +-
2: while en/ + (1 E)m Vo
3:
g 8/2
4: end while
5: En/ + (1 =)m
6: m - N N f(X;) -17 -
7: 6 - (f(X:) - -m)
8: if then
9: E min(e, u/6)
10: end if
11: return
C Extra details on numerical experiments
C.1 Runtime analysis
All the experiments were conducted using Python 3.13, jax 0.5 with GPU support, Cuda 12.5, and using float64. The hardware specifications CPU AMD EPYC 7702 64-Core Processor and GPU NVIDIA A100-PCIE-40GB except for SONAR, Census and MNIST datasets where EPYC 7713 and NVIDIA A100-PCIE-80GY were used. See Table
Table 1: For all conducted experiments, runtimes and max memory usage, across 5 repetitions. T is the number of iterations and N the number of samples whenever applicable. LR Logistic regression, BSL= = Bayesian Synthetic likelihood, MF Gaussian = mean-field Gaussian, Gaussian full-covariance Gaussian.
Experiement
Runtime (seconds)
resident (uigabvies)
usage)
Caussian, Alg. I, UN.T) = (100.50) (AX) BSI. Truncited MF Goussian, Alg. (,000, (AX)
(299
71.5
77.5
1,07
137 (++.6)
137.5
1.38.7
1.06
BSL MCMC, Blockjax (IAX)
268.1 (+33.4
266.5
274.3
1.16
Variable Selection, Alg
15)
60.3
Variable Selection,
(+1.7
284.1
298
0.45
LR Gmassian, PIMA, Alg. 1 9ch. 3. (IAX)
(6(())4)
1.0
4.1
1.28
A NGD, (10', (AX)
(2:35))))
1.51
4.9
3.29
Alg.
(10", 100), (AAX)
(.377)
7.3
1.28
A Alz. 3 sch. 2. (10", 100), (IAX)
3.101.. 1)
3.5
3.6
1.28
# JPyMC (f' 104)
(:-03))
5.5
6.4
0.84
GMMVIL
(TenseFFlow)
LR Ganssion, SOPAK,, Alg. 3cch, (10) 100, (AX)
5.1 (20.3)
5.l1
5.8
3.11
#t Alg. sch. (10?, 100), (IAX)
5.1 (+1.3)
4.4
7.4
3.11
(10')
9.6(11.))
LR Caussaan, MNNS, Alg. scch. (10), 10), (AX)
19.1 (2:10))
15.7
21.5
2,300
# Alg. 2.sch. 5, (10". 500), (IAX)
(0+4(+
10.3
10.6
2.36
(10' 500), (IAX)
(8.9)+
22.1
18.5
2.36
Blickjax ATMV,
(10', (10'
(+..0+
18.8
22.2
3.34
A HNG sM.), 5, (10", 500) (AXX)
25.3(+2.2)
29.2
24.1
4.38
LR MP Caussaan, 100)
Alg
(AX)
(2.00002)
12.5
13.3
3.30
(1.005)
# NGD (MF), 5, (10*, 10', 10") (AX))
(+1.77
69.0
73.1
3.33
C.2 Logistic regression
Data The Sonar (CC BY 4.0 License) and the Census Income (CC BY 4.0 License) datasets are available in the UCI repository while the Pima dataset (CCO: Public Domain License) in the example datasets of Python package particles (License MIT v0.4, [53, Ch. 1]) and MNIST (CC BY- 3.0 License) is available https: https: //tp:t//////b.t.// httm://ppeedddeeettsesstt following standard [e.g., 43] pre-processing strategy for Pima, Sonar and Census-Income datasets: we add an intercept, and we rescale the covariates so that non-binary predictors are centred with standard deviation 0.5, and the binary predictors centred 0 and range 1. For the third dataset
We use
17

---

(MNIST dataset), we restrict ourselves to the binary classification problem by selecting pictures labelled 0 or 8. The gray-scale features which range between 0 and 255 are normalised to between 0 and 1. No intercept is added. For the Census Income dataset, the categorical variables are mapped using one-hot encoding.
Table 2:
Logistic regression example: summary of datasets and approximation families, in paren-
theses the batch-size
Dataset
Gaussian family
d
72
Pima
full-covariance
768
Sonar
full-covariance
62
128
Census (subsampling)
mean-field
48 49 000 (1000)
MNIST
mean-held
784
11,774
Prior For all datasets except MNIST, the prior ((3) is zero-mean Gaussian distribution with diagonal covariance matrix, and the covariances are set to 25 for all the other covariates, except for the intercept, for which is set to 400. For the MNIST dataset, the prior Gaussian distribution with zero-mean and covariance matrix 25I,.
Initialisations, schedules and number of samples The initialisation distributions for all datasets except MNIST are standard normal distributions. The initialisation for the MNIST dataset is e( The learning schedules (Ex) are obtained via Algorithm 4 with specific inputs (u² st) summarised in Table 3 along with the number of samples N.
Table 3: Logistic regression setup. Left: Inputs to Algorithm 4 by dataset. Right: Schedule index
reference.
Dataset
Algeritim
Scheclule input (M".
Samples N
Schedule input (u², Ee)
Pima
Alg.
(0,
(10, 1)
Pina
NGD
(0c,
1/10
10a
Pima
Alg.
[10. 1).
Do, 17(2+
10"
2
(o0, 1/(t + 1))
Alg.
1),
(2c,, 1/(2
(20, 1)
MNIST
Alg.
(10,1), (cc, 10 (10,10
10a
(1,
MNIST
Blackjax (mnannicedd vi). (MF)
(20,
10
10'
5
[00, 10-3)
[10)
"). [oc,
(10, 10-3) 3)
Census
NGD
(00, 10
10
246
LSVI-FC (sch.. 4)
per iteration in seconds
1.4
LSVI FC Esch. 2, without Alp
ADVI (pyMC3, default indlmmentation)
2.0
ADVI (pyMC, cetault sch.
Alg.
0.6-
0.2
1+2 163
Se3
16t
Number samples
Iteration
Figure 4: Logistic regression posterior, Sonar data, full-covariance approximation, LSVI-FC and ADVI implementations. Left: average cost per iteration seconds as a function of the number of samples N, mean over 5 repetitions with 2 std interval. Right: KL divergence (up to unknown constant) between current Gaussian variational approximation and the posterior, as a function of 1, mean over 100 repetitions with standard deviation interval.
MNIST The PyMC3 (License Apache 2.0 v. 5.22, [6]) implementation fails in this context, and we resort to the stochastic gradient descent (SGD) implementation Blackjax (License Apache 2.0 1.2.55, [7]) of the mean-field ADVI Algorithm. For SGD, set the learning rate to 0.001 and the
18

---

number of samples for the Monte Carlo gradient estimates to 10 See Figure 5 for the average cost per iteration in seconds, and the same plot as Figure 2 with respect to elapsed time.
ADVI, sch, (Blcckjax, meanfeld vil NSD, scth
sch. (Blackja))
per iteration seconds
2000
NGD, sch.
LSVI4F, sch. (with Alg.41
LSVIME sch, 5. Alg.
0.30
ISVIMF, (with 41g:41
1500
LSVI-NI, sch. 6. with Alg.
0.20-
1000
0.10
0.05
0.00
163 Tea
5ed
105
Number of samples
Iteration
Figure 5: Logistic regression posterior, MNIST data, diagonal covariance approximation, LSVI- MF, NGD (JAX) and Blackjax (meanfield_y _i) implementations. Left: average cost per iteration in seconds as a function the number of samples N, mean over 5 repetitions with 2 std interval. Right: KL divergence (up to unknown constant) between current Gaussian variational approximation and the posterior, as function of I, mean over 100 repetitions with one standard deviation interval.
In addition, we provide missclassification rate for the logistic regression model using the mean (of the Gaussian approximation) as the regression parameter, see Figure 6.
Blacojan (msantistd vi (ich, 2. #g
Blacqje imeanfisdd vi sch. with 41 LSvi wibuut
6 Aly
w.whh
10-1
Ves Peschustcccttoo rate
10'
Iteucios
Figure 6: Logistic regression posterior, MNIST data, diagonal covariance approximation, LSVI- MF and Blackjax (meanfield_gg __i implementations. Top: Misclassification rate as function of the iterations, mean over 100 repetitions with standard deviation. Bottom: same in log-log axis.
Subsampling (Census dataset) placement from the dataset:
At each iteration t, new batch is sampled uniformly with
f(3) log ((3) logp(5) Σ log plyu, Izu, 3) + Σ log p(xu,),
(21)
where U1,..., Up u(1,..., n). A new batch is drawn at each iteration. The batch size is P 104. We also use f for evaluating the KL loss. See Figure
C.3 Variable selection
Dataset The Concrete Compressive Strength dataset [54] is made of 1030 observations and initial predictors denoted by C, W, CA. FA, BLAST, FASH, PLAST, and A. We enrich the dataset by adding predictors computed from the existing predictors. new predictors, LG_C, LG LG_CA, LG_FA, LG_A where LG X stands for the logarithm of the corresponding feature X. The cross- product of the predictors is also added, resulting in new predictors. Finally, we add an intercept. The total of possible predictors is d = 92.
19

---

1200
NGD, Sch.
LSVI-MF,
LSVI-ME, sch. (with 4)
divergence
400
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
Epachs (1000/49000
Figure 7: Logistic regression posterior. KI. loss for the Census-Income dataset (mean-field, with subsampling), mean over 100 repetitions with 1 standard deviation.
Prior
The hierarchical prior on 8, 0,02,7 is given by
TB (),,))= N o2 diag(?)) <(σ²) InvGamma(u/2, Aw/ T(7) 1((0,
We follow the recommendations of [55] by setting the hyperparameters w 4.0, of and 22 = 10/X, where of is the variance estimate of the residuals for the saturated linear model = (I)
Close-form expression for (>) For model y E {0, 10 let Z, /Zi]k/7:=1 be the selected covariates and let by Z, y. Consider the Cholesky decomposition C,,,C Z +o I and define the least squares estimate for the residuals based on the model given by 7, 02,,0 (y (C) ( (,), Then, the log-posterior for up the log-partition constant given by
1r
(T.) - I|rll log(v) + log( (w// o² log c,
log (r D)
SMC, extra numerical results As benchmark, we compute the posterior marginal probabilities of inclusion using waste-free variant of the tempering SMC algorithm of [48] with chain length P = 10 and N 10 particles.
Given any probability vector p E [0, 1jd, plot the histogram of the variable (1)(m"())//() p)) with y g( p) (Bernoulli product). The pendant for the SMC discrete measure is obtained by replacing q with the SMC empirical measure ** In Figure 8 plot the histograms when y is distributed according to the SMC empirical distribution 4" and when is distributed according to three different mean-field Bernoulli distributions 7 q(- p): i) p 21....... 1), i... the intercept is always included and the other coordinates has 0.5 probability to be included, ii) the L.SVI estimates, and iii) the marginal posterior probabilities estimated via SMC.
C.4 BSL and toads displacement model
Model The model assumes that M toads move along one-dimensional axis during D days. For any day
< D, the toad labelled by < i < M, has observed position yi,r- During the night of day L 7 1, the toad moves according to an overnight displacement, dyj, which is assumed to be a Lévy-alpha stable distribution with stability parameter a and scale parameter With probability Po. the toad takes refuge at Vi.t dy,,e- With probability Po. the toad moves back one the
20

---

7-0
0.4
7m 1 (aa
7 7~0( (011)
0.3
Density
0.0
-5600
-5750
-5700
-5630
-5800
-5550
Scores log(s"(>)/g(x)))
Figure 8: Variable selection example: distribution of scores log a (1)//(()) when y ~ #*, when y ~ 9 = q(( p(t) with p(t) given either by i), ii) or iii)).
previously explored sites you with e' chosen uniformly in I,...,t. Finally, for any day 1 D the observed position is
11,++1 (y,, ōyr,e) (1 B,,
(22)
with Bi Ber(po), 6' u(, and dyi, Lévy-alpha(a,, 8), all variables being mutually independent. The initial position yt is set to dy; ~ Lévy-alpha(a, 6). The model is parametrised by (a, 6, Po) [1,2] 100) [0, 0.9] O. Simulating from the previous model yields the observed data Y = (L
Summary statistic The summary statistic is the concatenation of 4 sets of statistics of size 12 resulting in total statistic of dimension 48. Each subset is computed from the displacement in- formation of lag for / E {1, 4, 8), denoted by Yi = (1 Vi.++1 (.l==))is:sDDLLIsss displacement from É day t + 1 of the toad y((id)) |yK. Vi,c+1| is less than 10. it assumed the toad has not moved. The first statistic is the number of pairs (i, () such that y(i.) < 10. We then compute the median displacement and the log difference between adjacent p-quantiles with p 0, 0.1........ for all the displacements greater than 10.
If the
Truncated Gaussian distributions approximation The dataset Y is generated with (M, = (66, 63) and underlying 6* (1.7,36 35, 0.6). The mean and covariance estimates are obtained with P= 100 samples for each evaluation of the synthetic likelihood. We follow the methodology of [49] and use [56] shrinkage covariance estimate given by Σ = D2/2() D where is the estimated correlation matrices and = 0.5 is the regularization parameter. The prior distribution the uniform distribution over O The variational family is the set of truncated Gaussian distributions over θ with diagonal covariances. The initial distribution has mean μ = (1.5, 50, 0.5) and diagonal covariances d² (0.05, 10, 0.01). We run Algorithm with N 100 samples and T 50 iterations, the step sizes are obtained by Algorithm with u = and linearly decreasing step sizes.
Full-covariance Gaussian distributions on transformed parameters To constrain the param- eters 8, we perform inference on the transformed parameters g(0) = logit(g; (#:)) with 9: (6,) (0, - aa)//b, with a; such that gi scales 0, to [o, 1]. The prior distribution on unconstrained parameters 0' is 1 g -1(0)) (('). The variational family is the set of full-covariance
Gaussian distributions. The initial distribution for 0' has mean n² = (0, O, 0) and covariance
matrix E' = diag(0.1, 0.1, 0.1). The benchmark is obtained via MCMC with random walk step
21

---

(0,00.133) the acceptance rate over the chain of length 10 is roughly 31%, excluding the first 10 states.
D Proofs
D.I
First order condition and critical points of the ukL objective
Proof of Proposition 2.3. Injecting we obtain
exp (f) and yn exp (n 8) into the objective function,
mz)= = /n1 s- - f) + --- exp((1 s).
uKL(gn π)
(23)
Using (23), we have
V, uKL. - /ssss 9n sfq--
(24)
Writing the first-order optimality condition for the following minimisation problem
rf" E argmin uKL((gn An (T)
(25)
and applying (24), then dividing by Zn < co, yield Es, [1] Ey [ss. s ( be some fixed statistic with first component 1. Assume that = (1)²), .7 ) E V is critical point, i.e., v (KL( π) = O. We have
=)= f sqn -f fan
000 uKL(qn
-
(26)
Injecting nTx s=== 10) + rj á into (26), setting V yy uKL(an () = 0 and normalising by Zn yields ( Eny [$ - 17 $].
(27)
from the definition of the KI. divergence, we deduce
()) KL F) + log
<(x)\/ /
exp (TT
(28)
x
Proof of Proposition 2.2. We have
f) KL π)) ZTj) / cn'((n f) log Un log Zx.
Computing the gradient of the KL requires computing the gradients of Zn, log Zuy and Z We have VnZin J Te s, and V, log(Zz) Z VvZ Es Similarly, Vn
Zn 2 VnZY Zn = Z (y Then, using the previous equalities, the gradient of the KL is
v KL(Āg x) V Z c "T 3 f)
V log log(Zg)
NØZØ / e9 " 3 + - (for " 5) V, log( (Zp)
(29)
- s| (7 5 fl + 1, |s π ff E, [s]
= E, [1] T - fl + E, 555 T ザ 3f
Now, let compute the gradient of the uKL objective with respect to 17 = (17')),1) TT. Using
ss = 3 33"
(30)
to expand (24) yields
1y uKL(gn x²) = // / (0) + ST 1)) + f fan:
(31)
Vn uKL(gn | x) /3(())))))
(5) + $3 m) - 55qn.
22

---

Assume that v,, (4 T) = O, then from (31), we obtain E, (7)(0) + 58 n)
[[ = 0.
Reinjecting the previous inequality into the gradient of the KL (29) yields
V ( F) = -,,,, v - fl - E, (3)
(32)
Injecting the expression for (7))) (27) into (32) yields V KL .) = O. Conversely, the previous
computations show that if V, KL л) 0 and 0, uKT.(4y x)) = 0, then V,, uKT.(g, .))
0.
D.2 The exact LSVI is natural gradient descent
Proof of Proposition 2.4. Assumption 2.1 ensures that for any 7) E V, F. is invertible (minimality assumption), and ensures the differentiability of all the involved functions (regularity). us denote by v, the Jacobian of 1j > ((7). Let V. using (24), we have
Vi, Z,,(F, zn);
(33)
where Z, (gy, is the normalisation constant of 9t- Let (nt) be the sequence obtained via natural gradient descent given by (6). Then, by (6) and (33), we have
Me+1 1t:
Zm "
F
V (Tte
1)t Ee Fr n
( (F 1)| 2m
(34)
= (1 - Ex) + E F 2n1
= (1 - Ex) M) + c,p(N))
Thus, the LSVI iteration with learning schedule (ex) given by (5) is the natural gradient descent (m)) with learning schedule (Ez/Z,(9) given by (6). Let now prove (7). We have
Z
/
88 On dp
(35)
2g For
By the chain rule and (35), the Jacobian of 17 + ((11) is
V V Vnw Val
(36)
= Za Fy x Val.
Finally, injecting (36) into (6) yields
M++1 7h - -, V
(37)
which is (7). This shows the first equivalence. Let wo E W, and define (wx) as given by (8). The first order condition on the minimisation problem (8) yields
V (W+(1) V.. - e Val(wr),
(38)
but (w*(w,) r7(we) The thus (38) is exactly (7).
D.3 The exact LSVI mapping for mean-field Gaussian distributions
Let s(x) =
where x (x (#1...., .a) and x (zf, .e......2))) 22). The set of admissi- ble natural parameters given by V R x Rd (R- (00))< C RT.
2d + Let (n) (2)T n (2))TT] V. The natural mapping from 7 to (2, is given by T(r) = " (- (2)-1 - (2)-1) where is the componentwise product and 1,2))-1 is the compo- nentwise inverse r)
Lemma D.1 (Reparametrissaioon of the regression in the mean-field case). Let N(H, o2) be mean-field Gaussian distribution with μl, Rdd and a; for all i E (1.....d). Let nj (n))) (1)). (n()) (7) R, ()) Rd. (2 E Rdd be the natural parameter associated with X for the statistic x E IRd (1, x 2) TR. 1R1+2d, Let be given by (42). For any E Rdd let
23

---

a(2) = = if (., then x(Z) ~ N(u, o²). Letr = ,, (T) E p2d+1
7(0) E E.. 700 E R ((2) E IR' be defined componentwise by
(1 + ()) p + ( + tz)) (2)) 2 o + (7) (2) @ uu@0
(39)
2m) x 0²
Then, for any 2 E Rd
7
(40)
Proof. Let identify y such that (40) is satisfied. Suppose that for all 2 Rdd, we have (40), then
"f s(x) - 11 + n (2), x + +) T2)
(0) n (1)T p 17 (or X 2) 1/2 (w))s
+ 2n (2) イ + 1) (d)
= n 11)T (7) To + n 22) (n3))s
lerms group
temm in group
term in group
(41)
21) (2) o
T
term in group 2
term in group
+ (n
8 1)s
term group
By identifying the factors in front of 1 (group 1), the 2j (group 2), and the 2f's (group 3), we obtain (39). By injecting (39) into (40), the equality is satisfied.
Li
Theorem D.2 (LSVI mapping & for the mean-fiee field Gaussian distributions). Let X ~ NH, and n E V be the corresponding natural parameter and let t be given by
t(2) 1,
27 -
(42)
V2
V2
Then, the LSVI mapping B = ((n) is desined recursively bottom top by
3())TT μ -- (3-()))TT (4² 102+0% +02)
B =
2 o -1 229(2) @ p
(43)
7(2 8 (V202)
and = E E((Z)f(u+o + & 2)}, with subcomponents = (1),, 7(2) E Rdd In addition, if J admits second-order derivatives such that Ex < cc, [VIl < 0o, and < Ex [Diag(V?f)]. then <in)) defines a Gaussian distribution with mean and variance given by
700) E R.
(n' (u',
(1= (E V ((xi)) & E [VX()) (R [diag( 4 ((X)))=
(44)
Proof. We know that (nn) realises the minimum of the OLS objective (4), i.e.,
(nn) argmin Ex~N(u00) (3 s(X) f(x))²]
(45)
Using Lemma D.I, we can rewrite the regression objective with covariates given by into a regres- sion with covariates given by t. Using the notations of Lemma D.I. we let 7 be given such that y t(2) 3 s(x(2)) for all Rd., and where 3 <(1) is the unique minimizer of the OLS objective (45). Then,
7 = argmin Ez [()) flu 8))
= (E2 [tt T (2))) -1 Ez ( Z)
(46)
Ez ((Z))(H o 2)]
24

---

since Ez (tt(()) = Im- Inverting the relation (39) given by Lemma D.1 between and 8, which possible since all the og 's are strictly positive, we obtain
(0)
3(D, ( (n² + x²))
② oo-oa 2g(22 @
(47)
(2) (v20²)
But 8 - ((n), this proves the first statement (43) of Theorem D.2. For the second statement, assume that f admits second-order derivatives. Using Stein's 's Lemma and (46), we obtain,
? = Ez-N(o, [t()) (( to Z)]
Ez (f(u + & 2))
(48)
Ez 2 (u to Z)
(02 2 E diag
V2 (u 2))
Injecting (48) into (43), we obtain for
<(1)) Ex Ex(V)) - y E diag V (X))
<(n) -0
(49)
2 (Ediag(V2f)
Using the natural mapping T(n) =
(- (-
1 ()))) we obtain (44).
D.4 The exact LSVI mapping for Gaussian distributions
Lemma D.3 (Reparametrisation of the regression in the full-covariance case). Let X ( be Gaussian distribution with E Rdd and Σ > O. Let 1j E V be the natural purameter associated E Rd + (1, (vee XX) Rl+d4d Let t be given by (16).
with X for the statistic
For any E Rd, let æzz) p + Cz with C' E Raxd such that CC = =. If Z N Ia), then xzZ) ~ (H, E), and for any Rd
7 =T=())=))) s(x(z)), (()()))
(50)
with 7 (7(0), 7 , TT E BB+d+ddd41)22 Furthermore, the components of 7, (000 E R,
y(L) Rdd, 7(2) Rdd(+1)/2 given by
( + (1)T n 17 12) voe voceupT + +E Pes
7
CT (n)a)) + 2 2 C) 172)) 7(2)
(51)
where
== unvec (C08))T1())
(52)
and
(2) =(V2T1,1,2
21 v 2 21 2.3,....... 2.d. v2 a.
(53)
Proof. The proof is similar to the proof of Lemma D.1. Let us rewrite the regression with respect to Z. Let = (n)(),, (0)),,())))))))))))) ().T, TT())TTT= E R 1ddddd with n7 e R, n7()) Rd, 17() E Rd Let
X = + CZ with C such that CCT = E. Rewriting the linear regression on s(X) with s(Z), we have
CZ + T s(x))= (0) ( p + 17 (TTCZ
(2)T vec 141
+ ()) vee uz CT (2) vec CZu
(54)
+ 2. vee CZZ TT
= (Z),
where 3 (4(0), &(),T, TT Il are left to be identified. By identifying the quadratic
terms in (54), we have for (2)
m (2). vec CzZ T T ()) (2 & C) vec ZZT =((2).T vec ZZ
(55)
25

---

where we used vec ABC = (CT @ A) vec B. Thus,
n) - (C @ C) Tg()) -CC- @ (-
(56)
where we used (A B) Aー @ B-1, For (1), expanding the linear term in (54), we have
11)) Cz2 + +1 12) vec Czu + +1(2) vec u TC = =n C2 + 21) vec Czp
= -nnOx C 2) (1 2 ())
(57)
=
i.e.
300) 10 + 2 (1 C)" 17(2)
(58)
Regrouping all the constants in (54), we obtain for (1))
70000 1イ (0) + 1) (1), 17 (2). T vec up T.
(59)
Now, we want to rewrite the regression on 2Z) in terms of t(2) where
((Z) (1,2T, 27
43 1
23 - 1
V2 2, Zida.
L2Z3.....
(60)
V2
V2
which satisfies 18 [tt [ Lmm with m' d d(d + 1)/2 + 1. We do that in two steps, let
tz(Z)= = 1, -1
Z2 -
Z3 - 1
v2 Z2 ........ Zd Z, Z2,
ZzZ3.
(61)
v22
V2
Let 5- (710),001,3 T) R Id of be such that
マ (z) sZ),
(62)
i.e, keeping only the constant terms and the terms quadratic in Z,
d
d-1
7 1 1/2
7(00) + (vec ZZT); = =7(0) +
Σ
(ve ZZT)
);-
(+d++1b
1=1
#+++(++1)
(63)
We set, for any & 2 0, (d++ = (2)) (1++ /2, and 7(0) = 7(00 +> (2)
(2]
(1+(++1)8) 1)k) Then, (62)
and (63) are satisfied. To go from ( to 1, we need to get rid of the coordinates ((Z); = Z:Zj
for some i > j. i.e., k E Idp 1, (d + 1)p] for some integer p. Let r unvec(f(2)), and let (2) E Rd(d+1)/2 be defined by
7000 = (V2 20,,,,,,1110 2....... 2D i,a 22 2T,,3...... 2T ...,, v2[
(64)
Then 7 = -7110, q(1)T, ,,,,,,00 (2). mmm satisfies
? =
(X)
(65)
All the previous computations give the expression of 7 as a function of r7.
We now turn to prove Theorem 4.1 using the previous Lemma.
Proof of Theorem 4.1. As in the proof of Theorem D.2, the least squares regression on s(X) can be rewritten in terms of t()). Then, by applying Lemma D.3, we can map the regressor y with respect to I, to the regressor with respect to 8, given 3 = Φ(1).. Since Ez(uT] = I. the OLS simplifies to 7 E (z)) CZ)). By Lemma D.3, the mapping from 7 to 3 is given by
Pe. ((1).T 4 3(2).T vec pi
B
C-T (1) - 21 3(2)
(66)
vec (C-11C-1)
26

---

where
1i 1/2
12)
(2)
12
(2)
12
7d+1 1/22
(3)
72d 1/2
(2)
(67)
(2) /2
Iald 1)/2 V2
(2)
or componentwise 1 1 01 (i-1)) 1) 2. r = 111/22224 -) 1)( for S and S < d- and P..j = 3 for < Regarding the complexity, the computation of the Cholesky matrix and its inverse requires (d operations; consequently, the computation of y can also be performed in (d9) operations. Using (66) to map to 1,, involves computing vec (C-ITC-T) and C-Tx(1) both of which can be performed in (dd))
(2)
D.5
Concentration bounds for the Fisher matrix in the compact case
We now prove Lemma to control the bias induced by inverting the estimated FIM, conditioned on the event that the estimated FIM is celliconditioned, which happens with high-probability given the number of samples N is large enough. We first prove version of the Lemma when is bounded (Lemma D.4), and then tackle the case where 8 is unbounded but with bounded-moments (Lemmaa D.5).
Lemma D.4 (Mean error bound for the inverse of F when 8 is uniformly bounded). Let & E (0, 1), N B(4/3r 3r 2 log (2mmmm))))))))) "), w E W, and A(w) = [I] - Full 11-'1 Then, under Assumptions 2.1, 3.2, and | I B. A(w) occurs with probability at least Furthermore,
-' - FF=' [(w) -ON-
(68)
where the constant in the big-O lerm can be chosen independently of w.
Proof of Lemma D.4 (exponential tail bound). Fix w E W. For the sake of notation, we drop the subscript in w but indicate the dependency in N. For any N 2 1, let EN = N Σ 88 (X,)
with X XN ~ Conditionally on A = [|| - Fl < | FN = F(1 -
iid
F-1 FN)) is invertible because F is invertible thanks to Assumption 2.1 and 0 < 1 - |
< Fx
F FNll
|| F- Fx)]l. Using the Neumann series, have
F-1-F F -1- -- F= (1))) FT.
(69)
Since I 1-F-1 FN is an average of i..dd random matrices with uniformly bounded second moments thanks to Assumptions 3.1 and 3.2, (VN(I - F-1 Fx),, is uniformly integrable (both in N and w). By the strong law of large numbers, have FN "3 F. thus 1[4(N)] 3 1. Since (VN(I- - F- FN)(,) is uniformly integrable, the central limit also holds for the sequence of conditional random matrix components VN((I F- Fx)es A Then for any S í,j S m, VN((I M) A(N) converges in law with uniformly bounded variance. Thus, conditioned on A(,) VN((1 F) Op()), which implies that vNiII F-1IF'le Op(1). Thus N] F- ((N) Op()), and taking the expectation in (69) yields 1Cf--1-F
A|]]] O(N -1), again using Assumption 3.2.
By [57, Th. 1.62] with uniform bound [[(ss T(X.) - F)/N < 2B/N and variance
D> [B[((ss] (X₂) F)/n)3II < B||II/N,, and the definition of r (Assumption 3.2), we have
P[[ - Fl| 2 | < P( - Fl| 2 r)
2
< emexp
B(4/3r + 2
(70)
Nr
< 2mexp
B(4/3r + +2B)
27

---

where to go from the second to the third line, we use [| S B. Setting N 2 B(4/3r + 2 log(mmo-) yields P( Fll > ||--)|1- < 6, i.e., P(A) 2 - d. The bound is independent of w, and true for any w W. finally yielding the result.
Lemma D.5. Under Assumptions 2.1, 3.1 and 3.2, for VN 2 rr20 (v lo + Sep log(m)), P(A(w)) 2 1 and (68) holds.
Proof of Lemma D.5 (polynomial tail bound). We follow the proof of D.4. The CLT for F still holds thanks to Assumption 3.1, and by the same argument as in D.4, the conditional CLT is still valid. Thus, ||B[/- F-1F | Alll O(N-)), and the constant inside the big-O notation is also inde- pendent on w using the uniform bounds on the fourth-moment Bienaymé-Tchebee 's inequality, and [38, Theorem 3.1]. we have
s. By the definition of r, the
P(IEE Fll 2 | < P([ P'l 2 r)
S E
FN
F12//2
log(m)
log(m)
<
&e
Nav + 8epu
(71)
N
N
Be log(m)uuu + Se log(m)}
where last line is solely here to simplify the requirement on N. Setting VN > r r-28-1 (v8e log(m) Seet log(m)) yields P([]Fx-FIl) 2 1 The bound is independent of and true for any w E W.
i,i.e Piche..Pl P(A) 2 1-
D.6 Convergence analysis of the stochastic LSVI algorithm
The proof of Theorem 3.4 relies on Lemmas D.6, D.7, D.9, and D.10. Lemma D.6 states the equiv- alence between stochastic mirror descent and stochastic natural gradient descent, the proof is very similar to the non-stochastic case (Proposition 2.4). Lemma D.7 gives the general convergence rate for stochastic mirror descent with the presence of additional bias under the assumption the bias has bounded variance. This is a generalisation of Hanzely and Richtárik [30, Th. 4.5). Both Lemmas D.9, D.10 are required handle the two first moments of the bias induced by inverting the FIM estimate. This analysis requires conditioning on the event that the estimated FIMs are well-conditioned, which happens with cech-rrbbbbiltty (Lemma D.5). Theorem 3.4 follows by suc- cessively applying Lemma D.6 and Lemma D.7, the latter requires Lemmas D.9 and D.10.
Lemma D.6 (Equivalence between stochastic mirror descent and stochastic natural gradient de- scent). Define the stochastic gradient Val by
Val w + n(ww) fo aws
(72)
given that Py is invertible, Then (10) is equivalent to
f = i - E, V M(se,, where wi w(fi). Farthermore, the previous dynamic is equivalent
(73)
{
Vi + Dz- (w, we))
@+t1 = -argmin
(74)
with fp+1 = 11(
Proof. The first equivalence follows from the same computations in Proposition 2.4. Let show that iteration (10) can be recovered as the dual in the natural parameter space of a stochastic mirror descent, i.e. that 1 n((++1) ) with (w)) given by (74) recovers (10). The first order condition on (74) gives
VZ7(4)11) = VZ*(2,) - - V.l.a.
(75)
However, since (2*(4)+) = 11(++11) 1 the desired equivalence between the two dynamics follows.
28

---

Lemma D.7 (General convergence for biased stochastic mirror descent). Let us define the bias B of the stochastic gradient at iteration by
B₁ [(V - vw((Wt) Wr],
(76)
given that Fa, is invertible and let us denote by m((s)) = Wh+1,- the exact mirror-descent iterate starting from We i.e.,
W1+1,+ = argmin
{ (i + Ef Dz- (wo wr)}-
(77)
Assume there exists 0² > (to be specified later) such that for any
E [[ (Wi+1,+ w(+1) [+)] S d²EL.
(78)
Let &, S + A for all f 2 0, let e = c2-167-1 - 11) for 2 1, and let co = 1. Let
Ck -t=1 C-1 for K 2 1. Then, under Assumptions 2.1, 2.5, and the additional bounded-noise assumption (78),
Ct- [((()))) - [(w0)] (En - (4) uKL( 9mm)
0
k-1
GEt
Ck
C₁
(79)
Ca B (w" - @(+1)).
1=0
Proof. Assumption 2.1 allows us to define (77). Under Assumption 2.5 and the boundedness of the gradient estimate (78), derive slightly modified version of the descent lemma [30, Lemma 4.3] which accounts for the presence of the bias. Next line follows from the calculations done the proof of Hanzely and Richtárik [30, Lemma 4.3]:
1-
1 DE (w* w(+1) [v)]
B((((++1) - I(w²) [e)] s
E1 -M)
Dz- Et
Et (1-1) EE (we+1, wx) wel + B, (w- we+1,a)
(80)
+ 5,0²
- E (#0+1 w ) We]l
where w" argmin (ww). Since E, 2 L, the fourth term negative. Therefore, (80) becomes
1
1 B[Dz- (00+))))) wxy]
E[((@++1) - I(w") |
11 Dz- Dz-(w",, (w",ar)
Et
(81)
2,0² 2 BT (w- w(+1,r) - E[B, (Wi41 wo-1,+) wr).
Taking the expectation of (81) gives
E[((0++1) - [[w*]] $
E1 -1 E[D2- (w= ,0,, - E[D2-(w** (w=,(+)))
(82)
+ +,,02 + E[B. (w- w.+1)..
Let c. ( -n)-1 for / > 1, and let co = 1. Let k: 2 I and define Cs = -ti C-11-
Since e, < f. we have et 2 Multiply by a (82) and sum for t k). then divide by Ck.
E[((())) -(w [(w"))
(5 1) (w" Wb)
#-1
Σ
GE
+
k-1
ct ElB (w) w (83)
C
C
C
C
t=
We essentially recover Hanzely and Richtárik [30, Th. 4.5], but with the additional bias terms. Finally, (79) follows from (83) and Dz (w" (w)) = uKT 9.
Remark D.8. The previous lemma requires boundedness assumption on the gradient estimate given by (78). This assumption is typically required for proving such descent lemmas, see [20, 30, 36, 37, 58]. In particular, this assumption is implied, by Cauchy-Schwand inequality, if the gradient estimate V.ll mapping given by (72) has bounded bias and variance, or directly if the gradient estimate unbiased. We will prove is satisfied under Assumptions 3.1 and 3.3 in Lemma D.10.
29

---

Lemma D.9 (Controlling the bias terms in (79)), Let 2 O, and let A₁ = - with A(w) =
Fl Full < || ||="] then under Assumptions 2.1, 3.1, 3.2 and 3.3, for any E [0, Ai]
ElB (w w++1) ) As] S ((N=+) m2m?/22u p4 x (B[lkw* - w 112 Aal) 1/2 (84)
Proof. By Cauchy Schwarz inequality, for any t E [0, R), B (w* w(+1) Ass & (Bil|Bll? A& - 2 AAll)²²
(85)
Conditionally Ak. Fa, is invertible, and by Lemma D.5, there exists C > 0 such that for N large enough, |I|E Fz.1 AR. will = C. Consequently.
N2 Bil2
2n = NN
20.
Ax Well?
Wi
1/51
Aus [ /حح 1|?
Fiol
Ass w1,]1 11
(86)
- |22,,,11 11 I C m2mpu
where we used go from the first the second line, the independency of zs. with Fs, conditioned on We and E[2e, [r] 20., and to go from the third to the fourth line, we use Lemma D.5, requiring Assumptions 2.1, 3.1 and 3.2, and to go from the fourth to the fifth line, we use bounds on the moment of and f² given by Assumptions 3.1, 3.3:
||ا||| || S (B||s)3) (B8)3))))
1/2
< m sup max E|s|? NEW Isism
X m2
(87)
S m 1/2 pam2
Taking the expectation of (86) conditioned on A₁, and plugging it into (85) yields (84).
Lemma D.10 EHigh-probabilie uniform bound for the variance of the gradient estimate). Let >> O. For any w E W, let M (w: argminn (V..((( WWW(Ww)' +D (w/, ,w) be the exact mirror-descent iterate starting from w with step size E. Similarly, let (w) be the mirror using gradient estimate Vil given by (72). Let o²(w) be defined by o²(w) E 'E((E(V))W))
(M)W)))(((()))))) - Mw))),,
(88)
where all the expectations are taken conditioned on A(w). Under Assumptions 2.1, 3.1, and 3.3, there exists some constant C O, such that for N large enough (see Lemma D.5), and any w E W.
o2(w)) S N [C + (1)),
(89)
for some constant C > O, and the little-o term is independent of W.
Proof. To obtain the uniform bound on (88), we independently bound both terms in the scalar prod-
uct. Let w E W, and let 7 be the corresponding natural parameter, 1/ = n7(w). Conditionally on
A(w), by the computations done in the proof of Lemmas D.4 and D.5 there exists a constant C > 0 such that for N large enough, [[[[]] F-1 11 < C, i.e., the mean-square error is O and thus the variance F - E| "1 O(1/N).
Using the definition of the stochastic gradient İi given by (72), and the previous bound, the first term is bounded by:
VE(]EE(__)(W)) +-]]]]= NE
[[(E|| P )m f (z a.m. ||]
< 2N
(| '1 - Ar )zm 11 + ETILE (aw 81))111
< 2N [|===]]² x [8I
アー ||
(90)
+ E[|l A (2w - f.m )1
< 2N (mpîm3 C/N + El| - (zw - -=mm(1 }=))
30

---

where we used (a + b)2 < 2a² + 2b from the first to the second line. from the second to the third line the uub-muliipliaatiiity of Il-Il. from the third to the fourth || I mylmg and the (1/N)
bound on the variance of _ Furthermore, we can bound the last term in (90) by
NE - 2 ) < NEI " - [1]]
[N - F - + F 11 2ull
(91)
2N 2 B||| 2ي
C²/N
where we used that r r sup 11? and E|| 2 (r
with Assump-
tions 3.1 and 3.3 to state that there exists C" >> 0 independent on such that R|| ~ - لي - C²/N. Gathering (90) and (91) yields for the first term of the scalar product:
NT which in turn can be bounded by some constant C" > 0 independent of w. Let us tackle the second term inside the scalar product. By Proposition 2.4,
- V. (13 S 2Cmymm3 + 4C' (r C/NN,
(92)
((w) = w(1 e [w)((w))
(93)
and similarly for M,
(w)
(94)
Under Assumption 2.1, the mapping 1 E V + w(n) is differentiable with Vupo F, ss Gn- see (35). Let us denote by H; D2 the Hessian of the i-th component application of w for any 1 S i < m, which is Rm*" matrix given by D2w(() = s;ss T 4n- and let D2 =
(H1, H2,.. ................................ Hm) be the collection of the Hessian matrices. For any E E., let us denote by DRwith, (h, hi] [h]² (h Hihh.... mmtl) E R". A Taylor expansion with Lagrange remainder yields.
M(w)- Mw) = w(n - =Va((w)) - (Vw(w)
= ((V.(w) - Vollw)) + [ (1 (D - (V))
(95)
- (w))) dt
Let R be the 5² remainder term in (95). Then, R/Zw is a R vector whose norm can be uniformly bounded using the uniform bounds on the fourth fourth-moment moment of using similar techniques as for the bound on || (see below), we omit the details. This implies that c² R Zwo)s) with constant in the little-o terms independent on w. Consequently,
/ (V((BI) M(
Val(w)) 2=)
< 2 - V=((w))111 + Z()
< (2 o())
(96)
< (m + +o1)
S (C + =())),
for some constant (44) > 0, and where we used || |F < m2p44:
| < 11
Σ
Es
(97)
15tJ5m
3mpp1 p4
using the definition of the Frobesnius norm and then Cauchy Schwarz inequality bound compo- nentwise Fu. By Cauchy Schwarz inequality, (92) and (96), for N large enough,
No²(w) < N( < (C (O" Z())
- (( (1)
(98)
This concludes the proof.
31

---

With previous Lemmas D.6, D.7, D.9, and D.10 in hand, we can prove the main result.
Proof of Theorem 3.4. Define wi as given in the theorem. By convexity of I,
((*x. - ((w')( - Σ
Ca-(((())) - ((w)))
(99)
C2 ==1
Combining Lemma D.7 with Lemma D.9 to control the bias terms, find that the expectation of the RHS (99) is upper bounded by
(00 11) uKI,(uw- Iwo
B(((a)) I(w") And
Ck
(100)
2
CEEL OONN x Sk.N.
1=
where Sk,N == mmm mM Σizo (Bilww" 0 1| Ael)) 1/2 where the big- term is independent
Ct
of l since it is independent of wo = wo W1 ... We and o² some upper bound of supk21 (x)) with (k) satisfying the assumption of Lemma D.7, for all t t k:
(E(V))(0)) We, Asl) (Wi+1, @ct1) cin] S o(((
(101)
By Lemma D010, we can set
o² (k) = -N-CC
max 0stsk-1
Zee
(102)
for some constant C' independent on N and the sequence Wo...., 0k-1
Let us tackle the terms which depend both upon N and k via the sequence Wo 21 .i.............................. Wy. By the law of large numbers, N >> oc, Fu.. Fwg and 2ws Awg almost surely. Then, by the continuous mapping theorem. (ul(w))) where wi is the first mirror-descent iterate. By induction, we obtain that for any k 2 1, We - wf a-s for all t E [1, kt.. i.e., the finite sequence (wp. ... w converges to the exact mirror-descent sequence (Wo, ........... ] We deduce that, almost surely, for all t 2 1, || - -W:ll - || w/ I| since the countable intersection of almost sure events is an almost sure event. By Frank-ska et al. [Thh 4 31], we know that the Mirror-Descent sequence l(w;) converges to I(w"). Since / is strongly-convex, I(w])) + I(w") implies that ||W> - will + 0 as goes to xo Combining with the previous almost-sure convergence, we obtain that for any & 2 1, the following equality holds almost surely,
(ad())) almost surely, and thus W1 - W1,+ wi almost surely,
lim max ||0" - well max ||0" - will == DE: -0000 IStck
(103)
Istsk
with supk21 Dr < 0o. For t 2 1, let U* C W be the closed-ball of center wt and of radius DR, let Un {wo), and let U be the reunion of Un and the ball centered at wt with radius suPR>1 Dr < xx. For each fixed l > 1, almost surely, when N + co, for any 0 < < k. ai E Ur. and therefore
wo, 11, wy E UEZOUR CU.
(104)
Since + Zu is continuous and U is compact, we have supweu Zn < x0. Then almost surely, as N + 0o, sup (k) + N- N''Csupwer Zu == o², and therefore almost surely, as N + co.
sup ( < 2N C sup Zu = 20² < 0o. A21
(105)
WEU
Almost surely, when N 0o, for any 0 < t 2 k, 1|0* well? < 2Dp. which implies that
1/2 < cc. Using Σi= C = Cz. and supostsk-1 E|||w* @++1ll 2| / < (2 suspp21 D22
bounding uniformly the summands of S yields
SK,N & m2m 1/2 p2 sup D&y 1/2 k21
(106)
32

---

Finally, plugging (105) and (106) into (100) yields the uniform bound over k:
B[((@x)) - ((W") As -
(50
(wa) + O(N- (4) u
Σ 1=0
GEt
O(N-)).
(107)
Ck
All the constants in the big-O terms can be chosen independently on the sequence of w.
Using Proposition D.5 with (/ + 1) and union bound, have P(ntzo chosen N.
ô for the
Finally, let prove the explicit convergence rates for linearly increasing stepsizes E (L + at) t 2 O. Similarly to Hanzely and Richtárik [Lemma 4.8 30J, we distinguish three cases depending on compared to If a > then C₂ = 12kP//a) and Σi=0 Gg = ((1) which yields Okk + ONN-1) for the RHS of (12). If μ4, then Ck (kk) and Di Ggε Ollog())). If p. then C = ((k4/a)) and k-1 Ge = O(ka/a-1).
33

---

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]
The main contributions of the paper are as follows: (i) KL minimisation within exponential families be performed via successive linear regressions under tem- pered variational approximations. This approach is equivalent to natural gradient descent (NGD) and mirror descent (MD) but avoids the need for explicit gradient- based procedures. This is detailed in Section 2 and Proposition 2.4, with relevant prior work cited. (ii) In the Gaussian variational family, exact LSVI can be tailored to eliminate the need invert the Fisher information matrix. The resulting procedures have computational complexity (dd) in the full-covariance case and (d) the mean-field case, as shown in Section 4, specifically Theorems 4.1 and D.2. (iii) Under standard optimization assumptions, LSVI converges at explicit rates, established in Theorem 3.4 (Section 3). (iv) Empirical results demonstrate that LSVI achieves performance comparable to state-of-the-art variational in- ference methods and remains effective for non-differentiable target densities, as shown in Section 5 and Appendix C.
Justification:
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 in 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: (i) The approach is currently restricted to exponential families. Extending LSVI to more general variational families, such as mixtures of exponential families, is a promising direction for future work. (ii) Convergence guarantees for the stochastic versions LSVI, namely, LSVI-MF and LSVI-FC, are not yet established. (ii) A comprehensive theoretical analysis of LSVI's convergence when the target distribution æ is replaced by unbiased estimator ă (e.g., via subsampling) remains an open problem. (iii) A comprehen- sive study of the constants involved in the convergence rates, in particular with respect to the smallest singular value of the FIM r, the latent dimension d and the dimension of the statistic is left for future work. We believe most of the proofs can be adapted but the analysis would be more involved. See Section 6.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The au- thors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g. if the approach was only tested on a few datasets or with few runs. general, empirical results often depend on implicit assumptions, which should be articulated.
34

---

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-tex system might not be used reliably provide closed captions for online lectures because 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 acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: All assumptions are provided in the core manuscript, see Theorem 2.1, The- orem 2.5, Theorem 3.1, Theorem 3.2 and Theorem 3.3. Furthermore, the assumptions are discussed in the core of the paper along with references mentioning existing similar as- sumptions in the VI literature. The proofs are deferred to the supplementary materials and are divided several comprehensive steps, including lemmas in order to make the proof- reading procedure easier.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theo- rems.
The proofs 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 proof's provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
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)?
Question:
Answer: [Yes]
Justification: Section 5 and Appendix C along with the pseudo-code Algorithms given in Section sufficient to reproduce all the experimental results. In particular, all input parameters are provided in Table 3.
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 not.
If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible verifiable.
35

---

Depending on the contribution, reproducibility can be aecomplished in various ways. For example, if the contribution is a novel architecture, deseribing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate the research performed.
While NeurlPS 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) the contribution is primarily a 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) the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to produce the model (e.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 case thors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source source models, it may be that access to the model is limited in some way (c.g., to registered users), but it should be possible for other researchers have some path reproducing 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?
Answer: [Yes]
Justification: The paper provides a Python (JAX) package that includes all discussed Al- gorithms (LSVI Algorithm 1, MF-LSVI Algorithm 2, FC-LSVI Algorithm 3, Variance control for the stepsizes Algorithm 4, NGD with details provided in Section C) as well as scripts to reproduce all the listed experiments. The package is explicitly divided into two parts, variational contains the generic implementations while experiments contains three sub-folders for the three distinct variational problems (logistic regression, variable selection and Bayesian synthetic likelihood). Full-pipeline for the experiments is provided (download and pre-processing of the datasets, inference procedures and post-processing scripts).
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/ SruuddesCooeeuumissttts
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for 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: / httpsscccppplccgggeeeeeeee htallii/gguiesseeee Sump://iidSSSSSeseeeeeee asionPolicy) 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. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
36

---

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 is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The paper provides all the necessary details reproduce the experiments, including the hyperparameters (the number of samples N, the number iterations T, initialisation distributions and the schedules) which are given in Section C. Different sched- ules have been considered to demonstrate robustness of the proposed methods while the number of samples is set to obtain reasonable numerical stability.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, 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: All experiments were conducted using multiple trials as indicated in the figure labels and the Appendix C. For the logistic regression examples, one -eviaaion confidence intervals are provided over 100 independent realisations. For the variable selec- tion problem, the means and the min-max intervals for the posterior marginal probabilities obtained via LSVI over 100 independent realisations. No statistical assumption is made for uncertainty measurement.
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, statistical significance tests, 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, or 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).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should prefer- ably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
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
37

---

Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory. time of execution) needed reproduce the experiments?
Answer: [Yes]
Justification: The full hardware and software specifications provided in Appendix C, specifically in Table 1, along with Figures 4 and 5, which report experiments runtime and computed using independent realisations for
memory usage. All performance statistics improved robustness. In addition, scripts for measuring the runtime and memory usage of the algorithms can be found in the package (/ //expermmenss/_ .J/time . pyy). All experiments were successfully performed and are reported in Section 5 and Appendix C. There is the exception of the applicability of ADVI (PyMC3, [6]) the MNIST dataset, which is explicitly stated Appendix C. Instead, ADVI as provided by Blackjax [7] was used as replacement PyMC3.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g.. preliminary or failed experiments that didn't '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 ttttps///eerrr wwwpeuuppp..ppppl////eeeecceeeee
Answer: [Yes]
Justification: We carefully read through the NeurIPS Code of Ethics, and we see viola- tion of any guideline.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS 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).
Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
The paper focuses on VI methods, emphasizing theoretical analysis and algo- rithmic implementability. As such, the work foundational in nature and does not directly pertain to real- world applications or deployments. Given its abstract and theoretical scope, it does not present identifiable positive or negative societal impacts, including concerns related to fairness, privacy, misuse.
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 socictal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact spe- cific groups), privacy considerations, and security considerations.
38

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or 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, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of MI...
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g.. pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
We no risk in the application of variational inference procedures.
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 quiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do require this, but encourage authors to take this into account and make best faith effort.
Licenses for existing assets
Question: Are the creators original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: All used Python packages are ppen-source, have permissive licenses, and are explicitly mentioned both in the manuscript and the code (pyproject. toml with com- plete dependency specifications). Specifically, Blackjax, PyMC, and JAX are mentioned in Section 1 and Section 5. Details on the datasets used, licenses and download links, provided in Appendix C.
Guidelines:
The answer NA means that the paper does not 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, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
license, copyright information, and terms of use the pack- age should be provided. For popular datasets, paperswi paporawithcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the li- cense of dataset.
If assets are released,
thcode
39

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and the documenta- tion provided alongside the assets?
Answer: [Yes]
Justification:
The provided Python (JAX) package for LSVI is well documented and in- cludes a README file with instructions for installation and usage. The license is also included in the package (Apache License 2.0). In addition, provide usage examples and accompanying commentaries.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/coodemmodl part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
For crowdsourcing experiments and research with human subjects, does the pa- per include the full text instructions given to participants and screenshots, applicable, well as details about compensation (if any)?
Question:
Answer: [NA]
Justification:
No experiment involving human subjects conducted.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing 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 in the main paper.
According to the NeurIPS Code Ethics, workers involved in data collection, cura- tion, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
No experiment involving human subjects were conducted.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equiva- lent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
40

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an 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 core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification:
There is no mention of LLMs in the manuscript, and LLM was used.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy pattps://neurihe cc/Conferences//2055 for what should or should be described.
41


Reviews:
[{"review_id": "yRgacAh80C", "reviewer": "Reviewer_X8ht", "summary": "This paper proposes an new variational inference method called least square variational inference.\nA connection has been drawn to natural gradient descent and mirror descent, which is used to prove a convergence rate of the algorithm.\nExperiments demonstrate marginal improvement against natural gradient descent.", "questions": "1. In the last line Algorithm 3, why there is an \"unvec\"?\nIf I understand correctly, \\\\(\\hat\\eta_{2, t + 1}\\\\) is already a matrix.\nSo I am not sure what \"unvec\" does here.\n\n1. Can we implement Eq (11) without explicitly inverting the Fisher information matrix for non-Gaussian variational distributions?\nThis is possible in natural gradient descent.", "limitations": "NA", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "1. The writing is not good.\nThe paper is notation dense (which is fine), but the authors did not do a good job walking readers through the material.\n    1. The matrix \\\\(F_\\eta\\\\) is the Fisher information matrix, but this is never mentioned in the paper.\n    1. I am very confused about Section 4 from Line 195 to Line 220.\n    I think this section is intended to derive an efficient implementation for Eq (11) in the special case of Gaussian variational distributions.\n    However, why does Theorem 4.1 have the name \"exact LSVI mapping \\\\(\\phi\\\\) for full-covariance Gaussian distributions\"?\n    If you have an exact form of \\\\(\\phi\\\\), then you are able to implement the fixed-point equation  (6) with zero Monte Carlo error.\n    Due to the connection between LSVI and NGD, this implies that natural gradient descent can be implemented with zero Monte Carlo error even if the likelihood is non-conjugate??\n    1. Not sure what \\\\(t(z)\\\\) is in Eq (16).\n\n1. The theory is weak and uninteresting.\nI feel like Assumption 2.5 basically does most of heavy lifting---it is unclear whether real VI objectives in practice are smooth or strongly convex.\nIn fact, there is evidence that the VI objective is not even convex for non-conjugate likelihoods, let alone strongly convex (Wu et al., 2024).\n\n1. In Eq (10), why do we use Monte Carlo samples to approximate the Fisher information matrix? The Fisher information matrix is already available in a closed-form for (at least many of) exponential families.\n\n1. Ultimately, Algorithm 3 is quite similar to natural gradient descent.\nAlthough the authors present the paper as a new method LSVI, I think the proposed method actually boils down to an (perhaps) improved stochastic gradient estimator for stochastic natural gradient descent.\nThis is just my gut feeling though.\n\nWu, K., & Gardner, J. R. (2024). Understanding stochastic natural gradient variational inference. arXiv preprint arXiv:2406.01870.", "was_revised": true, "final_justification": "I initialled gave a 2. In the rebuttal, the author has cleared some of my confusions, e.g., why do we estimate the Fisher information matrix stochastically even though a closed form exists. However, it also strengthens my believe that the proposed method ultimately boils down to an (perhaps) improved stochastic gradient estimator for stochastic natural gradient descent. I'll willing to raise my score to borderline reject. But I still think there are two major weaknesses remaining: 1) the author didn't do a good job in presentation; b) the novelty is limited.", "rebuttal": "We sincerely thank the reviewer for their feedback, which has helped us clarify and improve the paper. Below, we address\neach point in detail.\n\n### 1.1 Clarification on the Fisher information matrix (FIM)\n\nAgreed, we now mention that $F_\\eta$ is the Fisher information matrix as soon\nas it is defined.\n\n### 1.2. Exact LSVI for Gaussian Variational Families\n\nThe term \"exact\" in \"Exact LSVI Mapping for Full-Covariance Gaussian Distributions\" (Th. 4.1) refers to\nthe *oracle setting* where expectations (e.g., $ \\mathbb{E}_q[\\cdot] $) are computed exactly, as introduced in Section 2. \nThis does *not* imply zero Monte Carlo error in practice. More precisely, Theorem 4.1 gives the LSVI mapping in term\nof $\\gamma$ (Line 210), which is an intractable expectation (in most cases), and needs to be estimated via Monte Carlo.\nWe will clarify this distinction in the text to avoid confusion.\n\nRegarding the connection to NGD: Indeed, LSVI and NGD coincide in the exact setting (with oracle expectations, Proposition\n2.4). Our key contribution is a *gradient-free* stochastic approximation that avoids evaluating gradients of $ \\log \\pi $,\nplus matrix inversion of the FIM when the family is Gaussian (Algorithms 2, 3) .\n\n### 1.3. Reparameterization and the role of $t(Z)$ in Eq. (16)\n\nIn a nutshell, the point of $t(Z)$ is to define a simple, one-to-one\ntransformation such that the output vector has un-correlated components:\n$\\mathbb{E}[t(Z) t^\\top(Z)] = I$. That makes possible to perform linear\nregression without inverting the FIM.\n\nWe explicitly give the mapping from $\\gamma$ to $\\eta$ depending on $(\\mu,\n\\Sigma)$ (Theorem 4.1), and we derive the corresponding tailored stochastic\nschemes in 4.1. We now more thoroughly emphasize the role of $t$ in the text.\n\n\n### 2. Convexity assumption\n\nWe agree the strong-convexity/convexity conditions are rarely verified in\npractice, however, such assumptions are standard for analyzing convergence of\noptimization algorithms (including SGD and NGD), to ensure a unique minimizer\nand tractable rates (see Remark 2.6 and, e.g. Chap. 5 in the book of Bach,\n2024).   While non-conjugate VI objectives may not be globally convex, *local*\nconvexity near optima often suffices in practice. We emphasize that our theory\nprovides guarantees under idealized conditions (as is common in the optimization\nliterature), and empirical results demonstrate robustness beyond these\nsettings. We will further clarify that:\n\n- Convergence guarantees for mirror descent can be derived under non-strongly convex (but convex) assumption (see Lines\n  115-120).\n- Local convexity near optima often suffices in practice, and local convergence to local minimisers still hold in that\n  case.\n- Convexity holds when the variational family contains the target (e.g., conjugate models, Wu et al., 2024).\n\n### 3. Monte Carlo approximation of the FIM (Eq. 10) and variance reduction\n\nWe agree that closed-form FIMs exist for several exponential families, however we explicitly use a Monte Carlo estimate of\nthe FIM rather that its close-form expression. \nThis is because our method treats the KL minimisation problem *as a linear regression task*: the Monte Carlo\nestimator $ \\hat F_{\\eta} $, (Eq. 10), is the *empirical FIM matrix*, and $ \\hat z_{\\eta} $ is the empirical moment\nvector (Eq. 10). As a consequence, the variance-optimal estimator for the least-squares problem is $ \\hat F_{\\eta} \\hat z_\n{\\eta} $, and not $ F_{\\eta} \\hat z_{\\eta} $ where $ F_{\\eta} $ is the closed-form FIM (i.e., no Monte Carlo error).\nTo sum up:\n\n- Using the exact FIM $ F_{\\eta} $ and the Monte Carlo moment vector $ \\hat z_{\\eta} $ to compute the regressor is not\n  optimal with respect to the variance.\n- Using the same samples for $ \\hat F_{\\eta} $ and $ \\hat z_{\\eta} $ reduces variance (as in standard linear\n  regression, see Remark 3.6).\n\nIn the Gaussian case, we implicitly compute the exact inverse FIM via the reparameterization introduced in Section 4,\nspecifically, Theorem 4.1. The resulting Algorithms 2, 3, are faster (i.e., smaller CPU cost per iteration), but noisier.\n\n### 4. Relationship to natural gradient descent (NGD)\n\nWe agree with you, as we state in Section 2.1, exact LSVI recovers NGD on the *unormalised* KL objective. Our\nfocus is on *practical, gradient-free approximations* via a linear regression formulation:\n\n- Algorithm 1 does not require computing gradients of the log-target distribution, contrary to NGD, which enables us to\ntackle\npreviously inaccessible models.\n- Our algorithms yield very stable numerical results, and are easy to implement, and to tune.\n- All our algorithms can leverage highly optimized routines for performing the linear regression.\n\n**Questions**:\n\n1.`unvec` in Algorithm 3: The natural parameter $ \\eta $ is a vector for general exponential families (i.e., the variational distribution is $\nq_\\eta(x)\\propto\\exp(\\eta^\\top s(x)) $). For Gaussians, $ \\eta_2 = -\\frac{1}{2}\\textup{vec}(\\Sigma^{-1}) $ (Line 200), and the `unvec` reshapes the\nvector into the new precision matrix. Sorry for the confusion, we now clarify this.\n\n2.Avoiding FIM Inversion for Non-Gaussian Families: For general exponential families, FIM inversion can only be avoided\nif gradients with respect to the moment parameters $ \\omega $ are available (since $ \\nabla_\\omega \\ell = F_\n\\eta^{-1} \\nabla_\\eta \\ell $). Without such gradients, explicit inversion is required. \n\n**References**\n\n- Bach. Learning Theory from First Principles, MIT Press, 2024.\n- Kaiwen Wu and Jacob R. Gardner. Understanding stochastic natural gradient variational inference, ICML 2024."}, {"review_id": "3vckxYdbh2", "reviewer": "Reviewer_7qp3", "summary": "This work proposes a method of iteratively running least squares regressions to fit the parameters for variational inference over exponential families. In particular, the special case of fitting an exponential variational family can be shown to be equivalent to solving a fixed point problem. With this starting point, the authors of this paper then proceeded to demonstrate how such a fixed point method could be solved, first exactly assuming oracle access to expectations and then with practical Monte Carlo approximations. They further demonstrated that this MC estimate has convergence guarantees on the recovery of the true optimal variational parameters as would be produced under exact LSVI. Even this practical MC algorithm, however, involves an expensive step of inverting the Fisher Information Matrix (FIM), rendering this procedure more costly than would likely be desirable. In particular cases, specifically concentrating on the Gaussian case in this paper, such an inversion can be circumvented. Finally, the method was robustly demonstrated to perform well across many empirical settings.", "questions": "1. Mostly out of curiousity, where did the inspiration of introducing t(z) to compute the FIM inverse emerge from? Is this common in some other subfield?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "-------\nQuality\n-------\nThere were no technical flaws in either the theorem statements or proofs that I was able to identify. The proof assumptions were all stated precisely and explicitly referenced in the statements of the theorems. Additionally, the empirical findings were properly framed against alternate variational inference estimation techniques (i.e. in Figures 1 and 2).\n\n-------\nClarity\n-------\nThe paper was very clearly presented. As stated in the previous discussion, all the theorem statements were precisely stated and the assumptions similarly explicitly provided. There were a few spots that would have benefitted from some exposition to give clarity on significance, specifically on the convergence guarantees of Theorem 3.4, where the decomposition of the three terms (i.e. first term for initialization error, second for Monte Carlo error, and third for iteration convergence) that would have made parsing the significance of some of the results and where they arose from from the proofs more direct. Nonetheless, the authors provided comprehensive discussions around the difference between the provided results in the broader context of other works in this space, such as in the discussion of Remark 3.5, which made identifying contributions of the work quite immediate. \n\nIn addition to the individual theorem statements, the flow of the paper was also very well done, where complexities were gradually layered into the final formulation, starting with the exact solution to the fixed point problem, followed by a relaxation to the exactness of the expectations, and finally to the special case of Gaussians. This made following the progression of complexity much easier to follow along. Additionally, the text in the main body was further well-encapsulated in the algorithms, which also used consistent notation as the main body, making it straightforward to follow along with. Finally, the plots and presented results of the experiments section were also cleanly presented and easy to interpret, even in isolation of the main text.\n\n------------\nSignificance\n------------\nThis paper contributes significantly to the area of variational inference, both in likelihood-based and likelihood-free settings. As discussed, the method proposed is ultimately most useful when fitting Gaussian variational approximations due to the tractability of estimating the FIM inverse in the manner outlined in the paper. Nonethless, there are many contributions that all seem very worthwhile. The first is simply the insight that variational inference can be computed in this iterative least squares fashion: this is a very elegant and much more generally applicable insight than just in the special case of Gaussianity. Other researchers can, therefore, leverage this reframing for the development of similarly tractable estimation strategies to layer atop this initial insight of reframing. Additionally, the estimation of the inverse FIM matrix for Gaussian families was a very elegant trick to allow for the reduction from the naive $O(d^6)$ cost to $O(d^3)$ via backsubstitution. It is not exactly clear to me if this type trick has been done in other literature, but I have not seen it referenced in the space of the fitting of variational families, which alone makes its presentation to this community a worthwhile contribution. \n\n-----------\nOriginality\n-----------\nAs discussed above, this paper provides many novel insights, both in the original framing of the variational inference problem as a regression problem and in more granular technical details around the FIM inversion (i.e. in the reparameterization with Eqns 16 and 17). These are novel and worthwhile contributions to the variational inference and also likelihood-free inference communities, whether or not such insights were present in orthogonal communities of statistics.", "was_revised": true, "final_justification": "I maintain my original stance and believe this work is a technically sound and interesting contribution to the VI space. I, therefore, stand by my original score for acceptance.", "rebuttal": "Many thanks for your supportive comments and constructive feedback.\n\nRegarding the readability of the theoretical results, in particular the\nconvergence guarantees Theorem 4.1, we made it easier for the reader to parse\nand understand the practical consequences of this result (i.e., convergence in\naverage of the average of the iterates). We also added a remark on the\ndecomposition of the convergence error (i.e., deterministic error due to\nmirror and the initilisation, mirror descent error coupled to Monte Carlo error, and the vanishing Monte Carlo\nbias). We also added explanations between the steps of the proofs to\nenhance fluidity of the manuscript.\n\n**Questions**:\nThe intuition for the reparameterization comes from nonparametric regression,\ni.e. recovering a function $f$ from data $(X_i, Y_i)$ such that $Y_i = f(X_i) +\n\\epsilon_i$, for some zero-mean i.i.d random variables $\\epsilon_i$. A common strategy is to decompose $f$ on a functional orthonormal basis.\nThen the predictors (the elements of the basis) are independent, and it is\npossible to perform $p$ independent linear regressions (based on a single\npredictor) instead of one (based on all the predictors). In our case, the\nstatistic $t$ is constructed as the normalised Hermite polynomials of order at\nmost two, which form an orthogonal basis for the vector space of polynomials of\norder two under the Gaussian measure."}, {"review_id": "yamuZ39RKi", "reviewer": "Reviewer_WRbV", "summary": "In this work, the authors propose a gradient-free and Monte Carlo-based scheme version of the vanilla variational autoencoder, named the  Least Squares Variational Inference (LSVI). Remarkably, this framework proposes a gradient-free variational inference procedure that reframes the traditional KL divergence minimization target as a sequence of ordinary least squares (OLS) regression problems. This method enables optimization in exponential-family distributions that does not requiring gradients of the log-target density. This framework can be applied to various tasks, like the logistic regression, discrete variable selection, and Bayesian synthetic likelihood.", "questions": "1. Could you provide also a literature review of the paper in this VAE with more theoretical guarantees domain? Like, [1], [2].\n2. Could you apply this method to any real-word datasets and get the experiment results?\n\n\n\n\n[1] Least square variational Bayesian autoencoder with regularization. Ramachandra, et, al. 2017.\n\n[2] Forward \\chi^2 Divergence Based Variational Importance Sampling. Li, et, al. 2023.", "limitations": "Please refer to my 'Strengths And Weaknesses' and 'Questions' section.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "As for the strengthes of the method: The math part of the papers are rigorous and well-written. It's clear to me and the Propositions and the exact LSVI scheme Section 2.1 makes sense to me. Moreover, this work connect the LSVI theoretically to normal gradient descent and mirror descent, and shows the proved convergence rates under mild assumptions. Empirical results on the numerical experiments are widely conducted on the logistic regression, discrete variable selection, and Bayesian synthetic likelihood tasks, which I think is solid.\n\nFor the weaknesses, I think that in the Introduction section, the motivation of the work is not too clear and also the Eq. (1) does not look like the traditional VAE formula showed in its original paper [1]. Moreover, as the VAE framework is proposed for introducing the power of deep neural networks into the variational bayes framework, then adding the constraints listed in this work is actually a trade-off and the impact of this idea would be a bit limited. Thirdly, the method tends to provide the point estimates rather than with UQ.\n\n[1] Auto-Encoding Variational Bayes. Kingma, et, al. 2014.", "was_revised": true, "final_justification": "See my full comments in the review and below. I thank the authors for the clarifications and I will keep my rating at 4.", "rebuttal": "We are happy to cite more extensively the VAE litterature, including the paper\nyou mentioned. Please note however that our focus is on VB (Variational Bayes),\nthat is, a parametric approximation of a *fixed* posterior distribution, which\nis a bit different from VAEs. In a VAE, the variational approximation is only\none component of the whole procedure, which also comprises maximisation with\nrespect to both $\\theta$ (a model hyper-parameter) and $\\phi$ (parameter of the\napproximation family), using the notations of Kingma and Welling (2014). In our\ncase, we do not have a $\\theta$. In addition it does not look straightforward\nto plug neural networks in our approach, as done in VAE (although it may be an\ninteresting line of research in the future). In fact, the neural networks in\nVAE are introduced to map an input  $x$ to a Gaussian approximation of\n$p(z|x)$, based on a large sample of $x_i$ (hence one has as many targets as\nthere are $x_i$), but again in our case we have a single target."}, {"review_id": "nlYzZQJoOT", "reviewer": "Reviewer_aZSa", "summary": "This paper introduces Least Squares Variational Inference (LSVI), a novel gradient-free method for variational inference in exponential families. The core idea is to reformulate the variational fixed-point equation as a sequence of least squares regressions on tempered log-target evaluations. The authors show LSVI corresponds to a biased stochastic natural gradient descent and provide convergence guarantees under mild assumptions. They also develop efficient algorithms for Gaussian variational families, avoiding Fisher matrix inversion and enabling linear or cubic complexity depending on the covariance structure. The method is validated on tasks including logistic regression, discrete variable selection, and likelihood-free inference, demonstrating strong performance even when gradients are unavailable.", "questions": "1. Could LSVI be extended to handle mixtures or flows beyond exponential families?\n2. How stable is the method in higher-dimensional settings where Monte Carlo estimates of the FIM become noisy?\n3. Could the authors clarify under what conditions the regression operator stays within the variational family?\n4. How does LSVI compare in practice with Black-Box VI methods on models with intractable likelihoods?", "limitations": "- Currently restricted to exponential-family approximations.\n- No discussion on combining LSVI with amortized inference or VI for deep generative models.\n- Adaptive step size scheme is proposed but the method is quite sensitive to this point. A sensitivity and ablation study could be useful here.\n- The computational gain over NGD could be better quantified in real high-dimensional applications.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 2, "strengths_and_weaknesses": "**Strengths**\n- Proposes a technically novel and practically useful gradient-free VI method that works even in non-differentiable or discrete settings.\n- Theoretical analysis is solid, with convergence guarantees for both the exact and stochastic schemes under reasonable assumptions.\n- Offers efficient, tailored implementations for Gaussian families, with competitive time complexity.\n- Demonstrates flexibility across a range of problems including variable selection and Bayesian synthetic likelihood.\n- Writing is rigorous, with many technical details included in the appendices. Notation is generally clear, though some could be improved for clarity.\n\n**Weaknesses**\n- While well-motivated, the method is limited to exponential families, which constrains broader applicability.\n- Performance comparisons against deep VI models or more flexible approximation families are missing.\n- The connection to prior works, such as stochastic linear regression VI (e.g. Salimans & Knowles, 2013) could be more explicitly elaborated.\n- The choice of referring to the appendix to clarify the notation is not a good idea, as it makes the main text harder to follow. I suggest the authors make the notation clearer in the main text, or at least provide a summary of the notation there.\n- The practical choice of step size remains somewhat heuristic; adaptive methods are proposed but not systematically evaluated.\n\n**Minor**\n- The GitHub link could be made more visible, moved to a footnote or explicitly listed in the conclusion for better visibility.", "was_revised": true, "final_justification": "The comments by the authors have clarified some of my doubts, so I am more confident in the original evaluation of the paper.", "rebuttal": "Many thanks for your feedback. We address below the mentioned weaknesses.\n\n1. We agree that our method is currently limited to exponential families. Here are some ways of\n  lifting this restriction:\n\n- Mixtures of exponential families may be tackled by adapting the Expectation Maximimisation approach of\n  [Arenz et al., 2013].\n- For Gaussian approximations, if the posterior contains directions that are strongly non-Gaussian, then\n  conditional-Gaussian strategies like INLA [Rue et al., 2009] may be applied.\n- Another way to improve a Gaussian approximation is to use the skew-symmetric modification recently proposed \n  by [Pozza et al., 2024].\n- In discrete exponential families, independence can be lifted by considering tree-structured\n  dependencies [Wainwright & Jordan, 2008], which are quite flexible.\n\n  We will add this discussion to the conclusion.\n\n2. However, we currently do not have a way of directly using LSVI to fit a family of densities defined by deep models via\n   normalising flows. They can perhaps be fitted using layer-wise linearisation, but this is at this stage\n   just a direction for future research.\n3. Thank you for your remark, we agree, and we will elaborate on this in the final version of this manuscript.\n4. Thank you for this suggestion. We moved our notation table to the main text.\n5. Choosing the stepsizes for optimisation procedures in VI is recognised as a challenging problem [Welandawe et al., 2024]. We\n   do provide in the Appendix a comparison of the KL losses using different strategies for choosing the learning step,\n   including, adaptive rule (Alg. 4), and fixed step sizes, see Table 3, and Figures 5, 6, 7 in Appendix D. Any other\n   stepsize strategy could be adapted to our method, such as line search to ensure better convergence results,\n   however, at the cost of additional loss evaluations.  Additionally, as the proposed algorithms are stochastic, a\n   good strategy to stabilise the results is to average over the last iterations.\n\n**Questions**:\n\n1. The expectation-maximisation approach of Arenz, et al. could be extended to handle mixtures of exponential families.\n2. While a rigorous analysis of the convergence of LSVI with a control on the dimension $d$ is out of reach for now, we\n   found the sample complexities for the Fisher matrix in the mean-field case to be $O(d\\log(d)^2)$, and in the\n   full-covariance case $O(d\\log(d)^2)$.\n   Therefore, to maintain stable numerical results as the dimension increases, one needs to increase the number of\n   samples accordingly.\n3. When the target distribution $\\pi$ is in the variational family, the regression operator recovers the exact natural\n   parameter, in one iteration.\n   However, in general, the regression operator is not in the natural space, and therefore, we consider a momentum\n   iteration (Eq. 11), by taking a sufficiently small stepsize, which is always possible since the variational family is \n   an open set under Assumption 2.1.\n\n**References**\n\n- Arenz O., Neumann G., and Zhong M., Efficient gradient-free variational inference using policy search, PMLR, 2018.\n- Rue H., Martino S., and Chopin N., Approximate Bayesian inference for latent Gaussian models by using integrated\n  nested Laplace approximations, JRSSB, 2009.\n- Pozza F., Durante D. and Szabo B., Skew-symmetric approximations of posterior distributions, arXiv:2409.14167, 2024.\n- Martin J. Wainwright and Michael I. Jordan (2008), Graphical Models, Exponential Families, and Variational\n  Inference, Foundations and Trends® in Machine Learning: Vol. 1.\n- Welandawe M., Andersen M., Vehtari A., and Huggins J., A framework for improving the reliability of black-box\n  variational inference. Journal of Machine Learning Research."}, {"review_id": "CJ42knzfeF", "reviewer": "Reviewer_VdPC", "summary": "We have a target distribution of interest, $\\pi$. The goal is to design a VI method for exponential family variational distribution $q_\\eta$ in a gradient free way, i.e., we do not make use of $\\nabla \\log \\pi$. An unnormalised KL is introduced that shares the same optima as the standard KL. Then it is show that the stationary points of the unnormalised KL satisfy\n$$\nE_\\eta [s s^T] \\eta = E_\\eta [f s]\n$$\nwhere $s$ is the score, $f=\\log \\pi$. The connection to least squares is immediate. An iterative algorithm is outlined.\n\nFor moderate-dimensional exponential families, we form monte carlo estimates of the $E_\\eta$ above and directly invert the FIM. Some principled tricks are introduced in the case of high-dimensional Gaussian families that allows us to avoid this costly inversion.", "questions": "1. I believe the acronym MD, first used in Line 93 is never spelled out?\n2. In Prop 2.2, you write $\\eta=(\\eta_0, \\bar \\eta)$. Then in Algorithm 1, the notation $\\eta_0$ is double-loaded?\n3. The cost of inverting $F$ is mitigated only in the case of Gaussians. What about other exponential families where the computation time could still be formidable? Is this, at the moment, only a practical methodology for Gaussian variational families?\n4. Aren't there other gradient-free VI methods out there? Should these be discussed and used as empirical 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 technically solid and self-contained.\n* While the scope is restricted to exponential-family variational approximations, this focus is a strength, not a liability: it yields clean mathematics and clear algorithms.\n* Even ignoring the formal theory, the fixed-point (OLS) view and the Fisher-inversion-free Gaussian trick are valuable practical contributions to variational inference.\n\n**Weaknesses**\n\n* Experimental evidence is limited to small-scale benchmarks; larger or more diverse tasks would better demonstrate practical impact.\n* It remains unclear how often exponential-family variational forms are the best practical choice in modern applications—guidance or examples here would strengthen the paper.", "was_revised": true, "final_justification": "I maintain my very positive score. I disagree with the reviewer that gave a 3 who felt the paper was not well written or presented. This paper is very clear mathematically, far surpassing the typical ML conference submission. The same reviewer also found a lack of novelty. This is subjective of course but I disagree. Often with ML papers, the novelty is front-loaded and so much time is spent on narrative, and so little on the actual content. This paper is very clear in what it is doing, I understand what is new and what is not new just from the writing. I think what is new is sufficient novelty.", "rebuttal": "We thank the reviewer for his feedback and supportive comments. Below, we address the raised points.\n\n1. We agree that higher-dimensional benchmarks for our methods would further\n   strengthen the paper. Regarding the variety of the inference tasks: we\n   provide numerical examples on both discrete and continuous target\n   distributions, including case where the likelihood is intractable.  Another\n   potential example, where gradient-based procedures are not implementable, is\n   PAC-Bayes learning, in which the target distribution takes the form\n   $\\pi(\\theta)\\propto \\mathcal{N}(\\theta| \\theta_0, C) e^{-\\lambda\n   \\hat{R}(\\theta)}$, for some $\\lambda > 0$, and where $\\hat{R}$ is a empirical\n   risk function (e.g., the misclassification rate in a classification task). Indeed, since $\\hat{R}$ is typically piecewise constant, \n   its gradient is almost-everywhere zero.\n\n2. On the restriction to exponential families, we refer the reviewer to our\n   answer to Reviewer aZSa on the same point.\n\n**Questions**:\n\n1. We now spell out \"Mirror Descent (MD)\" the first time it occurs, thank you for catching this.\n2. We agree, and now use the notation $\\eta = (\\eta^{(0)}, \\bar{\\eta})$ in Proposition 2.2.\n3. When the variational family is the independent product of Bernoulli distributions indexed by\n   $p=(p_1, \\ldots, p_d)^{\\top}$, the structure of the FIM can be leveraged to propose a low-complexity scheme.\n   Indeed, in that case, $F$ is the sum of a diagonal matrix $D=-\\textup{diag}(p)$, and a rank-one\n   matrix $pp^{\\top} =(p_ip_j)$, and\n   therefore, the FIM can be inverted with little cost via Woodbury's matrix identity.\n   (This is something we have realised only recently, and have not yet\n   implemented.)\n4. To the best of our knowledge, only Arenz et al. (2018) propose a gradient-free method for VI.\n   We do mention it in the logistic regression example (Section 5.1), and compare it to our methods.\n\n**References**\n\n- Arenz O., Neumann G., and Zhong M., Efficient gradient-free variational inference using policy search, PMLR, 2018."}]

Author rebuttal:

Meta-review: This paper considers the setting of variational inference using exponential family distributions as the approximating family. They consider the problem of minimizing the KL divergence over this family, which they show becomes a least squares problem. A key advantage of this method is that it does not require differentiability of the target, which is often required in many black-box VI formulations.  The paper shows that the LSVI procedure is a biased stochastic natural gradient descent procedure, and use this equivalence to derive rates of convergence. The method is demonstrated empirically on several examples, including a likelihood free problem.

The majority of reviewers were in favor of accepting this paper. One reviewer, who initially was in favor of rejection, raised their score recommendation after the rebuttal phase, stating that the discussion helped clear up some misconceptions regarding the motivation behind the Fisher information matrix. The reviewer felt that the paper overall could improve in clarity, e.g. defining notation more explicitly, and other reviewers had clarifying questions that, while addressed in the discussion, should be incorporated into a future revision. To address concerns about novelty, the reviewers suggest providing a deeper discussion of context and related work, e.g., emphasizing / clarifying the gradient-free contribution when discussing both the connection with NGD and to expand the discussion of the broader VI literature, which has primarily been focused on the case where the scores of the target are available.  

I recommend accepting this paper, provided the authors address the remaining reviewers' comments.
Decision: Accept (poster)