Variational Bayes Gaussian Splatting

## Abstract

Abstract Recently, 3D Gaussian Splatting has emerged as a promising approach for modeling 3D scenes using mixtures of Gaussians. The predominant optimization method for these models relies on backpropagating gradients through a differentiable rendering pipeline, which struggles with catastrophic forgetting when dealing with continuous streams of data. To address this limitation, we propose Variational Bayes Gaussian Splatting (VBGS), a novel approach that frames training a Gaussian splat as variational inference over model parameters. By leveraging the conjugacy properties of multivariate Gaussians, we derive a closed-form variational update rule, allowing efficient updates from partial, sequential observations without the need for replay buffers. Our experiments show that VBGS not only matches state-of-the-art performance on static datasets, but also enables continual learning from sequentially streamed 2D and 3D data, drastically improving performance in this setting.

## 1 Introduction

Representing 3D scene information is a long-standing challenge for robotics and computer vision (Özyeşil et al., ). A recent breakthrough in this domain relies on representing the scene as a radiance field, i.e., using neural radiance fields (Mildenhall et al., ). Recently, Gaussian Splatting (Kerbl et al., ) has demonstrated the effectiveness of mixture models for 3D scene representation, as highlighted by a surge of subsequent research (see Chen and Wang () for a survey). This approach leverages the ability of Gaussians to represent physical space as a collection of ellipsoids.

The dominant method for optimizing these models involves backpropagating gradients through a differentiable renderer with respect to the parameters of the mixture model. However, in many real-world applications, such as autonomous navigation, data is continuously streamed and must be processed sequentially. Traditional backpropagation-based methods are prone to catastrophic forgetting (French, ), leading to performance degradation as new data overwrites old knowledge. To mitigate this, replay buffers are often employed to retain and retrain on older data (Matsuki et al., ), which can be computationally expensive and memory intensive.

In this paper, we propose Variational Bayes Gaussian Splatting (VBGS), casting Gaussian Splatting as a variational inference problem over the parameters of a generative mixture model, enabling a closed-form update rule (Blei et al., ). Our approach inherently supports continual learning through iterative updates that are naturally accumulative, eliminating the need for replay buffers.

Our contributions are summarized as follows:

- •
We propose a generative model of mixtures of Gaussians for representing spatial 2D (image) and 3D (point cloud) data and derive closed-form variational Bayes update rules.
- •
We show that VBGS achieves a good model fit with a single parameter update step and enables continual learning without catastrophic forgetting.
- •
We demonstrate and benchmark our approach on both 2D image (TinyImageNet (Le and Yang, )) and 3D point cloud (Blender 3D models (Mildenhall et al., ) and Habitat scenes (Savva et al., )) datasets.

Our approach offers a robust alternative to gradient-based optimization by allowing continual, efficient updates from sequential data, making it well-suited for real-world applications.

## 2 Related Work

3D Gaussian Splatting (3DGS): Recently, 3DGS (Kerbl et al., ) introduced a novel approach for learning the structure of the world directly from image data by representing the radiance field as a mixture of Gaussians. This method optimizes the parameters of the Gaussian mixture model by backpropagating gradients through a differentiable rendering pipeline. The Gaussians serve as ellipsoid shape primitives for scene reconstruction, with each component associated with deterministic, learned features such as opacity and color. This approach has triggered significant advances in novel view synthesis and sparked interest in many applications across robotics, virtual/augmented reality, and interactive media (Fei et al., ; Chen and Wang, ).

While our method also represents scenes as a mixture of Gaussians, it diverges from 3DGS by modeling a distribution over all features and parameters. Rather than relying on gradient-based optimization through differentiable rendering, we frame the learning process as variational inference. This approach allows us to perform continual learning through closed-form updates, enabling more efficient and robust handling of sequential data.

Image representations: Gaussian mixture models have also been used for representing image data, particularly for data efficiency and compression (Verhack et al., ). Steered Mixture-of-Experts (SMoE) models, for instance, typically estimate parameters using the Expectation-Maximization (EM) algorithm, although some work has explored gradient descent optimization (Bochinski et al., ).

Our approach can be seen as an SMoE model on 2D image data, but again, instead of relying on EM or gradient descent, we propose variational Bayes parameter updates for more robust and efficient parameter estimation.

Continual learning aims to develop models that can adapt to new tasks from a continuous datastream without forgetting previously learned knowledge (Wang et al., ). Backpropagation-based methods, such as 3DGS, face significant challenges in this area due to catastrophic forgetting (French, ), where the model’s performance on prior data deteriorates as it trains on new. This issue is particularly evident in real-world scenarios, such as localization and mapping, where data is streamed (Matsuki et al., ; Keetha et al., ). The common mitigation strategy involves maintaining a replay buffer of frames to revisit past data during updates (Sucar et al., ; Fu et al., ).

In contrast, Bayesian approaches offer a more elegant solution for online learning, especially when dealing with conjugate models, allowing for exact Bayesian inference in a sequential setting (Jones et al., ). We leverage this property to update Gaussian splats continuously as new data arrives, eliminating the need for replay buffers.

## 3 Method

VBGS relies on the conjugate properties of exponential family distributions. This is well suited for variational inference, as we can derive closed-form update rules for inferring the (approximate) posterior. We first write down the functional form of these distributions and then describe the particular generative model for representing space and color.

If the likelihood distribution is part of the exponential family, it can be written down as:

$$ $\displaystyle p({x}|{\theta})=\phi({x})\exp\big{(}{\theta}\cdot T({x})-A({ \theta})\big{)},$ (1) $$

where ${\theta}$ are the natural parameters of the distribution, $T({x})$ is the sufficient statistic, $A({\theta})$ is the log partition function, and $\phi({x})$ is the measure function. This likelihood has a conjugate prior of the following form:

$$ $\displaystyle p({\theta}|{\eta}_{0},\nu_{0})=\frac{1}{Z({\eta}_{0},\nu_{0})} \exp\left({\eta}_{0}\cdot{\theta}-\nu_{0}\cdot A({\theta})\right),$ (2) $$

where $Z({\eta}_{0},\nu_{0})$ is the normalizing term. This distribution is parameterized by its natural parameters $({\eta}_{0},\nu_{0})$. The posterior is then of the same functional form as the prior, whose natural parameters can be calculated as a function of the sufficient statistics of the data and the natural parameters of the prior, i.e. $({\eta}_{0}+\sum_{{x}}T({x}),\nu_{0}+\sum_{{x}}1)$ (Murphy, ).

### 3.1 The generative model

