Simulation has been an indispensable tool for the understanding and discovery of complex and open-ended phenomena in situ via the study of dynamic systems and processes in silico (Lavin et al., 2021). The studied phenomena run the gamut of scale and domain, from biology (Dada and Mendes, 2011) and climatology (Vautard et al., 2013) to economics and the social sciences (Elshafei et al., 2016). In simulation-based modeling, a parametric simulator is adopted as a mechanistic model of a physical system. Given specific parameter values, the simulator produces synthetic data. The general modeling principle is that parameters that lead to synthetic data close to the actual observations from the physical system are considered the most plausible ones to explain the measurements.

However, there are important challenges that have limited the adoption of simulators in many settings of scientific and engineering relevance. On the one hand, computational costs may motivate the imposition of simplifying assumptions, which may render the results unusable for reliable hypothesis testing. On the other hand, at a methodological level, simulators are poorly suited for statistical inference as they inherently provide only implicit access to the likelihood of an observation. In fact, simulators can sample from a distribution, but they cannot, typically, quantify the probability of a simulation output (Song et al., 2020; Simeone, 2022). These problems are currently being tackled using novel machine learning tools and probabilistic programming in the emerging field of simulation-based inference (SBI) (Cranmer et al., 2020).

In a frequentist setting, SBI produces point estimates for the simulator parameters, failing to capture epistemic uncertainty arising from the access to limited data from the physical system. Alternatively, adopting a Bayesian formulation, a distribution on the model parameter space can be optimized in order to reflect a probabilistic notion of uncertainty (Cranmer et al., 2020).

As illustrated in Figure 1, in this work, we study Bayesian SBI that leverages a parameterized quantum circuit (PQC) as the underlying simulator. PQCs are the subject of the field of quantum machine learning (Schuld and Petruccione, 2021). They consist of quantum circuits with a fixed ansatz whose parameters, typically rotation angles for some of the gates, are optimized using a classical computer. PQCs can be readily implemented on existing noisy intermediate scale quantum (NISQ) hardware, and are viewed as a potential means to demonstrate practical use cases for quantum computing.

The motivation for the proposed solution, integrating PQCs with SBI, is twofold. First, from the perspective of SBI, by leveraging quantum circuits as simulators, we follow the well-established principle that quantum computers are best suited for the simulation of certain physical phenomena, especially at the microscopic scale (Georgescu et al., 2014). Second, from the viewpoint of quantum machine learning, Bayesian learning methods have been argued to be potentially beneficial as they can be better account for uncertainty in the model space, enhancing test-time performance (Duffield et al., 2022). Our work thus contributes to the literature on quantum machine learning by moving beyond the likelihood-based method investigated in (Duffield et al., 2022) by leveraging state-of-the-art likelihood-free SBI methods.

This work explores the application of Bayesian SBI for the training of simulators implemented as PQCs. The main contributions are as follows.

• First, in a tutorial style, the article reviews two families of Bayesian SBI methods, namely, sampling-based and surrogate-based techniques. Sampling-based schemes aim at producing samples from the posterior distribution of the parameters, whilst surrogate-based methods estimate a surrogate for either the likelihood or directly for the posterior distribution in the model parameter space.

As shown in Figure 2, sampling-based Bayesian SBI, also known as ABC, models the data-generating mechanism via a hierarchical probability distribution. In it, the simulator p ( ⋅ | θ ) produces samples Z that are related to the true samples X by a mismatch model. This distribution posits an ancestral sampling procedure, whereby.

1. A model parameter is drawn from the prior θ ∼ p ( θ ) ;

The distribution p ( D | Z ) defining the mismatch model is subject to design, and it accounts for the fact that the simulator p ( ⋅ | θ ) is generally misspecified (Schmon et al., 2020). This is in the sense that there is typically no model parameter θ such that the simulator p ( ⋅ | θ ) matches exactly the data-generating distribution p * D = ∏ n = 1 N p * X n . (5) if no such mismatch is expected, one can set p D | Z = 1 D = Z , (6) where 1 ( ⋅ ) is the indicator function for discrete data and the Dirac impulse function for continuous-valued data.

