Shrinkage prior for implicit generative models

Editor:

Introduction

Implicit generative models are increasingly used to model complex random objects arising in scientific and real-world applications, including weather forecasting [46], population genetics [22] (Yelmen and Jay 2023), protein structure and molecular design [1], and image generation [18, 20]. Their ability to generate realistic samples without requiring an explicit likelihood has made them a powerful modelling paradigm. However, the lack of a tractable likelihood also creates significant challenges for statistical inference and model assessment.

In this work, we consider statistical inference for such models. Implicit generative models introduce a latent random variable $\mathbf{z} \in\mathbb{R}^{k} \sim q(\mathbf{z})$, where $q$ is a fixed, tractable distribution on $\mathbb{R}^{k}$, together with a deterministic nonlinear map $G_{\theta}: \mathbb{R}^{k} \to\mathbb{R}^{d}$, parameterized by $\theta\in\mathbb{R}^{p}$. The model generates samples through

with distribution

the pushforward distribution of $q$. While traditional approaches often focus on scientific models that incorporate prior knowledge of the underlying process, we consider high-dimensional implicit models where the data-generating mechanism is unknown. In particular, we focus on deep generative models, where $G_{\theta}$ is represented by a neural network (Mohamed and Lakshminarayanan [36]). For conditional tasks, where the objective is to model random objects given an observed variable $\mathbf{w} \in\mathbb{R}^{m}$, we instead consider the map $G_{\theta} : \mathbb{R}^{k} \times\mathbb{R}^{m} \to\mathbb{R}^{d}$, which generates

with conditional distribution $P_{\theta}(\cdot\mid\mathbf{w})$, as in conditional GANs (Mirza and Osindero [35]).

When $G_{\theta}$ is represented by a large neural network and trained on an observed dataset $\mathcal{D}_{n}=\{y_{1},\ldots,y_{n}\}$ with limited sample size, statistical inference becomes a large-$p$, small-$n$ problem ($p \gg n$). Bayesian methods based on shrinkage priors are well established for this setting in linear regression (Bhadra et al. [3]). However, these methods cannot be applied directly to implicit generative models because the probability density of $P_{\theta}$ is unavailable. Building on the generalized Bayesian framework of Pacchiardi et al. [43], which replaces the likelihood with a scoring rule, we develop a sparsity-inducing generalized Bayesian approach using the horseshoe prior (Carvalho et al. [6]). We further introduce a blocked Gibbs sampler based on preconditioned adjusted stochastic gradient Langevin dynamics (prec-adjSGLD) to efficiently sample from the resulting high-dimensional posterior. To the best of our knowledge, while Bayesian regularization has been extensively studied for deterministic neural networks ([13, 14, 37]), this is the first work to develop a sparsity-inducing Bayesian framework for neural implicit generative models.

We first study the large-$p$, small-$n$ sparse linear regression problem. Although standard Bayesian methods are available, Markov chain Monte Carlo algorithms are often prohibitively expensive in this setting (Johndrow et al. [24]). Our proposed method achieves strong sparsity recovery while providing the best predictive performance at a substantially lower computational cost. We also consider a sample-efficient training of a conditional implicit generative model for the Hopper, a reinforcement learning benchmark (from the OpenAI Gymnasium package (Kwiatkowski et al. [28])), which mimics a Markov decision process.

The rest of the chapter is structured as follows. § 2 defines a generalized posterior over the parameters of an implicit generative model using the adaptive shrinkage horseshoe prior. We then introduce an approximate blocked Gibbs sampling scheme to sample from the proposed posterior. § 3 presents simulation studies demonstrating the application and advantages of the proposed sampler over fixed priors. We consider two examples: high-dimensional sparse linear regression, and a conditional deep generative model for the reinforcement learning environment Hopper. Finally, § 4 concludes the chapter.

Generalized Bayesian regularization for implicit generative neural models

