Bernstein Flows for Flexible Posteriors in Variational Bayes

—Variational inference (VI) is a technique to approximate difficult to compute posteriors by optimization. In contrast to MCMC, VI scales to many observations. In the case of complex posteriors, however, state-of-the-art VI approaches often yield unsatisfactory posterior approximations. This paper presents Bernstein flow variational inference (BF-VI), a robust and easy-to-use method, flexible enough to approximate complex multivariate posteriors. BF-VI combines ideas from normalizing flows and Bernstein polynomial-based transformation models. In benchmark experiments, we compare BF-VI solutions with exact posteriors, MCMC solutions, and state-of-the-art VI methods including normalizing flow based VI. We show for low-dimensional models that BF-VI accurately approximates the true posterior; in higher-dimensional models, BF-VI outperforms other VI methods. Further, we develop with BF-VI a Bayesian model for the semi-structured Melanoma challenge data, combining a CNN model part for image data with an interpretable model part for tabular data, and demonstrate for the first time how the use of VI in semi-structured models.


I. INTRODUCTION
Uncertainty quantification is essential, especially if model predictions are used to support high-stakes decision-making.Quantifying uncertainty in statistical or machine learning models is often achieved by Bayesian approaches, where posterior distributions represent the uncertainty of the estimated model parameters.Determining the exact posterior distributions is often impossible when the posterior takes a complex shape and the model has many parameters.This is especially true for complex models such as Bayesian neural networks (NNs) or semi-structured models that combine an interpretable model part with deep NNs.Variational inference (VI) is a commonly used approach to approximate complex distributions through optimization [1,2].In VI, the complex posterior is approximated by a variational distribution by minimizing a divergence measure between the variational and the true posterior distribution.VI is currently a very active research field tackling different challenges, which can be categorized into the following groups: (1) constructing variational distributions that are flexible enough to match the true posterior distribution, (2) defining optimal variational objective for tuning the variational distribution, which boils down to finding the most suited divergence measure quantifying the difference between a variational distribution and posterior, and (3) developing robust and accurate stochastic optimization frameworks for the variational objective [3,4,5].Here, we focus on challenge (1) and propose a method to construct a variational distribution that is flexible enough to accurately and robustly approximate complex multidimensional posteriors.
To avoid model-specific calculations, we design our method as a Black Box VI (BBVI) approach [6].In BBVI, the approximative posterior is determined by stochastic gradient descent.The user simply defines the Bayesian model by specifying the likelihood and the prior, after which all subsequent calculations are carried out automatically.Due to its simplicity, BBVI is implemented in many packages for Bayesian modeling, like Stan [7] and Pyro [8] as an alternative to MCMC.Given BBVI's scalability to large datasets and its widespread applicability, it has emerged as the preferred technique in the field of machine learning [5].
Our approach uses transformation models (TMs) to construct complex posteriors.Transformation models (TMs) have been introduced for fitting potentially complex outcome distributions for probabilistic regression models [9].Since then, they have been mainly used to model different outcome types, such as ordinal [10,11], count [12], continuous [13], or time-to-event outcomes [14] based on tabular predictors.Moreover, TMs have been used to model multidimensional distributions [15].Neural networks can be used to extend TMs to model outcomes for unstructured predictors (e.g., images or text) or a combination of tabular and unstructured predictors [10,16,17,18].
The basic idea of TMs is to learn a flexible and monotone transformation function that transforms between a simple latent distribution and a potentially complex conditional outcome distribution.In TMs, the transformation function is parameterized as an expansion of basis functions.In the case of continuous target distributions, most often, Bernstein polynomials [19] are used because they can easily be constrained to be strictly monotone, and their flexibility can be tuned via the order M .A large order M ensures an accurate approximation of the distribution, which is robust against a further increase of M [20,21]; this is also demonstrated in our experiments for the BBVI setting.
Independently of TMs, normalizing flows (NFs) have been developed in the deep learning community.NFs and TMs rely on the same idea, but NFs usually construct the transformation by chaining many simple functions, while TMs construct one rather complex transformation function.In NFs, each simple function, such as shifting and scaling, incrementally adds to the complexity of the final transformation.Among the prominent NF implementations are RealNVP [22] and Masked Autoregressive Flow (MAF) [23].RealNVP stands out for its efficient, invertible transformations facilitated by a specialized neural network architecture.Its key advantage lies in the efficient computation of the Jacobian matrix's determinant, essential for direct density estimation in the change of variable function (see 6).This efficiency is achieved by iteratively splitting the components of the data into two parts.In each step, the first part of the components is used to train a neural network computing the scale and shift parameters of the transformation.This transformation is then applied to the other components, while the first part remains unaltered.This procedure is repeated multiple times with different partitioning, leading to a triangular Jacobian, thus enabling efficient and invertible transformations.In contrast, MAF adopts a fundamentally different approach to construct transformations [23].It utilizes a sequential (autoregressive) framework, facilitated by neural networks.In MAF each output component relies exclusively on its preceding components, a concept often referred to as causality in this context.This design also leads to a triagonal Jacobian matrix and thus a fast computation of the change of variable equation.The MAF ensures that the nth output of NN is solely dependent on the first n − 1 inputs, yielding an autoregressive model.However, some NF approaches use a single flexible transformation, such as sum-of-squares polynomials [24] or splines [25].Recently, also Bernstein-based polynomials have been used for modeling unconditional multivariate density distributions [21].
NFs were initially introduced for variational inference to approximate potentially complex distributions of latent variables in models such as variational autoencoders [26,27].In the past, oftenmembers from simple distribution families have been used to approximate the posterior in BBVI.In the "Bayes by Backprop" method, Blundell et al. used independent Gaussians to approximate the posterior of the weights in a Bayesian Neural Network (BNN).They determined the parameters of these Gausians using BBVI [28].This approach was made more flexible by using a multivariate Gaussians [29] as variational distribution.While it is clear that TMs or NFs have the potential to construct flexible variational distributions, the first attempts to use NF-based BBVI were proposed only recently [30].These NF-based BBVI approaches compare favorably against existing BBVI methods but require a complex training scheme, and sometimes exhibit pathological behavior [31].
Here, we introduce Bernstein flow variational inference (BF-VI), which, for the first time, uses TMs in BBVI.We use TMs based on Bernstein polynomials to construct a variational distribution that closely approximates a potentially complex posterior in Bayesian models.The proposed method is computationally efficient and applicable to typical statistical models.The proposed method yields superior results in our experiments compared to existing NF approaches [31].Using BF-VI we further demonstrate, for the first time, how VI can be used to fit Bayesian semi-structured models where interpretable statistical model parts (based on tabular data) and deep NN model parts (based on images) are jointly fitted.We define our method in section II-A for one-dimensional examples and generalize it to Bayesian models with multivariate posteriors in section II-B.In section III, we benchmark our BF-VI approach against exact Bayesian models, MCMC-Simulations, Gaussian-VI, and NF-based BBVI, showing accurate posterior approximations in low dimensions and superior approximations in higher dimensions when compared to NF-based BBVI and summarize in section IV.