More generally, the choice of the mismatch model must account for requirements of accuracy and efficiency, and it is typically specified as (Wilkinson, 2013; Schmon et al., 2020) p D | Z ∝ κ S D , S Z , (7) where κ ( ⋅ , ⋅ ) is a kernel function and S ( ⋅ ) ∈ R s is an s -dimensional summary statistic. The choice of the statistics S ( ⋅ ) is often done based on knowledge about the problem. For instance, for discrete-time sequences, one may use as statistics empirical transition rates (Sunnåker et al., 2013).

By (Equation 7), the data D are assumed to be more likely to correspond to samples Z if the data sets D and Z are “closer” in terms of the correlation between the statistics S ( D ) and S ( Z ) as measured by the kernel κ ( ⋅ , ⋅ ) . A typical choice for the kernel yields (Sunnåker et al., 2013) p D | Z ∝ 1 ρ S D , S Z ≤ ϵ , (8) where ρ ( ⋅ , ⋅ ) is an error measure and ϵ > 0 is tolerance level. Equivalently, the distribution (Equation 7) can be expressed as (Schmon et al., 2020) p D | Z ∝ e − ℓ S D , S Z , (9) where ℓ ( ⋅ , ⋅ ) is a loss function measuring the discrepancy between the statistics S ( D ) and S ( Z ) . Specifically, to match (Equation 7), one can set the loss as ℓ ( S , S ′ ) = − log ( κ ( S , S ′ ) ) .

By the model in Figure 2, the posterior distribution over model parameters and samples produced by the simulator given the data set D is given by p θ , Z | D ∝ p θ p Z | θ p D | Z . (10) by marginalizing out the simulator’s outputs Z , we obtain the model parameter posterior as p θ | D = ∫ p θ , Z | D d Z . (11) in the special case in which no mismatch is accounted for in the model, i.e., when p ( D | Z ) = 1 ( Z = D ) , then the distribution (Equation 11) evaluates as p θ | D = p θ , Z = D | D ∝ p θ p D | θ , (12) which corresponds to the conventional posterior distribution (Equation 1) in the ideal case of a well-specified model.

The goal of ABC is to produce samples θ of the model parameter vector that are approximately distributed according to the posterior distribution p ( θ | D ) in (Equation 11). This can be accomplished by generating samples ( θ , Z ) from the joint posterior p ( θ , Z | D ) in (Equation 10), and then discarding the simulator’s outputs Z . We next review two ABC methods of increasing complexity and efficacy.

Rejection-sampling ABC (RS-ABC) iteratively draws candidate samples θ ′ from the prior p ( θ ) . Each such sample is accepted with a probability that ensures that all accepted samples are drawn from the posterior p ( θ | D ) (Beaumont et al., 2002; Sunnåker et al., 2013).

To this end, for each candidate model parameter sample θ ′ , the simulator produces M samples Z ′ = { Z 1 ′ , Z 2 ′ , … , Z M ′ } distributed as p ( Z | θ ) = ∏ m = 1 M p ( Z m ′ | θ ′ ) . Therefore, RS-ABC produces the candidate pair ( θ ′ , Z ′ ) distributed as θ ′ , Z ′ ∼ p θ ∏ m = 1 M p Z m ′ | θ ′ . (13)

The sample ( θ ′ , Z ′ ) is selected with acceptance probability p ( a c c | θ ′ , Z ′ ) , where “ a c c ” denotes the event that a candidate sample is accepted. The acceptance probability generally depends on the model parameter θ ′ and on the generated samples Z ′ , as discussed next.

The distribution of an accepted sample can be computed as p θ , Z | a c c ∝ p θ p Z | θ p a c c | θ , Z . (14)

