Bayesian Online Natural Gradient (BONG)

Authors: matt  jones
peter  chang
kevin  murphy
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 07:53:08

Abstract: We propose a novel approach to  sequential Bayesian inference based on variational Bayes (VB).
  The key insight is that,
  in the online setting,
  we do not need to add the KL term to regularize to the prior (which comes from the posterior at the previous timestep);
  instead we can optimize just the expected log-likelihood,
  performing a single step of natural gradient descent
  starting at the prior predictive.
  We prove this method
 recovers exact Bayesian inference 
  if the model is conjugate.
 We also  show how to compute an
 efficient deterministic 
 approximation to the VB objective,
  as well as our simplified objective,
 when the variational distribution is
 Gaussian or a sub-family, including the case of
 a diagonal plus low-rank
precision matrix.
We show empirically that our
method outperforms other online VB methods
  in the  non-conjugate setting,
  such as online learning for neural networks,
  especially when controlling for computational costs.

Main paper: 
 introduction: 
 Bayesian methods for neural network (NN) training aim to minimize the Kullback-Leibler divergence between true and estimated posterior distributions. This is equivalent to minimizing the variational loss (or negative ELBO)EQUATIONHere θ are the network parameters, ψ are the variational parameters of the approximate posterior q ψ (θ), D is the training dataset, and p 0 (θ) is the prior. The two terms in the variational loss correspond to data fit and regularization to the prior, the latter being analogous to a regularizer r(θ) = -log p 0 (θ) in traditional point estimation methods like SGD.An important set of approaches learns the variational parameters by gradient descent on L(ψ) [Blundell et al., 2015] . More recently Khan and colleagues [Khan et al., 2018b , Khan and Rue, 2023 , Shen et al., 2024] have proposed using the natural gradient F -1 ψ ∇ ψ L(ψ) where F ψ is the Fisher information matrix of the variational family evaluated at q ψ . Natural gradient descent (NGD) is often more efficient than vanilla GD because it accounts for the intrinsic geometry of the variational family [Amari, 1998] . Khan and Rue [2023] call this approach the "Bayesian Learning Rule" or BLR. Using various choices for the variational distribution, generalized losses replacing negative log-likelihood, and other approximations, they reproduce many standard optimization methods such as Adam, and derive new ones.We study Bayesian NN optimization in online learning, where the data are observed in sequence, D t = {(x k , y k ) t k=1 }, and the algorithm maintains an approximate posterior q ψt (θ t ) ≈ p(θ t |D t ), which it updates at each step. Fast updates (in terms of both computational speed and statistical efficiency) are critical for many online learning applications [Zhang et al., 2024] . To allow for nonstationarity in the datastream, we include a time index on θ t , to represent that the parameters may change over time, as is standard for approaches based on state-space models and the extended Kalman filter (see e.g., [Sarkka and Svensson, 2023] ). The belief state is updated recursively using the prior q ψ t|t-1 derived from the previous time step so that the variational loss becomes L(ψ t ) = -E θt∼q ψ t [log p(y t |x t , θ t )] + D KL q ψt |q ψ t|t-1(2)One option for this online learning problem is to apply NGD on L(ψ t ) at each time step, iterating until ψ t converges before consuming the next observation. Our first contribution is a proposal for skipping this inner loop by (a) performing a single natural gradient step with unit learning rate and (b) omitting the D KL term in Eq. ( 2) so that learning is based only on expected loglikelihood:EQUATIONThese two modifications work together: instead of regularizing toward the prior explicitly using D KL q ψt |q ψ t|t-1 , we do so implicitly by using ψ t|t-1 as the starting point of our single natural gradient step. This may appear as a heuristic but we prove in Proposition 4.1 that it yields exact Bayesian inference when q ψ and p (y|x, θ) are conjugate and q ψ is an exponential family with natural parameter ψ. Thus our proposed update can be viewed as a relaxation of the Bayesian update to the non-conjugate variational case. As is common in work on variational inference, we view the result for the conjugate case as a motivating foundation that ensures our method is exact in certain simple settings. The experiments reported in Section 5 and Appendix B complement the theory by showing our method also works well in more general settings. We call Eq. ( 3) the Bayesian online natural gradient (BONG).Our second contribution concerns ways of computing the expectation in Eqs.(1) to (3). This is intractable for NNs, even for variational distributions that are easy to compute, since the likelihood takes the form p(y t |x t , θ t ) = p(y t |f (x t , θ t )) with f (x t , θ t ), representing the function computed by the network, is a complex, nonlinear function of θ t . Many previous approaches have approximated the expected loglikelihood by sampling methods which add variance and computation time depending on the number of samples [Blundell et al., 2015 , Shen et al., 2024] . We propose a deterministic, closed-form update that applies when the variational distribution is Gaussian (or a sub-family) and the likelihood is an exponential family with natural parameter f (x t , θ t ) and mean parameter h(x t , θ t ) (e.g., for classification, f returns the vector of class logits, h returns class probabilities, and h = softmax(f )). This update can be derived in two equivalent ways. First, we use a local linear approximation of the network h(x t , θ t ) ≈ ht (θ t ) [Immer et al., 2021a ] and a Gaussian approximation of the likelihood N (y t | ht (θ t ), R t ) [Ollivier, 2018 , Tronarp et al., 2018] . Under these assumptions the expectation in Eq. ( 3) can be calculated analytically. Alternatively, we use a different linear approximation f (x t , θ t ) ≈ ft (θ t ) and a delta approximation q ψ t|t-1 (θ t ) ≈ δ µ t|t-1 (θ t ) where µ t|t-1 = E q ψ t|t-1 [θ t ] is the prior mean, so that the expectation in Eq. ( 3) is replaced by a plugin prediction. The linear(h)-Gaussian approximation is previously known but the linear(f )-delta approximation is new, and we prove in Proposition 4.2 that they yield the same update, which we call linearized BONG, or BONG-LIN. Finally, we discuss different ways of approximating the Hessian of the objective, which is needed for NGD.Our BONG framework unifies several existing methods for Bayesian online learning, and it offers new algorithms based on alternative variational families or parameterizations. We define a large space of methods by combining 4 different update rules with 4 different ways of computing the relevant expected gradients and Hessians and 3 different variational families (Gaussians with full, diagonal, and diagonal-plus-low rank precision matrices). We conduct experiments systematically testing how these factors affect performance. We find support for all three principles of our approach-NGD, implicit regularization to the prior, and linearization-in terms of both statistical and computational efficiency. Code for our experiments is available at https://github.com/petergchang/bong/. 
related work: 
 Variational inference approximates the Bayesian posterior from within some suitable family in a way that bypasses the normalization term [Zellner, 1988 , Jordan et al., 1999] . A common choice for the variational family is a Gaussian. For online learning, the exact update equations for Gaussian variational filtering are given by the RVGA method of [Lambert et al., 2021] . This update is implicit but can be approximated by an explicit RVGA update which we show arises as a special case of BONG.Most applications of Gaussian VI use a mean-field approximation defined by diagonal covariance, which scales linearly with model size. More expressive but still linear in the model size are methods that express the covariance [Tomczak et al., 2020] or precision [Mishkin et al., 2018 , Lambert et al., 2023 , Chang et al., 2023] as a sum of diagonal and low rank matrices (DLR). In this paper, we consider variational families defined by full covariance, diagonal covariance, and DLR covariance.For NNs and other complicated models, even the variational approximation can be intractable, so methods have been developed to approximately minimize the VI loss. Bayes by backprop (BBB) [Blundell et al., 2015 ] learns a variational distribution on NN weights by iterated GD on the VI loss of Eq. ( 1). They focus on mean-field Gaussian approximations but the approach also applies to other variational families. Here we adapt BBB to online learning to compare to our methods.The Bayesian learning rule (BLR) replaces BBB's GD with NGD [Khan and Rue, 2023] . Several variants of BLR have been developed including VON and VOGN for a mean-field Gaussian prior [Khan et al., 2018b] and SLANG for a DLR Gaussian [Mishkin et al., 2018] . BLR has also been used to derive versions of many classic optimizers including SGD, RMSprop and Adam [Khan et al., 2018a , Khan and Rue, 2023 , Lin et al., 2024 , Shen et al., 2024] . Although BLR has been applied to online learning, we are particularly interested in Bayesian filtering including in nonstationary environments, where observations must be processed one at time and updates are based on the posterior from the previous step, often in conjunction with parameter dynamics. We therefore develop filtering versions of BLR to compare to BONG, some of which reduce to VON, VOGN and SLANG in the batch setting, while others are novel. We also note BLR is a mature theory including several clever tricks we have not yet incorporated into our framework. Khan and Rue [2023] observe that conjugate updating is equivalent to one step of BLR with learning rate 1. This is similar to our Proposition 4.1 except that BLR retains the KL term in the variational loss. BLR and BONG agree in this case because the gradient of the KL is zero on BLR's first iteration: ∇ ψ=ψ t|t-1 D KL q ψ |q ψ t|t-1 = 0. Therefore BONG can be seen as a special case of BLR with one update step per observation and learning rate 1. Our contribution is to recognize that doing a single update step allows the KL term to be dropped entirely, yielding a substantially simpler algorithm which our experiments show also performs better.While BLR allows alternative losses in place of the NLL in Eq. ( 2), we can also replacing the KL term with other divergences [Knoblauch et al., 2022] . Our approach fits within that "generalized VB" framework in that it drops the divergence altogether. Our approach of implicitly regularizing to the prior using a single NGD step is also similar to the implicit MAP filter of [Bencomo et al., 2023] which performs truncated GD from the prior mode. The principal difference is they perform GD on model parameters (θ t ) while we do NGD on the variational parameters (ψ t ). Thus BONG maintains a full prior and posterior while IMAP is more concerned with how the choice of optimizer can substitute for explicit tracking of covariance.We show two other ways to derive the BONG update in Appendix D, one of which is to replace the expected NLL in Eq. ( 2) with a linear approximation and solve the resulting equation exactly. Several past works have taken this approach, arriving at updates similar to ours. Chérief-Abdellatif et al.[2019] study streaming variational Bayes and propose solving Eq. ( 2) with a linearized expected NLL. When the variational family is an exponential family their update becomes NGD [Khan and Lin, 2017] and matches the BONG update. Hoeven et al. [2018] show how mirror descent can be derived as a special case of Exponential Weights [Littlestone and Warmuth, 1994] , which is closely related to Bayesian updating. The resulting algorithm is similar to BONG and follows from linearizing the NLL instead of expected NLL, with an additional delta assumption at the prior mean. Lyu and Tsang [2021] study relaxed block-box optimization where the objective is arg minψ E x∼q ψ [f (x)]for some target function f . They use a mirror descent formulation with linearized expected loss and KL regularizer and show the resulting update is NGD on expected loss, formally equivalent to our BONG update. From the perspective of this prior work, our contribution is to express the BONG update simply as NGD on the expected NLL, motivated by replacing the KL with implicit regularization, and to show how this yields a variety of known and novel algorithms for Bayesian filtering.EKF applications to NNs apply Bayesian filtering using a local linear approximation of the network, leading to simple closed form updates [Singhal and Wu, 1989, Puskorius and Feldkamp, 1991] .The classic EKF assumes a Gaussian observation distribution but it has been extended to other exponential families (e.g. for classification) by matching the mean and covariance in what we call the conditional moments EKF (CM-EKF) [Ollivier, 2018 , Tronarp et al., 2018] . Applying a KL projection to diagonal covariance yields the variational diagonal EKF (VD-EKF) [Chang et al., 2022] .Alternatively, projecting to diagonal plus low rank precision using SVD gives LO-FI [Chang et al., 2023] . We derive all these methods as special cases of BONG-LIN. Further developments in this direction include the method of [Titsias et al., 2024] which does Bayesian filtering on only the final weight layer, and WoLF [Duran-Martin et al., 2024] which achieves robustness to outliers through data-dependent weighting of the loglikelihood. 