Due to the intractability of the likelihood functions for the implicit generative neural models, they are trained by matching the pushforward distribution $P_\theta$ to the data distribution using adversarial objectives [17], $f$-divergences [39], proper scoring rules [19] (Pacchiardi and Dutta 2022) or using prequential scoring rules [42] for conditional generative models. Despite their empirical success, such training procedures are data-intensive and provide limited quantification of uncertainty. This is specially problematic in over-parameterized settings or when data are scarce, well known from many studies for $p \gg n$ sparse linear regression [11].

To account for parameter uncertainty in this setup, we propose to adopt a Bayesian approach. Let $\mathcal{G} = \{G_\theta: \theta\in\Theta\subseteq\mathbb{R}^p\}$ denote the model class and place a prior distribution $\pi(\theta)$ over parameters. Given observed data $\mathcal{D}_n$ and a likelihood function $l(\theta; \mathcal{D}_n)$ under the model $G_\theta$, Bayesian inference yields the posterior distribution

In the absence of a tractable likelihood $l(\theta; \mathcal{D}_n)$, we follow Pacchiardi et al. [42] and replace the likelihood with the exponential of the negative of a strictly proper scoring rule, yielding a scoring-rule posterior. To further regularize the model, we place a sparsity-inducing horseshoe prior [6, 14] on $\theta$. We discuss the details in the following sections.

Generalized Bayesian likelihood-free Inference

Likelihood-free Bayesian methods, such as approximate Bayesian computation (ABC) [2] and neural network–based surrogate approaches [31] often consider low-dimensional models with low dimensional parameters. However, they struggle to scale to the high-dimensional parameter spaces typical of deep generative models. Motivated by generalized Bayesian framework with loss functions [5, 23, 27], [43], we propose to define scoring rule posterior by choosing proper scoring rules [9] as the loss function. These scoring rules can be unbiasedly estimated from the easily available simulations of the implicit generative models, as shown in the previous chapters.

Scoring rules [16] were originally proposed as a diagnostic of a probabilistic forecast, which provides a penalty $(S(P,Y))$ to the forecast $(P)$ evaluated at an observed outcome $(Y)$. If $Y \sim Q$ denotes the true data-generating distribution, the expected score is defined as $S(P,Q) := \mathbb{E}_{Y\sim Q}[S(P,Y)]$. A scoring rule is proper with respect to a class of distributions $\mathcal{P}$ if

which ensures that truthful forecasting is optimal. A scoring rule is said to be strictly proper if $S(P,Q)$ is uniquely minimized when $P = Q$, the true data-generating distribution.

Strictly proper scoring rules induce statistical divergences between probability measures [9].

Given data $\mathcal{D}_n = \{y_1,\ldots,y_n\}$, an implicit generative model $G_\theta$, a prior $\pi(\theta)$, and a weight parameter $\omega> 0$, the scoring-rule posterior [43] becomes

where $P_\theta$ corresponds to the pushforward distribution induced by $G_\theta$ and, $\omega> 0$ is called the learning rate which determines the influence of the loss with respect to the prior. $\omega$ can be selected to achieve different goals (often, to match some frequentist property, see Lyddon et al. [32]; Syring and Martin (2019); Matsubara et al. [34], or to ensure two different scoring-rule posteriors have same scale [43]). However, we keep $\omega= 1$ fixed for our simulation studies.

When a proper scoring rule is used to define the posterior in Equation (1), the resulting distribution satisfies a Bernstein–von Mises–type result under mild regularity conditions [41], which we state here.

In particular, the scoring-rule posterior concentrates almost surely on the parameter value that minimizes the expected scoring rule, provided the minimizer is unique. In well-specified models with strictly proper scoring rules, this minimizer coincides with the true parameter value, yielding posterior consistency. Uniqueness of the minimizer, however, may also hold under model misspecification or for non-strictly proper scoring rules.

To facilitate computation in the likelihood-free scenario, we should choose a strictly proper scoring rule for which we can easily compute an unbiased estimate. From many alternative scoring rules [16], here we choose energy score [47] defined as