II. BERNSTEIN FLOW VARIATIONAL INFERENCE
In the following, we describe the Bernstein Flow-VI (BF-VI) approach, which we propose for accurately and robustly approximating potentially complex posteriors in Bayesian models.The main idea is to enable the VI procedure to approximate the joined posterior of the p model parameters by a flexible variational distribution.This is done by modeling the transformation function from a predefined simple latent distribution to a potentially complex variational distribution.The number of parameters p in the Bayesian model determines the dimension of both the latent and the variational distribution.
We first explain BF-VI for Bayesian models with a single parameter and hence a one-dimensional posterior and then generalize to models with multivariate posteriors.The code is publicly available on GitHub1 .

A. One-Dimensional Bernstein Flows
BF-VI approximates the bijective transformation function g : Z → θ between a latent variable Z ∈ R with predefined distribution F Z : R → [0, 1] with log-concave and continuous density f Z , and the model parameter θ ∈ R with a potentially complex distribution F θ : R → [0, 1] so that F Z (z) = F θ (g(z)).Fig. 1 visualizes this transformation on the scale of the densities, where f Z (z) = f θ (g(z)) | ∂g (z)  ∂z | according to the change-of-variable formula.
Hothorn et al. [20] give theoretical guarantees for the existence and uniqueness of g = F −1 θ • F Z .However, g cannot be computed directly if F θ is not known (in our application F θ is the unknown distribution of the posterior).The core of BF-VI is to approximate g, shown in Fig. 1, by f BP Bernstein polynomials (BP) 2 as with Be i (z) = Be (i+1,M −i+1) (z) being the density of a Beta distribution with parameters i + 1 and M − i + 1.To preserve the bijectivity of g, we use in BF-VI w.l.o.g. a strict monotone increasing BP to approximate g.With the approximation of the transformation function g by f BP , it holds that f Z (z) can be approximated by f θ (f BP (z)) | ∂fBP(z) ∂z |.Using a BP for approximating g gives the following theoretical guarantees [32] 1) With increasing order M of the BP, the approximation f BP to g gets arbitrarily close (the BP have been introduced for this very purpose in the constructive proof of the Weierstrass theorem by [19]); 2) the required strict monotonicity of the approximation f BP can be easily achieved by constraining the coefficients ϑ i of the BP to be increasing; 3) BPs are robust versus perturbations of the coefficients ϑ i ; 4) the approximation error decreases with 1/M (Voronovskaya Theorem).See [19,32] for more detailed discussions of the beneficial properties of BPs in general and [20,21] for transformation models.
While the output of f BP (z) is unrestricted, a BP requires a input z within [0, 1].We experimented with several approaches to ensure the restriction z ∈ [0, 1] resulting in slightly different behavior during the training (see Appendix B2).Based on these experiments, we decided to obtain z ∈ [0, 1] by sampling values from a standard normal distribution, z ′′ ∼ N (0, 1), then apply the affine transformation l(z ′′ ) = α • z ′′ + β, followed by a sigmoid σ(z ′ ) = 1/(1 + e −z ′ ).Altogether, we approximate 2 Some authors make a distinction between Bernstein polynomials in which ϑ i is fixed by the values of the function to be approximated and call expressions like expressions like in (1), where ϑ i is a fitting parameter, polynomials of Bernstein type.the transformation g by f : Z → θ via f = f BP • σ • l, which we call Bernstein flow.
To allow the application of unconstrained stochastic gradient descent optimization, which is typically used in the deep learning domain, we enforce the strict monotonicity of the flow f as follows: We optimize unrestricted parameters of f , i.e., ϑ ′ 0 , . . .ϑ ′ M , α ′ , β ′ , and apply the following transformations to determine the parameters of the bijective flow: ϑ 0 = ϑ ′ 0 , and ϑ i = ϑ i−1 + softplus(ϑ ′ i ) for i = 1, . . ., M for getting a strictly increasing BP and α = softplus(α ′ ), β = β ′ for getting an increasing affine transformation.
In the Appendix ??, we show that the resulting variational distribution is a tight approximation to the posterior in the sense that the KL divergence between q λ (θ) and p(θ | D) decreases with the order of the BP via 1/M .

