How to implement Bayesian variational inference of arbitrary functions like "blackbox"?

Hello everyone, I am currently doing Bayesian posterior parameter estimation. At first I wanted to use MCMC to do it, but due to the large amount of observation data, the efficiency is a bit low. So I want to try using Bayesian variational inference. Special thanks to Dr Matthew Pitkin’s tutorial for running well. (Using a “black box” likelihood function in PyMC3). Now I have a question, can the application of this blackbox also be applied to Bayesian variational inference? Does the formula also hold true in theory? I read some conventional variational inferences and derivations, but I didn’t find any information for this kind of blackbox. I am not good at deriving mathematical formulas, so do you have any good suggestions and recommendations? Thank you all and wish you all the best.

Can we customize the likelihood function and gradient function of ADVI?