Therefore, in order for the distribution (Equation 14) to match the desired posterior (Equation 10), one can set p a c c | θ , Z ∝ p D | Z . (15)

In particular, if p ( D | Z ) is selected as per the conventional choice in (Equation 8), then the acceptance step is simplified as accept candidate sample  θ ′ , Z ′ if ρ S D , S Z ′ ≤ ϵ . (16)

RS-ABC typically produces a low rate of acceptance of the generated samples, particularly when data are sufficiently high dimensional. To see this, consider the common case in which the prior p ( θ ) supports, with non-negligible probability, model parameters θ corresponding to simulators p ( ⋅ | θ ) that are very different from the ground-truth distribution p * ( ⋅ ) . Using the conventional acceptance rule (Equation 16), a sample ( θ ′ , Z ′ ) is retained only if the sufficient statistics S ( D ) and S ( Z ′ ) for simulator’s samples and data set are sufficiently close. Given that RS-ABC draws the model parameter θ from the prior as per (3.2), this acceptance event is quite unlikely, resulting in a low rate of acceptance of the generated samples.

To overcome this drawback, reference (Marjoram et al., 2003) proposed Metropolis-Hastings ABC (MH-ABC). MH-ABC proceeds to sample from Equation 10) in a sequential manner. To elaborate, let us denote as θ i − 1 the last sample accepted at the beginning of the i -th iteration. Furthermore, we introduce a transition probability distribution q ( ⋅ | θ i − 1 ) , which constitute a key design choice for MH-ABC. At the i -th iteration, a new sample ( θ ′ , Z ′ ) is drawn from the conditional proposal distribution q θ ′ , Z ′ | θ i − 1 = q θ ′ | θ i − 1 p Z ′ | θ ′ . (17) accordingly, a new parameter θ ′ is sampled from the Markov transition kernel q ( ⋅ | θ i − 1 ) dependent on the last accepted sample θ i − 1 , and the data set Z ′ is generated from the simulator p ( ⋅ | θ ′ ) . The rationale for this choice is that samples close to previously accepted samples, as dictated by the distribution q ( ⋅ | θ i − 1 ) may be more likely to have the desired distribution. The downside of the approach is that consecutive accepted samples are not independent, as in RS-ABC. Rather, the temporal correlation is determined by the Markov kernel q ( ⋅ | θ i − 1 ) .

Let p ( a c c | θ ′ , Z ′ ) be the acceptance probability of the proposed sample ( θ ′ , Z ′ ) . In order to ensure a stationary distribution given by the posterior p ( θ , Z | D ) , it is sufficient to impose detailed-balance condition (Hastings, 1970) p θ i − 1 , Z i − 1 | D q θ ′ , Z ′ | θ i − 1 p a c c | θ ′ , Z ′ = p θ ′ , Z ′ | D q θ i − 1 , Z i − 1 | θ ′ p a c c | θ i − 1 , Z i − 1 . (18) this equality can be ensured by setting p a c c | θ ′ , Z ′ = min 1 , p θ ′ , Z ′ | D p θ i − 1 , Z i − 1 | D ⋅ q θ i − 1 , Z i − 1 | θ ′ q θ ′ , Z ′ | θ i − 1 . (19) using Equations 10, 17, we finally get the acceptance probability adopted by MH-ABC as p a c c | θ ′ , Z ′ = min 1 , p θ ′ p Z ′ | θ ′ p D | Z ′ p θ i − 1 p Z i − 1 | θ i − 1 p D | Z i − 1 ⋅ q θ i − 1 | θ ′ p Z i − 1 | θ i − 1 q θ ′ | θ i − 1 p Z ′ | θ ′ ) = min 1 , p θ ′ p D | Z ′ p θ i − 1 p D | Z i − 1 ⋅ q θ i − 1 | θ ′ q θ ′ | θ i − 1 . (20)