Figure: Figure 1: The Generative Model: There are $K$ components in the mixture model, and $N$ observed data points, which are composed of a spatial vector ${s}$, and feature vector ${c}$. The parameters of the distributions of component $k$ ($\mu_{k}$, $\Sigma_{k}$) that generate ${s}$ and ${c}$ are random variables, which influence ${s}$ and ${c}$ respectively. $z$ is the associated mixture component for a given data point, dependent on the categorical parameters $\pi$. White and gray circles denote unobserved and observed variables respectively.

The generative model considered is a mixture model with $K$ components, each characterized by two conditionally independent modalities: spatial position (${s}$) and color (${c}$). For 2D images, ${s}$ represents the pixel location in terms of row and column coordinates (${s}\in\mathbb{R}^{2}$), while for 3D data, it corresponds to Cartesian coordinates (${s}\in\mathbb{R}^{3}$). The color component represents RGB values (${c}\in\mathbb{R}^{3}$) for both 2D and 3D data.
Both ${s}_{k}$, and ${c}_{k}$ are modeled as multivariate Normal distributions with parameters $(\mu_{{s},k},\Sigma_{{s},k})$, and $(\mu_{{c},k},\Sigma_{{c},k})$ respectively, and the components have a mixture weight $z$. The full generative model is visualized as a Bayesian network in Figure 1 and is factorized as:

$$ $\displaystyle p({s},{c},z,\mu_{{s}},\Sigma_{{s}},\mu_{{c}},\Sigma_{{c}},\pi)=$ $\displaystyle\bigg{(}\prod_{n=1}^{N}p({s}_{n}|z_{n},\mu_{{s}},\Sigma_{{s}})p({ c}_{n}|z_{n},\mu_{{c}},\Sigma_{{c}})p(z_{n}|\pi)\bigg{)}$ (3) $\displaystyle\bigg{(}\prod_{k=1}^{K}p(\mu_{k,{s}},\Sigma_{k,{s}})p(\mu_{k,{c}} ,\Sigma_{k,{c}})\bigg{)}p(\pi).$ (4) $$

The joint distribution in Equation (4) is factorized into two main components: one representing the likelihoods of $N$ data points $({s}_{n},{c}_{n})$ given their assignments, and another representing the priors over the mixture model parameters. The mixture components over space $p({s}|z,\mu_{{s}},\Sigma_{{s}})$ are conditionally independent from the component over color $p({c}|z,\mu_{{c}},\Sigma_{{c}})$, given mixture component $z$. The following distributions parameterize these random variables:

$$ $\displaystyle z_{n}$ $\displaystyle\sim\text{Cat}(\pi)$ (5) $\displaystyle{s}_{n}|z_{n}=k$ $\displaystyle\sim\text{MVN}(\mu_{k,{s}},\Sigma_{k,{s}})$ (6) $$

$$ $\displaystyle{c}_{n}|z_{n}=k$ $\displaystyle\sim\text{MVN}(\mu_{k,{c}},\Sigma_{k,{c}})$ (7) $$

The parameters of these distributions are also modeled as random variables. Treating the parameters as latent random variables allows us to cast learning as inference using the appropriate conjugate priors. The conjugate prior to a multivariate Normal (MVN) is a Normal Inverse Wishart (NIW) distribution, and the Dirichlet distribution is the conjugate prior to a categorical (Cat) distribution:

$$ $\displaystyle\mu_{k,{s}},\Sigma_{k,{s}}$ $\displaystyle\sim\text{NIW}(m_{0,{s}},\kappa_{0,{s}},V_{0,{s}},n_{0,{s}})$ (8) $\displaystyle\mu_{k,{c}},\Sigma_{k,{c}}$ $\displaystyle\sim\text{NIW}(m_{0,{c}},\kappa_{0,{c}},V_{0,{c}},n_{0,{c}})$ (9) $\displaystyle\pi$ $\displaystyle\sim\text{Dirichlet}(\alpha_{0}),$ (10) $$

See Appendix A.1 for a table with the values used for the hyperparameters.

To estimate the parameters, we infer their posterior distribution. However, as computing this posterior is intractable, we resort to variational inference (Jordan et al., ). We use a mean-field approximation to make inference tractable, which assumes that the variational posterior factorizes across the latent variables. This factorization allows for efficient coordinate ascent updates for each variable separately. Specifically, the variational distribution $q$ is decomposed as:

$$ $\displaystyle q(z,\mu_{{s}},\Sigma_{{s}},\mu_{{c}},\Sigma_{{c}},\pi)=$ $\displaystyle\bigg{(}\prod_{n=1}^{N}q(z_{n})\bigg{)}\bigg{(}\prod_{k=1}^{K}q( \mu_{k,{c}},\Sigma_{k,{c}})\bigg{)}\bigg{(}\prod_{k=1}^{K}q(\mu_{k,{s}},\Sigma _{k,{s}})\bigg{)}q(\pi),$ (11) $$

We select the distributions of the approximate posteriors to be from the same family as their corresponding priors. For the color variable ${c}$, we model the mean by a Normal distribution but keep the covariance fixed using a Delta distribution. This assures that the mixture components commit to a particular color and do not blend multiple neighboring colors in a single component.

$$ $\displaystyle q(z_{n})$ $\displaystyle=\text{Cat}(\gamma_{n})$ (12) $\displaystyle q(\mu_{k,{s}},\Sigma_{k,{s}})$ $\displaystyle=\text{NIW}(m_{t,{s}},\kappa_{t,{s}},V_{t,{s}},n_{t,{s}})$ (13) $\displaystyle q(\pi)$ $\displaystyle=\text{Dirichlet}(\alpha_{t})$ (14) $$

$$ $\displaystyle q(\Sigma_{k,{c}})$ $\displaystyle=\text{Delta}(\varepsilon I),$ (15) $\displaystyle q(\mu_{k,{c}})$ $\displaystyle=\text{Normal}(m_{t,{c}},\kappa_{t,{c}}^{-1}\varepsilon I)$ (16) $$

where the subscript $t$ indicates the parameters at timestep $t$, and $\varepsilon$ is a chosen hyperparameter.

Images can be generated using the mixture model by computing the expected value of the color, conditioned on a spatial coordinate ($\mathbb{E}_{p({c}|{s})}[{c}]$). For 3D rendering, we use the renderer from Kerbl et al. (), where the spatial component is first projected onto the image plane using the camera parameters, and the estimated depth is used to deal with occlusion. For more details, see Appendix B.

### 3.2 Coordinate ascent variational inference

