I don’t know if you are still interested in this issue. You can draw a bunch of samples at once, storing them in an extra dimension. You can then flatten this dimension and make the batch go through the decoder. You then reshape the tensor to reinsert the extra dimension (it should be of size L, following the notations in Kingma & Welling’s paper). You finally compute the likelihood and average across this dimension.
In principle, it sounds reasonable to expect that the running time goes up quite a bit. Imagine in each step the time it takes to perform the computations before the sampling is a, and the time after the sampling is b. For L = 1, the total time is a+b. If L>1, then total time is a+L*b. If b >> a, then total time is roughly multiplied by L, which means 1000% for L = 10.
By batching the computation as I described, you are able to benefit from the parallel matrix multiplications performed by the GPU. So if your batch size is small enough, the running time should not change.
Keep in mind I am not an expert.
I have written an implementation of this algorithm, but in my case setting L to anything other than 1 completely breaks the algorithm. It might be a bug, or a feature of the problem I am trying to solve that I don’t understand yet.
Keep me posted.