B. Multivariate Generalization
In the case of a Bayesian model with p parameters, θ 1 , θ 2 , . . ., θ p , the Bernstein flow bijectively maps a pdimensional Z ′ to a p-dimensional θ.We realize this flow by choosing p independent standard normal Gaussians as simple latent distribution for the p-dimensional Z ′ and apply on each component an affine transformation followed by a sigmoid function to achieve a [0,1] restricted Z.The possible complex dependencies in θ are modeled in the multivariate generalization f BP of the one-dimensional Bernstein Polynomial (see eq. 2 for the definition of the j-th component of f BP ).
To achieve an efficient computation, we use a triangular map for constructing coefficients ϑ j i j = 2, . . .p , i = 0, • • • , M from Z.This ensures that the j-th BP determining θ j only depends on the first j-1 components of Z (see eq. 2).It is known that bijective triangular maps with sufficient flexibility can map a simple p-dimensional distribution into arbitrary complex p-dimensional target distributions [33].We use a masked autoregressive flow (MAF) [23] to map Z to the BP coefficients The MAF architecture ensures that that ϑ j i depend only on those components of the latent variables z j ′ with j ′ ≤ j (as required in eq.2).Note that the first coefficients in all BPs ϑ 1 do not depend on z and are therefore not modeled via the MAF.Therefore, the Jacobian ∇f BP w.r.t.z is a triangular matrix, and hence det ∇f BP is given by the product of the diagonal elements of the Jacobian allowing for efficient computation of the resulting p-dimensional variational distribution q λ (θ) via the multivariate version of the change of variable formula (see eq. 6).The flexibility of such a p-dimensional bijective Bernstein flow is only limited by the order M of the Bernstein polynomial and the complexity of the MAF.In our experiments, we use an MAF with two hidden layers, each with 10 neurons.The weights w of the MAF are part of the variational parameters for f BP .In total, we have λ = (ϑ 1 , w, α, β) variational parameters.