To estimate the posterior over the model parameters, we maximize the Evidence Lower Bound (ELBO) $\mathcal{L}$ with respect to the variational parameters (Jordan et al., ):

$$ $\displaystyle\mathcal{L}=\mathbb{D}_{\text{KL}}[q(z,\mu_{{s}},\Sigma_{{s}},\mu _{{c}},\Sigma_{{c}},\pi)~{}||~{}p(z,\mu_{{s}},\Sigma_{{s}},\mu_{{c}},\Sigma_{{ c}},\pi|{s},{c})],$ (17) $$

This is done using coordinate ascent variational inference (CAVI) (Beal, ; Bishop, ; Blei et al., ) which contains two distinct steps, that are iteratively executed to optimize the model parameters. It parallels the well-known expectation maximization (EM) algorithm: the first step computes the expectation over assignments $q(z)$ for each data point.
Instead of computing the maximum likelihood estimate as is done in EM, in the second step we maximize the variational parameters of the posterior over model parameters. More in detail, in the first step the assignment for each data point $({s}_{n},{c}_{n})$ is computed by deriving the ELBO with respect to $q(z_{n})$.

$$ $\displaystyle\log q(z_{n}=k)=\log\gamma_{n}$ $\displaystyle=\mathbb{E}_{q(\mu_{k,{s}},\Sigma_{k,{s}})}[\log p({s}_{n}|\mu_{k ,{s}},\Sigma_{k,{s}})]$ (18) $\displaystyle\quad+\mathbb{E}_{q(\mu_{k,{c}},\Sigma_{k,{c}})}[\log p({c}_{n}| \mu_{k,{c}},\Sigma_{k,{c}})]$ (19) $\displaystyle\quad+\mathbb{E}_{q(\pi)}[\log p(\pi)]-\log Z_{n},$ (20) $$

where $Z_{n}$ is a normalizing term, i.e., if $\hat{\gamma}_{n}$ is the unnormalized logit, we can find the parameters of the categorical distribution as $\gamma_{n}=\frac{\hat{\gamma}_{n}}{\sum\hat{\gamma}_{n}}$.

In the second step, we compute the update of the approximate posteriors over model parameters. To this end, we derive the ELBO with respect to the natural parameters of the distributions. If the distribution is conjugate to the likelihood, we acquire updates of the following parametric form:
$\displaystyle\eta_{k}$
$\displaystyle=\eta_{0,k}+\sum_{x_{n}\in\mathcal{D}}\gamma_{k,n}T(x_{n})$
(21)
$\displaystyle\nu_{k}$
$\displaystyle=\nu_{0,k}+\sum_{x_{n}\in\mathcal{D}}\gamma_{k,n},$
(22)

where $\eta$ and $\nu$ are the natural parameters of the distribution over the likelihood’s natural parameters, and $T(x_{n})$ are the sufficient statistics of the data $x_{n}$.

For the approximate posterior $q(\mu_{{s},k},\Sigma_{{s},k})$ over the parameters of the spatial likelihood, the sufficient statistics are given by $T({s}_{n})=({s}_{n},{s}_{n}\cdot{s}_{n}^{T})$. The NIW conjugate prior consists of a Normal distribution over the mean and an inverse Wishart distribution over the covariance matrix. Hence, for each of the prior’s natural parameters, it has two values: $\eta_{0,{s}}=(\kappa_{0,{s}}\cdot m_{0,{s}},V_{0,{s}}+\kappa_{0,{s}}\cdot m_{0
,{s}}\cdot m_{0,{s}}^{T})$ and $\nu_{0,{s}}=(\kappa_{0,{s}},n_{0,{s}}+D_{{s}}+1)$. Here, $m_{0,{s}}$ is the mean of the Normal distribution over the mean, $\kappa_{0}$ is the concentration parameter over the mean, $n_{0,{s}}$ indicates the degrees of freedom, $V_{0,{s}}$ the inverse scale matrix of the Wishart distribution, and $D_{{s}}$ the dimensionality of the MVN.

For the approximate posterior $q(\mu_{{c},k},\Sigma_{{c},k})$ over the color likelihood, the sufficient statistics are again given by: $T({c}_{n})=({c}_{n},{c}_{n}\cdot{c}_{n}^{T})$. The prior is parameterized similarly as the NIW over ${s}$: $\eta_{0,{c}}=(\kappa_{0,{c}}\cdot m_{0,{c}},V_{0,{c}}+\kappa_{0,{c}}\cdot m_{0
,{c}}\cdot m_{0,{c}}^{T})$ and $\nu_{0,{c}}=(\kappa_{0,{c}},n_{0,{c}}+D_{{c}}+1)$. However, as we model the prior over $\Sigma_{k,{c}}$ as a delta distribution, we keep the values for $n_{k,{c}}$ and $V_{k,{c}}$ fixed.

Finally, the conjugate prior for the approximate posterior $q(\pi)$ over the component assignment likelihood $z$ is a Dirichlet distribution.
Here, sufficient statistics are given by $T(x)=1$, and the prior is parameterized by the natural parameter $\eta_{0,z}=\alpha$.

### 3.3 Continual updates

The proposed method supports continual learning because the parameters of each component are updated by aggregating the prior’s natural parameters with the data’s sufficient statistics via assignments $q(z)$. This iterative update process is order-invariant, allowing the model to adapt without forgetting past knowledge. Components without recent data assignments revert to their prior values, preserving flexibility in the model. Crucially, assignments $q(z)$ are always computed with respect to the initial posterior over parameters $q(\mu_{{s}},\Sigma_{{s}},\mu_{{c}},\Sigma_{{c}},\pi)$. This ensures that components without prior assignments can still be used to model the data.

Concretely, when data is processed sequentially, at each time point $t$, a batch of $\mathcal{D}_{t}$ of data points (${s}_{n},{c}_{n}$) is available, and the variational posterior over model parameters is updated. For each of these data points, the assignments $\gamma_{k,n}$ can be calculated similarly to Equation (18). We can rewrite the update steps from Equations (21) and (22), in a streaming way for $T$ timesteps:

$$ $\displaystyle\eta_{k}=\eta_{0,k}+\sum_{t=1}^{T}\sum_{x_{n}\in\mathcal{D}_{t}} \gamma_{k,n}T(x_{n})$ (23) $$

$$ $\displaystyle\nu_{k}$ $\displaystyle=\nu_{0,k}+\sum_{t=1}^{T}\sum_{x_{n}\in\mathcal{D}_{t}}\gamma_{k,n}$ (24) $$

Hence, we can write an iterative update rule for the natural parameters at timestep $t$ as a function of the natural parameters at $t-1$:

