The only reason Variational Inference (VI) requires knowing the density $q(z|\lambda)$ is if you seek to maximize the ELBO w.r.t. $\lambda$, which is equivalent to minimizing $\text{KL}(q(z|\lambda) \mid\mid p(z|\mathcal{D}, \theta))$ – a KL divergence to the true posterior.
If instead, you're interested in optimizing the ELBO w.r.t. model parameters $\theta$ (as in VAEs), then you only need samples. So let's focus on the posterior approximation part.
In general, the approximation's objective can be equivalently rewritten as a difference of a cross-entropy and entropy of $q$: $$ \text{KL}(q(z|\lambda) \mid\mid p(z|\mathcal{D})) = \text{CE}[q(z|\lambda), p(z|\mathcal{D})] -\text{H}[q(z|\lambda)] $$ The density $q(z|\lambda)$ is needed for the latter. In fact, it's known that estimating the entropy from the samples only is very hard, things are much easier when the density is known.
However, you might know something about the density $q(z|\lambda)$. For example, if $q(z|\lambda)$ is a marginal distribution of some hierarchical process $q(z|\lambda) = \int q(z, \varepsilon |\lambda) d\varepsilon$, then an efficient variational bound can be given on the entropy.
A black-box alternative is to learn the density of $q(z|\lambda)$. In practice, you won't capture all the quirks of $q(z|\lambda)$ exactly, so this method lacks theoretical guarantees.
But is KL-divergence the only way to measure discrepancy between the true posterior and its approximation? Certainly, no, and other options might be more convenient. For example, Ruiz and Titsias have leveraged the fact that MCMC is guranteed to improve upon proposal to cancel out chain's marginal distribution.