C. Variational Inference Procedure
In VI the variational parameters λ are tuned such that the resulting variational distribution q λ (θ) is as close to the posterior p(θ | D) as possible.Here, we do this by minimizing the KL-divergence between the variational distribution and the (unknown) posterior: The KL-divergence is commonly used in VI, and a recent study showed that it is easier to train than other divergences and applicable to higher-dimensional distributions [31].Instead of minimizing (3) usually only the evidence lower bound (ELBO) is maximized [34] which consists of the expected value of the log-likelihood, E θ∼q λ (log(p(D | θ))), minus the KL-divergence between the variational distribution q λ (θ) and the prior p(θ).Note that the ELBO does not explicitly contain the unknown posterior.In practice, we minimize the negative ELBO using stochastic gradient descent facilitated by automatic differentiation.For consistency with [31], we use TensorFlow's RMSprop in all our experiments, configured with the default settings.We follow the BBVI approach and approximate the expected log-likelihood by averaging over S samples θ s ∼ q λ (θ) via Hereby, we also assume the usual independence of the i = 1, . . .N training data points D i .To get these samples θ s , we use S samples z ′ s from the latent distribution, apply the transformation f = f BP • σ • l, and then compute the corresponding parameter samples via θ s = f (z s ).We use the same samples θ s ∼ q λ (θ) to approximate the Kullback-Leibler divergence between the variational distribution q λ (θ) and the prior p(θ) via: where the probability density q λ (θ s ) can be calculated, from the samples θ s using the change of variable function as:

D. Evaluation
Evaluating the quality of the fitted variational distributions requires a comparison to the true posterior.In the case of lowdimensional problems, the two distributions can be compared visually.In the case of higher-dimensional problems, this is not possible anymore.While the Evidence Lower Bound (ELBO) is a valuable metric for optimizing the parameters in VI, it is less helpful in comparing different approximations because it depends on the specific parametrization of the model [35].Therefore, [35] introduced k as a more suited approach for comparison, which since then has been used in other studies like [31] to which we compare.The computation of k is based on the importance ratios which are defined as If the variational distribution q λ (θ) would be a perfect approximation of the posterior p(θ | D) ∝ p(D | θ)p(θ), then important ratios r s would be constant.However, because of the asymmetry of the KL-divergence used in the optimization objective (see Eq. 3), the fitted q λ (θ) tends to have lighter tails than p(θ | D), with the effect that the distribution of r s are heavily right-tailed.To quantify the severity of the underestimated tails, a generalized Pareto distribution is fitted to the right tail of the r s .The estimated shape parameter k of the Pareto distribution can be used as a diagnostic tool.A large k indicates a pronounced tail in the r s distribution and, hence, a bad posterior approximation.According to Yao et al., [35] values of k < 0.5 indicates that the variational approximation q λ is good.Values of 0.5 < k < 0.7 indicate the variational approximation q λ is not perfect but still useful.

III. EXPERIMENTS
We performed several experiments to benchmark our BF-VI approach versus exact Bayesian solutions, Gaussian-VI, and recent NF-VI approaches.All experiments were conducted using five repetitions, with the observed stability of Evidence Lower Bound (ELBO) optimization generally appearing independent of the randomly chosen starting values.Table I shows an overview of the fitted models.The complete model definitions in Stan, along with the code for all experiments, can be found on GitHub.

A. Models with a Single Parameter
First, we demonstrate with two single-parameter experiments that BF-VI can accurately approximate a skewed or bimodal posterior, which is impossible with Gaussian-VI.To obtain complex posterior shapes, we work with small data sets.
a) Bernoulli Experiment: We first look at an unconditional Bayesian model for a random variable Y following a Bernoulli distribution Y ∼ Ber(π) which we fit based on data D consisting only of two samples (y 1 = 1, y 2 = 1).In this simple Bernoulli model, it is possible to determine the solution for the posterior analytically when using a beta distribution as prior.We choose p(π) = Be (1.1,1.1)(π) which leads to the conjugated posterior p(π | D) = Be (α+ yi,β+n− yi) (π) (see analytical posterior in Fig. 2).We now use BF-VI to approximate the posterior.To ensure that the modeled variational distribution for π is restricted to the support of π ∈ [0, 1], we pipe the result of the flow through an additional sigmoid transformation.Fig. 2 shows the achieved variational distributions after minimizing the negative ELBO and demonstrates the robustness of BF-VI: When increasing the order M of the BP, the resulting variational distribution gets closer to the posterior up to a certain value of M and then does not deteriorate when further  increasing M .The left part of Fig. 2 indicates a convergence order of M , which can also be proven for the one-dimensional case (see Appendix ??).As expected, the Gaussian-VI does not have enough flexibility to approximate the analytical posterior nicely (see Fig. 2).
b) Cauchy Experiment: Here, we follow an example from [36] and fit an unconditional Cauchy model Y ∼ Cauchy(ξ, γ) to six samples which we have drawn from a mixture of two Cauchy distributions Y ∼ Cauchy(ξ 1 = −2.5, γ = 0.5) + Cauchy(ξ 2 = 2.5, γ = 0.5).Due to the misspecification of the model, the true posterior of the parameter ξ has a bimodal shape which we have determined via MCMC (see Fig. 3).We use BF-VI and Gaussian-VI to approximate the posterior of the Cauchy parameter ξ by a variational distribution.As in the Bernoulli experiment, also here, BF-VI has enough flexibility to accurately approximate the complex shape of the posterior when M is chosen large enough.Further increasing M does not deteriorate the approximation.Gauss-VI fails as expected.