$$ $\displaystyle\eta_{t,k}=\eta_{t-1,k}+\sum_{x_{n}\in\mathcal{D}_{t}}\gamma_{k,n }T(x_{n})$ (25) $$

$$ $\displaystyle\nu_{t,k}$ $\displaystyle=\nu_{t-1,k}+\sum_{x_{n}\in\mathcal{D}_{t}}\gamma_{k,n}$ (26) $$

Here, $x_{n}$ is a placeholder for the particular data, e.g., when updating $q(\mu_{{s},k},\Sigma_{{s},k})$, this would be ${s}_{n}$. At $t=0$, the prior values for the natural parameters are used.

Note that when $\gamma_{k,n}$ is calculated using the initial parameterization of the variational posterior over parameters, applying these continual updates is identical to processing all the data in a single batch, avoiding the problem of catastrophic forgetting.

### 3.4 Component Reassignments

In many continual learning settings, the data statistics are not known in advance. The uniform initialization of the model might not adequately cover the data space, leading to a few components aggregating all the assignments with many components remaining unused. For example, consider a room where the largest object density is at the walls, while the center mainly consists of empty space. To mitigate this problem, we introduce a heuristic for reassigning the initial location of a components to data points which the model does not explain well.

The means of $n$ components for both space and color are replaced by the location and color of $n$ data points, sampled proportional to the negative ELBO under the current model. We choose $n$ as a fraction (5%) of these unused components. A component is unused if the concentration parameters ($\alpha_{k}$) of the prior over the mixture weights are equal to their prior value. By choosing $n$ as a fraction of the available components, more reassignments occur for the first frames when the object/scene is not yet known, and fewer near the end. Additionally, this ensures we do not reassign components already used at a different location. After updating the component means for these “initial” components, we can run CAVI as described in Section 3.2.

## 4 Results

We evaluate our approach by benchmarking it against backpropagating gradients through a differentiable renderer. In particular, we compare the following models on the Tiny ImageNet (Le and Yang, ), Blender 3D (Mildenhall et al., ) and Habitat rooms (Savva et al., ) datasets:

VBGS (Ours): We consider a variant of VBGS where the means of the initial posteriors ($q(\mu_{{s}},\Sigma_{{s}},\mu_{{c}},\Sigma_{{c}})$) are initialized on sampled points from the normalized data set (Data Init), as well as randomly initialized ($m_{k,{s}}\sim\mathcal{U}[-1,1]$, $m_{k,{c}}\sim\delta(0)$) (Random Init). For models with random initialization, data normalization is performed using estimated statistics; see Appendix C for further details.

Gradient: In 3DGS (Kerbl et al., ), the parameters are directly optimized using stochastic gradient descent on a weighted image reconstruction loss ($(1-\lambda)\cdot\text{MSE}+\lambda\cdot\text{SSIM}$). In the case of image data, $\lambda$ is set to 0 0; in the case of 3D data $\lambda$, it is set to $0.2$. We use spherical harmonics with no degrees of freedom, i.e., the specular reflections are not modeled. In order to be able to compare performance w.r.t. model size, we also do not do densification or shrinking (Kerbl et al., ) (for these results see Appendix D.2). When optimizing for images, we use a fixed camera pose at identity and keep the z-coordinate of the Gaussians fixed at a value of 1. Similar to the VBGS approach, we also consider a variant where the means of the Gaussian components are randomly initialized (Random Init), and on sampled points from the dataset (Data Init).

Figure: (a)
Refer to caption: extracted/5902389/Figures/image-psnr/boxplot_psnr_img.png

### 4.1 Images

We first evaluate performance on the Tiny ImageNet test set consisting of 10k images. Reconstruction accuracy, measured in PSNR (dB), is evaluated across varying numbers of components (Figure 2). Our results indicate that VBGS achieves reconstruction errors comparable to the gradient-based approach. We then evaluate the models in a continual learning setting, where data patches are sequentially observed and processed (Figure 2). While VBGS maintains consistent reconstruction quality across all observed patches, the gradient-based method disproportionately focuses on the most recent patch. The PSNR over timesteps is shown in Figure 3. VBGS converges to the same accuracy as training on the static dataset, whereas Gradient degrades due to catastrophic forgetting.

We further compared the computational efficiency by measuring the wall-clock time required for the Gradient approach to reach the performance level of VBGS after a single update step. This was computed over all images of the Tiny ImageNet validation set. We observe that VBGS is significantly (t-test, $p=0$) faster in wall clock time ($0.03\pm 0.03$ seconds) compared to Gradient ($0.05\pm 0.02$ seconds).

To assess performance in a continual learning setting, we divide each image into 8x8 patches, feeding the model one patch at a time (Figure 2). VBGS performs a single update per patch, whereas the Gradient approach performs 100 training steps per patch with a learning rate of 0.1. Figure 3 shows how the reconstruction PSNR evolves at each time step over the tiny ImageNet test set for the model with capacity 10K components. It’s important to note that the reconstruction error of VBGS after observing all the patches converges to the same value as observing the whole image at once, as the inferred posterior ends up being identical. Even though the gradient approach, when observing all data at once, achieves much higher PSNR, when the data is fed in continuously, it does not achieve that performance in the continual setting, as it always focuses on the latest observed patch (see Figure 2).

Figure: (a)
Refer to caption: extracted/5902389/Figures/continual/image_psnr.png

### 4.2 Blender 3D Objects

Next, we evaluate VBGS on 3D objects from the Blender dataset (Mildenhall et al., ). Model parameters are inferred using 200 frames which include depth information. VBGS is trained on the 3D point cloud, which is acquired by transforming the RGBD frame to a shared reference frame. In contrast, the gradient-based approach is optimized using multi-view image reconstruction. We evaluate reconstruction performance on 100 frames from the validation set. The results, measured as PSNR, for a model with a capacity of 100K components are shown in Table 1. Our results show improved performance for both methods when components are initialized from data rather than randomly. VBGS achieves performance similar to the gradient-based approach under Data Init conditions. In random initialization, our approach outperforms the gradient-based approach except for the “ship” object. For a more in-depth analysis where we vary the number of components, see Appendix D.1.