background: 
 We study online supervised learning where the agent receives input x t ∈ R D and observation y t ∈ R C on each time step, which it aims to model with a function f t (θ t ) = f (x t , θ t ) such as a NN with weights θ t ∈ R P . The predictions for y t are given by some observation distribution p(y t |f t (θ t )).For example, f may compute the mean for regression or the class logits for classification.We work in a Bayesian framework where the agent maintains an approximate posterior distribution over θ t after observing dataD t = {(x k , y k ) t k=1 }.The filtering posterior q ψt (θ t ) ≈ p(θ t |D t ) is approximated within some parametric family indexed by the variational parameter ψ t . We allow for nonstationarity by assuming θ changes over time according to some dynamic model p(θ t |θ t-1 ). By pushing the posterior from step t -1 through the dynamics we obtain a prior for step t given by q ψ t|t-1 (θ t ) ≈ p(θ t |D t-1 ). For example suppose the variational posterior from the previous step is Gaussian, q ψt-1 (θ t-1 ) = N (θ t-1 |µ t-1 , Σ t-1 ), and the dynamics model is an Ornstein-Uhlenbeck process, as proposed in prior work [Kurle et al., 2020 , Titsias et al., 2024] to handle non-stationarity, i.e., the dynamics model has the formθ t ∼ N (γ t θ t-1 + (1 -γ t )µ 0 , Q t ), where Q t = (1 -γ 2t )Σ 0 is the covariance of the noise process, 0 ≤ γ t ≤ 1 is the degree of drift, and p(θ 0 ) = N (µ 0 , Σ 0 ) is the prior. In this case, the parameters of the prior predictive distribution areµ t|t-1 = γ t µ t-1 +(1-γ t )µ 0 and Σ t|t-1 = γ 2 t Σ t-1 + Q t .In general the predict step may require approximation to stay in the variational family (e.g., if the dynamics are nonlinear). In this paper, our focus is the update step from ψ t|t-1 to ψ t upon observing (x t , y t ), so for simplicity we assume constant (static) parameters, i.e., p(θ t |θ -1 ) = δ(θ t -θ t-1 ) (equivalently γ t = 1), so ψ t|t-1 = ψ t-1 ; however, our method can trivially handle non-stationary parameters.Variational inference seeks an approximate posterior that minimizes the KL divergence from the exact Bayesian update from the prior. In the online setting this becomesEQUATIONwhere L t is the online VI loss defined in Eq. ( 2), and the normalization term Z t (which depends on x t ) drops out as an additive constant. Our goal is an efficient approximate solution to this variational optimization problem.We will sometimes assume the variational posterior q ψ is an exponential family distribution with natural parameter ψ so that q ψt (θ t ) = exp (ψ t T (θ t ) -Φ(ψ t )), with log-partition function Φ and sufficient statistics T (θ t ). Assuming Φ is strictly convex (which holds in the cases we study) there is a bijection between ψ t and the dual (or expectation) parameterρ t = E θt∼q ψ t [T (θ t )].Classical thermodynamic identities imply that the Fisher information matrix has the form F ψt = ∂ρ t /∂ψ t . This has important implications for NGD on exponential families [Khan and Rue, 2023] because it implies that for any function defined on the variational parameter space the natural gradient wrt natural parameters ψ t is the regular gradient wrt the dual parameters ρ t , i.e., F -1 ψt ∇ ψt = ∇ ρt . 
methods: 
 We propose to approximate the variational optimization problem in Eq. (4) using the BONG update in Eq. (3). When q ψ is an exponential family, the fact that the natural gradient wrt the natural parameters ψ t is the regular gradient wrt the dual parameters ρ t implies an equivalent mirror descent form (see Appendix D for further analysis of BONG from the MD perspective):EQUATIONThis is NGD with unit learning rate on the variational loss in Eq. ( 2) but ignoring the D KL q ψ |q ψ t|t-1 term. In this section we first prove this method is optimal when the model is conjugate and then describe extensions to more complex cases of practical interest.Our approach is motivated by the following result which states that BONG matches exact Bayesian inference when the variational distribution and the likelihood are conjugate exponential families: Proposition 4.1. Let the observation distribution (likelihood) be an exponential family with natural parameter θ t (where T l (y t ) = y t is the sufficient statistics for the likelihood and A(θ t ) is the log-partition function) p t (y t |θ t ) = exp (θ t y t -A(θ t ) -b(y t )) (6) and let the prior be the conjugate exponential familyq ψ t|t-1 (θ t ) = exp ψ t|t-1 T (θ t ) -Φ(ψ t|t-1 ) (7) with T (θ t ) = [θ t ; -A(θ t )].Then the exact Bayesian update agrees with Eq. (5).The proof is in Appendix C. Writing the natural parameters of the prior asψ t|t-1 = [χ t|t-1 ; ν t|t-1 ],we show the Bayesian update and BONG both yield χ t = χ t|t-1 + y t and ν t = ν t|t-1 + 1. Intuitively, we are just accumulating a sum of the observed sufficient statistics, and a counter of the sample size (number of observations seen so far).In practical settings the conjugacy assumption of Proposition 4.1 will not be met, so Eqs. ( 3) and ( 5) will only approximate the Bayesian update. In this paper we restrict to Gaussian variational families. We refer to the unrestricted case as FC (full covariance), defined by the variational distributionEQUATION) where Σ t|t-1 can be any positive semi-definite (PSD) matrix. The natural and dual parameters are ψ = (Σ -1 µ, -1 2 vec(Σ -1 )) and ρ = (µ, vec(µµ + Σ)). Appendix E.1.1 shows that Eq. ( 5) translated back to (µ, Σ) gives the following BONG update for the FC case:EQUATIONEQUATION)which matches the explicit update in the RVGA method of [Lambert et al., 2021] .The integrals over the prior q ψ t|t-1 in Eqs. ( 9) and ( 10) are generally intractable and must be approximated. One option is to use Monte Carlo, in what we call BONG-MC. Given M independent samples θ(m) t ∼ q ψ t|t-1 , we estimate the expected gradientg t = E θt∼q ψ t|t-1 [∇ θt log p(y t |f t (θ t ))]and expected Hessian G t = E θt∼q ψ t|t-1 ∇ 2 θt log p (y t |f t (θ t )) as the empirical meansEQUATIONWe use G MC-HESS only for small models. Otherwise we use empirical Fisher (Section 4.5).As an alternative to BONG-MC, we propose a linear approximation we call BONG-LIN that yields a deterministic and closed-form update. Assume the likelihood is an exponential family as in Proposition 4.1 but with natural parameter predicted by some functionEQUATIONWe also define the dual (moment) parameter of the likelihood as h t (θ t ) = E [y t |f t (θ t )]. In a NN, f t and h t are related by the final response layer. For example in classification f t and h t give the class logits and probabilities, with h t (θ t ) = softmax(f t (θ t )), with y t being the one-hot encoding.We now define two methods for approximating the expected gradient g t and expected Hessian G t , based on linearizing the predictive model at the prior mean µ t|t-1 in terms of either f t (θ t ) or h t (θ t ), and then prove their equivalence.The linear(h)-Gaussian approximation [Ollivier, 2018 , Tronarp et al., 2018] EQUATIONEQUATIONand approximates the likelihood by a Gaussian with variance based at µ t|t-1EQUATIONThe linear(f )-delta approximation linearizes f t (θ t ) and maintains the original exponential family likelihood distribution in Eq. ( 13)EQUATIONEQUATIONIt also uses a plug-in approximation that replaces q ψ t|t-1 (θ t ) with a point mass δ µ t|t-1 (θ t ) so that the expected gradient and Hessian are approximated by their values at the prior mean, i.e., ∇ θt=µ t|t-1 log pLD t (y t |θ t ) and ∇ 2 θt=µ t|t-1 log pLD t (y t |θ t ), rather than being sampled. Proposition 4.2. Under a Gaussian variational distribution, the linear(h)-Gaussian and linear(f )delta approximations yield the same values for the expected gradient and HessianEQUATIONSee Appendix C for the proof. The main idea for the g LIN t part is that the linear-Gaussian assumptions make the gradient linear in θ t so the expected gradient equals the gradient at the mean. The main idea for the G LIN-HESS t part is that eliminating the Hessian of the NN requires different linearizing assumptions for the Gaussian and delta approximations, and the remaining nonlinear terms (from the log-likelihood in Eq. ( 13)) agree because of the property of exponential families that the Hessian of the log-partition A equals the conditional variance R t .Applying Proposition 4.2 to Eqs. ( 9) and ( 10) gives the BONG-LIN update for a FC Gaussian priorEQUATIONwhere K t is the Kalman gain matrix (see Appendix E.1.2). This matches the CM-EKF [Tronarp et al., 2018 , Ollivier, 2018] . The methods in Sections 4.3 and 4.4 require explicitly computing the Hessian of the loss (MC-HESS) or the Jacobian of the network (LIN-HESS). These are too expensive for large models or high-dimensional observations. Instead we can use an empirical Fisher approximation that replaces the Hessian with the outer product of the gradient (see e.g, [Martens, 2020] ).For the MC-EF variant, we make the following approximation:G MC-EF t = - 1 M Ĝ(1:M) t Ĝ(1:M) t (26) where Ĝ(1:M) t = [ĝ (1) t , . . . , ĝ(M) t] is the P × M matrix of gradients from the MC samples.We can also consider a similar approach for the LIN-EF variant that is Jacobian-free and sampling-free. Note that if ŷt were the true value of E [y t |x t ] (i.e., if the model were correct) then we would haveE [(y t -ŷt )(y t -ŷt ) ] = R t , implying E [g LIN t (g LIN t ) ] = -G LIN-HESS t. This suggests usingEQUATIONwhere Eq. ( 29) is the EF approximation to Eq. ( 22).A more accurate EF approximation is possible by sampling virtual observations ỹt from p(•|f t ( θt (m) )) or p(•|f t (µ t|t-1 )) and using them for the gradients in Eq. ( 26) or Eq. ( 29) (respectively) [Martens, 2020 , Kunstner et al., 2020] . However, in our experiments we use the actual observations y t which is faster and follows previous work (e.g., [Khan et al., 2018b] ).Loss UpdateBONG E [NLL] NGD(I = 1) BOG E [NLL] GD(I = 1) BLR VI NGD(I ≥ 1) BBB VI GD(I ≥ 1)Table 2 : The 4 update algorithms.In addition to the four ways of approximating the expected Hessian (summarized in Table 1 ), we also consider four variants of BONG, based on what kind of loss we optimize and what kind of update we perform, as we describe below. See Table 2 for a summary.BONG (Bayesian online natural gradient) performs one step of NGD on the expected log-likelihood. We set learning rate to α t = 1 since this is optimal for conjugate models. The update (for an exponential variational family) is as in Eq. ( 5):EQUATIONBOG (Bayesian online gradient) performs one step of GD (instead of NGD) on the expected loglikelihood. We include a learning rate α because GD does not have the scale-invariance of NGD:EQUATIONBLR (Bayesian learning rule, [Khan and Rue, 2023] ) uses NGD (like BONG) but optimizes the VI loss using multiple iterations, instead of optimizing the expected NLL with a single step. When modified to the online setting, BLR starts an inner loop at each time step with ψ t,0 = ψ t|t-1 and iteratesEQUATIONFor an exponential variational family this can be written in mirror descent formEQUATIONBBB (Bayes By Backprop, [Blundell et al., 2015] ) is like BLR but uses GD instead of NGD. When adapted to online learning, it starts each time step at ψ t,0 = ψ t|t-1 and iterates with GD:ψ t,i = ψ t,i-1 + α t ∇ ψt,i-1 E θt∼q ψ t,i-1 [log p(y t |f t (θ t ))] -D KL q ψt,i-1 |q ψ t|t-1 (34) Variational Family Update Hessian Full Diag DLR BONG MC-EF O(M P 2 ) [RVGA ] O(M P ) O((R + M ) 2 P ) BLR MC-EF O(IP 3 ) O(IM P ) [VON ] O(I(R + M ) 2 P ) [SLANG ] BOG MC-EF O(P 3 ) O(M P ) O(RM P ) BBB MC-EF O(IP 3 ) O(IM P ) [BBB ] O(IR(R + M )P ) BONG LIN-HESS O(CP 2 ) [CM-EKF ] O(C 2 P ) [VD-EKF ] O((R + C) 2 P ) [LO-FI ] BLR LIN-HESS O(IP 3 ) O(IC 2 P ) O(I(2R + C) 2 P ) BOG LIN-HESS O(P 3 ) O(C 2 P ) O(C(C + R)P ) BBB LIN-HESS O(IP 3 ) O(IC 2 P ) O(I(C + R)RP ) BONG LIN-EF O(P 2 ) O(P ) O(R 2 P ) BLR LIN-EF O(IP 3 ) O(IP ) O(IR 2 P ) BOG LIN-EF O(P 3 ) O(P ) O(RP ) BBB LIN-EF O(IP 3 ) O(IP ) O(IR 2 P )Table 3 : Time complexity of the algorithms. P : params, C: observation dim, M : MC samples, I: iterations, R: DLR rank. We assume P {I, C, R, M } so display only the terms of leading order in P . Time complexities for MC-HESS algorithms (not shown) are always at least as great as for the corresponding MC-EF. Full (full covariance) and Diag (diagonal covariance) columns indicate natural parameters; corresponding algorithms using moment parameters have the same complexities except BOG-FC_MOM which is O(M P 2 ) for MC-EF, O(CP 2 ) for LIN-HESS, and O(P 2 ) for LIN-EF.[Names] correspond to the following existing methods (or variants thereof) in the literature: RVGA: [Lambert et al., 2021] (explicit update version); VON: [Khan et al., 2018b] (modified for online); SLANG: [Mishkin et al., 2018] (modified for online); BBB: [Blundell et al., 2015] (modified for online and uses moment parameters); CM-EKF: [Ollivier, 2018 , Tronarp et al., 2018] ; VD-EKF: [Chang et al., 2022] ; LO-FI: [Chang et al., 2023] .We investigate five variational families for the posterior distribution: (1) FC Gaussian using natural parameters ψ = (Σ -1 µ, -1 2 Σ -1 ), (2) FC Gaussian using central moment parameters ψ = (µ, Σ), (3) diagonal Gaussian using natural parameters ψ = (σ -2 µ, -1 2 σ -2 ) (using elementwise exponents and products), (4) diagonal Gaussian using central moment parameters ψ = (µ, σ 2 ), and (5) DLR Gaussian with parameters ψ = (µ, Υ, W) and precision Σ -1 = Υ + WW where Υ ∈ R P ×P is diagonal and W ∈ R P ×R with R P . The moment parameterizations are included to test the importance of using natural parameters per Proposition 4.1. The diagonal family allows learning of large models because it scales linearly in the model size P . DLR also scales linearly but is more expressive than diagonal, maintaining some of the correlation information between parameters that is lost in the mean field (diagonal) approximation [Lambert et al., 2023 , Mishkin et al., 2018 , Chang et al., 2023] .Optimizing the BONG objective wrt (µ, Υ, W) using NGD methods is challenging because the Fisher information matrix in this parameterization cannot be efficiently inverted. Instead we first derive the update wrt the FC natural parameters (leveraging the fact that the prior Σ -1 t|t-1 is DLR to make this efficient), and then use SVD to project the posterior precision back to low-rank form, following our prior LO-FI work [Chang et al., 2023] . However, if we omit the Fisher preconditioning matrix and use GD as in BOG and BBB, we can directly optimize the objective wrt (µ, Υ, W) (see Appendix E.5).Crossing the four algorithms in Table 2 , the four methods of approximating the Hessian in Table 1 , and the five variational families yields 80 algorithms. Table 3 shows the 36 based on the three tractable Hessian approximations and the three variational families that use natural parameters. Update equations for all the algorithms are derived in Appendix E. Pseudocode is given in Appendix A. 