Ratio estimation (RE) applies contrastive learning to estimate the ratio between the likelihood p ( D | θ ) , which is not available in the likelihood-free setting of interest, and the data marginal p ( D ) , i.e. (Thomas et al., 2022), r D , θ = p D | θ p D . (21) given an estimate r ̂ ( D , θ ) of the ratio, the likelihood can be in principle estimated as p ̂ ( D | θ ) ∝ r ̂ ( D , θ ) , and the posterior distribution as p ̂ θ | D ∝ p θ p ̂ D | θ . (22)

In practice, the unnormalized posterior in (Equation 22) can be used, without the need for an explicit normalization, to obtain samples θ ∼ p ̂ ( θ | D ) with the aid of Markov chain Monte Carlo (MCMC) techniques (Cranmer et al., 2020; Simeone, 2022).

RE methods train a binary classifier to distinguish between data sets generated according to the distributions at the numerator and denominator of the ratio (Equation 21). To this end, for any fixed value θ , the simulator is first leveraged to generate two classes of data sets, with each data set containing a number N of examples. The first class of data sets contains data sets X drawn according to the distribution p ( X | θ ) = ∏ m = 1 M p ( X m | θ ) ; while the second class contains data sets drawn from the marginal p ( X ) = ∫ p ( X | θ ) p ( θ ) d θ .

To generate data sets in the first class, one directly runs the simulator with the given value θ . For the second class, one first samples θ ′ ∼ p ( θ ) from the prior, and then a data set X = { X m ∼ p ( X | θ ′ ) } m = 1 M , discarding the sample θ ′ . We assign label t = 1 to all data sets in the first class, and the label t = 0 to all data sets in the second class.

The binary classifier takes as input a data set X of M examples, computes a fixed function S ( X ) , and outputs a probability distribution p ( t | S ( X ) , ϕ ) quantifying the confidence of the predictor in X belonging to either class. The classifier depends on a model parameter vector ϕ . If the classifier is well trained, the probability p ( t | S ( X ) , ϕ ) provides a good approximation of the true posterior distribution p ( t | S ( X ) ) , as discussed next.

By the construction of the data set, a data set X is conditionally distributed as. p X | t = 1 = p X | θ (23) p X | t = 0 = p X . (24) furthermore, the posterior distribution is p t = 1 | X = p t = 1 p X | t = 1 p t = 0 p X | t = 0 + p t = 1 p X | t = 1 = p t = 1 p X | θ p t = 0 p X + p t = 1 p X | θ . (25) writing p ( t = 1 ) / p ( t = 0 ) = α we get p t = 1 | X = α p X | θ p X + α p X | θ (26) and p t = 0 | X = 1 − p t = 1 | X = p X p X + α p X | θ . (27)

Making the approximation p ( t | S ( X ) , ϕ ) ≈ p ( t | X ) , we can finally estimate the ratio of likelihood and data marginal as r ̂ D , θ ≈ p t = 1 | S D , ϕ p t = 0 | S D , ϕ . (28)

To improve the performance of RE, amortization techniques can be used, whereby the classifier is amortized using the parameters from the simulator. To explain, note that the true ratio can be equivalently expressed as r D , θ = p X | θ p X = p X , θ p X p θ . (29) This modification suggests a way to train the binary classifier to distinguish between dependent sample-parameter pairs ( X , θ ) ∼ p ( X , θ ) , which are assigned class label t = 1 , from independent sample-parameter pairs ( X , θ ) ∼ p ( X ) p ( θ ) , which are assigned class label t = 0 . To do so, samples from the first class are generated by running the simulator with the given value θ , and concatenating the sampled sequence to the vector θ . For the second class, one again first samples θ ′ ∼ p ( θ ) from the prior along with a data set X = { X m ∼ p ( X | θ ′ ) } m = 1 M , discarding the sample θ ′ ; and then generates new, independent, θ ′ ∼ p ( θ ) to which one concatenates the output of the simulator. This method is referred to as amortized RE (Hermans et al., 2020).