**Table 1: Prediction performance for the 3D dataset. Measured as PSNR (dB). Values ($\mu\pm\sigma$) are computed over 100 validation frames for each of the 8 blender objects. All models in this table have 100K components. The best performance for each column is marked in bold.**
|  | chair | drums | ficus | hotdog | lego | materials | mic | ship |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| VBGS<br>${}^{\text{(Data Init)}}$ | $22.82$<br>^±0.94 | $\mathbf{19.50}$<br>^±0.47 | $\mathbf{22.06}$<br>^±0.79 | $\mathbf{23.62}$<br>^±1.23 | $\mathbf{22.53}$<br>^±0.89 | $\mathbf{20.55}$<br>^±1.64 | $24.78$<br>^±0.60 | $21.23$<br>^±0.64 |
| VBGS<br>${}^{\text{(Random Init)}}$ | $21.35$<br>^±0.69 | $18.71$<br>^±0.43 | $21.49$<br>^±0.75 | $22.24$<br>^±0.97 | $20.59$<br>^±0.86 | $20.47$<br>^±1.36 | $24.42$<br>^±0.58 | $20.80$<br>^±0.90 |
| Gradient<br>${}^{\text{(Data Init)}}$ | $\mathbf{22.98}$<br>^±1.23 | $19.05$<br>^±0.57 | $21.08$<br>^±0.92 | $21.47$<br>^±1.67 | $19.97$<br>^±2.24 | $20.53$<br>^±1.64 | $\mathbf{25.25}$<br>^±0.90 | $\mathbf{23.55}$<br>^±1.20 |
| Gradient<br>${}^{\text{(Random Init)}}$ | $20.59$<br>^±0.85 | $15.04$<br>^±1.13 | $19.41$<br>^±0.90 | $19.81$<br>^±1.86 | $19.10$<br>^±1.02 | $16.11$<br>^±1.45 | $23.03$<br>^±0.72 | $21.15$<br>^±0.82 |

Figure: (a)
Refer to caption: extracted/5902389/Figures/volume-psnr/qualitative_3D_100k.png

Novel view predictions for the eight blender objects are shown in Figure 4. These renders are generated using a VBGS with 100K components, observed from a camera pose selected from the validation set. Note that these are rendered on a white background, and only the 3D object is modeled by the VBGS. Figure 4 shows the reconstruction performance as a function of the number of available components. It can be observed that for lower component regimes, VBGS renders patches of ellipsoids, while the gradient approach fills the areas more easily. We attribute this to a strong prior over the covariance shape encoded in the Wishart hyperparameters.

Finally, we also conduct the continual learning experiment for 3D and observe that the same properties from the 2D experiment hold, reaching an average reconstruction error over all objects of $11.19\pm 3.53$ dB for VBGS (Random Init) and $21.26\pm 1.76$ dB for Gradient (Random Init). Crucially, in the continual learning setting, the model does not have access to the data for initialization and has to be initialized randomly.

For the continual setting, 200 images are streamed in a continuous stream to the model, similar to the patches in the image experiment. For each observed (RGBD) image, VBGS applies a single update, and the gradient baseline applies 100 gradient steps with a learning rate of 0.1. We evaluate performance as the measured PSNR on novel-view prediction. Figure 3 shows the evolution of performance as a function of observed data. Note how the gradient approach’s performance deteriorates after observing more frames. This figure also indicates how VBGS converges quickly to an adequate level of performance after 50 steps.

### 4.3 Habitat Rooms

Figure: (a)
Refer to caption: extracted/5902389/Figures/habitat/rooms_qual.png

Finally, we evaluate on a dataset of 3D rooms, which is relevant for e.g. autonomous navigation for robotics in unseen environments. We consider three rooms from the Habitat test suite (“Van Gogh Room”, “Apartment 1” and “Skokloster Castle”) (Savva et al., ). The model parameters are inferred using 200 randomly sampled views from each room using the preprocessing pipeline of (Wang et al., ). On top of the models from the previous sections, we now also evaluate the model with component reassignments as detailed in Section 3.4.

Figure 5 shows the qualitative performance for the various considered models. Similar to previous experiments, we observe that initializing using the data yields the best performance. In contrast to the earlier experiments, random initialization does not capture well the structure of the room. This is mainly due to the randomly sampled initial means of the mixture model not covering the data well. For the Gradient approach, there are too few components that provide a gradient signal to optimize the observed view. In contrast, for VBGS, a single component might always be closer to the data than others, aggregating all the assignments. This is further evidenced by the noticeable performance improvement when using the reassign mechanism on top of VBGS. Quantitative results for all models can be found in appendix D.3

We also evaluate the models in a continual setting in Figure 5, where 200 frames (RBGD + pose) of the environment are observed sequentially. Note that these frames are randomly sampled from the environment and are not captured using a trajectory through the room. We observe that the gradient approach does not integrate the information well, even though the views are randomly sampled and cover the room. VBGS, without reassignments, does integrate the information and reaches performance on par with the Gradient trained in the non-continual setting. Finally, we observe that adding reassignments drastically improves performance on all three rooms.

## 5 Discussion and Conclusion

In this paper, we introduced a novel approach to optimizing Gaussian splats using variational Bayes. Our approach achieves comparable performance to backpropagation-based methods on 2D and 3D datasets, while offering key advantages for scenarios involving continual learning from streaming data.

One limitation of VBGS compared to 3D Gaussian Splatting (3DGS) is its reliance on RGBD data, as opposed to optimizing directly on RGB projections. For many use cases we envision, such as robot navigation, depth information is often readily available from stereo vision, lidar or other sensing technologies. Furthermore, most 3DGS approaches still rely on pre-computed camera parameters obtained from structure-from-motion algorithms (Özyeşil et al., ), which implicitly also provide depth information through triangulation.
One could also use a pretrained neural network that predicts depth from monocular RGB data (Yang et al., ), as evidenced by Fu et al. ().

Additionally, 3DGS dynamically adjusts the model size by growing and shrinking as needed, a feature that we have not yet incorporated. In future work, we aim to explore principled approaches for dynamic model sizing, leveraging model evidence (Friston et al., ) to guide this process.

While VBGS requires only a single update step per observation, this step is computationally more expensive than a single gradient descent iteration in backpropagation. For extremely large datasets that do not fit into memory, VBGS may need to process data in distinct update steps, where stochastic gradient descent (SGD) might offer better efficiency. However, in streaming or partial data settings, VBGS excels by allowing updates without replay buffers, opening avenues for active data selection, a potential future research direction. Investigating active learning mechanisms, particularly in robotic SLAM or autonomous systems, could optimize data usage, ensuring that only the most informative observations are processed.

The continual learning capabilities of VBGS make it especially suited for applications requiring real-time adaptation and learning, such as autonomous navigation, augmented reality, and robotics. By continuously integrating new information without the need for replay buffers, VBGS reduces the computational burden typically associated with processing large-scale data streams. Furthermore, the variational posterior over model parameters enables parameter-based exploration (Schwartenbeck et al., ), allowing agents to efficiently explore and build a model of their environment in real-time. This suggests exciting opportunities for VBGS to contribute to adaptive, real-time learning systems in various real-world domains.

