This problem is solved in our 2007 paper "Convergence of adaptive mixtures of importance sampling schemes" that considers the optimisation of $\mathbb{E}_{f}\left[ \log f_{\theta}(X) \right]$ using an integrated EM algorithm, that is an EM algorithm with no data. In short, you can run an EM algorithm with an extra integral over the (virtual) data.
If we consider the target $$\arg\max_{\alpha,\theta}\int \log\left( \sum_{d=1}^D\alpha_d q_d(x; \theta_d) \right) \pi(x)\,\text{d}x $$ the EM algorithm relies on introducing the latent variable $z$ such that the joint density $f$ of $x$ and $z$ satisfies $$ f(z)=\alpha_z \quad \text{and} \quad f(x|z)=q_z(x;\theta_z) \, . $$ The expectation corresponding to the E step of the EM algorithm is the expected complete log-likelihood, namely, at iteration $t$ of the algorithm, $$ \mathbb{E}_{\pi}^X \left[ \mathbb{E}_{(\alpha^t,\theta^t)}^Z \left\{ \log\left( \alpha_Z q_Z(X; \theta_Z) \right) | X \right\} \right]\,, $$ where the inner expectation is computed under the conditional distribution of $Z$ in the mixture model given the current value $(\alpha^t,\theta^t)$ of the parameters, i.e. $$ f(z|x)=\alpha^t_z q_z(x;\theta_z^t) \bigg/ \sum_{d=1}^D \alpha^t_d q_d(x;\theta_d^t)\,, $$ while the outer expectation is under the distribution $X \sim \pi$.
The corresponding M step amounts to setting the new parameters $(\alpha^{t+1},\theta^{t+1})$ equal to $$ (\alpha^{t+1},\theta^{t+1})=\arg\max_{(\alpha,\theta)}\,%\left( \mathbb{E}_\pi^X\left[\mathbb{E}_{(\alpha^t,\theta^t)}^Z\left\{\log(\alpha_Z q_Z(X; \theta_Z))|X\right\}\right]\,,%\right)\,. $$ It is straightforward to check that the convexity argument used for the EM algorithm also applies in this setup and, hence, that the Kullback-Leibler distance $(\mathfrak{K}(\pi,q_{(\alpha^t,\theta^t)}))_{t\geq 1}$ is a non-decreasing sequence. Setting $$ \rho_d(X;\alpha,\theta) = {\alpha_d q_d(X;\theta_d)}\bigg/{\sum_{\ell=1}^D \alpha_\ell q_\ell(X;\theta_\ell)} \, , $$ the maximisation program reduces to \begin{align*} & \alpha^{t+1}=\arg\max_{\alpha}\mathbb{E}_\pi^X\left[\sum_{d=1}^D \rho_d(X;\alpha^t,\theta^t)\log(\alpha_d)\right]\,, \\ & \theta^{t+1}=\arg\max_{\theta}\mathbb{E}_\pi^X\left[\sum_{d=1}^D \rho_d(X;\alpha^t,\theta^t)\log(q_d(X;\theta_d))\right]\,. \end{align*} where the first maximisation to be carried out under the constraint that $$\sum_{d=1}^D \alpha^{t+1}_d = 1$$ Hence, $$\begin{align*} & \alpha_d^{t+1}=\mathbb{E}_\pi^X\left[\rho_d(X;\alpha^t,\theta^t) \right]\,, \\ & \theta_d^{t+1}=\arg\max_{\theta_d}\mathbb{E}_\pi^X\left[ \rho_d(X;\alpha^t,\theta^t)\log(q_d(X;\theta_d)) \right]\,. \end{align*}$$