experiments: 
 This section presents our primary experimental results. These are based on MNIST (D = 784, N train = 60k, N test = 10k, C = 10 classes) [LeCun et al., 2010] . See Appendix B for more details on these experiments, and more results on MNIST and other datasets. We focus on training on a prefix of the first T = 2000 examples from each dataset, since our main interest is in online learning from potentially nonstationary distributions, where rapid adaptation of a model in response to a small number of new data points is critical.Our primary evaluation objective is the negative log predictive density (NLPD) of the test set as a function of the number of training points observed so far.1 It is defined as NLPD t = -1 Ntest i∈D test log p(y i |f (x i , θ t ))q ψt (θ t )dθ t . We approximate this integral in two main ways:(1) using Monte Carlo sampling2 ,or (2) using a plugin approximation, where we replace the posterior q ψt (θ t ) with a delta function centered at the mean, δ(θ t -µ t ).For methods that require a learning rate (i.e., all methods except BONG), we optimize it wrt mid-way or final performance on a holdout validation set, using Bayesian optimization on NLL. All methods require specifying the prior belief state, p(θ 0 ) = N (µ 0 , Σ 0 = σ 2 0 I). We optimize over σ 2 0 and sample µ 0 from a standard NN initializer. As Hessian approximations, we use MC-EF with M = 100 samples, as well as the deterministic approximations LIN-HESS and LIN-EF. (In the appendix, we also study MC-HESS but find that it works very poorly, even with M = 1000.)In Fig. 1 we compare the 4 main algorithms using DLR family with rank R = 10. We apply these to a CNN with two convolutional layers (each with 16 features and a (5,5) kernel), followed by two linear layers (one with 64 features and the final one with 10 features), for a total of 57,722 parameters. Shaded regions in these and all other plots indicate ±1 SE, based on 5 independent trials randomly varying in the prior mean µ 0 , data ordering, and MC sampling. From this figure (and additional results in Appendix B) we conclude • Linearization helps: LIN-HESS and LIN-EF both outperform the MC variants.• NGD helps: BONG outperforms BOG.• Implicit regularization helps: BONG outperforms BLR.• The LIN-HESS approximation outperforms LIN-EF, at least for BONG.• BBB generally does poorly.• The BONG posterior predictive (using LIN-HESS) is slightly better calibrated than BOG, and both are much better than BLR and BBB, at least for small sample sizes, as shown in Fig. 7 . • The plugin posterior predictive is similar to the Lin-MC predictive (see Footnote 2), and both are generally much better than the simple MC predictive, as shown in Fig. 5 .In Fig. 2 we compare BONG using different variational families, and conclude• DLR-10 outperforms DLR-1, which is similar to diagonal (except when using BONG-LIN-EF, where DLR-1 is worse than diagonal). Also, we find (in results not reported here) that rank 5-10 often gives results as good as FC, but is much cheaper. • Both natural and moment parameterizations for the diagonal representation perform comparably to each other, although with LIN-EF, the moment parameterization can be numerically unstable. 3Finally, in Fig. 3 we report the runtimes from these experiments and conclude• One-step methods (BONG and BOG) are faster than iterative methods (BLR and BBB), as expected.• Linearized methods (LIN-HESS and LIN-EF) are faster than MC methods (MC-EF).• EF methods (LIN-EF) are a bit faster than methods that compute the Hessian exactly (LIN-HESS), especially for diagonal family. (This speedup is larger when the output dimensionality C is big.) -dlr10-hess_lin [sec:2.7]  bong-dlr10-hess_lin [sec:5.2]  bbb-dlr10-ef_lin [sec:31.3]  blr-dlr10-ef_lin [sec:28.6]  bog-dlr10-ef_lin [sec:2.7]  bong-dlr10-ef_lin [sec:3.8] Figure 1 : Performance on MNIST using Lin-MC posterior predictive, where the posterior is computed using BONG, BOG, BBB and BLR and the 3 tractable Hessian approximations with DLR-10 variational family. 0-diag-ef_mc100 [sec:3.2] bong-diag-hess_lin [sec:3.3] bong-diag-ef_lin [sec:1.7] bong-diag_mom-ef_mc100 [sec:3.1] bong-diag_mom-hess_lin [sec:2.2] bong-diag_mom-ef_lin [sec:1.7] bong-dlr1-ef_mc100 [sec:14.5] bong-dlr10-ef_mc100 [sec:13.1] bong-dlr1-hess_lin [sec:4.8] bong-dlr1-ef_lin [sec:3.1] bong-dlr10-hess_lin [sec:5.0] bong-dlr10-ef_lin [sec:3.8](lin-MC) bong-diag-ef_mc100 [sec:3.2] bong-diag-hess_lin [sec:3.3] bong-diag-ef_lin [sec:1.7] bong-diag_mom-ef_mc100 [sec:3.1] bong-diag_mom-hess_lin [sec:2.2] bong-diag_mom-ef_lin [sec:1.7] bong-dlr1-ef_mc100 [sec:14.5] bong-dlr10-ef_mc100 [sec:13.1] bong-dlr1-hess_lin [sec:4.8] bong-dlr1-ef_lin [sec:3.1] bong-dlr10-hess_lin [sec:5.0] bong-dlr10-ef_lin [sec:3.8]Figure 2 : Performance on MNIST using Lin-MC posterior predictive, where the posterior is computed using BONG with different variational families, namely diagonal (natural and moment), DLR-1, DLR-10. 