## 6 Code availability

The code for using VBGS is available in the vbgs repository, which can be found at the following link: [https://github.com/VersesTech/vbgs](https://github.com/VersesTech/vbgs).

## 7 Acknowledgments and Disclosure of Funding

The authors would like to thank the members of the VERSES AI Research lab for critical discussions and feedback that improved the quality of this work, with special thanks to Jeff Beck, Tommaso Salvatori, and Conor Heins.

## Appendix A Hyperparameters

### A.1 Prior parameters

The considered conjugate priors over the likelihood parameters are parameterized by the canonical parameters shown in Table 2. Some of these values are a function of the number of available components, indicated by nc.

**Table 2: Parameters of the conjugate priors over likelihood parameters. Some values are a function of the number of components, indicated by nc, $I$ represents the identity matrix of size the multivariate dimension $D$. Parameters are in the canonical form of the corresponding distribution.**
|  |  | 2D | 3D |
| --- | --- | --- | --- |
| $p(\mu_{k,{s}},\Sigma_{k,{s}})$ | $m_{{s},k}$ | $\mathbf{0}$ | $\mathbf{0}$ |
| $\kappa_{{s},k}$ | $10^{-2}\cdot\mathbf{1}$ | $10^{-2}\cdot\mathbf{1}$ |  |
| $V_{{s},k}$ | $2.25\cdot 10^{4}\cdot\text{nc}\cdot I$ | $2.25\cdot 10^{6}\cdot\text{nc}\cdot I$ |  |
| $n_{{s},k}$ | 4 | 5 |  |
| $p(\mu_{k,{c}},\Sigma_{k,{c}})$ | $m_{{c},k}$ | $\mathbf{0}$ | $\mathbf{0}$ |
| $\kappa_{{c},k}$ | $10^{-2}\cdot\mathbf{1}$ | $10^{-2}\cdot\mathbf{1}$ |  |
| $V_{{c},k}$ | $10^{6}\cdot I$ | $10^{8}\cdot I$ |  |
| $n_{{c},k}$ | 5 | 5 |  |
| $p(\pi)$ | $\alpha_{k}$ | $\frac{1}{\text{nc}}$ | $\frac{1}{\text{nc}}$ |

### A.2 Initial parameters

In the first step of the coordinate ascent algorithm, $q(z)$ is inferred with an initial configuration of $q(z,\mu_{{s}},\Sigma_{{s}},\mu_{{c}},\Sigma_{{c}},\pi)$. Crucially, this configuration is distinct from the prior described in Appendix A.1. The initial canonical parameters are shown in Table 3.

The values of ${m}_{{s},k\text{init}}$ and ${m}_{{c},k\text{init}}$ vary for the two considered cases. When the means are initialized on the data point (Data Init), ${m}_{{s},k\text{init}}$ and ${m}_{{c},k\text{init}}$ are set to $K$ values $({s}_{n},{c}_{n})$ sampled from the data $\mathcal{D}$. When the means are randomly initialized (Random Init), then $m_{k,{s},\text{init}}\sim\mathcal{U}[-1,1]$, $m_{k,{c},\text{init}}\sim\delta(0)$).

**Table 3: Parameters of the initial approximate posteriors over likelihood parameters. Some values are a function of the number of components, indicated by nc, $I$ represents the identity matrix of size the multivariate dimension $D$. Parameters are in the canonical form of the corresponding distribution.**
|  |  | 2D | 3D |
| --- | --- | --- | --- |
| $q(\mu_{k,{s}},\Sigma_{k,{s}})$ | $m_{{s},k}$ | $\mathbf{m}_{{s},k,\text{init}}$ | $\mathbf{m}_{{s},k,\text{init}}$ |
| $\kappa_{{s},k}$ | $10^{-5}\cdot\mathbf{1}$ | $10^{-6}\cdot\mathbf{1}$ |  |
| $V_{{s},k}$ | $2.25\cdot 10^{4}\cdot\text{nc}\cdot I$ | $2.25\cdot 10^{6}\cdot\text{nc}\cdot I$ |  |
| $n_{{s},k}$ | 4 | 5 |  |
| $q(\mu_{k,{c}},\Sigma_{k,{c}})$ | $m_{{s},k}$ | $\mathbf{m}_{{c},k,\text{init}}$ | $\mathbf{m}_{{c},k,\text{init}}$ |
| $\kappa_{{s},k}$ | $10^{-2}\cdot\mathbf{1}$ | $10^{-2}\cdot\mathbf{1}$ |  |
| $V_{{c},k}$ | $10^{6}\cdot I$ | $10^{8}\cdot I$ |  |
| $n_{{c},k}$ | 5 | 5 |  |
| $q(\pi)$ | $\alpha_{k}$ | $\frac{1}{\text{nc}}$ | $\frac{1}{\text{nc}}$ |

## Appendix B Image rendering

Rendering is the process of generating an image, given an internal representation. Typically, this refers to the process of projecting from a 3D representation to the image plane. In our generative model, this process boils down to evaluating the expected color value for each considered pixel, i.e. $\mathbb{E}_{p({c}|{s})}[{c}]$. This is straightforward to compute for image data:

$$ $\displaystyle\mathbb{E}_{p({c}|{s})}[{c}]$ $\displaystyle=\sum_{k}\bigg{(}p(z_{n}=k|{s})\sum_{{c}_{i}}{c}_{i}\underbrace{p ({c}_{i}|z_{n}=k)}_{\approx\delta({c}_{k})}\bigg{)}\approx\sum_{k}p(z_{n}=k|{s }){c}_{k},$ (27) $$

where we approximate the distribution over the color features by a delta distribution positioned at the mean of $q(\mu_{{c}})$.

In equation (27), the distribution is conditioned on the spatial location of the pixel. For rendering in 3D, we leverage the computationally efficient 3D renderer designed by Kerbl et al. []. Here, the 3D Gaussians are first projected to the image plane, and by using alpha blending along a casted ray, color values are combined into a pixel color. As we do not optimize on image reconstruction, our approach does not have sensible alpha blending. We, therefore, set the alpha value for all Gaussians at 1, i.e., all components are opaque.

## Appendix C Data Normalization

Before training, the data is normalized to have zero mean and a standard deviation of one. This is not strictly necessary but ensures that we can use the same hyperparameters and initial parameters for all models. When all data is readily available, the data is simply normalized using the statistics from the data itself. Note that each dimension is considered individually.