In cases when the divergence between the densities is large, the classifier can obtain almost perfect accuracy with a relatively poor estimate of the density ratio. This failure mode is known as the density-chasm problem, and can be overcome by transporting samples from one distribution to the other, creating a chain of intermediate data sets. The density-ratio between consecutive datasets along this chain can be then accurately estimated via classification. The chained ratios are then combined via a telescoping product to obtain an estimate of the original density-ratio. This method is referred to as telescopic amortized RE (Montel et al., 2023).

Finally, as practical note, we emphasize that, both in the amortized and non-amortized settings, to avoid numerical errors one can extract the logit, log ( r ̂ ( D , θ ) ) , from the classifier before applying the activation in the output layer. This choice also mitigates vanishing gradient issues.

The proposed quantum SBI solution aims at developing simulators for the generation of a quantity of interest X that can take values in a set of 2 d elements for some integer d . Note that this requires the quantity X to be either discrete to start with, or to be quantized as finely as allowed by a resolution of d bits. To this end, we propose to implement a PQC that acts on a register of d qubits. Accordingly, the allowed resolution of quantity X increases exponentially with the physical dimension of the qubit register. In such a setting, one can assume, without loss of generality, that the quantity X – or its quantized version–assumes values in the set of integers { 0,1 , … , 2 d − 1 } , or equivalently in the set of all binary strings of d bits.

As reviewed in (Schuld and Petruccione, 2021; Simeone et al., 2022), PQCs implement a parameterized unitary transformation U ( θ ) , a 2 d × 2 d complex-valued matrix, on a register of a given number, d , of qubits. The unitary transformation U ( θ ) is described by a quantum circuit that is specified by quantum gates placed according to a predefined arrangement. The arrangement is referred to as the ansatz of the PQC. Some of the quantum gates in the quantum circuit can be controlled by selecting real-valued parameters, which are collectively denoted as vector θ , see Figure 3.

Initializing the register of d qubits in a reference state | 0 〉 ⊗ d , the PQC produces the output state | ψ θ 〉 = U θ | 0 〉 ⊗ d . (30) Furthermore, by Born’s rule, the probability distribution p ( X | θ ) is given by p X | θ = | ⟨ X | ψ θ ⟩ | 2 , (31) where | X 〉 represents the state in the computational basis corresponding to integer X .

In general, choosing a good ansatz for the quantum simulator entails a difficult trade-off between adherence to the physics of the problem and complexity of implementation. In particular, if prior knowledge about the structure of the data are available, this may be encoded as an inductive bias into the choice of the quantum circuit architecture, assuming that the complexity of the implementation allows it.

The typical way to encode structure into the ansatz is to leverage symmetries in the data. Symmetries refer to transformations of the data that leave it invariant or change it in a predictable, equivariant manner. For example, the binding energy of a molecule does not change by permuting the order of the atoms, and a picture of a cat still depicts a cat regardless of the position of the cat within the image. This prior knowledge can be encoded into the simulator ansatz as a geometric prior. Notable examples include quantum graph neural networks (QGNNs) (Verdon et al., 2019; Mernyei et al., 2022) and quantum convolutional neural network (QCNNs) (Cong et al., 2019), which preserve equivariance to permutations and shifts, respectively. Other examples include quantum recurrent neural networks for time series processing (Nikoloska et al., 2023).

In the absence of prior knowledge, or when the practitioner is concerned with efficient hardware implementation, they may choose to use a hardware-efficient architecture (HEA). Such architectures use only single qubit and two qubit gates, placed along the existing connectivity of the quantum computer, which are easily implemented on both gate-based or pulse-based NISQ machines (Zulehner et al., 2018; Gyongyosi and Imre, 2021).

We consider four circuit architectures. All of the considered architectures are comprised of four qubits and two layers.