B. Models with Multiple Parameters
The following experiments use BF-VI in multi-parameter Bayesian models and benchmark the achieved solutions versus MCMC or published state-of-the-art VI approximations (see Table I).In the following experiments, we did not tune the flexibility of our BF-VI approach but allowed it to be relatively high (M = 50) since BF-VI does not suffer from being too flexible.Further, in all experiments in this section, we trained for 10 5 epochs with 5 repetitions and set the number of samples for MC estimation to S = 10 to be comparable with [31].From the repetitions and the posterior samples, we estimated k and the 90 % confidence intervals using the Rubins rule for the BF-VI and Gaussian-VI methods, with S = 50 ′ 000 samples.Though runtime is not a consideration in this study, to provide context, a 2023 MacBook Pro's CPU processes approximately 100 epochs per second.
a) Toy Linear Regression Experiment: To investigate if dependencies between model parameters are correctly captured, we use a simulated toy data set with two predictors and six data points to which we fit a Bayesian linear regression modeling the conditional outcome distribution Fig. 4 gives a visual impression of the joint true posterior of the four model parameters (µ 0 , β 1 , β 2 , σ) determined via MCMC samples (red) and its variational approximation (blue) achieved via BF-VI.The strong correlation between the regression coefficients (β 1 , β 2 ) is nicely captured by the BF-VI approximation.Further, the skewness of posterior marginals involving sigma is similar.However, we can see that BF-VI slightly underestimates the long tails of the posterior, confirming the known shortcoming of using the asymmetric KL-divergence in the objective function [2].BF-VI ( k = 0.68(0.51,0.86)) is superior to MF-Gaussian-VI ( k = 0.90(0.77,1.01)), which can not (by construction) capture the dependencies (MF) or the non-Gaussian shapes (see Fig. S4).
b) Diamond: Linear Regression Experiment: The Diamond linear regression benchmark example (Y | x) ∼ N (µ x = µ 0 + x ⊤ β, σ) has 26 model parameters and 5000 data points (see [37] for reference MCMC samples and Stan code for a complete model definition).Since we have much more data than parameters, the posterior is expected to be a narrow Gaussian around the maximum-likelihood solution, which is indeed seen in the MCMC solution (see Fig. S5).In this setting, BF-VI or NF-VI cannot profit from their ability to fit complex distributions.Still, [31] used this data set for benchmarking different VI methods, e.g., planar NF (PL-NF-VI), non-volume-preserving NF (NVP-NF-VI), and MF-Gaussian-VI.They achieved the best approximation via the simple Gaussian-VI ( k = 1.2).The posterior approximation via PL-NF-VI and NVP-NF-VI have both been unsatisfactory ( k = ∞).We use the same amount of sampling (S = 10) and achieve with BF-VI a better approximation of the posterior k = 5.34(−2.52,13.20) but is still worse than the Gaussian-VI.
A large spread in k indicates an unstable training procedure.This dataset is also quite challenging for MCMC simulations; we did not get satisfactory MCMC samples and took the reference posterior samples from the posteriorDB3 .See Fig. S5 for a comparison of BF-VI and MCMC.c) 8schools: Hierarchical Model Experiment: The 8schools data set is a benchmark data set for fitting a Bayesian hierarchical model and is known to be challenging for VI approaches [35], [38].It has 8 data points, corresponding to eight schools that have conducted independent coaching programs to enhance the SAT (Scholastic Assessment Test) scores of their students.There are two commonly used parameterizations of the model: centered parameterization (CP) and non-centered parameterization (NCP).NCP uses a transformed parameter to facilitate the MCMC sampling.See Table I for more details of these parametrizations and [37] for complete model definitions in Stan.In both parametrizations, the model has 10 parameters.In [31], this benchmark data set was fitted with two NF-based methods and MF-Gaussian-VI.For the CP version [31] kCP = 1.3, 1.1, 0.9 was reported for PL-NF-VI, NVP-NF-VI, MF-Gaussian-VI respectively, which are all outperformed by our BF-VI method with kCP = 0.53(0.11,0.95).For the NCP version [31] kNCP = 1.2, 0.7, 0.7 was reported (same order), and again BF-VI yields a superior kNCP = 0.36(0.17,0.55).A visual inspection of the true MCMC posterior and its variational approximation again shows the underestimated distribution tails (see Fig. S6).For 8Schools and Diamond, a comparison with state-of-the-art models from the literature is summarized in Table II.
d) NN-Based Non-linear Regression Experiment: For this experiment, we use a small Bayesian NN for nonlinear regression fitted on 9 data points.The model for the conditional outcome distribution is (Y | x) ∼ N (µ(x), σ = 0.2).The small size of the used BNN, with only one hidden layer comprising 3 neurons and one neuron in the output  BNN.Because the weights in a BNN with hidden layers are not directly interpretable, they are not of direct interest, and therefore, the fit of a BNN is commonly assessed on the level of the posterior predictive distribution (see Fig. 5).In this example, the more flexible BF-VI shows a slight improvement in approximating the true posterior predictive distribution when compared to the less complex approach with MF-Gaussian-VI, especially inside the regions where there is data (around x = 0 in Fig. 5).