conclusions, limitations and future work: 
 Our experiment results show benefits of BONG's three main principles: NGD, implicit regularization to the prior, and linearization. The clear winner across datasets and variational families is BONG-LIN-HESS, which embodies all three principles. BLR-LIN-HESS nearly matches its performance but is much slower. Several of the best-performing algorithms are previously known (notably CM-EKF and LO-FI) but we explain these results within a systematic theory that also offers new methods (including BLR-LIN-HESS).BONG is motivated by Proposition 4.1 which applies only in the idealized setting of a conjugate prior. Nevertheless we find it performs well in non-conjugate settings. On the other hand our experiments are based on relatively small models and datasets. It will be important to test how our methods scale up, especially using the promising DLR representation. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We report running times. Experiments all run on any standard GPU/TPU so further details are not necessary.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 as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: We have read the Code of Ethics and affirm our research adheres to all points therein.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [No] Justification: We believe there are no direct societal impacts (positive or negative) beyond the generic impacts of foundational work on training neural networks.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: The paper poses no such risks. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere 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 not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: We provide citations and URLs for the two datasets we use in Section 5. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a 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 (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If 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 is the documentation provided alongside the assets? Answer: [NA] Justification: There are no new data. As noted above, we will post our model code after the blind review period. Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used.  
 
 Appendix: Algorithms 1 to 7 give pseudocode for applying the methods we study. For the predict step, we assume the dynamics model has the form θ t ∼ N (F t θ t-1 + b t , Q t ). The update step in Algorithm 3 calls one of four grad functions (Algorithms 4 to 7) that estimate the expected gradient and Hessian using either MC or linearization combined with either the direct Hessian (or Jacobian and observation covariance) or EF. The update step also calls an inner step function that implements BONG, BLR, BOG or BBB on some variational family corresponding to the encoding of ψ (not shown). In practice the grad-fn and inner-step-fn are not as cleanly separated because the full matrix Ḡt,i is not passed between them except when the variational family is FC. When the family is diagonal, grad-fn only needs to return diag( Ḡt,i ). When grad-fn uses EF, it only needs to return Ĝ(1:M)t,i(grad-MC-EF) or g LIN t,i (grad-LIN-EF) and inner-step-fn will implicitly use the outer product of this output as Ḡt,1 . Finally, note the expressions for g LIN t,i in Algorithms 6 and 7 are equivalent ways of computing the same quantity, as explained after Eq. ( 28).for t = 1 : ∞ do ψ t|t-1 = predict(ψ t-1 ) ψ t = update(ψ t|t-1 , x t , y t ) end Algorithm 1: Main loop. def predict(ψ t-1 = (µ t-1 , Σ t-1 )): µ t|t-1 = F t µ t-1 + b t Σ t|t-1 = F t Σ t-1 F t + Q t Return ψ t|t-1 = (µ t|t-1 , Σ t|t-1 ) Algorithm 2: Predict step.def update(ψ t|t-1 , x t , y t , f (), A(), grad-fn, inner-step-fn, α, I, M ):ψ t,0 = ψ t|t-1 f t (θ t ) = f (x t , θ t ) h t (θ t ) = E[y t |x t , θ t ] = ∇ η=ft(θt) A(η) V t (θ t ) = Cov[y t |x t , θ t ] = ∇ 2 η=ft(θt) A(η) t (θ t ) = log p(y t |x t , θ t ) = log p(y t |f t (θ t )) for i = 1 : I do (ḡ t,i , Ḡt,i ) = grad-fn(ψ t,i-1 , t , h t , V t , M ) ψ t,i = inner-step-fn(ψ t|t-1 , ψ t,i-1 , ḡt,i , Ḡt,i , α) end Return ψ t,IAlgorithm 3: Update step. The inner-step-fn is BONG, BLR, BOG or BBB (not shown).In this section, we give a more thorough set of experimental results.The running times of the methods for the experiments in Figs. 1 and 2 , where we fit a CNN to MNIST, are shown in Fig. 3 .The running times of the methods for the FC and DLR case, where we fit an MLP to a synthetic regression dataset, are shown in Fig. 4 . The slower speed of BLR (even with I = 1) relative to BONG is at least partly attributable to the fact that BLR must compute the SVD of a larger matrix (see Appendices E.5.3 and E.5.4).def grad-MC-HESS(ψt,i-1 , t , h t = [], V t = [], M ): for m = 1 : M do θ(m) t,i ∼ q ψt,i-1 (θ) ĝ(m) t,i = ∇ θt= θ(m) t,i t (θ t ) Ĝ(m) t,i = ∇ 2 θt= θ(m) t,i t (θ t ) end g MC t,i = 1 M M m=1 ĝ(m) t,i G MC-HESS t,i = 1 M M m=1 Ĝ(m) t,i Return (g MC t,i , G MC-HESS t,i ) Algorithm 4: MC gradient/Hessian estimator def grad-MC-EF(ψ t,i-1 , t , h t = [], V t = [], M ): for m = 1 : M do θ(m) t,i ∼ q ψt,i-1 (θ) ĝ(m) t,i = ∇ θt= θ(m) t,i t (θ t ) end g MC t,i = 1 M M m=1 ĝ(m) t,i Ĝ(1:M) t,i = [ĝ (1) t,i , . . . , ĝ(M) t,i ] G MC-EF t,i = -1 M Ĝ(1:M) t,i Ĝ(1:M) t,i Return (g MC t,i , G MC-EF t,i) Algorithm 5: MC gradient/Hessian estimator with Empirical Fisher def grad-LIN-HESS(ψ t,i-1, t = [], h t , V t , M = []): µ t,i-1 = E[θ t |ψ t,i-1 ] ŷt,i = h t (µ t,i-1 ) H t,i = ∂ht ∂θt |θ=µt,i-1 R t,i = V t (µ t,i-1 ) g LIN t,i = H t,i R -1 t,i (y t -ŷt,i ) G LIN-HESS t,i = -H t,i R -1 t,i H t,i Return (g LIN t,i , G LIN-HESS t,i ) Algorithm 6: Linearized gradient/Hessian estimator def grad-LIN-EF(ψ t,i-1 , t = [], h t , V t , M = []): µ t,i-1 = E[θ t |ψ t,i-1 ] R t,i = V t (µ t,i-1 ) g LIN t,i = ∇ θt=µt,i-1 -1 2 (y t -h t (θ t )) R -1 t,i (y t -h t (θ t )) G LIN-EF t,i = -g LIN t,i g LIN t,i Return (g LIN t,i , G LIN-EF t,i) Algorithm 7: Linearized gradient/Hessian estimator with empirical Fisher -MC100-I1-LR0-EF0-Lin0 bong_fc-MC100-I1-LR0-EF0-Lin1 bong_fc-MC100-I1-LR0-EF1-Lin0 bog_fc-MC100-I1-LR0_001-EF0-Lin0 bog_fc-MC100-I1-LR0_001-EF0-Lin1 bog_fc-MC100-I1-LR0_001-EF1-Lin0 blr_fc-MC100-I10-LR0_001-EF0-Lin0 blr_fc-MC100-I10-LR0_001-EF0-Lin1 blr_fc-MC100-I10-LR0_001-EF1-Lin0 bbb_fc-MC100-I10-LR0_001-EF0-Lin1 bbb_fc-MC100-I10-LR0_001-EF1-Lin0-dlr10-EF1-MC100 bong-dlr10-lin bog-dlr10-EF1-MC100-LR0_001 bog-dlr10-lin-LR0_001 blr-dlr10-EF1-MC100-I1-LR0_001 blr-dlr10-lin-I1-LR0_001 bbb-dlr10-EF1-MC100-I1-LR0_001 bbb-dlr10-lin-I1-LR0_001 (a) (b)Here we report further metrics for the experiments in Figs. 1 and 2 . We show 3 approximations to the NLPD: plugin, MC, and Linearized MC. 4 For each of these approximations to the posterior predictive, we also measure the corresponding misclassification rate based on picking the most probable predicted class. Results are shown in Figs. 5 and 6 . We see that the plugin and lin-MC approximations are similar, and both are generally much better than standard MC.Finally, in Fig. 7 and Fig. 8 method is the most well-calibrated (in addition to exhibiting the strongest plugin and linearized-MC NLPD results), when compared to other DLR-10 methods as well as other BONG variants. We see that BONG and BOG are both well calibrated, as least when combined with LIN-HESS, with BONG have a slight edge, especially for small sample sizes, where there is more posterior uncertainty.In addition to MNIST, we report experiments on the SARCOS regression dataset (D = 22, N train = 44,484, N test = 4449, C = 1). This dataset derives from an inverse dynamics problem for a seven degrees-of-freedom SARCOS anthropomorphic robot arm. The task is to map from a 21-dimensional input space (7 joint positions, 7 joint velocities, 7 joint accelerations) to the corresponding 7 joint torques. Following Rasmussen and Williams [2006] , we pick a single target output dimension, so C = 1. The data is from https://gaussianprocess.org/gpml/data/.We use a small MLP of size 21-20-20-1, so there are P = 881 parameters. For optimizing learning rates for SARCOS, we use grid search on NLPD-PI. We fix the variance of the prior belief state to σ 2 0 = 1.0, which represents a mild degree of regularization. 5 We fix the observation variance to R t = 0.1 R, where R = Var(y 1:T ) is the maximum likelihood estimate based on the training sequence; we can think of this as a simple empirical Bayes approximation, and the factor of 0.1 accounts for the fact that the variance of the residuals from the learned model will be smaller than from the unconditional baseline. We focus on DLR approximation of rank 10. This gives similar results to full covariance, but is much faster to compute. We also focus on the plugin approximation to NLPD, since the MC approximation gives much worse results (not shown).In Fig. 9 we show the results of using the LIN-HESS approximation. For 1 iteration per step, we see that BONG and BLR are indistinguishable in performance, and BBB and BOG are also indistinguishable, but much worse. For 10 iterations per step, we see that BBB improves significantly, and approaches BONG and BLR. However, BLR and BBB are now about 10 times slower. (In practice, the slowdown is less than 10, due to constant factors of the implementation.) (Note that BONG and BOG always use a single iteration, so their performance does not change.)In Fig. 10 we show the results of using the MC-EF approximation with MC = 100 samples. The trends are similar to the LIN-HESS case. In particular, for I = 1, BONG and BLR are similar, with BONG having a slight edge; and for I = 10, BBB catches up with both BONG and BLR, with BOG always in last place. Finally, we see that the performance of MC-EF is slightly worse than LIN-HESS when I = 1, but catches up with I = 10. However, in larger scale experiments, we usually find that LIN-HESS is significantly better than MC-EF, even with I = 10. . Data: sarcos. Expt: bakeoff-it1 bbb-dlr10-LinHess-It1-LR0_05 [sec:5.9+-0.1 nlpd-pi@T=2000:3.50+-0.09] blr-dlr10-LinHess-It1-LR0_5 [sec:6.5+-0.0 nlpd-pi@T=2000:3.28+-0.21] bog-dlr10-LinHess-LR0_05 [sec:5.0+-0.0 nlpd-pi@T=2000:3.50+-0.09] bong-dlr10-LinHess [sec:6.1+-0.1 nlpd-pi@T=2000:3.32+-0.21] (a) 1 iteration. . Data: sarcos. Expt: bakeoff-it10 bbb-dlr10-LinHess-It10-LR0_01 [sec:14.6+-0.1 nlpd-pi@T=2000:3.60+-0.25] blr-dlr10-LinHess-It10-LR0_1 [sec:21.2+-0.2 nlpd-pi@T=2000:3.29+-0.14] bog-dlr10-LinHess-LR0_05 [sec:5.0+-0.0 nlpd-pi@T=2000:3.50+-0.09] bong-dlr10-LinHess [sec:6.1+-0.1 nlpd-pi@T=2000:3.32+-0.21] (b) 10 iterations. In Fig. 11 we show the test set performance for BLR (with LIN-HESS approximation) for 5 different learning rates (namely 5 × 10 -3 , 1 × 10 -2 , 5 × 10 -2 , 1 × 10 -1 , and 5 × 10 -1 ).When using 1 iteration per step, the best learning rate is α = 0.5, which is also the value chosen based on validation set performance. With this value, BLR matches BONG. For other learning rates, BLR performance is much worse. When using 10 iterations per step, there are several learning rates all of which give performance as good as BONG. . Data: sarcos. Expt: bakeoff-it1 bbb-dlr10-MCEF100-It1-LR0_05 [sec:6.9+-0.1 nlpd-pi@T=2000:5.29+-0.59] blr-dlr10-MCEF100-It1-LR0_5 [sec:13.4+-0.0 nlpd-pi@T=2000:4.45+-0.24] bog-dlr10-MCEF100-LR0_1 [sec:5.4+-0.1 nlpd-pi@T=2000:6.19+-0.90] bong-dlr10-MCEF100 [sec:12.3+-0.1 nlpd-pi@T=2000:4.12+-0.03] (a) 1 iteration. . Data: sarcos. Expt: bakeoff-it10 bbb-dlr10-MCEF100-It10-LR0_01 [sec:15.8+-0.1 nlpd-pi@T=2000:3.78+-0.10] blr-dlr10-MCEF100-It10-LR0_1 [sec:78.1+-0.5 nlpd-pi@T=2000:3.45+-0.05] bog-dlr10-MCEF100-LR0_1 [sec:5.4+-0.1 nlpd-pi@T=2000:6.19+-0.90] bong-dlr10-MCEF100 [sec:12.3+-0.1 nlpd-pi@T=2000:4.12+-0.03] (b) 10 iterations. . Data: sarcos. Expt: blr1-lr-bong bong-dlr10-LinHess [sec:6.1+-0.1 nlpd-pi@T=2000:3.32+-0.21] blr-dlr10-LinHess-It1-LR0_005 [sec:6.2+-0.1 nlpd-pi@T=2000:8.13+-0.22] blr-dlr10-LinHess-It1-LR0_01 [sec:6.2+-0.0 nlpd-pi@T=2000:6.12+-0.47] blr-dlr10-LinHess-It1-LR0_05 [sec:6.3+-0.0 nlpd-pi@T=2000:3.68+-0.70] blr-dlr10-LinHess-It1-LR0_1 [sec:6.4+-0.1 nlpd-pi@T=2000:3.51+-0.49] blr-dlr10-LinHess-It1-LR0_5 [sec:6.5+-0.0 nlpd-pi@T=2000:3.28+-0.21] (a) 1 iteration. . Data: sarcos. Expt: blr10-lr-bong bong-dlr10-LinHess [sec:6.1+-0.1 nlpd-pi@T=2000:3.32+-0.21] blr-dlr10-LinHess-It10-LR0_005 [sec:20.0+-0.1 nlpd-pi@T=2000:3.66+-0.68] blr-dlr10-LinHess-It10-LR0_01 [sec:20.5+-0.2 nlpd-pi@T=2000:3.50+-0.46] blr-dlr10-LinHess-It10-LR0_05 [sec:20.5+-0.1 nlpd-pi@T=2000:3.30+-0.25] blr-dlr10-LinHess-It10-LR0_1 [sec:21.2+-0.2 nlpd-pi@T=2000:3.29+-0.14] blr-dlr10-LinHess-It10-LR0_5 [sec:22.9+-0.3 nlpd-pi@T=2000:3.36+-0.17] (b) 10 iterations.Figure 11 : Same setup as Fig. 9 , except now we plot performance for BLR for 5 different learning rates. We also show BONG as a baseline, which uses a fixed learning rate step size of 1.0. . Data: sarcos. Expt: new-bbb1-lr bong-dlr10-LinHess [sec:6.1+-0.1 nlpd-pi@T=2000:3.32+-0.21] bbb-dlr10-LinHess-It1-LR0_005 [sec:5.9+-0.0 nlpd-pi@T=2000:8.13+-0.22] bbb-dlr10-LinHess-It1-LR0_01 [sec:5.9+-0.0 nlpd-pi@T=2000:5.90+-0.55] bbb-dlr10-LinHess-It1-LR0_05 [sec:5.9+-0.1 nlpd-pi@T=2000:3.50+-0.09] bbb-dlr10-LinHess-It1-LR0_1 [sec:16.7+-8.5 nlpd-pi@T=161:inf+-nan] bbb-dlr10-LinHess-It1-LR0_5 [sec:27.7+-0.3 nlpd-pi@T=54:1.09e+28+-inf] (a) 1 iteration. . Data: sarcos. Expt: new-bbb10-lr bong-dlr10-LinHess [sec:6.1+-0.1 nlpd-pi@T=2000:3.32+-0.21] bbb-dlr10-LinHess-It10-LR0_005 [sec:14.6+-0.1 nlpd-pi@T=2000:3.75+-0.52] bbb-dlr10-LinHess-It10-LR0_01 [sec:14.6+-0.1 nlpd-pi@T=2000:3.60+-0.25] bbb-dlr10-LinHess-It10-LR0_05 [sec:74.1+-33.4 nlpd-pi@T=516:5.10+-1.35] bbb-dlr10-LinHess-It10-LR0_1 [sec:132.5+-1.0 nlpd-pi@T=92:9.17+-2.51] (b) 10 iterations.Figure 12 : Same setup as Fig. 9 , except now we plot performance for BBB for 5 different learning rates. We also show BONG as a baseline.In Fig. 12 , we show the analogous plot for BBB. When using 1 iteration per step, all learning rates result in poor performance, with many resulting in NaNs. When using 10 iterations per step, there are some learning rates that enable BBB to get close to (but still not match) the performance of BONG. . Data: sarcos. Expt: new-bog-lr bong-dlr10-LinHess [sec:6.1+-0.1 nlpd-pi@T=2000:3.32+-0.21] bog-dlr10-LinHess-LR0_005 [sec:4.9+-0.1 nlpd-pi@T=2000:8.13+-0.22] bog-dlr10-LinHess-LR0_01 [sec:4.9+-0.0 nlpd-pi@T=2000:5.90+-0.55] bog-dlr10-LinHess-LR0_05 [sec:5.0+-0.0 nlpd-pi@T=2000:3.50+-0.09] bog-dlr10-LinHess-LR0_1 [sec:16.6+-8.7 nlpd-pi@T=161:inf+-nan] bog-dlr10-LinHess-LR0_5 [sec:26.8+-0.1 nlpd-pi@T=54:1.08e+28+-inf] (a) BOG with LIN-HESS. . Data: sarcos. Expt: new-bog-lr-ef bong-dlr10-MCEF100 [sec:12.3+-0.1 nlpd-pi@T=2000:4.12+-0.03] bog-dlr10-MCEF100-LR0_005 [sec:5.4+-0.1 nlpd-pi@T=2000:8.17+-0.24] bog-dlr10-MCEF100-LR0_01 [sec:5.3+-0.0 nlpd-pi@T=2000:8.18+-0.24] bog-dlr10-MCEF100-LR0_05 [sec:5.3+-0.1 nlpd-pi@T=2000:5.24+-0.45] bog-dlr10-MCEF100-LR0_1 [sec:5.4+-0.1 nlpd-pi@T=2000:6.19+-0.90] bog-dlr10-MCEF100-LR0_5 [sec:11.9+-9.3 nlpd-pi@T=270:3.39e+26+-inf] (b) BOG with MC-EF. Finally, in Fig. 13a , we show the analogous plot for BOG with LIN-HESS, and in Fig. 13b with MC-EF, where results are much worse.Overall we conclude that all the methods (except BONG) are quite sensitive to the learning rate. In our experiments, we pick a value based on performance on a validation set, but in the truly online setting, where there is just a single data stream, picking an optimal learning rate is difficult, which is an additional advantage of BONG.Proposition 4.1. To ease notation we write the natural parameters of the prior as ψ t|t-1 = [χ t|t-1 ; ν t|t-1 ], which can be interpreted as the prior sufficient statistics and prior sample size. Note that x t can be omitted as a constant. Based on the prior q ψ t|t-1 (θ t ) the exact posterior isEQUATION∝ q ψt (θ t ) (37)ψ t = χ t|t-1 + y t ν t|t-1 + 1 (38)For BONG, we first note the dual parameter is given byEQUATIONTherefore the natural gradient in Eq. ( 5) is∇ ρ t|t-1 E θt∼q ψ t|t-1 [log p (y t |θ t )] = ∇ ρ t|t-1 E θt∼q ψ t|t-1 θ t y t -A (θ t ) (41) = ∇ ρ t|t-1 ρ t|t-1 y t 1 (42) = y t 1 (43)Therefore the BONG update yields ψ t = [χ t|t-1 + y t ; ν t|t-1 + 1] in agreement with Eq. ( 38).Proposition 4.2. The intuition behind this proof is as follows. For the mean update in Eq. ( 9) ∇ θt t is linear in θ t so the expectation equals the value at the mean. For the covariance update in Eq. ( 10)∇ 2 θt t is independent of θ t so we can drop the expectation operator. The tricky part is why we need different linearizations for the Hessians to agree. It has to do with making the Hessian of the NN disappear (as in GGN). In the Gaussian approximation this happens when the predicted mean (ŷ t = f t (θ t )) is linear in θ t . In the plugin approximation it happens when the outcome-dependent part of the loglikelihood (h t (θ t ) y t ) is linear in θ t . In the latter case the only nonlinear term remaining in the log-likelihood is the log-partition A, and the two methods end up agreeing because of the property that the Hessian of the log-partition equals the conditional variance R t .Formally, under the linear(h)-Gaussian approximation in Eqs. ( 14) and ( 17) the expected gradient and Hessian can be calculated directly:EQUATION)For the linear(f )-delta approximation we use the properties of exponential families that (1) the gradient of the log-partition A with respect to the natural parameter f t (θ t ) equals the expectation parameter h t (θ t ), (2) the Hessian of the log-partition with respect to the natural parameter equals the conditional variance, and consequently (3) the Jacobian of the expectation parameter with respect to the natural parameter equals the conditional variance R t :∇ η=ft(µ t|t-1 ) A(η) = h t (µ t|t-1 ) = ŷ (48) ∇ 2 η=ft(µ t|t-1 ) A(η) = V y t |θ t = µ t|t-1 = R t (49) ∂h t (θ t ) ∂f t (θ t ) |θt=µ t|t-1 = R t (50)The last of these implies F t = R -1 t H t . Therefore the expected gradient and Hessian can be calculated asE θt∼δµ t|t-1 ∇ θt log pLD t (y t |θ t ) = ∇ θt=µ t|t-1 log pLD t (y t |θ t ) (51) = F t y t -F t ∇ η=ft(µ t|t-1 ) A(η) (52) = F t (y t -ŷt ) (53) = H t R -1 t (y t -ŷt ) (54) E θt∼δµ t|t-1 ∇ 2 θt log pLD t (y t |θ t ) = ∇ 2 θt=µ t|t-1 log pLD t (y t |θ t ) (55) = -F t ∇ 2 η=ft(µ t|t-1 ) A(η) F t (56) = -F t R t F t (57) = -H t R -1 t H t (58)In this section we give a more detailed derivation of BONG as mirror descent and use this to give two alternative interpretations of how BONG approximates exact VB: (1) by approximating the expected NLL as linear in the expectation parameter ρ, or (2) by replacing an implicit update with an explicit one.Assume the variational family is an exponential one as introduced at the end of Section 3, with natural and dual parameters ψ and ρ, sufficient statistics T (θ), and log-partition Φ(ψ):EQUATIONWe first review how NGD on an exponential family is a special case of mirror descent [Khan and Rue, 2023, Martens, 2020] . The mirror map is the gradient of the log-partition, which satisfies the thermodynamic identityρ = ∇Φ(ψ) (61)This is a bijection when Φ is convex (which includes the cases we study), so we can implicitly treat ψ and ρ as functions of each other. Given a loss function L(ψ), MD iteratively solves the local optimization problemEQUATIONThe first term is a linear (in ρ) approximation of L about the previous iteration ψ i and the second term is the Bregman divergenceEQUATIONThe Bregman divergence acts as a regularizer toward ψ i and captures the intrinsic geometry of the parameter space because of its equivalence with the (reverse) KL divergenceD KL q ψi+1 |q ψi = E θ∼q ψ i+1 [(ψ i+1 -ψ i ) T (θ) + Φ(ψ i ) -Φ(ψ i+1 )] (64) = D Φ (ψ i , ψ i+1 ) (65)Importantly, this recursive regularizer is not part of the loss and serves only to define an iterated algorithm that converges to a local minimum of L. Solving Eq. ( 62) by differentiating by ρ yields the MD updateEQUATION) Because the Fisher matrix for an exponential family is F ψ = ∂ρ/∂ψ, this is equivalent to NGD with respect to ψ. Khan and Rue [2023] offer this as a derivation of the BLR, when L(ψ) is taken to be the variational loss from Eq. ( 1).By applying this analysis to the online setting, our approach can be seen to follow from two insights. First, the online variational loss in Eq. ( 2) already includes KL divergence from the previous step, so we do not need the artificial regularizer in Eq. ( 62). That is, if we start from the online variational problem in Eq. ( 4) and define L t (ψ) as the expected NLL,EQUATIONthen replacing L t (ψ) with a linear approximation based at ψ t|t-1 and applying Eq. ( 65) leads toEQUATIONBy comparing to Eq. ( 62) we see this defines an MD algorithm with unit learning rate that works in a single step rather than by iterating. Paralleling the derivation of Eq. ( 66) from Eq. ( 62) we getEQUATIONwhich matches the BONG update in Eq. ( 5). Thus BONG can be seen as an approximate solution of the online variational problem in Eq. ( 4) based on linearizing the expected NLL wrt ρ. (Note this is different from the assumption underlying BONG-LIN that f t (θ t ) or h t (θ t ) is linear in θ t .)Second, in the conjugate case, this linearity assumption is true: L t is linear in ρ (see proof of Proposition 4.1). Therefore 68 is equivalent to solving Eq. ( 4) exactly:EQUATIONThis recapitulates Proposition 4.1 that BONG is Bayes optimal in the conjugate case. In general the exact solution to Eq. ( 70) isEQUATIONThis is an implicit update because the gradient is evaluated at the (unknown) posterior, whereas Eq. ( 69) is an explicit update because it evaluates the gradient at the prior. (In the Gaussian case these can be shown to match the implicit and explicit RVGA updates of [Lambert et al., 2021] .) Therefore BONG can be also interpreted as an approximation of exact VB that replaces the implicit update, Eq. ( 71), with an explicit update, Eq. ( 69).This section derives the update equations for all 80 algorithms we investigate (Table 3 plus the MC-HESS and LIN-EF variants). In Appendix E.6 we also translate the BLR algorithms from our online setting back to the batch setting used in Khan and Rue [2023] .For an exponential variational family with natural parameters ψ and dual parameters ρ, we can derive all 16 methods (BONG, BLR, BOG, BBB under all four Hessian approximations) from four quantities:EQUATIONThe NGD methods (BONG and BLR) use gradients with respect to ρ t,t-i while the GD methods (BOG and BBB) use gradients with respect to ψ t,i-1 . For BONG and BOG the D KL term is not relevant, and there is no inner loop so ψ t,i-1 = ψ t|t-1 and gt,i = g t , G t,i = G t .When ψ is not the natural parameter of an exponential family we must explicitly compute the inverse-Fisher preconditioner for the NGD methods. Therefore the updates can be derived from these three quantities:EQUATIONWe will frequently use Bonnet's and Price's theorems [Bonnet, 1964 , Price, 1958 ]EQUATIONEQUATIONFor diagonal Gaussians with covariance Diag σ 2 , Price's theorem also implies6 EQUATIONUpdate equations for MC-HESS, MC-EF and LIN-EF methods are displayed together in the subsections that follow, because for the most part they differ only in the choice of G t = -1 M Ĝ(1:M) t Ĝ(1:M) t or G LIN-EF t = -g LIN t (g LIN t )allows a more efficient update. We derive updates for BONG-LIN-HESS and BOG-LIN-HESS from the corresponding BONG-MC-HESS and BOG-MC-HESS updates using Proposition 4.2 which entails substitutingEQUATION)For the algorithms with inner loops (BLR and BBB) we adapt the notation of Section 4.4 as follows:EQUATIONThis corresponds to basing the linear(f )-Gaussian and linear(h)-delta approximations at µ t,i-1 instead of µ t|t-1 . Thus the updates for BLR-LIN-HESS and BBB-LIN-HESS are obtained by substitutingEQUATIONE.1 Full covariance Gaussian, natural parametersThe natural and dual parameters for a general Gaussian are given by ψ (1)t,i-1 = Σ -1 t,i-1 µ t,i-1 ρ (1) t,i-1 = µ t,i-1 (93) ψ (2) t,i-1 = -1 2 Σ -1 t,i-1 ρ (2) t,i-1 = µ t,i-1 µ t,i-1 + Σ t,i-1 (94) Inverting these relationships gives µ t,i-1 = -1 2 ψ (2)-1 t,i-1 ψ (1) t,i-1 = ρ (1) t,i-1 (95) Σ t,i-1 = -1 2 ψ (2)-1 t,i-1 = ρ (2) t,i-1 -ρ (1) t,i-1 ρ (1) t,i-1 (96)The KL divergence in the VI loss isEQUATIONFollowing Appendix C of [Khan et al., 2018a] , for any scalar function the chain rule givesEQUATIONand ∇ ρ (1)EQUATIONFollowing the same approach for ψ givesEQUATIONEQUATIONand∇ ψ (1) t,i-1 D KL q ψt,i-1 |q ψ t|t-1 = Σ t,i-1 Σ -1 t|t-1 µ t,i-1 -µ t|t-1 (116) ∇ ψ (2) t,i-1 D KL q ψt,i-1 |q ψ t|t-1 = 2Σ t,i-1 Σ -1 t|t-1 µ t,i-1 -µ t|t-1 µ t,i-1 + Σ t,i-1 Σ -1 t|t-1 -Σ -1 t,i-1 Σ t,i-1 (117) E.1.1 BONG FC (explicit RVGA)Substituting Eqs. ( 104) and ( 105) into Eq. ( 5) givesEQUATIONTranslating to (µ t , Σ t ) gives the BONG-FC updateEQUATIONEQUATIONThis is equivalent to the explicit update form of RVGA [Lambert et al., 2021] . Using G MC-HESS t this update takes O(P 3 ) because of the matrix inversion. Using G MC-EF t and the Woodbury matrix identity we can write the update in a form that takes O(M P 2 + M 3 ):EQUATIONEQUATIONLikewise using G LIN-EF t takes O(P 2 ):EQUATIONEQUATIONK t = Σ t|t-1 g LIN t 1 + (g LIN t ) Σ t|t-1 g LIN t (127) E.1.2 BONG-LIN-HESS FC (CM-EKF)Applying Proposition 4.2 to Eqs. ( 120) and ( 121) gives the BONG-LIN-HESS-FC updateEQUATIONEQUATIONThis is equivalent to CM-EKF [Tronarp et al., 2018 , Ollivier, 2018] and can be rewritten using the Woodbury identity in a form that takes O(CP 2 + C 3 ):EQUATIONEQUATIONE.1.3 BLR FC Substituting Eqs. ( 104) to (107) into Eq. ( 33) givesEQUATIONTranslating to (µ t,i , Σ t,i ) gives the BLR-FC updateEQUATIONEQUATIONThis update takes O(P 3 ) per iteration because of the matrix inversion. In Appendix E.1.1 we were able to use the Woodbury identity to exploit the low rank of G MC-EF t and G LIN-EF t and obtain BONG updates with complexity quadratic in P . This does not appear possible with Eq. ( 136) because of the extra precision term on the RHS (applying Woodbury would require inverting (1 -α)Σ -1 t,i-1 + αΣ -1 t|t-1 ). Therefore unlike BONG-FC, BLR-FC requires time cubic in the model size, for reasons that can be traced back to the KL term in Eq. ( 33).Applying Proposition 4.2 to Eqs. ( 135) and ( 136) gives the BLR-LIN-HESS-FC updateEQUATIONEQUATIONThis update takes O(P 3 ) per iteration because of the matrix inversion.Substituting Eqs. ( 114) and (115) into Eq. (31) givesEQUATIONTranslating to (µ t , Σ t ) gives the BOG-FC updateEQUATIONThis update takes O(P 3 ) because of the matrix inversion. The greater cost of the BOG-FC update relative to BONG-FC can be traced to the difference between GD and NGD: the NLL gradients wrt ψ t|t-1 in Eqs. ( 114) and ( 115) are more complicated than the gradients wrt ρ t|t-1 in Eqs. ( 104) and ( 105).Applying Proposition 4.2 to Eqs. ( 141) and ( 142) gives the BOG-LIN-HESS-FC updateEQUATIONEQUATIONThis update takes O(P 3 ) because of the matrix inversion.Substituting Eqs. ( 114) to (117) into Eq. (34) givesψ (1) t,i = ψ (1) t,i-1 + αΣ t,i-1 g t,i -αΣ t,i-1 Σ -1 t|t-1 µ t,i-1 -µ t|t-1 (145) ψ (2) t = ψ (2) t,i-1 + 2αΣ t,i-1 g t,i µ t,i-1 + αΣ t|t-1 G t Σ t|t-1 -2αΣ t,i-1 Σ -1 t|t-1 µ t,i-1 -µ t|t-1 µ t,i-1 -αΣ t,i-1 Σ -1 t|t-1 -Σ -1 t,i-1 Σ t,i-1 (146)Translating to (µ t,i , Σ t,i ) gives the BBB-FC updateEQUATIONThis update takes O(P 3 ) per iteration because of the matrix inversion.Applying Proposition 4.2 to Eqs. ( 147) and ( 148) gives the BBB-LIN-HESS-FC updateEQUATIONEQUATIONThis update takes O(P 3 ) per iteration because of the matrix inversion.The Bonnet and Price theorems giveEQUATIONEQUATIONFrom Appendix E.1 we haveEQUATIONWe write the Fisher with respect to the moment parameters ψ = (µ, vec(Σ)) as a block matrix:EQUATIONThe blocks can be calculated by the second-order Fisher formulaEQUATIONIn the final line we usedEQUATIONwith ij and k treated as composite indices in [P 2 ]. Therefore the preconditioner in the NGD methods isEQUATIONSubstituting Eqs. ( 151), ( 152) and ( 168) into Eq. ( 3) gives the BONG-FC_MOM updateEQUATIONThis update takes O(P 3 ) using G MC-HESS t, O(M P 2 ) using G MC-EF t , and O(P 2 ) using G LIN-EF t . The update is similar to the BONG-FC update except that Eq. ( 170) ignores the Ĝ(1:M)t Σ t|t-1 Ĝ(1:M) tterm in Eq. ( 124) or the (g LIN t ) Σ t|t-1 g LIN t term in Eq. ( 127) which estimate the epistemic part of predictive uncertainty.Applying Proposition 4.2 to Eqs. ( 169) and ( 170) gives the BONG-LIN-HESS-FC_MOM updateEQUATIONEQUATIONThis update takes O(CP 2 ).Substituting Eqs. ( 151) to ( 154) and (168) into Eq. ( 32) gives the BLR-FC_MOM updateEQUATIONEQUATIONThis update takes O(P 3 ) per iteration because of the matrix inversion. UsingG MC-EF t = -1 M Ĝ(1:M) t Ĝ(1:M) t or G LIN-EF = -g LIN t (g LIN t )allows the BONG-FC_MOM covariance update in Eq. ( 170) to scale quadratically, but this does not help here. Instead the BLR-FC_MOM update scales cubically because of the additional Σ -1 t|t-1 term that comes from the KL divergence in the VI objective.Applying Proposition 4.2 to Eqs. ( 173) and ( 174) gives the BLR-LIN-HESS-FC_MOM updateEQUATIONThis update takes O(P 3 ) per iteration because of the matrix inversion in the Σ -1 t|t-1 term that comes from the KL divergence in the VI objective.Substituting Eqs. ( 151) and ( 152) into Eq. ( 31) gives the BOG-FC_MOM updateEQUATION)Σ t = Σ t|t-1 + α 2 G t (178)Note the mean update is vanilla online gradient descent (OGD) and does not depend on the covariance. This update takes O(M P 2 ) using G MC-HESS t or G MC-EF t and O(P 2 ) using G LIN-EF t .Applying Proposition 4.2 to Eqs. ( 177) and ( 178) gives the BOG-LIN-HESS-FC_MOM updateEQUATIONEQUATIONThis update takes O(CP 2 ).Substituting Eqs. ( 151) to (154) into Eq. ( 34) gives the BBB-FC_MOMEQUATIONEQUATIONThis update takes O(P 3 ) per iteration because of the matrix inversion, which traces back to the VI objective. Comparing to the BOG-FC_MOM update in Eqs. ( 177) and ( 178) (which has quadratic complexity in P ), the extra terms here come from the KL part of Eq. ( 34).Applying Proposition 4.2 to Eqs. ( 181) and ( 182) gives the BBB-LIN-HESS-FC_MOM updateEQUATIONEQUATIONThis update takes O(P 3 ) per iteration because of the matrix inversion. Comparing to the BOG-LIN-HESS-FC_MOM update in Eqs. ( 179) and ( 180) (which has quadratic complexity in P ), the extra terms here come from the KL part of Eq. ( 34).Throughout this subsection, vector multiplication and exponents are elementwise.The natural and dual parameters for a diagonal Gaussian are given by ψ (1)EQUATIONInverting these relationships givesµ t,i-1 = -1 2 ψ (2) t,i-1 -1 ψ (1) t,i-1 = ρ (1) t,i-1 (187) σ 2 t,i-1 = -1 2 ψ (2) t,i-1 -1 = ρ (2) t,i-1 -ρ (1) t,i-1 2 (188)The KL divergence in the VI loss isEQUATIONFor any scalar function the chain rule givesEQUATIONThereforeEQUATIONand∇ ρ (1)t,i-1EQUATIONFollowing the same approach for ψ givesEQUATION∇ ψ (2)EQUATIONandEQUATIONOur implementations often make use of the following trick: Suppose A ∈ R n×m and B ∈ R m×n .Then we can efficiently compute diag(AB) in O(mn) time using (AB) ii = M j=1 A ij B ji .For MC-HESS methods, we approximate the diagonal of the Hessian for each MC sample θ(m) t using Hutchinson's trace estimation method [Hutchinson, 1989] which has been used in other DNN optimization papers such as adahessian Yao et al. [2021] . This involves an extra inner loop with size denoted N .Substituting Eqs. ( 196) and (197) into Eq. ( 5) givesEQUATIONTranslating to µ t , σ 2 t gives the BONG-DIAG updateEQUATIONThis update takes O(M P ) to estimate G t using G MC-EF t , O(N M P ) using G MC-HESS t and Hutchinson's method, and O(P ) using G LIN-EF t .Applying Proposition 4.2 to Eqs. ( 210) and ( 211) gives the BONG-LIN-HESS-DIAG updateEQUATION)This update is equivalent to VD-EKF [Chang et al., 2022] and takes O(C 2 P ).Substituting Eqs. ( 196) to (199) into Eq. ( 33) givesEQUATIONTranslating to µ t,i , σ 2 t,i gives the BLR-DIAG updateEQUATIONEQUATIONThis update takes O(M P ) per iteration to estimate G t using G MC-EF t , O(N M P ) per iteration using G MC-HESS t and Hutchinson's method, and O(P ) per iteration using G LIN-EF t .The MC-HESS and MC-EF versions of this update are respectively equivalent to VON and VOGN [Khan et al., 2018b] in the batch setting where we replace q ψ t|t-1 with a spherical prior N (0, λ -1 I P ) (see Appendix E.6).with W ∈ R P ×R and diagonal Υ t|t-1 ∈ R P ×P . We sometimes abuse notation by writing Υ t|t-1 for the vector diag Υ t|t-1 when the meaning is clear from context.Substituting the DLR form in the gradients for the KL divergence derived in Appendix E.2 givesEQUATIONFor any function the chain rule givesEQUATIONTherefore the gradients we need areEQUATIONEQUATIONandEQUATIONThe Fisher matrix can be decomposed as a block-diagonal with blocks for µ t,i and for (Υ t,i , W t,i ), but (in contrast to the FC Gaussian case in Eq. ( 168)) we have not found an efficient way to analytically invert the latter block, which has size P + RP . To avoid the O R 3 P 3 cost of brute force inversion, we use a different strategy for BONG and BLR of performing the update on the natural parameters of the FC Gaussian as in Appendix E.1 and projecting the result back to rank R using SVD. Specifically, assume the updated precision from applying Eq. ( 121) for BONG or Eq. ( 136) for BLR can be written asEQUATIONand let the SVD of Wt,i be Wt,i = U t,i Λ t,i V t,i (275) with U t,i , V t,i orthogonal and Λ t,i rectangular-diagonal. Following [Mishkin et al., 2018 , Chang et al., 2023] we define the update W t,i = U t,i [:, : R] Λ t,i [: R, : R] (276)EQUATIONso that W t,i contains the top R singular vectors and values of the FC posterior and the diagonal is preserved: diag Υ t,i + W t,i W t,i = diag Σ-1 t,i . This approach works for MC-EF and LIN-EF methods but not MC-HESS, which we omit.Finally, in the MC-EF methods we sample from the DLR prior using the routine in [Mishkin et al., 2018 , Lambert et al., 2023] which takes O(R(R + M )P ).Substituting the DLR prior from Eq. ( 264) into the FC precision update from Eq. ( 121) and using the MC-EF approximation yields the posterior precisionEQUATIONNote Wt ∈ R P ×(R+M ) . Using this posterior precision in the mean update from Eq. ( 120) yieldsEQUATIONwhere the second line comes from the Woodbury identity and can be computed in O((R + M ) 2 P + (R + M ) 3 ). Applying the SVD projection gives W t = U t [:, : R] Λ t [: R, : R] (284)EQUATIONwhich takes O((R + M ) 2 P ) for the SVD. Therefore the BONG-MC-EF-DLR update is defined by Eqs. ( 281) and ( 283) to ( 286) and takes O((R + M ) 2 P + (R + M ) 3 ).The BONG-LIN-EF-DLR update comes from replacing 1 √ M Ĝ(1:M) t with g LIN t in Eq. ( 281) and replacing I R+M with I R+1 in Eq. ( 283). This update takes O((R + 1) 2 P + (R + 1) 3 ).Applying Proposition 4.2 to Eqs. ( 281) and ( 283) gives the BONG-LIN-HESS-DLR update:µ t = µ t|t-1 + Υ -1 t|t-1 -Υ -1 t|t-1 Wt I R+C + W t Υ -1 t|t-1 Wt -1 W t Υ -1 t|t-1 × H t R -1 t (y t -ŷt ) (287) W t = U t [:, : R] Λ t [: R, : R] (288) Υ t = Υ t|t-1 + diag Wt W t -W t W t (289) (U t , Λ t , _) = SVD Wt (290) Wt = W t|t-1 , H t A t (291) A t = chol R -1 t (292)This is equivalent to LO-FI [Chang et al., 2023] and takes O((R + C) 2 P + (R + C) 3 ).Substituting DLR forms for q ψ t|t-1 and q ψt,i-1 into the FC precision update from Eq. ( 136) and using the MC-EF approximation yields the posterior precision Σ-1EQUATIONEQUATIONNote Wt ∈ R P ×(2R+M ) . Using this posterior precision in the mean update from Eq. (135) yieldsEQUATIONEQUATIONwhere the second line comes from the Woodbury identity and can be computed in O((2R + M ) 2 P + (2R + M ) 3 ). Applying the SVD projection gives W t,i = U t,i [:, : R] Λ t,i [: R, : R] (299)EQUATIONEQUATIONwhich takes O((2R + M ) 2 P ) for the SVD. Therefore the BLR-MC-EF-DLR update is defined by Eqs. ( 295), ( 296) and ( 298) to (301) and takes O((2R + M ) 2 P + (2R + M ) 3 ) per iteration. Notice Wt has larger rank for BLR than for BONG (2R + M vs. R + M ) because of the extra √ αW t|t-1 term that originates in the KL part of the VI loss. This difference will slow BLR-MC-EF-DLR relative to BONG-MC-EF-DLR especially when R is not small relative to M . This method closely resembles SLANG Mishkin et al. [2018] in the batch setting where we replace q ψ t|t-1 with a spherical prior N (0, λ -1 I P ) (see Appendix E.6).We can also define a BLR-LIN-EF-DLR method by replacingα M Ĝ(1:M) t with √ αg LIN tin Eq. ( 296) and I 2R+M with I 2R+1 in Eq. ( 298). This update takes O((2R + 1) 2 P + (2R + 1) 3 ) per iteration.Applying Proposition 4.2 to Eqs. ( 296) and ( 298) gives the BLR-LIN-HESS-DLR update:EQUATIONEQUATIONEQUATIONEQUATIONThis update takes O((2R + C) 2 P + (2R + C) 3 ) per iteration. As with the EF versions of BLR-DLR, Wt has larger rank for BLR-LIN-HESS-DLR than for BONG-LIN-HESS-DLR (2R + C vs. R + C). This difference will slow BLR-LIN-HESS-DLR relative to BONG-LIN-HESS-DLR especially when R is not small relative to C.Substituting Eqs. ( 269) to (271) into Eq. ( 31) gives the BOG-DLR updateEQUATIONEQUATIONUsing the EF approximation and Woodbury, the BOG-MC-EF-DLR update can be rewritten asEQUATIONEQUATIONB t = Υ -1 t|t-1 -Υ -1 t|t-1 W t|t-1 I R + W t|t-1 Υ -1 t|t-1 W t|t-1 -1 W t|t-1 Υ -1 t|t-1 Ĝ(1:M) t (314) which takes O(RM P + M R 2 + R 3 ).The BOG-LIN-EF-DLR update comes from replacing Ĝ(1:M) t with g LIN t in Eq. ( 314) and dropping the M -1 factors in Eqs. ( 312) and ( 313). This update takes O(RP + R 3 ).Applying Proposition 4.2 to Eqs. ( 309) and ( 312) to (314) gives the BOG-LIN-HESS-DLR updateEQUATIONEQUATIONEQUATIONEQUATIONThis update takes O(C(C + R)P + CR 2 + R 3 ).Substituting Eqs. ( 265) and ( 269) to (273) into Eq. ( 34) gives the BBB-DLR updateµ t,i = µ t,i-1 + α Υ t|t-1 + W t|t-1 W t|t-1 µ t|t-1 -µ t,i-1 + αg t (319) Υ t,i = Υ t,i-1 + α 2 diag Σ t,i-1 Υ t|t-1 -Υ t,i-1 + W t|t-1 W t|t-1 -W t,i-1 W t,i-1 -G t,i Σ t,i-1 (320) W t,i = W t,i-1 + αΣ t,i-1 Υ t|t-1 -Υ t,i-1 + W t|t-1 W t|t-1 -W t,i-1 W t,i-1 -G t,i Σ t,i-1 W t,i-1 (321)The previous covariance can be written using Woodbury asEQUATIONThe BBB-MC-EF-DLR update can be computed efficiently by expanding terms in Eqs. (320) to (322). For example the terms involving G t,i can be calculated asEQUATIONUsing this strategy the update takes O((R + M )RP + M R 2 + R 3 ).The BBB-LIN-EF-DLR update comes from replacing Ĝ(1:M) t with g LIN t and dropping the M -1 factors in Eqs. ( 323) to (325). This update takes O(R 2 P + R 3 ).Applying Proposition 4.2 to Eqs. ( 319) to (321) gives the BBB-LIN-HESS-DLR updateEQUATIONThis can be computed in O((C + R) 2 P + CR 2 + R 3 ) using Eq. ( 322) and following a computational approach similar to Eqs. ( 323) to ( 325).It is interesting to translate the BLR updates derived here back to the batch setting where BLR was developed [Khan and Rue, 2023] , by replacing N µ t|t-1 , Σ t|t-1 with a centered spherical prior N 0, λ -1 I P .The batch BLR-FC update becomesEQUATIONEQUATIONThe batch BLR-LIN-HESS-FC update becomesEQUATIONThe batch BLR-FC_MOM update becomesEQUATIONThe batch BLR-LIN-HESS-DLR update becomesEQUATIONEQUATIONThis algorithm could be called SLANG-LIN-HESS and would be deterministic and faster than SLANG since it does not need MC sampling. We can also define SLANG-LIN-EF which would be even faster, by replacingα M Ĝ(1:M) i with √ αg LIN iin Eq. ( 349) and I R+M with I R in Eq. ( 345).• 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 or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.We assume the training and test sets are drawn from the same static distribution. Alternatively, if there is only one stream of data coming from a potential notstationary source, we can use the prequential or one-step-ahead log predictive densityGama et al. [2013]. We leave studying the non-stationary case to future work.That is, we compute S = 100 posterior samples θ s t ∼ p(θt|D train 1:t ) and then use p(y|x) ≈ 1 S S s=1 p(y|x, θ s t ). For efficiently sampling from a DLR Gaussian we followMishkin et al. [2018]. FollowingImmer et al. [2021b], we find it better to linearize the model (i.e., replacing h(xi, θt) with h(xi, θt) defined in Eq. (14)) before pushing posterior samples through. We call this the Lin-MC approximation.This is likely due to the fact that in the moment update in Eq. (249), we add diag(Gt) to the variance, whereas in the natural update in Eq. (213), we subtract diag(Gt) from the precision; the latter is more stable since Gt is negative semi-definite for allHessian approximations.Lin-MC is defined in Footnote 2. The motivation for this approximation (fromImmer et al. [2021b]) is the following: If we push posterior samples through a nonlinear predictive model, the results can be poor if the samples are far from the mean, but if we linearize the predictive model, extrapolations away from the mean are more sensible. This is true even if the posterior was not explicitly computed using a linear approximation.This value was based on a small amount of trial and error. Using a smaller value of σ0 results in underfitting relative to a linear least squares baseline, and using a much larger value results in unstable posterior covariances, causing the NLPD-MC samples to result in NaNs after a few hundred steps.We use diag(A) to denote the vector of diagonal elements of matrix A and Diag(v) to denote the matrix whose diagonal entries are v and off-diagonal entries are 0. 