As a representative of sampling-based schemes, we implement RS-SBI with the classical kernel (Equation 8) with summary statistics given by the histogram of the generated samples Z . We set ϵ = 0.3 , and we draw M = 1000 examples for each draw from the prior distribution p ( θ ) . For surrogate-based schemes, we adopt the amortized RE technique whereby the surrogate model is implemented as three-member ensemble in which each member is comprised of a transformer layer with three attention heads followed by a linear layer with ReLU activations (Vaswani et al., 2017). The outputs of the ensemble members are averaged to obtain the logit. We use dropout with rate 0.1, and the Adam optimiser with learning rate 0.001.

We are interested in evaluating the adherence of the distribution of the samples produced by the simulator to the ground-truth data-generating distribution. To this end, for any fixed simulator parameters θ , we use the simulator to generate a large number of samples, namely, 1,000, from the corresponding model distribution p ( X | θ ) . The probability distribution p ( X | θ ) is estimated using the histogram of the generated samples. The quality of the samples is then quantified via the total variation distance (TVD) D ( p ( X | θ ) , p * ( X ) ) , where p * ( X ) is the ground-truth distribution, which in the examples at hand can be obtained from the training set.

In Bayesian SBI, the model parameter θ is drawn from the learnt posterior distribution, which represents the uncertainty of the learner on the optimal parameters of the simulator. In the examples at hand, the training data sets are sufficiently informative to fully describe the data-generating distribution p * ( X ) . However, epistemic uncertainty remains, owing to the unknown likelihood. In fact, the lack of access to the likelihood limits the information that the learning algorithm can extract from the simulator to the M samples drawn to evaluate the statistics S ( ⋅ ) , here the histogram.

Each draw of the model parameter vector θ yields a generally different TVD D ( p ( X | θ ) , p * ( X ) ) between the distribution produced by the simulator, p ( X | θ ) , and the ground-truth distribution, p * ( X ) . In the following, we evaluate the epistemic uncertainty produced by Bayesian SBI by plotting an estimate of the distribution of the TVD produced due to the randomness on the model parameter vector θ . The estimate is obtained via a kernel density estimator (KDE) with bandwidth equal to 0.9.

As a benchmark learning algorithm, we also show the performance of a scheme that produces a point estimate for the parameters θ . This strategy may be considered as a representative of frequentist learning methods that do not attempt to characterize epistemic uncertainty. Specifically, we implement a simple approach that looks for the value of parameters θ that maximizes the likelihood estimated via the TVD between the histogram of the M samples generated by the simulator and the ground-truth distribution. To this end, we retain the parameter vector θ that yields the minimum mentioned TVD across all samples θ generated by the considered Bayesian SBI schemes.

The distributions of the TVD for both sampling- and surrogate-based Bayesian SBI schemes are shown in Figure 4 for the two tasks under study. For this figure, we adopt the best-performing ansatz for each task, namely, the QCNN and QGNN, respectively. It is observed that, by accounting for the uncertainty on the likelihood, Bayesian SBI schemes can outperform conventional frequentist techniques. In fact, samples produced from the posterior distribution can yield significantly lower TVD values, which indicate a closer match of the ground-truth distribution. Furthermore, the spread of the distribution produced by Bayesian SBI strategies is task-dependent. Similarly, the choice between sampling-based and surrogate-based schemes is also seen to depend on the task, with the latter having a clear advantage in the BAS task.

We now analyze the impact of different ansatzes by showing in Figure 5 distributions of the TVD for various architectures of the quantum simulator. Whilst we do not claim that quantum circuits are provably better than classical counterparts for the problem at hand, the separable circuits is observed to result in a very large TVD. In contrast, for both tasks, the symmetry-preserving simulators–QCNN for the BAS task and QGNN for the molecular topology task–result in the smallest TVD between the generated samples and the true distribution, suggesting that encoding inductive-biases in the simulator is indeed helpful for SBI.