C. Melanoma: Semi-structured NN Experiment
In this experiment, we use BF-VI for semi-structured transformation models [10] (see Fig. 6), where complex data like images can be modeled by deep NNs and tabular data by interpretable model components.Please note that here, both the conditional distribution of the outcome (y | B, x) and the unconditional posterior of the parameters are modeled by transformation models.Because of the deep NN model components involved, MCMC is not feasible anymore to determine the posterior.As a dataset, we use the SIIM-ISIC Melanoma Classification Challenge4 data.The data comes from 33126 patients (6626 as test set, 21200 train, and 5300 validation set) with a confirmed diagnosis of their skin lesions, which is in ≈ 98% benign (y = 0) and in ≈ 2% malignant (y = 1).The provided data D = (B, x) is semi-structured since it comprises (unstructured) image data B from the patient's lesion along with (structured) tabular data x, i.e., the patient's age.
We fit the conditional outcome distribution (Y | D) ∼ Ber π D ) by modeling the probability for a lesion to be malignant π D = p(y = 1 | D) = σ(h) applying the sigmoid function σ(•) to a fitted transformation function h : Y → Z.We study three models for h depending on x alone, B alone, and in combination B and x: M1 (DL-Model) h = µ(B): As a baseline, we use a deep convolutional neural network (CNN) based on the melanoma image data (see Fig. 6 c) with a total of 419,489 weights to take advantage of the predictive power of DL on complex image data.For this DL model, we use deep ensembling [39] by fitting three CNN models with different random initializations and averaging the predicted probabilities.The achieved test predictive performance and its comparison to other models are discussed in the last paragraph of this section.M2 (Logistic Regression) h = µ 0 + β 1 • x: When using only tabular features x, interpretable models can be built.We consider a Bayesian logistic regression with age as the only explanatory variable x and use a BNN without a hidden layer to set up the model (see Fig. 6 b with only one input feature x).In logistic regression, a latent variable is modeled by a linear predictor h = µ 0 + β 1 • x, which determines the probability for a lesion to be malignant via π x = σ µ 0 + β 1 • x allowing to inter et e β1 as the odds-ratio, i.e., the factor by which the odds for lesions to be malignant changes when increasing the predictor x by one unit.In Fig. 7, we compare the exact MCMC posterior of β 1 with the BF-VI approximation, demonstrating that BF-VI accurately approximates the posterior.
M3 (semi-structured) h = µ(B) + β 1 • x: This model integrates image and tabular data and combines the predictive power of M1 with the interpretability of M2.We use a (non-Bayesian) CNN that determines µ(B) and BF-VI for the NN without a hidden layer that determines β 1 (see Fig. The resulting posterior for β 1 differs from the simple logistic regression (see Fig. 7), indicating a diminished effect of age after including the image.Again, e β1 can be interpreted as the factor by which the odds for a lesion to be malignant change when increasing the predictor age by one unit and holding the image constant.
While the main interest of our study is on the posteriors, we also determine the predictive performance on the test set.To quantify and compare the test prediction performances, we look at the achieved log scores (M1: -0.076, M2: -0.085, M3: -0.076) and the AUCs with 95% CI (M1: 0.83(0.79,0.86), M2: 0.66(0.61,0.71), M3: 0.82(0.79,0.85)).For both measures, higher is better.Interestingly, the image-based models (M1, M3) have higher predictive power than M2, which only uses tabular data.The semi-structured model M3, including tabular and image information, has a similar predictive power compared to M1, which only uses images.The benefit of the semi-structured model here is that it provides interpretable parameters for the tabular data along with uncertainty quantification without losing predictive performance.

IV. SUMMARY AND OUTLOOK
The proposed BF-VI is flexible enough to approximate any posterior in principle without being restricted to variational distributions with known parametric distribution families like Gaussians.In benchmark experiments, BF-VI accurately fits non-trivial posteriors in low-dimensional problems in a BBVI setting.For higher dimensional models, BF-VI outperforms published results from other NF-VI methods [31] on the studied benchmark data sets.Still, we observe that the posterior can not be fitted perfectly in high dimensions by BF-VI, especially since the tails of the approximation are too short.We attribute this limitation to known difficulties in the optimization process and the asymmetry of the KL divergence.These challenges of VI were not in the focus of our study, and we leave it to further research.
To the best of our knowledge, we are the first to demonstrate how BBVI can be used in semi-structured models.We used BF-VI on the public melanoma challenge data set, integrating image data and tabular data by combining a deep CNN and an interpretable model part.We see a valuable application of BF-VI in models with interpretable parameters, i.e., statistical or semi-structured models where we can model complex posterior distributions of the interpretable parameters.Especially in semistructured models with deep NN components that cannot be fitted with MCMC, BF-VI allows determining the variational distribution for the interpretable model parts.Moreover, efficient SGD optimizers can be used in BF-VI to fit all model parts jointly.We plan to extend our research on BF-VI for semi-structured models in the future and investigate the quality of the posterior approximations.
This shows that the KL-Divergence in (3) asymptotically converges to zero as 1/M .
In the following, we give additional details for the experiments.For the examples using MCMC simulations as a benchmark, the Stan-code can be found at https://github.com/tensorchiefs/bfvipaper/tree/main/mcmc.If not otherwise stated, the training has been as described in the main text.If not given below, the used data can be found at: https://github.com/tensorchiefs/bfvipaper/tree/main/data.

A. Bernoulli Example
The Data Generation Process: Two samples are drawn from a Bernoulli distribution.
The Actual Data: The realizations are y 1 = 1 and y 2 = 2.

Details of the Training Procedure:
To demonstrate the possibility that the described solution converges in principle to the exact solution, we used a large number of MC samples S = 2500 and epochs=2500 for the training.1) Additional Results (Tightness of the ELBO): For this example, we can estimate the KL divergence via The sum is calculated from the samples θ s of the trained approximative distribution, which can be obtained via the trained transformation function.The likelihood p(D | θ s ) and prior p(θ s ) for those samples can be easily calculated and probabilities q λ (θ s ) can be calculated via (6).In this simple case, the posterior can be calculated analytically to p(θ s | D) = Be (α+ yi,β+n− yi) (θ s ) = Be (3.1,1.1)(θ s ).
In Fig. 2, we show a direct comparison of the posterior and its approximation of the dependence of the tightness of the ELBO on M , importantly increasing M and thus the flexibility does not deteriorate the approximation.
Details of the Model: We use an unconditional Cauchy model y ∼ Cauchy(ξ, γ = 0.5) were γ = 0.5 is given.Note that this model is misspecified and gives rise to a multimodal posterior.The location parameter ξ is assigned a standard normal prior ξ ∼ N (0, 1).