where $X \perp X' \sim P_\theta$ and $\beta\in(0,2)$. In this work, we will keep $\beta$ fixed at 1. It is strictly proper for distributions satisfying $\mathbb{E}_{X\sim P_\theta}\lVert X\rVert^\beta< \infty$. Given samples $x_m=\{x_j\}_{j=1}^{m}$ with $x_j \sim P_\theta$, the energy score can be unbiasedly estimated as,

Sparsity inducing horseshoe prior

Most works applying Bayesian inference for deterministic neural networks [12], adopt simple isotropic Gaussian priors or heavy-tailed t-distribution. However, the “Bad Prior Hypothesis” (Wenzel et al., 2020) highlights the limitations of such priors in over-parameterized models, which can be especially hazardous in data-scarce regimes. Gaussian priors implicitly correspond to $\ell_2$-type regularization [3], which may be insufficient to control model complexity in high-dimensional settings.

To address this, shrinkage priors have been proposed as alternatives. One class consists of discrete mixture priors, such as spike-and-slab [25, 7], which combine a point mass at zero with a diffuse continuous component. However, spike-and-slab priors pose substantial computational challenges, as posterior inference involves a combinatorial model space with discrete inclusion indicators and prevents efficient block updates of model parameters. An alternative class comprises continuous global-local shrinkage priors [45], which adaptively shrink small signals while preserving large ones. Among these, the horseshoe prior [6] has emerged as a popular and effective choice, enabling efficient block sampling strategies. In the context of Bayesian neural networks, Ghosh et al. [14] demonstrated that horseshoe priors can induce sparsity at the level of network nodes, enabling automatic model selection and improved interpretability. Motivated by these results, we employ horseshoe priors over the parameters of generative neural networks to encourage sparse representation.

Let $\theta$ denote the vector of all weights and biases in the network. The horseshoe prior is specified hierarchically as

where $\theta=(\theta_1,\ldots,\theta_p)$. Here, $\lambda_j$ are parameter-specific local shrinkage parameters for all $j=1,\ldots,p$; and $\tau$ is a global shrinkage parameter that controls the overall degree of shrinkage. $C^{+}$ denotes the half-Cauchy distribution, and $b_0$ and $b_g$ are hyperparameters governing the prior distributions of the local and global shrinkage parameters, respectively.

Blocked Gibbs sampler

Combining (1) and (2), we can define the generalized posterior over the hierarchical model parameters as following,

where $\lambda_j \sim C^{+}(0,b_0)$ $\forall j = 1,\ldots,p$ and $\tau\sim C^{+}(0,b_g)$. Here, exact block sampling would require conjugacy between prior and likelihood, which breaks down for parameters $\theta$ as we have replaced log-likelihood with proper scoring rules. The posterior on the hyperparameters, however, still remains conjugate. Thus, we propose an approximate sampling scheme based on stochastic gradient Monte Carlo Markov chain (SGMCMC) [38] for $\theta$ conditional on hyperparameters and an exact conditional update for the hyperparameters conditional on $\theta$. This is an approximate blocked SGMCMC within the Gibbs algorithm.

When $\theta= (\theta_j)_j$ are the network parameters and $\alpha= ((\lambda_j^2)_j,\tau^2)$ are the hyperparameters, to perform the Gibbs sampling we update $\theta$ and $\alpha$ in an iterative manner using the following update equations.

Block 1: Update parameter $\theta$ conditional on hyperparameters $\alpha$ and data $\mathcal{D}_n$

We update $\theta$ by using a variant of stochastic gradient Monte Carlo Markov chain algorithm (SGMCMC), specifically preconditioned adjusted stochastic gradient Langevin dynamics (prec-adjSGLD) described in § 2.4