Decision: This paper studies the problem of online Bayesian inference in Bayesian neural networks using variational inference. The paper considers an online natural gradient update rule given by only considering the expected log-likelihood term and using a unit learning rate. This Bayesian online natural gradient update rule is studied theoretically under a conjugate prior setting. Next, the paper proposes a deterministic way of efficiently computing the expectations in the update rule for Gaussian variational families. The method is then evaluated on a number of data sets, including the MNIST data set. 

Overall, the reviewers found the contribution to be interesting to the NeurIPS community. Reviewer 6NTx commented on a lack of experiments evaluating the resulting model uncertainty; these concerns were addressed during the rebuttal and discussion phase. Several reviewers brought up issues of novelty and wanted to see additional discussion that much more explicitly discusses the contributions of this paper in the context of related update rules in the literature. Please revise the paper to include these discussions in the final version.

Based on the reviewer comments, I recommend accepting the paper.

Reviews: review: paper_summary: paper_summary: This paper proposes a generalized framework for variational sequential inference of Bayesian neural networks with Gaussian prior distributions by further approximating the KL-divergence and expectations. It presents a very elegant unification of a wide range of (approximate) variational inference algorithms for Bayesian neural networks developed over the past 10 years and contains an extensive empirical evaluation on subsets of the MNIST dataset.
strength_weakness: strength_weakness: The space of 9 pages is way to short to do this full justice and I often had to read into the 30 pages (!) appendix. I recommend the authors to further considering a journal submission.
questions: questions: Page 4: The abbreviation PSD as "positive-semi-definite" has never been introduced.
Page 5: jac has not been introduced
Page 7: In line 239, the component-wise product of \sigma^{-2} and \mu should be indicated somewhere
limitations: limitations: This is not applicable.
paper_summary: paper_summary: This paper proposes a new approximate Bayesian technique specifically for online learning whereby posterior distributions over modeling parameters at time $t$ can be achieved with a single natural gradient step evaluated on a model using the previous posterior distribution at time $t-1$ as the prior distribution. This method is rigorously derived in a conjugate setting with an exponential family and extended approximately with various ablations to more expressive black-box settings. Extensive empirical studies are conducted validating the performance of the method.
strength_weakness: strength_weakness: The main weakness of the method is the lack of experiments evaluating the resulting model uncertainty achieved after training. From what I read, it seems to only be predictive performance via misclassification and negative log predictive density. Granted, the latter does have some ablations where it is integrated over posterior samples; however, the presentation of these results are a bit difficult to interpret other than just methods with lower values being better. It would be interesting to investigate the predictive uncertainty via expected calibration error, or some other calibration metric.
questions: questions: No other questions, please address the main weaknesses. Please let me know if I missed or misunderstood anything; I am more than happy to be convinced otherwise.
limitations: limitations: The authors adequately addressed the limitations of the proposed method.
paper_summary: paper_summary: The paper proposed a novel parameter update rule. The rule is natural gradient descent on the variational inference loss, where the prior term is dropped. Such rule is motivated in theory under the assumption of likelihood and variational distribution being conjugate, and in practice though experiments on MNIST.
Several variants are considered: various loss, various hessian approximations and various way of dealing with the intractable expectation.
strength_weakness: strength_weakness: The theoretical motivation is quite weak:\
-Proposition 4.1 is nice, but the assumptions are overly restrictive. It feels like this is such a special case that is has no implication on real case scenarios \
-The claimed contribution (in Line 44) is a very weak contribution. The closed from expression with linearization+gaussian likelihood was already observed by Immer, and I don't think it is fair to claim it as a novel contribution of the paper. Specifically because there is no motivation to such approximations.

