The E-step in the variational EM algorithm computes approximate posterior distributions of the latent variables, which are subsequently used to compute the expectations in Eq. REF . Particularly, there are two parameters to be estimated for each latent variable \(\mathbf {z}_{k,i}\) , which are the mean vector \(\mathbf {m}_{k,i}\)  and the covariance matrix \(\mathbf {S}_{k,i}\) . It is straightforward to maximize over these parameters by using the completing-the-square approach [1]} (See Appendix ). The terms that quadratically depend on \(\mathbf {z}_{k,i}\)  in the joint log-likelihood yield the posterior covariance update:
\( \mathbf {S}_{k,i} = \big [ \Sigma _k^{-1} + \sum _{l=1}^{L} N_{kl,i} \Theta _{kl}^T \mathbf {A}_{l} \Theta _{kl} \big ]^{-1},\) 