Details of the Training Procedure:
To demonstrate the possibility that the described solution converges in principle to the solution, we used a rather large number of MC samples S = 1000 and epochs=1000 for the training.p(D|θs)p(θs) dθ analytically anymore.However, we can rewrite (A1) to: Numerical integration obtains the log-evidence log(P (D)) = −21.43069for this example.In Fig. S3, we show a direct comparison of posteriors.The dependence of the tightness of the ELBO on M calculated in the lower part is calculated via (S4).
2) Additional Results: Effect of Pre-Processing: In Fig. S3, the effect of different transformations before the Bernstein polynomials is studied.
There are different options to ensure that the input to BP fulfills the restriction z ∈ [0, 1].One option is to use a latent distribution with unrestricted.A random variable z ′ from that distribution is then "piped" through a sigmoid function σ : [−∞, ∞] → [0, 1].Note that the total transformation f BP • σ stays monotonically increasing.This construction allows us to choose the Standard-Normal for F Z ′ .We have observed in low-dimensional examples a more efficient training and better fitting with fewer parameters M when additionally prepending an affine transformation function l = α • z + β before σ so that the total transformation is f = f BP • σ • l (see section B).We also experimented with an additional affine function after the Bernstein polynomials but did not find a beneficial effect.We use f = f BP • σ • l in the following, except in B where we studied other possibilities to restrict the input to f BP to [0, 1].Moreover, a distribution F Z with support [0, 1] can be used.We experimented with a truncated Gaussian.
While a for large number M of the performance for the different methods become similar (see Fig. S3), the transformation f 1 = f BP • σ • l works better for smaller of M , we thus choose this transformation for all experiments.

C. Toy Linear Regression Example
In this example, a posterior is generated, which displays a strong correlation between the results.