On the other hand, as the authors also correctly point out in Line 307, the empirical evaluation is based on MNIST only. A very small dataset which doesn't really motivate the superiority of the proposed algorithm.

Overall, I don't think the authors provided enough evidence (either theoretical or empirical) of the good performance of the proposed update rule.

TYPOS: \
-Line 126. I think the $\theta_{-1}$ was supposed to be a $\theta{t-1}$ \
-There is a duplicate reference: Line 336 and Line 340 are the same
questions: questions: The notation $f_t$ that appears in Equations (9) and (10) was never introduced. I guess it refers to the evaluation of $f$ on $x_t$, is it correct?
limitations: limitations: Limitations are well discussed. Good job on this point
paper_summary: paper_summary: The authors explore sequential Bayesian inference using variational inference (VI). They propose to remove the regularizing KL term, performing a single natural gradient step using the expected log-likelihood only, and provide relevant ways to approximate the hard-to-compute expectations. They support their approach with detailed experiments.
strength_weakness: strength_weakness: My main concern is about the claimed novelty of the approach ("removing the KL" and "implicit regularization towards the prior"). The idea of using a one-step natural gradient descent on the expected log-loss directly (without the KL term) has already been proposed in the literature and analyzed; for example, [Lyu & Tsang, 2021] studied this algorithm in a model-free setting where a general function f is considered instead of the log-loss. This one-step natural gradient update is in fact simply the solution of Equation (2) of the present paper, but with a linearized expected log-loss. This dual interpretation resonates with the dual interpretation of standard Euclidean gradient descent (single update or minimization of a first-order development of the loss to be minimized). 

