Implementing VI is difficult in the likelihood-free setting. In the context of ABC, due to the intractability of ABC posteriors, traditional VI algorithms fail to work. Some adaptations have been made to the likelihood-free setting. EP-ABC, first introduced in [1]}, is an adaption of EP to the ABC setting. However, EP-ABC works only for likelihoods with a special structure, and is not designed for multimodal distributions. As shown in [1]}, when the posterior is multimodal, EP-ABC tends to produce an approximation covering all the modes, hence leading to a larger support. A variational Bayes with intractable likelihood (VBIL) method is developed in [3]}, which can be applied to commonly used statistical models without requiring an analytical solution to model-based expectations. However, VBIL aims at minimizing the KL of densities on an extended space \(\mathrm {KL}(q(\theta ,z) \Vert \pi (\theta ,z))\) , where \(\pi (\theta ,z)\)  admits the posterior \(\pi (\theta )\)  as its marginal and \(z\)  as a latent variable. This differs from the usual way of minimizing \(\mathrm {KL}(q(\theta ) \Vert \pi (\theta ))\) . To recover VI based on \(\mathrm {KL}(q(\theta ) \Vert \pi (\theta ))\)  in the likelihood-free setting, [4]} uses multilevel Monte Carlo [5]}, [6]} to derive unbiased estimates of the gradient of the KL.
On the other hand, [7]} modified the VBIL method to work with unbiased log-likelihood estimates in the synthetic likelihood framework, resulting in the variational Bayes synthetic likelihood (VBSL) method. A limitation of VBSL is its assumption on the normality of summary statistics.
Many variational methods including VBIL and VBSL are based on stochastic gradient decent (SGD) algorithms, which require an unbiased estimator of the gradient of KL. A difficulty with SGD is that plain Monte Carlo sampling to estimate the gradient can be error prone or inefficient. Some variance reduction methods have been adopted to improve SGD [8]}, [9]}.