When the data is not available, we assume that the random variable is distributed uniformly within a certain range: $x\sim\mathcal{U}(r_{\text{min}},r_{\text{max}})$. The normalized value is then calculated as:

$$ $\displaystyle\hat{x}=\frac{x-\mathbb{E}[x]}{\sqrt{\text{Var}[x]}},$ (28) $$

for which the ranges for each of the parameters are displayed in Table 4.

**Table 4: Range for data normalization. The subscript i indicates a single dimension of the vector.**
|  | Range 2D | Range 3D Objects | Range 3D Rooms |
| --- | --- | --- | --- |
| ${s}_{i}$ | $[0,64]$ | $[-1,1]$ | $[-5,5]$ |
| ${c}_{i}$ | $[0,255]$ | $[0,255]$ | $[0,255]$ |

## Appendix D Additional results

### D.1 Reconstruction performance as a function of number of components

Figure: (a)
Refer to caption: extracted/5902389/Figures/image-psnr/qualitative.png

Figure 6 shows the reconstruction of an image from the TinyImageNet dataset using various model sizes for both VBGS and the gradient model. Notice how in low component regimes, the initialization on data yields smoother Gaussian components (e.g., in the 100 components regime); however when the capacity increases, this impact is drastically reduced. In high component regimes, the gradient approach yields a smoothed version of the surface, while the VBGS models better capture the high-frequency textures.

We evaluated reconstruction performance on the Blender 3D models for a variety of models, as a function of model size, measured as the amount of available components. The results are shown in Figure 6.

### D.2 Growing and shrinking of 3D Gaussian Splatting

**Table 5: Dynamic 3DGS. The first row shows reconstruction performance, measured in PSNR (dB), evaluated on the validation set of the eight objects of the Blender set. The second row for each model shows the resulting amount of components after optimization.**
|  | chair | drums | ficus | hotdog | lego | materials | mic | ship |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Gradient<br>${}^{\text{(Data Init)}}$ | 26.73<br>^±1.15 | 20.65<br>^±0.71 | 22.93<br>^±1.09 | 27.92<br>^±0.84 | 26.94<br>^±1.10 | 17.02<br>^±1.48 | 28.10<br>^±0.81 | 25.91<br>^±1.62 |
|  | 456K | 378K | 281K | 177K | 337K | 46K | 186K | 263K |
| Gradient<br>${}^{\text{(Random Init)}}$ | 26.72<br>^±1.16 | 20.65<br>^±0.69 | 22.94<br>^±1.10 | 28.06<br>^±0.80 | 26.98<br>^±1.11 | 17.03<br>^±1.48 | 28.26<br>^±0.79 | 25.81<br>^±1.59 |
|  | 455K | 384K | 282K | 177K | 339K | 53K | 182K | 262K |

3D Gaussian Splatting, in its default implementation, dynamically grows and shrinks the model. In this paper, we kept the number of components fixed to make a comparison as a function of model size. Here, we optimize 3D Gaussian Splats using the gradient-based approach with dynamic model sizes. The results are reported in Table 5. We can see that the reconstruction quality achieved is much higher, but they also require a larger number of components.

### D.3 Novel view prediction performance on Habitat rooms

The results for the various model configurations with 100K components on the rooms data are shown in Table 6. We observe that the performance of all models lies within the same range, with the gradient descent approach yielding better performance when initialized on the data. However, VBGS performs better when initialized randomly. Using component reassignments also pushes the performance of Random Init close to Data Init while outperforming it for “Apartment 1”.

**Table 6: Performance for novel view prediction on the Habitat rooms. Measured as PSNR (dB). Values ($\mu\pm\sigma$) are computed over 100 validation frames for each of the three rooms. All models in this table have 100K components. The best performance for each column is marked in bold.**
|  | Van Gogh Room | Skokloster Castle | Apartment 1 |
| --- | --- | --- | --- |
| VBGS<br>${}^{\text{(Data Init)}}$ | $19.83$<br>^±4.56 | $\mathbf{18.66}$<br>^±2.02 | $22.92$<br>^±3.33 |
| VBGS<br>${}^{\text{(Random Init - Reassign)}}$ | $19.78$<br>^±4.49 | $17.93$<br>^±2.52 | $22.98$<br>^±3.72 |
| VBGS<br>${}^{\text{(Random Init)}}$ | $18.98$<br>^±3.89 | $15.62$<br>^±2.27 | $20.38$<br>^±4.53 |
| Gradient<br>${}^{\text{(Data Init)}}$ | $\mathbf{20.97}$<br>^±4.29 | $18.45$<br>^±1.89 | $\mathbf{24.96}$<br>^±6.17 |
| Gradient<br>${}^{\text{(Random Init)}}$ | $20.55$<br>^±4.02 | $14.30$<br>^±4.36 | $19.18$<br>^±5.32 |

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: **Q1:** Why not evaluate using datasets and protocols from recent SLAM works such as (Matsuki et al., 2024; Keetha et al., 2024)?

**Q2.** When should one expect the mean-field approximation to hurt performance?
soundness: 3
strengths: **S1.** The method addresses an important limitation of 3DGS.

**S2.** The method is principled and presented with clarity.

**S3.** The experimental validation lends support to the practical usefulness of the method: in particular, in the online setting 3DGS is clearly impacted by catastrophic forgetting while VBGS is not.
summary: This paper proposes Variational Bayes Gaussian Splatting (VBGS) for modeling 3D scenes by framing training of Gaussian Splats as variational inference. The formulation enables closed-form variational update of model parameters which is particularly beneficial in online settings. In particular, updating the model from a sequence of observations is equivalent to processing all data in a single batch. Thus, catastrophic forgetting is circumvented.

Experimental results indicate VBGS matches the performance of 3D Gaussian Splatting (3DGS) on static datasets while enabling continual learning with streamed 3D data.
weaknesses: **W1.** The method is motivated by the application to continual learning scenarios such as SLAM. However, the experimental validation does not compare against recent SLAM methods or datasets, e.g., (Matsuki et al., 2024; Keetha et al., 2024).

**W2.** The significance of the mean-field approximation is not discussed or investigated experimentally.

**W3.** The method's reliance on depth maps is, as acknowledged by the authors, a limitation of the method. The author’s propose inferring depth maps from RGB data using a pretrained model. Experiments along these lines would have been interesting to include.

---

confidence: 5
contribution: 3
initial_rating: 5
presentation: 2
questions: These are minor concerns that are not counted in my final scores.