Y. Lyu and I. W. Tsang: Black-box optimizer with stochastic implicit natural gradient. In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 217–232, 2021.
questions: questions: NA
limitations: limitations: NA
paper_summary: paper_summary: The paper introduces a framework for variational online learning, in which a method is defined by three ingredients: (1) implicit or explicit regularisation (2) type of approximation of the expected gradient and hessian (3) geometry in which the gradient is taken (ie employing natural gradients or standard gradients). In implicit regularisation, one learning step is taken for each datapoint, with no regularisation; explicit method is the more standard approach of optimising with multiple steps and regularisation to the posterior at the previous step. A case study on different combinations of the ingredients is conducted; some of the methods appeared in previous papers. The study finds the best performing method is the one with implicit regularisation, linearisation of the hessian, and natural gradient.
strength_weakness: strength_weakness: Given that the main strength of the paper is the unifying recipe, and the comparison of using different ingredients, a more exhaustive case study would strengthen the work (and the less essential comments, like the equivalence to the full Bayesian update in the conjugate case, in my opinion can be placed in the appendix and referenced). Nevertheless, I find the contribution of the paper warrants it being accepted.
questions: questions: * Given the well documented efficiency of natural gradient steps, is there a danger of the method with no regularisation being more prone to "forgetting"?
limitations: limitations: N/A