Block 2: Update hyperparameters conditional on $\theta$ Half-Cauchy priors on $\lambda_j$ and $\tau$ break conjugacy. However, a scale-mixture representation of the Half-Cauchy prior proposed by [33], can be used to arrive at the following update equations using the auxiliary variables $\nu_0,\nu_1,\ldots,\nu_p \sim\mathcal{IG}(1/2,1)$,

for all $j = 1,2,\ldots,p$, where $\mathcal{IG}(a,b)$ is the inverse-gamma distribution with shape parameter $a > 0$ and scale parameter $b > 0$. Hence, we can update hyperparameters $\alpha$ exactly, using the above update equations.

Prec-adjSGLD to update Block 1

We propose a variant of the SGMCMC algorithm to update the model parameters in Block 1. SGMCMC methods are based on continuous-time diffusion processes using gradients of the negative log target distribution (or the potential function $U(\theta)$); whose stationary distribution matches the target distribution. In our case, $U(\theta) = \sum_{i=1}^{n} S_E(P_\theta, y_i) + \sum_{j=1}^{p} \frac{1}{2\lambda_j^2\tau_j^2}\theta_j^2$, for which the gradient is not analytically tractable.

But we can compute an unbiased estimate of the gradient, by first noticing that $S_E(P_\theta,y) = \mathbb{E}_{X,X'\sim P_\theta} f(X,X',y)$ for a smooth differentiable function $f$. Thus, the scoring rule is defined through an expectation over (possibly multiple) samples from $P_\theta$. Hence the gradient of the scoring rule reduces to,

