In recent years, deep generative models for multimodal data have gained significant attention. Among these, multimodal variational autoencoders (VAEs) have emerged as a promising approach, aiming to capture a shared latent representation by integrating information across different modalities through their inference models. A primary challenge for multimodal VAEs is accurately inferring representations from arbitrary subsets of modalities after learning a multimodal inference model. Naively, this would require training \( 2^M \) different inference networks ( \(M\) is # of modalities) to handle every possible combination of modalities, which is infeasible for a large number of modalities. Mixture-based models address this challenge by requiring only as many inference models as there are modalities, aggregating unimodal inferences to perform multimodal inference. However, when modalities are missing, these models suffer from information loss, particularly of modality-specific information, leading to deteriorated inference performance. Alternatively, alignment-based multimodal VAEs aim to align unimodal inference models with a multimodal inference model by minimizing the Kullback–Leibler (KL) divergence between them. Yet, the multimodal amortized inference, which is alignment source in these models inherently suffers from amortization gaps, preventing it from perfectly approximating the true inference and compromising the accuracy of unimodal inference. To address both issues, we introduce an iterative amortized inference mechanism within the multimodal VAE framework, termed multimodal iterative amortized inference. By iteratively refining the multimodal inference using all modalities, this method overcomes the information loss due to missing modalities in mixture-based models and minimizes the amortization gap in alignment-based models. Furthermore, by aligning the unimodal inference to approximate this refined multimodal posterior, we obtain unimodal inferences that effectively incorporate multimodal information while requiring only unimodal inputs at inference time. Experimental results on two benchmark datasets demonstrate that the proposed method improves the performance of the inference itself, suggested by higher linear classification accuracy and cosine similarity, and that the learned representations effectively capture the distributions of other modalities, as indicated by lower Fréchet Inception Distance (FID) scores in cross-modal generation. This indicates that the proposed approach significantly enhances the inferred representations from unimodal inputs.