D. Diamond
The diamond data set is modeled with a linear model.To be in accordance with [37], we choose standard normal priors, except for σ and the intercept, for which we choose Student's t-distributions as priors.Specifically, the intercept is assigned a Student's t-distribution with 3 degrees of freedom, a location parameter of 8, and a scale parameter of 10, while σ is modeled using a half-Student's t-distribution with 3 degrees of freedom, a location of 0, and a scale parameter of 10.

E. 8School
The 8-Schools example comes in NCP and CP parameterizations.See Table I for the model definition, including prios. 1) Additional Results: In Fig. S6, we show samples from the marginals of the 8-dimensional parameter θ (right side, centered parameterization) and θ (left side, non-centered) parameterization.For the MCMC samples, the easier centered parameterization yielding θ is used for both cases, and θ is calculated from it.We see an underestimation of the true posterior.

Fig. 1 .
Fig.1.Overview of the transformation model.A: shows the bijective transformation function g : Z → θ (or its approximation f BP ) mapping between B: a predefined latent density f Z and C: a potentially complex posterior (or its variational distribution).

Fig. 2 .
Fig. 2. Bernoulli experiment.Left panel: Comparison of the analytical posterior for the parameter π in the Bernoulli model Y ∼ Ber(π) with variational distributions achieved via Gaussian-VI and BF-VI with BP order M = 1, 10, 50.Right panel: The dependence of the divergence KL(q λ (w) ∥ p(w | D)) on M for 20 runs.

Fig. 4 .
Fig. 4. Toy linear regression example visualization of the posterior.The model has four parameters: two regression coefficients β 1 and β 2 , the intercept µ 0 , and the standard derivation σ.Samples from the true posterior resulting from MCMC (red) are overlayed with samples from the BF-VI approximation (blue).

Fig. 6 .
Fig. 6.The architecture of the used NN models or model parts.a) Dense NN with one hidden layer to model non-linear dependencies from the input (used in the NN-based non-linear regression example).b) Dense NN without hidden layer to model linear dependencies from tabular input data (used in M1 and M3 of the Melanoma experiment).c) CNN to model non-linear dependencies from the image input (used in M1 and M3).

Fig. S3 .
Fig.S3.Effect of the different transformations before the f BP on the performance.In f 1 , f 2 z ∼ N (0, 1) and the transformation aref 1 = f BP • σ • l and f 2 = f BP • σ.For f 3 ,z is sampled from a truncated normal with shape 0.5, scale 0.15, and with bound 10 −4 to 1 − 10 −4 and transformed via f BP .

0
Fig. S4.Mean-Field Gaussian results for the linear regression posteriors.

1 ) 15 MCMCFig. S5 .
Fig. S5.Comparison of samples from the marginals of MCMC samples and samples from the approximative posterior for the diamond data set.Here, the MCMC solution is very narrow, and the BF-VI (as other approximations) fails to converge to the narrow posterior.The red crosses indicate the maximum likelihood solution.

F
Fig. S6.Marginal samples from MCMC and BF-VI for the 8 schools examples.For the NCP-parametrization θi (i = 1, . . ., 8) is shown on left two plots and θ right two plots (CP-Parametrization).

Fig
Fig. S7.Loss curve for the different experiments.The loss is scaled by subtracting the minimum loss and dividing by the absolute value of the minimum loss.

TABLE I OVERVIEW
OF THE FITTED BAYESIAN MODELS IN THE BENCHMARK EXPERIMENTS AND THE METHODS USED TO GET THE POSTERIORS.MCMC OR ANALYTICAL SOLUTIONS ARE USED AS GROUND TRUTH, AGAINST WHICH THE QUALITY OF DIFFERENT VI APPROXIMATIONS IS COMPARED (BOLD-FACED METHODS OUTPERFORM ALL OTHERS).THE DIFFERENT VI APPROXIMATIONS ARE EITHER COMPARED BY THE ACHIEVED k (SEE SECTION II-D) OR VISUALLY BY INSPECTING THEIR DEVIATION FROM THE GROUND TRUTH.

TABLE II COMPARISON
[31]ERIOR APPROXIMATION WITH RESULTS FROM THE LITERATURE FOR THE DIAMOND AND 8SCHOOLS DATASET IN CP AND NCPPARAMETRIZATION.SHOWN IS k (THE LOWER, THE BETTER) FOR DIFFERENT APPROXIMATIONS: MEAN-FIELD GAUSSIAN (MF-G) AND STUDENT-T (MF-T); NVP (NVP-VI) AND PLANAR (P-VI) FLOW FROM[31], AND BF-VI.