Comments: comment: We fully agree with your suggestions about presentation. We will put more emphasis on the idea of implicit regularization while also acknowledging that the resulting update rule matches what others have proposed based on different motivations. Thanks again for your thoughtful feedback!
I think that your statement "our primary proposal of dropping the KL term and doing one NGD step still appears to be new" in the first paragraph of your last reply is only partly true: the motivation of dropping the KL is new, but the one NGD step is not. It's exactly what I meant in my first review by "the idea of using a one-step natural gradient descent on the expected log-loss directly (without the KL term) has already been proposed in the literature and analyzed". You're right saying that linearizing the NLL in the full ELBO is an alternative means of deriving your update, but to be fair it exactly results in the same update. This is either your Equation (5) or Equation (4) of Lyu & Tsang (2021). In a sense, the KL is dropped in both cases, even though the initial motivation may be different.

I'm not sure I fully understand why Khan told you there is no close connection to his works, because unless I'm mistaken, the SVB Equation (7) in Cherief-Abdellatif, Alquier & Khan (2019) is exactly the update rule (2) of Lyu & Tsang (i.e. linearization of the NLL in the full ELBO), which is then reformulated as a single-step update in Equation (10) of Cherief-Abdellatif, Alquier & Khan. This is not the same rule as your Equation (5) or Equation (4) of Lyu & Tsang because of a different choice of the variational parameter, but this still provides a one-step GD-like update based on the NLL only? 

Once again, I want to insist on the fact that I enjoyed the paper. I am convinced by the comments of the authors and will raise my score to 6 accordingly. However, I strongly advise the authors to work on the presentation of their paper, and recommend to put emphasis on the idea rather than on the rule (5) that has already been proposed and studied in the literature (in the exact same form by Lyu & Yang, and under a slightly different variant by Cherief-Abdellatif, Alquier & Khan?). A particular point worth highlighting is that their approach is actually framed in its most general form in Equation (3), which is broader than the specific one-step NGD rule (5) that has already been derived from a different perspective for exponential families.
Thank you for the additional comments.

To clarify our previous reply, linearizing the loss$\^*$ in the full ELBO is an alternative means of deriving our update which we explain in our new appendix (to be added in the revision). The papers you reference are all good examples of this approach and we will cite them there. Our primary proposal of dropping the KL term and doing one NGD step still appears to be new.

The algorithms in Chérief-Abdellatif, Alquier & Khan (2019) all use linearization of the loss, and we don't see any suggestions there to perform a single gradient descent step for the expected loss without the KL. Their SVB algorithm (eq 7) removes the KL divergence from the time 0 prior but replaces it with KL from the previous time step, which is the same KL we start with (before dropping it). Please let us know if we're missing something.

We have discussed our work with Khan since our submission, and the closest his papers have come to the BONG update is a passing comment in section 5.1 of Khan and Rue (2023) that conjugate updating is equivalent to running the BLR for one step with learning rate 1. This is close to our proposition 4.1 except that BLR includes the KL term that BONG drops. The reason BLR and BONG agree in this case is that the gradient of the KL is zero on the first iteration: $\\nabla\_{\\boldsymbol{\\psi}=\\boldsymbol{\\psi}\_{t-1}} KL(q\_{\\boldsymbol{\\psi}} | q_{\\boldsymbol{\\psi}\_{t-1}}) = \\boldsymbol{0}$. When BLR is run for multiple iterations (as it normally is) the gradient of the KL contributes to the update. In our revision we will note that BONG can be obtained as a special case of BLR with $\\alpha=I=1$. Nevertheless we believe that omitting the KL entirely because its regularizing role can be replaced by truncated (one step) updating is a substantially new insight relative to noting that mathematically its gradient drops out on the first iteration.

We appreciate your careful thought about our work. Hopefully this discussion increases your estimation of the novelty of our primary proposal, in addition to the other contributions of our paper.

—

$\^*$Just to avoid any possible confusion: Linearizing the loss should not be confused with linearizing the model predictions $f\_t(\\theta)$ or $h\_t(\\theta)$ as in our proposition 4.2. The latter is an additional trick for approximating the expected Hessian.
Thank you for the additional experiments conducted. These do indeed resolve the concerns I had regarding the method. 

I will maintain my initial score of a 7 with the hope of this paper being accepted.
The other is the linear($f$)-delta approximation.

The linear($h$)-Gaussian approximation linearizes the mean parameters of the observation distribution (e.g., the predicted class probabilities) and replaces the likelihood with a moment-matched Gaussian. This is the method of Ollivier and Tronarp et al. The linear($f$)-delta approximation linearizes the natural parameters of the observation distribution (e.g., the predicted class logits) and replaces the expectations in (9) and (10) with plugin values at the prior mean. This is new.

Proposition 4.2 shows the linear($f$)-delta and linear($h$)-Gaussian approximations yield the same expressions for the expected gradient $\\boldsymbol{g}_t$ and expected Hessian $\\boldsymbol{G}_t$ and consequently the same updates.

We see that our initial rebuttal wrote linear($f$) in two places where we meant linear($h$). Apologies for the confusion. There is a similar error at lines 51-52 in the paper which we will correct.
Thank you very much for your reply and the detailed explanations. However, I remain unconvinced about the novelty of "connecting all these ideas to the ELBO", as well as the claim that the BONG update's "derivation and application to variational inference are novel".

While it is true that Lyu & Tsang's paper "goal is optimization, not inference or Bayesian updating", it was merely an example of a paper that provides a general analysis of this specific algorithm. The idea has been explored in various papers with different approaches and motivations, such as [van der Hoeven et al., 2018; Cherief-Abdellatif et al., 2019; Khan & Rue, 2023]. In particular, the authors of [Khan & Rue, 2023] discuss the online VI case at the end of Section 5 and mention the linearization of the expected loss in the ELBO, referring to [van der Hoeven et al., 2018], which delves into these questions in detail (though not through the lens of VI). From a different perspective, [Cherief-Abdellatif et al., 2019] establish a connection between VI and online learning by proposing to perform a single gradient descent step for the expected loss without the KL, which they show to be equivalent to minimizing the ELBO with a linearized expected loss.

Therefore, while I do believe this is a nice piece of work, I am still not fully convinced of its novelty. And although I recognize that the paper offers more than just this point, I feel that it is the aspect most emphasized (particularly in the abstract).


[van der Hoeven et al., 2018]: D. van der Hoeven, T. van Erven, and W. Kotłowski. The many faces of exponential weights in online learning. In Conference On Learning Theory, pages 2067–2092. PMLR, 2018.

[Cherief-Abdellatif et al., 2019]: B.E. Cherief-Abdellatif, P. Alquier, and M.E. Khan. A generalization bound for online variational inference. In Asian Conference on Machine Learning, pages 662–677. PMLR, 2019.

[Khan & Rue, 2023]: M.E. Khan and H. Rue. The Bayesian Learning Rule. In Journal of Machine Learning Research, 2023.
"...follows from two different approximations. One is known (and we credit Immer, Ollivier, and Tronarp et al.) while the other is new" 

What do you refer to with "the other"?