where $q$ stands for an auxiliary distribution independent of $\theta$. In the latter equality, the exchange between expectation and gradient is not a trivial step. Luckily, Theorem 5 in (Bińkowski et al., [4] proved the above step to be valid almost surely with respect to a measure on $\theta$, under mild conditions on $G_\theta$ (such conditions are satisfied if the function is differentiable, which holds for all the maps defined by neural networks considered in this work). Based on this, we can unbiasedly estimate the gradient of the scoring rule as follow:

where $Z_i \perp\!\!\!\perp Z'_j \sim q$. In practice, this can be easily obtained by implementing the function $G_\theta$ using automatic-differentiation libraries ([44]).

The bias due to this approximate stochastic gradient of $U(\theta)$ in SGMCMC algorithms can be mitigated by incorporating momentum ($p$) and friction ($\xi$) terms, resulting in the adjusted stochastic gradient Langevin dynamics (adjSGLD) ([29]). adjSGLD operates in an extended state space consisting of the model parameters, and the auxiliary momentum and friction variables. This has been successfully applied in Bayesian inference ([10]) as well.

To improve the efficiency of SGMCMC in the high-dimensional parameter space of $\theta$, we can use preconditioning strategies ([15, 8]) using the squared gradient of $U(\theta)$, similar to adaptive optimizers like RMSprop ((Tieleman, 2012)) and Adam ([26]). Here we use a preconditioned version of the adjSGLD scheme (prec-adjSGLD) to sample $\theta$, for which the stochastic differential equation (SDE) defining the diffusion becomes,

where $W$ denotes a standard Brownian motion and $H(\theta)$ encodes geometric information of the potential energy $U(\theta)$, called the Riemannian metric ([15]), which is commonly defined by the Fisher information matrix. Theoretical properties and convergence of adSGLD and prec-adjSGLD algorithm have been studied in [29, 30] and [10] correspondingly.

The other source of error in SGMCMC arises due to the time discretization of the above SDE. For example, when the first order Euler–Maruyama discretization scheme is used. To further reduce these numerical errors in SDE integration, we use a symmetric splitting scheme (SSS) [8], described in Appendix A. This 2nd-order integrator provides higher accuracy with minimal extra computation, and we arrive at Algorithm 1 presented below.

Algorithm 1: Preconditioned adjusted stochastic gradient Langevin dynamics

Input: Learning rate $\eta$, parameters $(\lambda_0 = 0.99,\lambda_1 = 10^{-6})$, noise $\{\zeta_t \sim\mathcal{N}(0,I_p)\}$

Initialize: $\theta_0$, $u_0 = \sqrt{\eta}\cdot\mathcal{N}(0,I_p)$, $a_0 = \sqrt{\eta}$, $v_0 = 0$

for $t = 1,2,\ldots$ do

Compute stochastic gradient: $\tilde{f}_t \leftarrow\nabla_\theta\widetilde{U}(\theta_{t-1})$

$v_t \leftarrow\lambda_0v_{t-1} + \frac{1-\lambda_0}{\eta^2}\tilde{f}_t\odot\tilde{f}_t$

$g_t \leftarrow1\oslash\sqrt{\lambda_1+v_t}$

$\theta_t \leftarrow\theta_{t-1}+\frac{1}{2}g_t\odot u_{t-1}$

$a_t \leftarrow a_{t-1}+\frac{1}{2}(u_{t-1}\odot u_{t-1}-\eta)$

$u_t \leftarrow\exp(-a_t/2)\odot u_{t-1}$

$u_t \leftarrow u_t-\eta\cdot g_t\odot\tilde{f}_t+\sqrt{2\cdot g_{t-1}\eta^{3/2}}\odot\zeta_t$

$u_t \leftarrow\exp(-a_t/2)\odot u_t$

$a_t \leftarrow a_t+\frac{1}{2}(u_t\odot u_t-\eta)$

$\theta_t \leftarrow\theta_t+\frac{1}{2}g_t\odot u_t$

end

Simulation studies

To evaluate our method, we first consider a sparse linear regression task with a small number of observations and varying numbers of regressors. In this setting, we assume Gaussian noise, which allows access to the exact likelihood. We compare the performance of the (likelihood-based) exact Gibbs sampler with our proposed (likelihood-free) approximate sampler under the horseshoe prior. Finally, we consider a Markov decision process, called Hopper. Here, the simulator model is a conditional generative model and we analyze the effect of different priors on the model parameters compared to the horseshoe prior.

Bayesian sparse high-dimensional linear regression

Let us consider a high-dimensional linear regression problem,

where $\boldsymbol{\epsilon}\sim\mathcal{N}(0,\sigma^2 I_n)$, $I_n$ is $n\times n$ identity matrix and $p\gg n$. Inference for $\theta$ under sparsity has been widely addressed using penalized estimation methods. Classical examples include ridge regression [21], LASSO (Tibshirani, 1996), and elastic net (Zou and Hastie, 2005), among many other regularization approaches. Similarly, Bayesian inference for $\theta$ can be performed by using a shrinkage-inducing prior, like a horseshoe prior.

Although the true likelihood function is available in this case, the traditionally used Gibbs sampler is prohibitively expensive, and approximations to sample from the posterior is an active research area [3, 24].

Considering sparse linear regression as an implicit generative model, we demonstrate inference on $\theta$ by sampling from a scoring-rule posterior using our proposed methodology. It is scalable and achieves superior predictive performance. In our experiments, we fix the sample size at $n = 50$ and consider three settings with $p = 500$, $1000$, and $3000$.

Our experimental setup follows that of Johndrow et al. [24]. We set $\sigma^2 = 4$, and generate $W$ from a standard Gaussian distribution. We define the true sparse parameter vector $\theta$ as

We compare the proposed sampler with an exact Gibbs sampler, which is feasible in this setting since the likelihood is available in closed form. For both methods, we use $b_0 = b_g = 1$ to define the horseshoe prior. The sensitivity to these hyperparameters is investigated in Appendix B, where we also consider $(b_0, b_g) = (10, 10)$ and $(50, 50)$. We set the learning rate to $\eta= 10^{-8}$ in Algorithm 1 and we use $m = 10$ simulations to estimate the gradient in Equation (6) for the approximate sampler. Both the samplers are run under matched time budgets. In Figure 1, we report the log likelihood of the posterior samples.

Figure 1 illustrates that our method exhibits improved sample efficiency as the dimensionality $p$ increases, while the exact Gibbs sampler struggles in higher-dimensional settings. All experiments are repeated across ten independently generated datasets, and the results are averaged.

For comparison, we also implement ridge regression, LASSO, and elastic net for the linear regression model. For these methods, we evaluate the log-likelihood at the estimated parameter. For the Bayesian approach with the horseshoe prior, we compute the mean log-likelihood over posterior samples obtained from both the exact and the proposed approximate samplers. To ensure a fair comparison, we discard the initial samples corresponding to the first 25, 75, and 100 seconds of runtime for $p = 500$, $1000$, and $3000$, respectively.

The mean log-likelihood is then computed using the remaining samples. Table 1 reports the results averaged over the 10 simulated datasets.

Method	$p = 500$	$p = 1000$	$p = 3000$
Ridge	$-30.79$ (45.52)	$-36.62$ (65.32)	$-127.20$ (165.12)
LASSO	$-34.70$ (36.70)	$-98.18$ (80.73)	$-111.54$ (116.90)
Elastic Net	$-34.34$ (37.00)	$-88.61$ (84.56)	$-92.38$ (124.10)
Horseshoe (exact Gibbs)	$-1.616$ (0.001)	$-1.830$ (0.032)	$-23.290$ (1.169)
Horseshoe (approximate, ours)	$-2.805$ (1.486)	$-1.626$ (0.035)	$-1.612$ ($2 \times10^{-5}$)

The Bayesian methods consistently outperform the three frequentist regularization methods in this example. Moreover, the approximate sampler performs better in larger dimensions, consistent with the trends observed in Figure 1.

Next, we report the runtime comparison of each iteration of the exact and approximate samplers in Table 2. The runtime of the exact sampler blows up with dimension, whereas the proposed approximate sampler does not have that trend.

	$p = 500$	$p = 1000$	$p = 3000$
Exact Sampler	0.0987 (0.0019)	0.6736 (0.0218)	21.9139 (1.8124)
Approximate Sampler	0.0309 (0.0025)	0.0310 (0.0015)	0.0345 (0.0005)

Neural generative Markov processes: Hopper

Next, we consider a model-based reinforcement learning (RL) task. In model-based RL, the environment dynamics are first learned from interactions with the environment. The learned model is then used to train a sample-efficient policy. Therefore, improving model learning with fewer interactions can directly benefit the final policy performance. This motivates the use of shrinkage priors, which are known to perform well in data-scarce settings.

We consider the ‘Hopper’ environment from the OpenAI Gymnasium package [28]. The task is to move a two-dimensional, single-legged structure forward by applying forces to three joints connecting its four body parts. The action space consists of torques applied to the three joints, with values ranging from $-1$ Nm to $1$ Nm. The state is represented by a 12-dimensional real-valued vector containing joint angles, angular velocities, and the positions and velocities of the body parts. The reward function encourages forward movement while penalizing large torques that may destabilize the system.

We generate a trajectory of state-action-next state from the above system using a random policy. We wish to model this Markov process by a deep generative model following the setup of Sinha-Roy et al. (2024). Here, we consider the following conditional generative model,

where $z$ is a Gaussian noise, $a_t$ and $s_t$ are action and state at time $t$ and $f_{\theta}$ is fully connected feedforward neural network with four hidden layers of 100 neurons each and swish activation function, which amounts to a total number of parameters approximately 33,000.

To infer the model parameters $\theta$, we consider a scoring-rule prequential posterior (Sinha Roy et al., 2026) under three different priors: a horseshoe prior with hyperparameters $b_0 = b_g = 1$, an independent standard Gaussian prior, and an independent Student-$t$ prior with 10 degrees of freedom on each network parameter. We focus on the data-scarce regime with $n = 10, 20, 50,$ and $100$ observations, in contrast to the conventional setting where thousands of observations are used for training. For sampling, we set the learning rate to $\eta= 10^{-12}$ in Algorithm 1 and use $m = 10$ simulations to estimate the gradient in Equation (6). Each sampler is run for 10000 iterations, and the last 1000 samples are retained to evaluate predictive performance.

Figure 2 reports the NRMSE, coefficient of determination $R^2$, and average energy score on an independent test set of 1000 observations. We also consider an ensemble model in which the network parameters are initialized from the horseshoe prior and optimized using Adam. The objective function is estimated using $m = 10$ simulations. The horseshoe prior consistently outperforms the Gaussian and Student-$t$ priors, especially when the number of training observations is small. The ensemble model also performs well for small values of $n$, suggesting that the horseshoe prior helps recover the sparse structure of the dynamics model. In contrast, the Gaussian and Student-$t$ priors produce substantially larger prediction errors, indicating poorer calibration when only fewer observations are available.

We next study the sensitivity of the horseshoe hyperparameters for additional cases when $n = 500$ and 1000. We consider five independently generated trajectories and compare the standard Gaussian prior with horseshoe priors using five different choices of $(b_0, b_g)$. Table 3 reports the results for the Gaussian prior, while Table 4 reports the results for the horseshoe prior with $b_0 = b_g = 1$. The remaining results are provided in Appendix C. Overall, the predictive performance is mostly insensitive to the choice of horseshoe hyperparameters when $b_0, b_g \ge1$. Finally, Table 5 reports the average runtime per iteration for the Gaussian and horseshoe priors.

$n$	NRMSE	$R^2$	Energy score
50	16998.45 (3058.44)	-2.31e10 (8.50e9)	78153.29 (16198.17)
100	17301.19 (2740.90)	-2.58e10 (6.68e9)	80490.82 (13381.43)
500	16442.49 (3102.66)	-2.74e10 (1.05e10)	63382.63 (9684.82)
1000	18561.19 (2940.66)	-3.27e10 (1.11e10)	70151.97 (12845.99)

$n$	NRMSE	$R^2$	Energy score
50	0.212 (0.070)	-3.070 (3.260)	3.001 (0.079)
100	0.215 (0.074)	-3.120 (3.750)	2.984 (0.065)
500	0.152 (0.073)	-1.550 (3.200)	2.934 (0.046)
1000	0.115 (0.055)	-0.470 (1.890)	2.914 (0.031)

Method	$n = 50$	$n = 100$	$n = 500$	$n = 1000$
Gaussian prior	0.062 (0.002)	0.105 (0.006)	0.487 (0.022)	1.077 (0.006)
Horseshoe prior	0.118 (0.016)	0.153 (0.012)	0.517 (0.028)	1.208 (0.050)

Table 5 reports the runtime of the approximate sampler using the Gaussian and horseshoe priors. The horseshoe prior requires slightly more computational time due to the additional componentwise updates of the shrinkage parameters. However, this increase in runtime is small compared to the substantial predictive improvements obtained. These results highlight the advantage of using adaptive shrinkage priors, such as the horseshoe prior, in severely data-poor settings.

Discussion and conclusion

In recent years, neural implicit generative models have been rapidly adopted across many application domains. This work extends their applicability to settings with limited training data, which may arise due to data collection constraints or high acquisition costs. We achieve this by incorporating shrinkage priors within the scoring-rule posterior framework.

First, we demonstrate the strong predictive performance of the proposed methodology on a high-dimensional sparse linear regression. Next, we implement our method on a neural generative model, motivated by applications in meteorology and reinforcement learning. Beyond enabling scalable training of larger neural architectures with fewer observations, our approach also helps identify sparse structures within generative neural networks due to the shrinkage prior. Further investigation in this direction may lead to more meaningful and potentially interpretable neural network architectures for specific applications.

In this work, we fix the scaling parameter of the scoring-rule posterior, $\omega$. Future work could investigate the use of multivariate scaling parameters or data-driven approaches, such as cross-validation, to select $\omega$ and obtain more informative scoring rules. Another interesting direction is to develop model selection approaches for implicit generative models, or for combinations of models and scoring rules viewed as approximate likelihoods, through comparisons of their marginal scoring-rule posteriors.

Acknowledgment

We thank Sherman Khoo for his initial investigations into the horseshoe prior and for many fruitful discussions.

Splitting scheme for preconditioned adjusted stochastic gradient Langevin dynamics

To generate more accurate samples, we resort to an easily implemented 2nd-order symmetric splitting scheme (SSS). The idea of SSS is to split an infeasible SDE into several sub-SDEs, with each part analytically solvable; approximate samples are then generated by sequentially evolving parameters via these sub-SDEs. Specifically, we split the SDE in Equation (7) into the following three sub-SDEs:

We then use the order ABOBA to solve the above sub-SDEs to generate approximate samples. Specifically, based on the splitting, parameters first evolve following $A$ for time $\eta/2$, then $B$ for time $\eta/2$, $O$ for time $\eta$, $B$ again for time $\eta/2$, and finally $A$ again for time $\eta/2$, resulting in the update equations in Algorithm 1. The resulting parameters then serve as approximate samples from the target distribution.

Sensitivity analysis of prior choices for sparse linear regression

We present the impact of the choice of the hyperparameters $(b_0,b_g)$ on the exact and the approximate sampler for linear regression in Figure 3. The choice of hyperparameters has very little impact on the exact sampler; however, the approximate sampler performs the best when $b_0,b_g = 1$.

Sensitivity analysis of prior choices for Hopper

For the Hopper experiment we provide the ablation study on $b_0,b_g$ in Table 6. We observe similar predictive performance via horseshoe prior when $b_0,b_g \ge1$. However, the performance slightly deteriorates when the hyperparameters are set to values smaller than 1.

$n$	NRMSE	$R^2$	Energy score
	$\mathbf{b}_0, \mathbf{b}_g = 1/p$
50	0.543 (0.155)	-27.810 (23.520)	3.295 (0.155)
100	0.626 (0.157)	-38.750 (28.170)	3.407 (0.172)
500	0.674 (0.189)	-45.410 (30.210)	3.470 (0.150)
1000	0.651 (0.219)	-46.080 (31.280)	3.431 (0.148)
	$\mathbf{b}_0, \mathbf{b}_g = 0.01$
50	0.319 (0.136)	-9.580 (11.300)	3.046 (0.129)
100	0.333 (0.135)	-10.990 (12.190)	3.051 (0.116)
500	0.281 (0.124)	-8.770 (9.170)	2.996 (0.079)
1000	0.227 (0.145)	-6.600 (9.880)	2.982 (0.105)
	$\mathbf{b}_0, \mathbf{b}_g = 1$
50	0.212 (0.070)	-3.070 (3.260)	3.001 (0.079)
100	0.215 (0.074)	-3.120 (3.750)	2.984 (0.065)
500	0.152 (0.073)	-1.550 (3.200)	2.934 (0.046)
1000	0.115 (0.055)	-0.470 (1.890)	2.914 (0.031)
	$\mathbf{b}_0, \mathbf{b}_g = 10$
50	0.204 (0.070)	-3.010 (3.380)	2.995 (0.071)
100	0.202 (0.071)	-2.930 (3.920)	2.979 (0.060)
500	0.135 (0.075)	-1.400 (3.740)	2.926 (0.035)
1000	0.098 (0.061)	-0.320 (2.130)	2.908 (0.023)
	$\mathbf{b}_0, \mathbf{b}_g = 100$
50	0.205 (0.071)	-3.060 (3.270)	2.995 (0.068)
100	0.204 (0.072)	-3.070 (3.690)	2.980 (0.055)
500	0.154 (0.082)	-1.810 (3.430)	2.935 (0.032)
1000	0.121 (0.079)	-0.740 (2.020)	2.916 (0.032)