- Having only three components in Figure 2.a without Gradient (Data init) seems incomplete. Please have it updated.
- As the authors have demonstrated with a continual learning framework, is it workable to scene in motion?
- Does this method scale up? Can this method scale up to tens of millions of Gaussian splatting for larger datasets (e.g., MipNeRF-360)?
soundness: 3
strengths: The main strengths I find in this paper lie in its novelty, clarity in writing, and reproducibility.

- Both relating 3D Gaussian splatting as Gaussian mixture model and variational inference is new to this field, and this novelty should be acknowledged.
- Application of continual learning well justifies the usage of variational inference framework.
- The training scheme is described in detail and the described motivation and methodology is quite clear.
- The supplementary material has the code implementation for reproducibility.
summary: The paper proposes a reformulation of 3D Gaussian splatting with variational inference framework. This allows a closed-form update rule and applicability to continual learning framework. The model is compared with gradient-based update rule in small-scale (~100k) 3DGS fitting to Habitat test suite and NeRF synthetic dataset and shows performance boost with better robustness to randomized initialization.
weaknesses: Although I find the method informative, there are several points to be clarified in order to improve the presentation.

- Common 3DGS involves several millions of Gaussians, but the presented experiments are only based on 100k Gaussians, which is very slim. This small parameter constraint makes the PSNR scores in Table 1 way below the numbers reported by other NeRF/3DGS papers which go high above 30 dB for the same NeRF synthetic dataset. Is there any specific reason (e.g., computational burden) that constraints the experiments? If so, this should be specified.
- Unlike the authors’ claim in Appendix D.3, Table 6 shows that the presented VBGS shows comparable or worse results than Gradient-based method, which needs further clarification.
- It seems like the method does not incorporate cloning/pruning/splitting of Gaussians which take crucial roles in the original 3DGS. By scaling up the Gaussians during the training, the performance can be boosted significantly, so unless there is a counterpart in VBGS for resizing the Gaussian pool, I believe there is a little advantage of using this framework in practical scenarios. I believe there should be an experiment with VBGS having initial Gaussians of the same size of the final Gaussians in Table 5, so we can compare fairly. For example, if we start from 456K Gaussians in VBGS and performs better than dynamically resized 3DGS in Table 5, we can clearly say that there is a reason to migrate from 3DGS to VBGS.
- I believe it is better to have a dedicated discussion on changes in computational demand relative to the size of Gaussians in different frameworks. Inference time should be the same, but training time and the memory requirement will be very informative to followers to this work.

In summary, unless there is a clear advantage of VBGS over standard size-varying gradient-based 3DGS, the method is not practically persuasive, yet. Furthermore, currently, there is no guarantee of scalability, which I believe important in learnable 3D representation. It is OK to have a slow and memory-intensive method if the performance (PSNR/SSIM) is stronger, but regarding Table 6, the presented method is not yet persuasive. Regarding these, the current version of the manuscript seems not ready to be published despite its novel and intuitive idea. Please note that I am flexible with my score, and looking forward to further discussion in the rebuttal phase.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: As stated in the weaknesses above, I have some additional questions as follows:
1. Why is it said that 3DGS struggles with catastrophic forgetting when dealing with continuous streams of data?

2. Could you provide numerical results on the HABITAT ROOMS dataset?
soundness: 2
strengths: 1. Effectively handles continuous data streams: By framing the training process as variational inference over model parameters, VBGS avoids the issue of catastrophic forgetting typically seen in continuous data streams with traditional methods.

2. No need for replay buffers: VBGS leverages the conjugacy properties of multivariate Gaussians, deriving a closed-form variational update rule that allows efficient incremental updates without requiring replay buffers.
summary: This paper introduces Variational Bayes Gaussian Splatting (VBGS), a new method for modeling 3D scenes with Gaussian splats that addresses the challenge of catastrophic forgetting in continuous data streams. Unlike traditional methods that rely on backpropagation through a differentiable renderer and require replay buffers, VBGS formulates training as variational inference over model parameters, enabling efficient updates through a closed-form variational update rule.
weaknesses: 1. Limited improvement. In Tab. 1, your method is comparable to the Gradient method (3DGS), which cannot prove the effectiveness of your method.

2. Unlike 3D Gaussian Splatting (3DGS), it relies on RGBD data. It means your method needs a depth camera and does not outperform 3DGS.

3. The experiments are not comprehensive, as they do not include evaluations on standard datasets like mip-NeRF 360, which could provide a more robust comparison with existing methods.

4. The comparison with baselines is limited; a more thorough evaluation would include comparisons with established methods such as Mip-NeRF 360, Instant-NGP, and NeRF. Including these methods as baselines would provide a clearer understanding of VBGS's performance relative to a broader range of state-of-the-art approaches.

---

confidence: 3
contribution: 1
initial_rating: 3
presentation: 2
questions: Please refer to the weaknesses.
soundness: 1
strengths: 1 The paper is well-structured.
2 The idea of introducing the conjugacy properties of multivariate Gaussians seems to make sense.
summary: This paper proposes Variational Bayes Gaussian Splatting (VBGS), which frames training a Gaussian splat as variational inference over model parameters.
weaknesses: 1 In the introduction, and background (e.g., 3D Gaussian Splatting) and research gap are not detailed enough, so the author's research motivation cannot be understood clearly.
2 The novelty is unclear and insufficient, it's confused why the proposed method can eliminate the need for replay buffers. More relevant discussions and technical details should be analyzed.
3 The authors claim that this work can "more robust and efficient parameter estimation." However, no computational efficiency analysis (e.g., computational time, memory usage) is provided.
4 In the experiments, the baseline only includes "Gradient". More discussions are need to explain why choosing"Gradient". Besides, the superiority of the proposed method over "Gradient" needs more analysis from both theoretical and experimental perspectives. Are there any more recent and advanced baselines?
5 The evaluation metrics should be introduced in detail, e.g., PSNR, to enhance the understandability and readability of the experimental part.

Meta-review: The paper proposes a variational Bayesian approach to learning Gaussian splats. The main claim is that the method matches state-of-the-art for static datasets but works much better for continual and streaming data. 

The work shows promising results, especially in mitigating the 'catastrophic forgetting' encountered when running traditional Gaussian Splatting methods in a streaming setting. 

However, the experimental part seems too preliminary to warrant publication at ICLR. Comparisons to baselines (as used in SLAM methods, which work in a streaming setting) are mostly missing, and the results are rather small scale. 

If these aspects are improved, this work will make a strong submission. However, at this stage, I recommend to reject the paper.

Decision: Reject