FuseMoE: Mixture-of-Experts Transformers for Fleximodal Fusion

Authors: xing  han
huy  nguyen
carl  harris
nhat  ho
suchi  saria
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-11 08:17:10

Abstract: As machine learning models in critical fields increasingly grapple with multimodal data, they face the dual challenges of handling a wide array of modalities, often incomplete due to missing elements, and the temporal irregularity and sparsity of collected samples. Successfully leveraging this complex data, while overcoming the scarcity of high-quality training samples, is key to improving these models' predictive performance. We introduce ``FuseMoE'', a mixture-of-experts framework incorporated with an innovative gating function. Designed to integrate a diverse number of modalities, FuseMoE is effective in managing scenarios with missing modalities and irregularly sampled data trajectories. Theoretically, our unique gating function contributes to enhanced convergence rates, leading to better performance in multiple downstream tasks. The practical utility of FuseMoE in the real world is validated by a diverse set of challenging prediction tasks.

Main paper: 
 introduction: 
 Multimodal fusion is a critical and extensively studied problem in many significant domains [78, 94, 90, 9] , such as sentiment analysis [26, 54] , image and video captioning [42, 41] , and medical prediction [33, 83] . Previous research has shown that embracing multimodality can improve predictive performance by capturing complementary information across modalities, outperforming single-modality approaches in similar tasks [71, 32] . However, an ongoing challenge lies in the creation of scalable frameworks for fusing multimodal data under a variety of conditions, and in creating reliable models that consistently surpass their single-modal counterparts.Table 1 : We evaluated the characteristics of FuseMoE against various benchmarks. The pipeline approach [83] relies on a simple feature extraction scheme for each modality, followed by concatenation and classification. It doesn't incorporate irregularities or missingness in its process, but its use of concatenation and zero-imputation for missing modalities allows it to be adapted to FlexiModal settings. Both [101] and [96] tackle multi-modality fusion, but as modalities increase, their method demands exponentially more cross-modal computations and significant model architecture modifications. Finally, [58] presents MoE for language-image alignment, yet it also requires substantial adjustments for the more intricate and universal FlexiModal context we explore.Theory FlexiModal Adaptive?Soenksen et al. [83] Data Pipeline ✗ ✗ ≥4 ✗ ✓ Zhang et al. [101] Modality Fusion ✓ ✗ 2 ✗ ✗ Zadeh et al. [96] Modality Fusion ✗ ✗ 3 ✗ ✗ Mustafa et al. [58] Multimodal MoE✗ ✗ 2 ✗ ✗ FuseMoE This Paper ✓ ✓ ≥4 ✓ AdaptedHandling a variable number of input modalities remains an open challenge in multimodal fusion, due to challenges with scalability and lack of unified approaches for addressing missing modalities. Many existing multimodal fusion methods are designed for only two modalities [25, 102, 101] , rely on costly pairwise comparisons between modalities [90] , or employ simple concatenation approaches [83] , rendering them unable to scale to settings with a large number of input modalities or adequately capture inter-modal interactions. Similarly, existing works are either unable to handle missing modalities entirely [101, 98] or use imputation approaches [89, 51, 83] of varying sophistication. The former methods restrict usage to cases where all modalities are completely observed, significantly diminishing their utility in settings where this is often not the case (such as in clinical applications); the latter can lead to suboptimal performance due to the inherent limitations of imputed data. In addition, the complex and irregular temporal dynamics present in multimodal data have often been overlooked [101, 88] , with existing methods often ignoring irregularity entirely [83] or relying on positional embedding schemes [90] that may not be appropriate when modalities display a varying degree of temporal irregularity. Consequently, there is a pressing need for more advanced and scalable multimodal fusion techniques that can efficiently handle a broader set of modalities, effectively manage missing and irregular data, and capture the nuanced inter-modal relationships necessary for robust and accurate prediction. We use the term FlexiModal Data to capture several of these key aspects, which haven't been well-addressed by prior works:"Flexi" suggests flexibility, indicating the possibility of having any combination of modalities, even with arbitrary missingness or irregularity.FlexiModal data is most evident in clinical scenarios, where extensive monitoring results in the accumulation of comprehensive electronic health records (EHRs) for each patient. A typical EHR encompasses diverse data types, including tabular (e.g., age, demographics, gender), images (Xrays, magnetic resonance imaging, and photographs), clinical notes, physiological time series (ECG and EEG), and vital signs (blood chemistry, heart rate). In this setting, we observe a variety of modalities, sampled with varying irregularity and a high degree of missingness and sparsity.Contributions In this paper, we introduce a novel mixture-of-experts (MoE) framework, which we call FuseMoE, specifically designed to enhance the multimodal fusion of FlexiModal data. FuseMoE incorporates sparsely gated MoE layers in its fusion component, which are adept at managing distinct tasks and learning optimal modality partitioning. In addition, FuseMoE surpasses previous cross-attention-based methods in scalability, accommodating an unlimited array of input modalities. Furthermore, FuseMoE routes each modality to designated experts that specialize in those specific data types. This allows FuseMoE to effectively handle scenarios with missing modalities by dynamically adjusting the influence of experts primarily responsible for the absent data, while still utilizing the available modalities. Lastly, FuseMoE integrates a novel Laplace gating function, which is theoretically proven to ensure better convergence rates compared to traditional Softmax functions, thereby enhancing predictive performance. We have conducted comprehensive empirical evaluations of FuseMoE across a range of application scenarios to validate its effectiveness. and irregular health status measurements over time; patients with milder conditions only require monitoring across fewer categories. FuseMoE is adept at handling inputs featuring any combination of modalities, including those with missing elements. It starts by encoding inputs using modality-specific feature extractors, followed by employing a multi-time attention mechanism [82] to address temporal irregularities. The core of FuseMoE lies the MoE Fusion Layer, where a routing mechanism is trained to categorize multimodal inputs and direct them to the appropriate combinations of MLPs. The outputs from these MLPs are weighted through a gating function, resulting in fused embeddings, which are subsequently utilized for further processing.2 FuseMoE: Enhance Predictive Performance for FlexiModal DataIn this section, we delve into the fundamental components of FuseMoE, illustrated in Figure 1 . We focus on two critical elements: the modality and irregularity encoder, and the MoE fusion layer. 
sparse moe backbone: 
 The main components of a sparse MoE layer are a network G as a sparse gate and an expert network E. [79] proposed a Top-K gating function that takes as an input a token representation x ∈ R D and then routes it to the Top-K experts out of the set {E i } S i=1 . The gating network parameter W ∈ R D×N produces logits h s (x) = Top K(x • W ), which are normalized via Softmax:EQUATIONEach expert network (E i : R D → R D ) contains a feed-forward layer (FFN) and its parameters are independent of other models. The final output of the expert network y is the linearly weighted combination of each expert's output on the token by the gate's output: y = S i=1 G(x) i E i (x). Gating Network Design The gating network's advantage lies in its capacity to be concurrently trained with FFNs, facilitating the learning of an optimal sparse combination of experts. Essentially, by evaluating the similarity between the input token and the experts, the gating network/router optimally matches the input partition with the most suitable experts. In many cases, variations in the routing mechanism can greatly influence performance across diverse applications [48] . The Softmax gating is the most widely adopted across domains [77, 79] . We introduce a novel Laplace gating function that offers enhanced convergence guarantees and delivers superior predictive performance, particularly in FlexiModal applications. The function is formulated as follows:h l (x) = Top K(-∥W -x∥ 2 ).(EQUATIONThe Laplace gating function, characterized by its Euclidean term exp(-∥W -x∥ 2 ), is less prone to converge towards extreme weight distributions due to the bounded nature of this term. In subsequent sections, we will illustrate how this gating function facilitates faster parameter estimation rates compared to Softmax gating. Moreover, our empirical findings indicate that the Laplace gating exhibits enhanced performance in managing FlexiModal data.To encode the irregularity of sampling in each modality, we utilize a discretized multi-time attention (mTAND) module [82] , which leverages a time attention mechanism [43, 92] to discretize irregularly sampled observations into discrete intervals. Specifically, given a set of l k continuous time points, an input scenario with three modalities: Time-Series (TS), Text (TXT), and images (IMG). (a) The joint router design utilizes a concatenated embedding of all modalities, directing this combined input to selected experts. (b) In the modality-specific router design, each modality's embedding is independently assigned to a shared pool of experts. (c) The third design variant also uses modality-specific routers but assigns each modality's embedding to separate pools of experts, each pool uniquely tailored to process a specific modality type.t ∈ R l k , corresponding to the k th dimensionality of a given modality, we employ H embedding functions ϕ h (τ ) to embed each τ k ∈ t k in a d h dimensional vector space (detailed definition and examples can be found in Appendix B and E). The i th dimension of the h th embedding is defined asϕ h (τ )[i] = w i τ k , if i = 1 sin(w i τ k + ϕ i ), if 1 < i ≤ d h , where {w i , ϕ i } d hi=1 are learnable parameters. By performing this for each continuous time point in t k , we create a d h dimensional representation of each time point in H different embedding spaces. We then leverage these embeddings to discretize the irregularly sampled observations into discretized bins. Specifically, we seek to discretize x k (with l k corresponding observation times t k ) into γ regularly sampled intervals γ. We do this via an attention mechanism, which, for each embedding function ϕ h (τ ), takes γ as queries, t k as keys, and x k as values and produces xk,h ∈ R γ embeddings for each sequence. Formally,xk,h = Softmax ϕ h (γ)Q h K ⊤ h ϕ h (t k ) ⊤ √ l k x k ,where Q h and K h are learnable parameters. This formulation allows us to discretize univariate observations x k into γ regularly-sampled bins. To model irregularity across a multivariate set of observations for a given modality with d m dimensions, we repeat this process for each dimension of the input. This allows us to obtain an interpolation matrix Xh = [x 1,h , x2,h , ..., xdm,h ] ∈ R γ×dm for each of the h embedding functions. We then concatenate the interpolation matrices across all H embedding functions (i.e., I = [ X1 , X2 , ..., Xh ] ∈ R γ×(H•dm) ) and employ a linear projection to achieve a final, discretized embedding for each modality, Z ∈ R γ×de , where d e denotes the desired dimensionality of each modality's representation. The discretization procedure offers a standardized approach to managing irregularly sampled time series across various input types; however, it can inevitably result in information loss. On the other hand, relying solely on the mTAND module may yield suboptimal performance due to the potentially varying sampling rates of different variables [31] , especially in scenarios where the sample sizes are small. To mitigate this, we combine discretized outputs with continuous representations learned through the mTAND module.Encoding Multiple Modalities The process described above allows us to discretize an arbitrarily long irregular, multivariate sequence into a regularly sampled, discretized embedding with length γ and dimensionality d e . We repeat this for each of the M modalities, to create M embeddings, {Z j } M j=1 , which are then combined to generate predictions.Router Design Study Upon obtaining embeddings from each of the j modalities, we propose multiple complementary approaches for processing multimodal inputs. Figure 2 illustrates a range of router design options. The most straightforward strategy involves employing a common router that handles the concatenated embeddings of all j modalities, without imposing any gating constraints.As the complexity increases with additional modalities, we consider more sophisticated alternatives: deploying separate routers for each modality's embedding and assigning these embeddings to a shared pool of experts. This allows for distinct processing while maintaining a unified expert framework. Additionally, we further segregate these common expert pools, allowing each router to direct its respective embedding to dedicated experts skilled in handling such specific inputs. These varied router design choices offer users enhanced flexibility, enabling more fine-grained control of both inter-modal and intra-modal relationships. Details of the respective advantages and challenges of these router design mechanisms can be found in Appendix C.We implement an entropy regularization loss to ensure balanced and stable expert utilization, a concept supported by various previous studies [58, 57, 22] . It maximizes the mutual information between modalities and experts and serves as an auxiliary loss function in addition to task-specific loss. Given a total of M modalities, and denoting H as the entropy, we define the loss function E asEQUATIONwhere pmj (E) is the distribution over the experts {E i } S i=1 for the j th modality. This distribution can be approximated by pmj (E) = 1l j l j i=1 p mj (E | x mj i ), where l j is the number of observations of the j th modality. Intuitively, we actively encourage the input embeddings to diminish the uncertainty in selecting experts. By incorporating the loss E, we aim to stabilize the experts' preferences within each modality, while promoting a diverse range of expert selections across different modalities.In scenarios where certain modalities are missing throughout the data trajectories, we substitute the original embedding Z missing with a learnable embedding Z, acting as a generic "missing indicator". This strategy is facilitated by employing per-modality routers, which, in conjunction with entropy regularization, guide Z predominantly toward a specific group of lessutilized experts. The new embeddings Z are dynamically adjusted throughout the model training process to minimize the task-specific loss and the entropy regularization loss. As a result, the router will assign lower weights to the experts responsible for processing these embeddings. 
theoretical contribution: 
 In this section, we provide a theoretical guarantee of the benefits of the Laplace gating over the standard Softmax gating in MoE. In particular, we conduct a convergence analysis for maximum likelihood estimation (MLE) under the Lapace gating Gaussian MoE, and demonstrate that the MLE under this model has better convergence behaviors than that under the softmax gating Gaussian MoE.Problem Setup. Since the convergence analysis of MLE under the Top-K sparse gating MoE has been studied in [63] , we will focus on examining the Laplace gating solely in the sequel. Assume thatEQUATIONwhere we define for any vectorsv = (v i ) k * i=1 that softmax(v i ) := exp(vi) k *j=1 exp(vj ) . Above, f (•|µ, ν) denotes a univariate Gaussian density function with mean µ and variance ν. For ease of the presentation, we denoteG * := k * i=1 exp(β * i )δ (W * i ,a * i ,b * i ,ν * i )as a true but unknown mixing measure associated with unknown parameters(β * i , W * i , a * i , b * i , ν * i ) for i ∈ {1, 2, . . . , k * }.In the paper, we specifically consider two settings of the true number of experts k * : (i) Exact-specified setting: when k * is known; (ii) Over-specified setting: when k * is unknown, and we over-specify the model in equation 4 by a Laplace gating MoE model with k > k * experts. However, due to the space limit, we present only the latter setting, and defer the former setting to Appendix J.Maximum Likelihood Estimation. We use the maximum likelihood method to estimate the unknown mixing measure G * [91] . In particular, the MLE is given byEQUATIONwhereG k (Θ) := {G = k ′ i=1 exp(β i )δ (Wi,ai,bi,νi) : 1 ≤ k ′ ≤ k, (W i , a i , b i , ν i ) ∈Θ} denotes the set of all mixing measures with at most k components. Given the MLE defined in equation 5, we are ready to present the main results. Before that, let us introduce some necessary notations for our analysis.Notations. We denote [n] := {1, 2, . . . , n} for any n ∈ N. For any vector v ∈ R d , ∥v∥ stands for its 2-norm value. Additionally, the notation |S| indicates the cardinality of a given set S, while δ denotes the Dirac delta measure. Finally, for any two probability densities p, q dominated by the Lebesgue measure µ, we denote V (p, q) = 1 2 |p -q|dµ as their Total Variation distance.Firstly, we demonstrate in Theorem 3.1 that the convergence rate of density estimation under the Laplace gating Gaussian MoE is parametric on the sample size n. Theorem 3.1 (Density estimation). The density estimation p Gn (Y |X) converges to the true density p G * (Y |X) under the Total Variation distance at the following rate:E X [V (p Gn (•|X), p G * (•|X))] = O( log(n)/n).Proof of Theorem 3.1 is in Appendix K.2. The parametric rate O( log(n)/n) of the conditional density function p Gn indicates that if there exists a loss function among parameters D such thatE X [V (p Gn (•|X), p G * (•|X))] ≳ D( G n , G * ), then we will achieve the parameter and expert estimation rates via the boundD( G n , G * ) = O( log(n)/n).Voronoi Loss. Following the above implication, we now define a loss function among parameters based on a notion of Voronoi cells as in [55] . Given some mixing measure G, we distribute its components θ i := (W i , a i , b i , ν i ) to the following Voronoi cells, which are generated by the components θ * j := (W * j , a * j , b * j , ν * j ) of the true mixing measure G * :EQUATIONfor any 1 ≤ j ≤ k * . Note that, the cardinality of the Voronoi cell A j is exactly the number of fitted components approximating θ * j . For ease of the presentation, let us denoteΦ ij (ρ 1 , ρ 2 , ρ 3 , ρ 4 ) := ∥W i -W * j ∥ ρ1 + ∥a i -a * j ∥ ρ2 + |b i -b * j | ρ3 + |ν i -ν * j | ρ4 , for any (ρ 1 , ρ 2 , ρ 3 , ρ 4 ) ∈ R 4 .Then, the Voronoi loss function D 2 (G, G * ) used for our analysis under the over-specified setting is given by:EQUATION)The notation r(|A j |) stands for the minimum value of r ∈ N such that the following system of polynomial equations does not have any non-trivial solutions for the unknown variables {(q 1i , q 2i , q 3i )}EQUATIONA solution to the above system is regarded as non-trivial if at least among variables q 1i is different from zero, whereas all the variables q 3i are non-zero. It is worth noting that the function r(•) was previously studied in [28] to characterize the convergence behavior of parameter estimation under the location-scale Gaussian mixture models. [28] also gave some specific values of that function, namely r(2) = 4 and r(3) = 6. Meanwhile, they claimed that it was non-trivial to determine the value of r(m) when m ≥ 4, and further techniques should be developed for that purpose. Since Gaussian MoE models are generalization of the Gaussian mixture models, we also involve the function r(•) in our convergence analysis. Now, we provide in the following theorem the convergence rate of parameter estimation under the over-specified setting of the Laplace gating Gaussian MoE model (see also Figure 3 for the empirical convergence rates justifying the theoretical rates in Theorem 3.2).Table 2 : Parameter estimation rates under the Softmax and Laplace gating Gaussian MoE models. The function r(•) represents the solvability of a system of polynomial equations considered in [64] while r(•) ≤ r(•) and r(2) = 4, r(3) = 6. Additionally, A n j := Aj( Gn) denotes a Voronoi cell defined in equation 6.Gates exp(β * j ) W * j a * j b * j ν * j Softmax [64] O(n -1/2 ) O(n -1/2 r(|A n j |) ) O(n -1/ r(|A n j |) ) O(n -1/2 r(|A n j |) ) O(n -1/ r(|A n j |) ) Laplace (Ours) O(n -1/2 ) O(n -1/4 ) O(n -1/4 ) O(n -1/2r(|A n j |) ) O(n -1/r(|A n j |) )Theorem 3.2 (Parameter Estimation). When k > k * becomes unknown, the following Total Variation bound holds true for any mixing measure G ∈ G k (Θ):E X [V (p G (•|X), p G * (•|X))] ≳ D 2 (G, G * ).Consequently, we obtain thatD 2 ( G n , G * ) = O( log(n)/n).Proof of Theorem 3.2 is in Appendix K.3. The results of Theorem 3.2 together with the formulation of the loss function D 2 in equation 7 reveal that (see also Table 2 ):(i) The parameters W * i , a * i , b * i , ν * i which are fitted by exactly one component, i.e. |A n i | := |A i ( G n )| = 1,enjoy the same estimation rate of order O(n -1/2 ) (up to some logarithmic factor), which match those in [64] .(ii) The rates for estimating the parameters W * i , a * i , b * i , ν * i which are fitted by more than one component, i.e. |A n i | > 1, are no longer homogeneous. On the one hand, the estimation rates for the parameters b * i and ν * i are of orders O(n -1/2r(|A n i |) ) and O(n -1/r(|A n i |) ), respectively, both of which are determined by the function r(•) and vary with the number of fitted components |A n i |. Those rates are comparable to their counterparts in [64] . On the other hand, the estimation rates for the gating parameters W * i and the expert parameters a * i are all of order O(n -1/4 ), which remains constant with respect to the number of fitted components. Meanwhile, those rates in [64] depend on a different system of polynomial equations from that in equation 8, which are significantly slower.Advantage of Laplace Gating on FlexiModal Setting In the standard Softmax gating [64] , the similarity score is computed as the inner product of a token's hidden representation and an expert embedding. However, this approach can lead to representation collapse [13, 69] , where a subset of experts dominates the decision-making process, resulting in the redundancy of other experts. This issue likely contributes to the slow rates of estimating expert parameters a * i in this setting (see Table 2 ). By contrast, the Laplace gating function partially alleviates this problem by computing the similarity score as the L 2 -distance between token representations and expert embeddings. This approach does not inherently favor any expert based on magnitude, unlike inner product which can be biased towards experts with larger norms. The Laplace gating ensures that all experts have a more balanced opportunity to be selected based on how close they are to the token representation. Therefore, Laplace gating is beneficial when dealing with heterogeneous inputs, such as multimodal data, where its feature distributions can be very different across modalities. This is because it can handle these differences without being overly sensitive to the scale and variance of the input features. In addition, it can gracefully degrade in the presence of missing data, rather than causing abrupt changes in gating probabilities that might occur with inner product-based measures. The improved estimation rates for expert parameters a * i under the Laplace gating Gaussian MoE, along with our empirical results on multiple large-scale datasets, substantiate these insights. 
experiments: 
 Overview We demonstrate that FuseMoE can provide accurate and efficient predictions when applied to the FlexiModal setting. We tested FuseMoE on a diverse set of benchmarks, including MIMIC-III [40] and MIMIC-IV [36] , CMU-MOSI and MOSEI [97] , the Physical Activity Monitoring (PAM) dataset [75] , and CIFAR-10 [46] . Compared to CMU-MOSI and MOSEI, the MIMIC ecosystem exhibits irregular and missing modality patterns and includes distinct modalities unlike PAM and CIFAR-10. Evaluating FuseMoE across these diverse datasets provides various empirical insights into critical aspects of our model's performance. Comprehensive details on the datasets, metrics, parameters, and additional results are thoroughly presented in the Appendices.We first apply our method to the CMU-MOSI and MOSEI datasets [97] , which utilize visual, acoustic, and textual data for sentiment analysis and emotion recognition tasks. Our methodology employs pre-trained T5 [73] for text encoding, librosa [56] for audio feature extraction, and EfficientNet [84] for video feature encoding. This improvement is attributed to the MoE's ability to effectively allocate specific components to handle distinct input modalities, thus better addressing both inter-and intra-modal relationships.CIFAR-10 Dataset Subsequently, we evaluate our method using the Vision-MoE framework [77] on the CIFAR-10 classification task [46] , with results illustrated in Figure 4 (a). In this experiment, we selectively replace the FFN layers with an even number in the Vision Transformer (ViT) models with MoE layers. These results, along with Table 3 on the CMU-MOSI and MOSEI datasets comparing Softmax-gating MoE, indicate that the Laplace gating function surpasses the standard Softmax gating function in performance. This outcome is consistent with our theoretical claims.We then conduct comprehensive evaluations of FuseMoE on MIMIC-IV [36] , and the Physical Activity Monitoring (PAM) dataset [75] . These datasets feature multiple input modalities, each characterized by varying degrees of irregular sampling or significant levels of missingness. Our tasks of interest for MIMIC datasets include the 48-hour in-hospital mortality prediction (48-IHM), 25-type phenotype classification (25-PHE), and length-of-stay (LOS) prediction. In addition to the previously mentioned baselines, we have incorporated the HAIM method [83] , a data pipeline specifically designed for integrating multimodal data from the MIMIC-IV dataset. We also include the cross-attention combined with irregular sequences modeling approach (MISTS) [101] . Table 4 shows the outcomes of combining irregular vital signs and clinical notes from the MIMIC-IV dataset. In addition to the commonly used Softmax gating function, we also evaluated the Gaussian gating function [93] as a comparative benchmark. The FuseMoE-based methods surpass baselines in most scenarios, often by a non-trivial margin. Furthermore, we observe that HAIM shows considerable efficacy in extracting features from time series, resulting in a strong performance in the 48-IHM and LOS tasks, which are heavily reliant on such data. However, its performance appears more moderate on the 25-PHE task. The PAM dataset captures daily living activities through 17 sensors, with data from each sensor treated as a separate modality. These modalities are individually processed through time-series and irregularity encoders before being integrated into the FuseMoE framework. Our baselines include the Transformer [92] , GRU-D [10] , SeFT [31] , a mTAND-only configuration, and IP-Net [81] . We use the Laplace gating and its joint experts & router structure in these experiments. The results in Figure 4 (b) have again shown the efficacy of integrating the irregularity encoder with the MoE fusion layer.Scalability of FuseMoE with Increasing Modalities Table 5 presents the revised outcomes of the MIMIC-IV dataset after integrating CXR and ECG of corresponding patients, employing the per-modality router and the entropy loss E within FuseMoE. This setup was chosen as it slightly outperformed the joint router with an increase in modalities. Relative to their two-modality versions, FuseMoE has effectively harnessed additional information (notably from CXR), resulting in a significant enhancement in performance. Conversely, the addition of new modalities did not benefit the HAIM method, possibly due to its reliance on vital signs and clinical notes without adequately addressing the dynamics between different modalities. Furthermore, HAIM's notably high F1 scores on the 25-PHE task can be attributed to XGBoost's proficiency in managing missing minority classes. Note that, except for HAIM, other baselines were not designed to be agnostic to the quantity and variety of input modalities. Therefore, adapting them to manage extra and missing modalities requires considerable model changes, which might compromise their performance.Missing Modalities Figure 4 (c) illustrates the effectiveness of utilizing per-modality routers and the entropy loss E in addressing missing modalities. Initially, we compare the performance of FuseMoE on patients with fully available modalities against those with missing components, employing a joint router mechanism with the importance loss function [79] , to ensure load balancing.The inclusion of datasets with missing modalities, while expanding the sample size, resulted in a decrease in performance due to the compromised data quality. However, a performance enhancement was observed upon integrating per-modality or disjoint routers with E. Notably, the outcomes for the 48-IHM and LOS tasks with missing modalities surpassed those obtained from datasets without any missingness. This is because the per-modality approach can better separate the present and missing modalities, reducing the influence of experts responsible for processing the absent inputs. Therefore, this leads to a more efficient exploitation of a broader array of samples. 
discussions and limitations: 
 In this paper, we introduced FuseMoE, a model adept at managing multimodal data characterized by random missingness or irregularity-a crucial yet relatively unexplored challenge. FuseMoE integrates MoE fusion layers with modal embeddings and offers multiple router configurations to adeptly handle multimodal inputs across different complexity levels. FuseMoE also employs an innovative Laplace gating function, which provides better theoretical results. Through empirical evaluation, FuseMoE has demonstrated superior performance across diverse scenarios. However, our current approach to encoding irregularities may potentially lead to over-parameterization when the input size is small. In our future work, we aim to identify simpler and more efficient methods to handle the irregularities of input samples while preserving the model's overall performance.Appendix for "FuseMoE: Mixture-of-Experts Transformers for Fleximodal Fusion" H. [7, 11, 70] , graphical models [59, 21, 76] , and neural networks [60, 20, 67] .With the diverse evolution of deep learning models, numerous advanced methods have now been employed in the fusion of multimodal data. In the realm of sentiment analysis, [96, 52] employ a low-rank Tensor Fusion method that leverages both language and video content. Attentiongating mechanisms are used by [74, 95] to generate displacement vectors through cross-modal self-attention, which are then added to the input vectors from the primary modality. [90] takes an alternative approach by integrating multiple layers of cross-modal attention blocks in a word-level vision/language/audio alignment task.In the context of clinical prediction, [45, 16] adopt a late fusion approach to combining vital sign and text data by concatenating embeddings from pre-trained feature extractors. [83] developed a generalizable data preprocessing and modeling pipeline for EHR encompassing four data modalities, albeit through a direct concatenation of existing feature embeddings for each modality followed by an XGBoost classifier [12] . Recently, [101] expanded on the work of [90] by introducing a discretized multi-time attention (mTAND) module [82] to encode temporal irregularities in time series and text data. Their fusion approach involves layering sets of self-and cross-modal attention blocks. However, this approach is limited to just two modalities and is not easily extendable to include additional modal components or handle missing modalities. To the best of our knowledge, existing works are tailored to application-specific settings that necessitate the computation of pairwise cross-modal relationships, which are not scalable to more general settings with arbitrary modalities. Moreover, these studies typically do not account for scenarios where modalities are missing, or rely on imputation approaches based on observed data.Mixture-of-Experts MoE [35, 93] has gained significant popularity for managing complex tasks since its introduction three decades ago. Unlike traditional models that reuse the same parameters for all inputs, MoE selects distinct parameters for each specific input. This results in a sparsely activated layer, enabling a substantial scaling of model capacity without a corresponding increase in computational cost. Recent studies have demonstrated the effectiveness of integrating MoE with cutting-edge models across a diverse range of tasks [79, 19, 103] . These works have also tackled key challenges such as accuracy and training instability [66, 104, 72] . Given its ability to assign input partitions to specialized experts, MoE naturally lends itself to multimodal applications. This approach has been explored in fields such as vision-language modeling [58, 80] and dynamic image fusion [9] . However, the application of MoE in complex real-world settings, such as those involving FlexiModal Data, remains largely unexplored. This gap presents an opportunity to leverage MoE's potential in handling its intricate and multifaceted nature such as multimodal EHR, where reliable multimodal integration is crucial.MoE Theory While MoE has been widely employed to scale up large models, its theoretical foundations have remained nascent. Recently, [64] provided convergence rates for both density and parameter estimation of Softmax gating Gaussian MoE. They connected these rates to the solvability of systems of polynomial equations under Voronoi-based loss functions. Later, [63] extended these theories to the top-K sparse softmax gating Gaussian MoE. Their theories further characterize the effect of the sparsity of gating functions on the behaviors of parameter estimation and verify the benefits of using top-1 sparse softmax gating MoE in practice. Other theoretical results for MoE include estimation rates of parameters and experts for multinomial logistic MoE [62] , for dense-tosparse gating MoE [61] , for Gaussian gating MoE [65] , and for input-independent gating MoE [29] .Definition of FlexiModal Data We provide a generic definition for the FlexiModal Data as we used throughout the paper. Let D = {(x m1 i , t m1 i ), (x m2 i , t m2 i ), . . . , (x mj i , t mj i ), y i } N i=1 to be the FlexiModal dataset with N units, where x mj i represents the input sequence from the i th unit of the j th modality, t mj i denotes the corresponding time points, and y i is the task-specific outcome. Take multimodal EHR as an example, each j th modality, which may vary from time-series data like heart rate, blood pressure, and glucose levels to high-dimensional inputs such as clinical notes and X-rays, contains l mj observations. Figure 5 is a more specific illustration of the FlexiModal example.Tasks In the ICU, where rapid and informed decisions are crucial, accurate mortality prediction is essential to provide clinicians with advanced warnings of patient deterioration, aiding in critical decision-making processes [5] . Similarly, the prediction of patient length-of-stay is indispensable for optimizing treatment plans, resource allocation, and discharge processes [6] . Further, phenotyping of critical care conditions is highly relevant to comorbidity detection and risk adjustment and presents a more challenging task than binary classification, due to the heterogeneous presentation of conditions and the larger number of prediction tasks [101] . We concentrate on three critical care tasks as highlighted in [44] , performing extensive empirical analysis on each building block of the proposed framework. care conditions [18, 53] (e.g., congestive heart failure, pneumonia, shock, etc.) at the end each each patient's ICU stay. Because the original task was designed for diagnoses based on ICD-9 codes, but MIMIC-IV includes both ICD-9 and ICD-10 codes, we map patients with diagnoses coded using ICD-10 using the conversion database provided by [8] . We implement an in-hospital mortality prediction (48-IHM) task to evaluate our method's ability to predict short-term patient deterioration. Similarly, an accurate determination of patient discharge times is crucial for optimizing patient outcomes and hospital resource allocation [6] , which motivates our length-of-stay (LOS) task. We frame 48-IHM and LOS as binary classification problems and use a 48-hour observation window (for patients who spent at least 48 hours in the ICU) to predict in-hospital mortality (48-IHM) and discharge (without expiration) within the 48 hours following the observation window (LOS). Lastly, identifying the presence of specific acute care conditions in patient records is essential for various clinical objectives, including the construction of cohorts for clinical studies and the detection of comorbidities [1] . Traditional methods, often reliant on manual chart reviews or simple billing code-based definitions, are increasingly being supplemented by machine learning techniques [27] ; automating this process requires high-fidelity classifications, motivating our 25-type phenotype classification (25-PHE) task. In this multilabel classification problem, we attempt to predict one of 25 acute care conditions using data from the entire ICU stay.Evaluation In our initial analysis, we focused on patients with no missing modalities, resulting in a dataset comprised of 8,770 ICU stays for the 48-IHM and LOS tasks, and 14,541 stays for the 25-PHE task. For our analyses with missing observations, we include a total of 35,129 stays for 48-IHM and LOS, and 71.173 for 25-PHE. To evaluate the single-label tasks, 48-IHM and LOS, we employ the F1-score and AUROC as our primary metrics. In line with previous studies [101, 50, 3] , we use macro-averaged F1-score and AUROC to assess the 25-PHE task.We leveraged data from MIMIC-IV [36] , a comprehensive database with records from nearly 300k patients admitted to a medical center from 2008 to 2019, focusing on the subset of 73,181 ICU stays. We were able to link core ICU records (containing lab results and vital signs) to corresponding chest X-rays [39] , radiological notes [38] , and electrocardiogram (ECG) data [24] taking place during a given ICU stay. We allocated 70 percent of the data for model training, with the remaining 30 percent evenly split between validation and testing.Missingness Rates The total number of samples for each of our three tasks (i.e., those in which at least one vital sign was recorded in the specified observation window), along with the total number of observations per-modality, are shown in Table 6 . Task We focus on the multimodal sentiment analysis (MSA) task which aims to predict sentiment polarity ∈ {positive, negative, and neutral} and sentiment intensity, which is a real number ranging from -3 to +3 under a multimodal setting.Evaluation Following previous work such as [26] , we adopt mean absolute error (MAE), Pearson correlation (Corr), binary classification accuracy, F1 score computed for non-negative/negative class as evaluation metrics. ,000 32x32 color images containing one of 10 object classes ("plane", "car", "bird", "cat", "deer", "dog", "frog", "horse", "ship", "truck"), with 6000 images per class.In this approach, a concatenated embedding of all modalities is created, and this combined input is directed to selected experts by the router. This method allows the model to capture interactions between modalities at the input level, as the concatenated embedding provides a unified representation that includes all modalities. The router and experts work with this comprehensive view, enabling the model to learn correlations and interactions directly from the fused data. However, this approach might not fully capture modality-specific nuances since the characteristics of each modality are blended into a single representation.1. Captures inter-modal relationships by considering all modalities together.2. Simplifies the routing mechanism by treating the concatenated embedding as a single input.1. May overlook modality-specific features due to the blending of all modalities into one representation.2. Could be less efficient if some modalities are irrelevant for certain tasks or experts.Each modality's embedding is independently assigned to a shared pool of experts by modalityspecific routers. This setup allows the model to maintain the distinctiveness of each modality while still leveraging a common pool of expertise. By doing so, it can better capture modality-specific nuances and how they contribute independently to the overall task. However, this approach might be less effective in capturing complex inter-modal interactions since the initial routing is done independently for each modality.1. Preserves modality-specific information by routing each modality independently.2. Flexible in directing modalities to the most relevant experts, potentially improving efficiency.3. Captures interactions between modalities to some extent, but may not be as effective as joint routing approaches.Challenge needs additional coordination between independent routes to leverage cross-modal insights.In this configuration, modality-specific routers assign each modality's embedding to separate pools of experts, with each pool uniquely tailored to process a specific modality type. This method maximizes the ability of the model to capture and exploit modality-specific features and relationships, as each pool of experts is optimized for a particular type of data. However, this setup might limit the model's ability to learn from the interactions between modalities, as each is processed in isolation.1. Allows for highly specialized processing of each modality, potentially improving performance on modality-specific tasks.2. Modality-specific experts can develop deeper insights into the characteristics and patterns within their designated data type.1. Inter-modal relationships might be underutilized due to the segregated processing of each modality.2. Requires additional coordination or subsequent integration stages to combine insights from different modality-specific experts.Each router type offers unique benefits and faces specific challenges in capturing the subtle relationships between modalities. The choice among them depends on the specific requirements of the application, including the importance of preserving modality-specific information versus capturing inter-modal interactions, and the computational efficiency of managing multiple experts and routers.For example, we found that modality-specific routers are more effective in ameliorating the effect of missing modality in our experiments.In the preprocessing stage, we focused on 30 pertinent lab and chart events from each patient's ICU record for vital sign measurements. For chest X-rays, we utilized a pre-trained DenseNet-121 model [14] , which was fine-tuned on the CheXpert dataset [34] , to extract 1024-dimensional image embeddings. For radiological notes, we obtained 768-dimensional embeddings using the BioClini-calBERT model [2] . ECG signals were processed using a convolutional autoencoder, adapted from [4] , to generate a 256-dimensional embedding for each ECG.We selected 30 time series events for inclusion, following [83] . Nine of these were vital signs: heart rate, mean/systolic/diastolic blood pressure, respiratory rate, oxygen saturation, and Glascow Coma Scale (GCS) verbal, eye, and motor response. We also included 21 lab values: potassium, sodium, chloride, creatinine, urea nitrogram, bicarbonate, anion gap, hemoglobin, hematocrit, magnesium, platlet count, phosphate, white blood cell count, total calcium, MCH, red blood cell count, MCHC, MCV, RDW, platlet count, neutrophil count, and vancomycin. We standard scale each time series value to have mean 0 and standard deviation 1, based on the values in the training set.Chest X-rays To incorporate a medical imaging modality into our analyses, we use the MIMIC-CXR-JPG [37] module available from Physionet [23] , which includes 377,110 JPG format images derived from the DICOM-based MIMIC-CXR database [39] . Following [83] , for each image, we resize each JPG image to 224 × 224 pixels and then extract embeddings from the last layer of the Densenet121 model. We identify X-rays taken while the patient was in the ICU by first matching subject IDs in MIMIC-CXR-JPG with the core MIMIC-IV database, then limiting these matched X-rays to those with a chart time occuring between an ICU admission and discharge.To incorporate text data, we use the MIMIC-IV-Note module [38] , which contains 2,321,355 deidentified radiology reports for 237,427 patients that can be matched with patients in the main MIMIC-IV via a similar approach to chest X-rays. We note that we were unable to obtain intermediate clinical notes (i.e., notes made by clinicians throughout a patient stay), as those have not yet been publicly released. We extract note embeddings using Bio-Clinical BERT [2] .Electrocardiograms (ECGs) To include ECGs as an additional modality in our models, we utilize the MIMIC-IV-ECG [24] module, which includes approximately 800,000 ECGs (10 seconds, sampled at 500 Hz) collected from nearly 160,000 unique patients. To transform the ECGs so that they are suitable for input to our model, we adopt a convolutional autoencoder approach, adapted from [4] , that compresses each ECG into a 256-dimensional vector. Specifically, each diagnostic ECG contains a 5000 × 12 dimensional vector (5000 time points × 12 ECG leads). To prepare the ECG for input to the autoencoder, we only include the first 4096 time points. We then train the autoencoder to compress the ECG into a 256-dimensional latent vector, and then reconstruct the original ECG using upsampling layers, using mean squared error as our loss function. The architecture is shown in Figure 6 . We train the autoencoder with 90% of the ECGs available in the MIMIC-IV-ECG projection and use the rest for validation. We selected a batch size of 2048, and reduced the learning rate by a factor of 0.5 if the validation loss had plateaued for 3 epochs. Training stopped if the validation loss had not decreased for 6 epochs. For our encoder, we use filter numbers of [16, 16, 32, 32, 64, 64] , kernel widths of [5, 5, 5, 3, 3, 3] and a dropout rate of 0.1. For the decoder, we use the same filter numbers and kernel widths in reverse, and maintain a dropout rate of 0.1.Figure 6 : CNN Autoencoder Architecture The encoder consists of 6 convolutional blocks (temporal convolution, batch normalization, dropout, and max pooling layers), followed by a dense layer that reduces the dimensionality of the representation of 256. The decoder reconstructs the input ECG (dimensionality 4096 × 12) from this latent vector via a dense layer, followed by 6 upsampling convolutional blocks (upsampling, convolutional, batch normalization, and dropout layers).We follow the preprocessing procedure from [100] as published from their official GitHub repository 1 .Unlike the embeddings in chest X-rays, clinical notes, and ECGs, vitals/lab/time-series values present temporal irregularity across dimensions. That is, for the former three modalities, each dimension of the corresponding is observed at each irregular time point τ . By contrast, the sampling for vitals/labs is irregular in both within and across dimensions. For example, we might observe heart rate values sampled at times τ HR = {0, 0.2, 0.8, 1.2, 2.8} and glucose values sampled at time τ Glu = {0.1, 0.7, 3.4}. Given this unique challenge present in vitals/labs, we adapt the Unified Temporal Discretization Embedding (UTDE) approach described in [101] , which combines the mTAND approach described in Section 2.2 with a simpler imputation-based discretization scheme. Specifically, given a set of t observations x ∈ R t observed at irregular times τ ∈ R t , we a simple imputation scheme to discretize x into target bins γ (e.g., γ = {0, 1, 2, ..., γ}. Specifically, given bin value γ i ∈ γ, we apply the following rules:• If there exists a previously observed value of x (i.e., ∃τ ∈ τ st. τ ≤ γ i ), we set the imputed value of x at time γ, xγi , to the closest previously observed value. • If no previously observed value exists, we set the value of xγi to the global mean of x.We do this for each possible vitals/lab, to generate a matrix of imputation embeddings I ∈ R γ×dvitals , were d vitals is the number of vitals/labs. We then input this embedding into a 1D causal convolutional layer with stride 1 to obtain our final imputation embeddings with hidden dimension d h , E Imp ∈ R γ×d h .We combined simple imputation and mTAND embeddings via a gating function g. Following [101] , we let E mTAND ∈ R γ×d h denote the mTAND embeddings for vitals/labs derived from the process described in Section 2.2 and let E Imp ∈ R γ×d h denote the simple imputations from the process described above. We use each of these discretization embeddings to derive a final set of embeddings for vitals/labs E vitals via a one-layer MLP gating function f . Specifically, we let g = f (E Imp ⊕ E mTAND ), where ⊕ denotes the concatenation operator. We then calculate E vitals asE vitals = g ⊙ E Imp + (1 -g) ⊙ E mTAND ∈ R γ×d h ,where ⊙ denotes point-wise multiplication.This approach, from [101] , casts time series and clinical notes as multivariate, irregularly-sampled time series (MISTS) and uses layers of self-and cross-attention to fuse modalities. The method uses a Time2Vec [43] encoding scheme to represent the irregularity of observation times. We use the same hyperparameters as in the original paper (e.g., 3 self-and cross-attention blocks, 128dimensional time embedding, etc.).This model from [90] relies on multiple stacks of pairwise and bidirectional cross-modal attention blocks (without a self-attention mechanism) to attend to low-level features. The results of crossmodal attention are then sent to modality-specific transformers, concatenated, and used to make predictions.This method introduces the Multimodal Adaptation Gate (MAG) as an extension to BERT and XL-Net, allowing these pre-trained models to incorporate visual and acoustic data during fine-tuning.By generating a modality-conditioned shift in their internal representations, MAG enables enhanced sentiment analysis performance on multimodal datasets, achieving human-level accuracy in the field [74] .The proposed Tensor Fusion Network approach (TFN) integrates three core components: Modality Embedding Subnetworks for generating rich embeddings from unimodal inputs, a Tensor Fusion Layer for capturing all levels of modality interactions through a 3-fold Cartesian product, and a Sentiment Inference Subnetwork tailored to perform sentiment analysis based on the fusion layer's output [96] .The multimodal fusion approach detailed by [83] extracts a single set of features for each ICU stay, and uses this to predict the outcome of interest (in-hospital mortality, etc.). For vitals/lab values, the authors extract a set of 11 generic time series features: signal length, maximum, minimum, mean, median, SD, variance, number of peaks, and average time-series slope and piece-wise change over time of these metrics. This is done independently for each of the 30 events, leading to 30×11 = 330 vital/lab features per ICU stay. To provide a fair comparison with our method, we only include the most recent five notes and 128 vitals measurements in calculating embeddings. We only include entries for which all modalities are observed. For note/X-ray/ECG embeddings, we compute the mean embedding across all observations occurring during the specified time frame (i.e., the first 48 hours of 48-IHM and LOS, the entire stay for PHE). As with our method, we standardize scale values based on the training set. [83] uses an XGBoost [12] classifier to predict the outcomes of interest. We follow the hyperparameter optimization approach described in the paper. Specifically, we conduct a grid search across the following sets of hyperparameters: max depth = {5, 6, 7, 8}, number of estimators = {200, 300}, learning rate = {0.3, 0.1, 0.05}. Hyperparameters are selected based on the maximum AU-ROC from five-fold cross-validation.We integrate F.1 through F.4 into our workflow using the implementation provided by [101] . For F.5, we adapt the time series (e.g., series variance, mean, etc.) feature extraction and model fitting code from the repository released by the corresponding paper. The original paper doesn't use ECG waveforms, so we adopt a similar approach to ECG embeddings as with image and note embeddings, and take the mean value of the latent vector across all included observations.We train models using a Lambda Workstation with four A550 GPUs with 24 GB of memory. We are able to train models using a single GPU. An analysis of computation time and memory requirements is shown in Figure 12 .The set of parameters we used for experiments can be found in Table 8 .We present additional results comparing FuseMoE to baselines using the MIMIC-III dataset, which includes only vital signs and clinical notes (Table 9 ), and the MIMIC-IV dataset, featuring vital signs and CXR (Table 10 ). All experiments utilize the "joint experts and router" configuration. In these settings, FuseMoE demonstrates noticeable advantages. In Figure 7 , we evaluate the impact of various irregularity encoders on the performance of the FuseMoE framework. Our baseline approaches include the following methods:1. employing only the imputation (discretization) module from the time-series irregularity encoder, as detailed in Appendix E 2. utilizing solely the mTAND module [82] within the time-series irregularity encoder 3. implementing the SeFT method [31] as an irregularity encoder 4. adopting the RAINDROP method [100] as an irregularity encoder In Figure 8 , we evaluate the impact of various time-series encoders on the performance of the Fuse-MoE framework. The original FuseMoE framework feeds time-series embeddings obtained from the irregularity encoder into the Transformer [92] and extracts the last hidden states of the Transformer output to pass through fully connected layers to make predictions. Our baseline approaches include CNN [47] and LSTM [30] to encode time-series embeddings from the irregularity encoder. The performance outcomes of these approaches are derived from averages over 5 random runs. We utilized the vital signs and clinical notes components of the MIMIC-IV dataset.In Figure 9 , we assess the effect of different CXR encoders on the FuseMoE framework. Currently, the FuseMoE framework incorporates DenseNet-121 as the feature extractor for CXR images before their integration into the mTAND module. This setup is compared with the application of the stateof-the-art vision transformer (ViT-B) [17] as an alternative CXR encoder.In Figure 10 , we evaluate the influence of text encoders on the FuseMoE framework. Currently, FuseMoE incorporates Clinical-Longformer [49] as the text encoder before integrating it into the mTAND module. This setup is compared with other state-of-the-art text encoders: GRU-D [10] , FT-LSTM [99] , and HierTrans [68] .Finally, in Figure 11 , we investigate the effect of the mTAND module on each modality, while we removed mTAND for a particular modality, the rest of FuseMoE's components remained constant.We then conducted ablation studies to explore the efficiency and effectiveness of MoE architecture on model performance. We mainly use MIMIC-IV as our test bed. the comparison. Despite the increase in model parameters due to the incorporation of the MoE layer, its sparse nature does not significantly escalate the computational load. Figure 12 (b) illustrates the correlation between the number of experts and task performance across different modalities. Generally, performance improves with the addition of more experts, plateauing once the count exceeds 16. To achieve a compromise between performance and computational expense, we opted to utilize the top 4 experts out of 16 in our experiments. Figure 12 (c) and Figure 13 study the influence of each modality on the top-k chosen experts. For every expert selected, we calculate the number of samples that include a specific modality, weighted by corresponding weight factors from the gating functions. The outcomes are subsequently normalized across modalities. The analysis of Figure 12 (c) reveals that predictions across all tasks heavily depend on vital signs and clinical notes. This reliance is attributed to the abundant samples in these two modalities. Despite the notably smaller quantity of CXR, they play more significant roles in the 25-PHE and 48-IHM tasks, which aligns with our findings in Table 5 . The results in Figure 13 demonstrate that the modality weight distribution in the MOSI and MOSEI datasets is more "spread out", with the audio component carrying a greater weight in the MOSEI dataset.We conduct multiple numerical experiments to illustrate the theoretical convergence rates of the MLE G n to the true mixing measure G * under both exact-specified and over-specified settings.Synthetic Data. Assume that the true mixing measureG * = k * i=1 exp(β * i )δ (W * i ,a * i ,b * i ,ν * i ) is of order k * = 2.The true parameters for the router, (W * i , β * i ) ∈ R d × R, are drawn independently from an isotropic Gaussian distribution with zero mean and variance σ 2 r = 0.01/d for 1 ≤ i ≤ 6, and otherwise are set to zero. Similarly, the true parameters of the experts, (a * i , b * i ) ∈ R d × R, are drawn independently of an isotropic Gaussian distribution with zero mean and variance σ 2 e = 1/d for all experts. For the variances ν * i , we also sample from the Gaussian distribution N (0, σ 2 e ), and then take the absolute value of the sample. These parameters remain unchanged for all experiments.Then, we generate i.i.d samples {(X i , Y i )} n i=1 by first sampling X i 's from the uniform distribution Uniform[0, 1] and then sampling Y i 's from the true conditional density g G * (Y |X) of the Laplace gating Gaussian mixture of experts (MoE) given in equation 4.A popular approach to determining the MLE G n for each set of samples is to use the EM algorithm [15] . However, since there are not any closedform expressions for updating the gating parameters β 0i , β 1i in the maximization steps, we have to leverage an EM-based numerical scheme, which was previously used in [63] . In particular, we utilize a simple coordinate gradient descent algorithm in the maximization steps. Additionally, we select the convergence criterion of ϵ = 10 -6 and run a maximum of 2000 EM iterations.Initialization. For each k ∈ {k * , k * + 1}, we randomly distribute elements of the set {1, 2, ..., k} into k * different Voronoi cells A 1 , A 2 , . . . , A k * , each contains at least one element. Moreover, we repeat this process for each replication. Subsequently, for each j ∈ [k * ], we initialize parameters W i by sampling from a Gaussian distribution centered around its true counterpart W * j with a small variance, where i ∈ A j . Other parameters β i , a i , b i , ν i are also initialized in a similar fashion.Under the exact-specified settings, we conduct 5 sample generations for each configuration, across a spectrum of 10 different sample sizes n ranging from 10 3 to 10 5 . It can be seen from Figure 3 (left) that the MLE G n empirically converges to the true mixing measure G * under the Voronoi metric D 1 at the rate of order O(n -0.49 ), which matches the theoretical parametric convergence rate established in Theorem J.2.Under the over-specified settings, we continue to generate 5 samples of size n for each setting, given 10 different choices of sample size n ∈ [10 3 , 10 5 ]. From Figure 3 (right), we observe that the MLE G n empirically converges to G * under the Voronoi metric D 2 at the rate of order O(n -0.47 ), which aligns with the theoretical parametric convergence rate established in Theorem 3.2.In this appendix, we study the theoretical behaviors of the MLE under the exact-specified setting, i.e., k = k * , of the Laplace gating Gaussian MoE. We demonstrate that under the exact-specified setting, the rate of estimated conditional density function p Gn to p G * is parametric On -1/2 ) (up to some logarithmic factor). Theorem J.1. The density estimation p Gn (Y |X) converges to the true density p G * (Y |X) under the Total Variation distance V at the following rate:E X [V (p Gn (•|X), p G * (•|X))] = O( log(n)/n).The proof of Theorem J.1 can be done similarly to that of Theorem 3.1 in Appendix K.2. The result of Theorem J.1 indicates that as long as we can establish the lower bound of the total variation distance between p Gn and p G * based on certain loss function between the MLE G n and the true mixing measure G * , we directly achieve the rate of the MLE under that loss function.Voronoi Loss We now define that loss function between the MLE and the true mixing measure for the exact-specified setting:EQUATIONAbove, for any(ρ 1 , ρ 2 , ρ 3 , ρ 4 ) ∈ R 4 , we define Φ ij (ρ 1 , ρ 2 , ρ 3 , ρ 4 ) = ∥W i -W * j ∥ ρ1 +∥a i -a * j ∥ ρ2 + |b i -b * j | ρ3 + |ν i -ν * j | ρ4 for any i ∈ A j and j ∈ [k * ].We demonstrate in the following theorem that the rate of MLE to the true mixing measure under the Voronoi loss function D 1 is O(n -1/2 ) (up to some logarithmic factor). Theorem J.2 (Exact-specified setting). When k = k * is known, the following Total Variation bound holds guarantetrue for any G ∈ G k (Θ):E X [V (p G (•|X), p G * (•|X))] ≳ D 1 (G, G * ).Therefore, we haveD 1 ( G n , G * ) = O( log(n)/n).Proof of Theorem J.2 is in Appendix K.1. The convergence rate of MLE under the Voronoi loss function D 1 implies that the rates of estimating the true parametersW * i , a * i , b * i , ν * i are also O(n -1/2) (up to logarithmic factors). These rates are comparable to those under the exact-specified setting of softmax gating Gaussian MoE (cf. Theorem 1 in [64] ).In this appendix, we provide proofs for all theoretical results in the paper. Throughout this appendix, for any vector v ∈ R d and α:= (α 1 , α 2 , . . . , α d ) ∈ N d , we denote v α = v α1 1 v α2 2 . . . v α d d , |v| := v 1 + v 2 + . . . + v d and α! := α 1 !α 2 ! . . . α d !. K.1 Proof of Theorem J.2First of all, we need to establish the following bound:E X [V (p G (•|X), p G * (•|X))] ≳ D 1 (G, G * ).For that sake, it is sufficient to demonstrate two following inequalities:• Inequality A. inf G∈G k * (Θ):D1(G,G * )≤ε ′ E X [V (p G (•|X), p G * (•|X))] D 1 (G, G * ) > 0; • Inequality B. inf G∈G k * (Θ):D1(G,G * )>ε ′ E X [V (p G (•|X), p G * (•|X))] D 1 (G, G * ) > 0,for some constant ε ′ > 0.The inequality A is equivalent tolim ε→0 inf G∈G k * (Θ):D1(G,G * )≤ε E X [V (p G (•|X), p G * (•|X))] D 1 (G, G * ) > 0.Assume that the above inequality is not true, then, there exists a sequence of mixing measureG n := k * i=1 exp(β n i )δ (W n i ,a n i ,b n i ,ν n i ) ∈ G k * (Θ) such that both D 1 (G n , G * ) and E X [V (p Gn (•|X), p G * (•|X))]/D 1 (G n , G * ) go to zero as n → ∞. Now, we define A n j = A j (G n ) := {i ∈ [k * ] : ∥θ n i -θ * j ∥ ≤ ∥θ n i -θ * τ ∥, ∀τ ̸ = j},for any j ∈ [k * ] as Voronoi cells with respect to the mixing measure G n , where we denoteθ n i := (W n i , a n i , b n i , ν n i ) and θ * j := (W * j , a * j , b * j , ν * j ).In this proof, since our arguments are assymptotic, we can assume without loss of generality (WLOG) that these Voronoi cells does not depend on n, that is, A j = A n j . Next, it follows from the hypothesis D 1n := D 1 (G n , G * ) → 0 as n → ∞ that each Voronoi cell contains only one element. Therefore, we may assume WLOG that A j = {j} for any j ∈ [k * ], which implies that (W n j , a n j , b n j , ν n j ) → (W * j , a * j , b * j , ν * j ) and exp(β n j ) → exp(β * j ) as n → ∞. Then, the loss function between G n and G * is given byD 1 (G n , G * ) = k * i=1 exp(β n i ) ∥∆W n i ∥ + ∥∆a n i ∥ + ∥∆b n i ∥ + ∥∆ν n i ∥ + exp(β n i ) -exp(β * i ) ,where we denote In this step, we aim to decompose the term∆β n 1i := β n 1i -β * 1i , ∆a n i := a n i -a * i , ∆b n i := b n i -b * i and ∆ν n i := ν n i -ν * i . Now,Q n := k * i=1 exp(-∥W * i -X∥ + β * i ) • [p Gn (Y |X) - p G * (Y |X)],which can be represented as follows:EQUATIONwhere we denoteF (Y |X; W, a, b, ν) := exp(-∥W -X∥)f (Y |a ⊤ X + b, ν) and H(Y |X; W ) = exp(-∥W -X∥)p Gn (Y |X). By applying the first-order Taylor expansion, we can rewrite A n asA n = k * i=1 |α|=1 exp(β n i ) α! • (∆W n i ) α1 (∆a n i ) α2 (∆b n i ) α3 (∆ν n i ) α4 × ∂ |α1|+|α2|+α3+α4 F ∂W α1 ∂a α2 ∂b α3 ∂ν α4 (Y |X; W * i , a * i , b * i , ν * i ) + R 1 (X, Y ) = k * i=1 |α|=1 exp(β n i ) α! • (∆W n i ) α1 (∆a n i ) α2 (∆b n i ) α3 (∆ν n i ) α4 × ∂ |α1| g ∂W α1 (X; W * i ) • ∂ |α2|+α3+α4 f ∂a α2 ∂b α3 ∂ν α4 (Y |(a * i ) ⊤ X + b * i , ν * i ) + R 1 (X, Y ),where R 1 (X, Y ) is a Taylor remainder that satisfies R 1 (X, Y )/D 1 (X, Y ) → 0 as n → ∞ and g(X, W ) := exp(∥W -X∥). Recall that f is the univariate Gaussian density, then by denoting h 1 (X; a, b) := a ⊤ X + b, we can verify that∂ α4 f ∂ν α4 (Y |(a * i ) ⊤ X + b * i , ν * i ) = 1 2 α4 • ∂ 2α4 f ∂h 2α4 1 (Y |(a * i ) ⊤ X + b * i , ν * i ).Consequently, we getEQUATIONwhere we denote η = α 3 + 2α 4 ∈ N.Subsequently, we also apply the first-order Taylor expansion to the term B n defined in equation 10 and get thatEQUATIONwhereR 2 (X, Y ) is a Taylor remainder such that R 2 (X, Y )/D 1 (G n , G * ) → 0 as n → ∞.From the above results, the term Q n can be rewritten asEQUATIONin which we respectively define for each i ∈ [k * ] thatS n i,α1,α2,η := α3+2α4=η, 0≤α3+α4≤1-|α1|-|α2| exp(β n i ) 2 α4 α! • (∆W n i ) α1 (∆a n i ) α2 (∆b n i ) α3 (∆ν n i ) α4 , T n i,γ := exp(β n i ) γ! (∆W n i ) γ , for any (α 1 , α 2 , η) ̸ = (0 d , 0 d , 0) and γ ̸ = 0 d . Otherwise, S n i,0 d ,0 d ,0 = T n i,0 d := exp(β n i )-exp(β * i ). Stage 2 -Non-vanishing coefficients:Moving to the second step, we will show that not all the ratios S n i,α1,α2,η /D 1 (G n , G * ) and T n i,γ /D 1 (G n , G * ) tend to zero as n → ∞. Assume by contrary that all of them approach zero when n → ∞, then for (α 1 , α 2 , η) = (0 d , 0 d , 0), it follows thatEQUATIONAdditionally, for tuples (α 1 , α 2 , η) where α 1 ∈ {e 1 , e 2 , . . . , e d } with e j := (0, . . . , 0, 1 j-th , 0, . . . , 0), α 2 = 0 d and η = 0, we get1 D 1 (G n , G * ) • k * i=1 exp(β n i )∥∆W n i ∥ 1 = k * i=1 |S n i,α1,α2,η | D 1 (G n , G * ) → 0.For (α 1 , α 2 , η) where α 1 = 0 d , α 2 ∈ {e 1 , e 2 , . . . , e d } and η = 0, we have1 D 1 (G n , G * ) • k * i=1 exp(β n i )∥∆a n i ∥ 1 = k * i=1 |S n i,α1,α2,η | D 1 (G n , G * ) → 0.For (α 1 , α 2 , η) where α 1 = α 2 = 0 d and η = 1, we have1 D 1 (G n , G * ) • k * i=1 exp(β n i )∥∆b n i ∥ 1 = k * i=1 |S n i,α1,α2,η | D 1 (G n , G * ) → 0.For (α 1 , α 2 , η) where α 1 = α 2 = 0 d and η = 2, we have1 D 1 (G n , G * ) • k * i=1 exp(β n i )∥∆ν n i ∥ 1 = k * i=1 |S n i,α1,α2,η | D 1 (G n , G * ) → 0.As a result, we achieve that1 D 1 (G n , G * ) • k * i=1 exp(β n i ) ∥∆W n i ∥ 1 + ∥∆a n i ∥ 1 + |∆b n i | + |∆ν n i | → 0.Due to the topological equivalence between norm-1 and norm-2, the above limit implies thatEQUATIONCombine equation 14 with equation 15, we deduce that D 1 (G n , G * )/D 1 (G n , G * ) → 0, which is a contradiction. Consequently, at least one among the ratios S n i,α1,α2,η /D 1 (G n , G * ) and T n i,γ /D 1 (G n , G * ) does not vanish as n tends to infinity. Stage 3 -Fatou's contradiction:In this step, we use the Fatou's lemma to point out a contradiction to the results achieved in Step 2.In particular, we denote by m n the maximum of the absolute values of S n i,α1,α2,η /D 1 (G n , G * ) andT n i,γ /D 1 (G n , G * ).Since at least one of the previous ratios does not converge to zero, we deduce that 1/m n ̸ → ∞.E X [V (p Gn (•|X), p G * (•|X))]/D 1 (G n , G * ) → 0 as n → ∞. Ac- cording to the Fatou's lemma, we have 0 = lim n→∞ E X [V (p Gn (•|X), p G * (•|X))] D 1 (G n , G * ) ≥ 1 2 • lim inf n→∞ |p Gn (Y |X) -p G * (Y |X)| D 1 (G n , G * ) d(X, Y ) ≥ 0. This result indicates that |p Gn (Y |X) -p G * (Y |X)|/D 1 (G n , G * )tends to zero as n goes to infinity for almost surely (X, Y ). As a result, it follows thatlim n→∞ Q n m n D 1 (G n , G * ) = lim n→∞ |p Gn (Y |X) -p G * (Y |X)| m n D 1 (G n , G * ) = 0.Next, let us denote S n i,α1,α2,η /[m n D 1 (G n , G * )] → ξ i,α1,α2,η and T n i,γ /[m n D 1 (G n , G * )] → κ i,γ with a note that at least one among them is non-zero. From the formulation of Q n in equation 13, we deduce thatEQUATIONfor almost surely (X, Y ). The above equation is equivalent toEQUATION2(1-|α1|-|α2|) η=0 ξ i,α1,α2,η • X α2 ∂ α2+η f ∂h α2+η 1 (Y |(a * i ) ⊤ X + b * i , ν * i ) + κ i,α1 p G * (Y |X)   × ∂ |α1| g ∂W α1 (X; W * i ) = 0,for almost surely (X, Y ). It is worth noting that parameters W * 1 , . . . , W * K are pair-wise distinct, thus, the set ∂ |α1| g ∂W α1 (X; W * i ) : i ∈ [k * ], 0 ≤ |α 1 | ≤ 1 is a linearly independent, which implies that 1-|α1| |α2|=0 2(1-|α1|-|α2|) η=0 ξ i,α1,α2,η • X α2 ∂ α2+η f ∂h α2+η 1 (Y |(a * i ) ⊤ X + b * i , ν * i ) + κ i,α1 p G * (Y |X) = 0, for any i ∈ [k * ], 0 ≤ |α 1 | ≤X α2 ∂ α2+η f ∂h α2+η 1 (Y |(a * i ) ⊤ X + b * i , ν * i ), p G * (Y |X) : 0 ≤ |α 2 | ≤ 1 -|α 1 |, 0 ≤ η ≤ 2(1 -|α 1 | -|α 2 |)is also linearly independent. Consequently, we obtain that ξ i,α1,α2,η = κ i,γ = 0 for any i∈ [k * ], 0 ≤ |α 1 | + α 2 ≤ 1, 0 ≤ η ≤ 2(1 -|α 1 | -|α 2 |) and 0 ≤ |γ| ≤ 1, which contradicts the fact that at least one among those terms is different from zero.Hence, we can find some constant ε ′ > 0 such that infG∈G k * (Θ):D1(G,G * )≤ε ′ E X [V (p G (•|X), p G * (•|X))] D 1 (G, G * ) > 0.Proof of inequality B: Assume by contrary that the inequality B does not hold, then there exists a sequence of mixing measuresG ′ n ∈ G k * (Θ) such that D 1 (G ′ n , G * ) > ε ′ and lim n→∞ E X [V (p G ′ n (•|X), p G * (•|X))] D 1 (G ′ n , G * ) = 0. This result leads to E X [V (p G ′ n (•|X), p G * (•|X))] → 0 as n → ∞.Recall that Ω is a compact set, therefore, we can replace the sequence G ′ n by one of its subsequences that converges to a mixing measureG ′ ∈ G k * (Θ). Since D 1 (G ′ n , G * ) > ε ′ , this result induces that D 1 (G ′ , G * ) > ε ′ .Subsequently, by means of the Fatou's lemma, we achieve that0 = lim n→∞ E X [2V (p G ′ n (•|X), p G * (•|X))] ≥ lim inf n→∞ p G ′ n (Y |X) -p G * (Y |X) d(X, Y ). It follows that p G ′ (Y |X) = p G * (Y |X)for almost surely (X, Y ). According to Lemma L.1, the noisy top-K sparse softmax gating Gaussian mixture of experts is identifiable, thus, we obtain that G ′ ≡ G * . As a consequence, we obtain that D 1 (G ′ , G * ) = 0, which contradicts to the fact thatD 1 (G ′ , G * ) > ε ′ > 0.Hence, the proof is completed.In this appendix, we employ results for M-estimators in [91] to establish the density estimation rate under the Laplace gating Gaussian mixture of experts (MoE).Firstly, we introduce some necessary notations and fundamental results. In particular, let P k (Θ) := {p G (Y |X) : G ∈ G k (Θ)}) := 1 2 ( √ p - √ q) 2 dµ 1/2for any probability densities p, q dominated by the Lebesgue measure µ. Then, we provide in the following lemma the upper bounds of those terms. Lemma K.1. If Θ is a bounded set, then the following inequalities hold for any 0 < η < 1/2:(i) log N (η, P k (Θ), ∥ • ∥ ∞ ) ≲ log(1/η); (ii) H B (η, P k (Θ), h) ≲ log(1/η).Proof of Lemma K.1 is in Appendix K.2.2. Subsequently, we denoteP k (Θ) := {p (G+G * )/2 (Y |X) : G ∈ G k (Θ)}; P 1/2 k (Θ) := {p 1/2 (G+G * )/2 (Y |X) : G ∈ G k (Θ)}.In addition, for each δ > 0, we define a Hellinger ball centered around the conditional density function p G * (Y |X) and intersected with the set P1/2 k (Θ) as P 1/2 k (Θ, δ) := {p 1/2 ∈ P 1/2 k (Θ) : h(p, p G * ) ≤ δ}.To capture the size of the above Hellinger ball, [91] suggest using the following quantity:EQUATIONwhere t∨δ := max{t, δ}. Given those notations, let us recall a standard result for density estimation in [91] .Lemma K.2 (Theorem 7.4, [91] ). Take Ψ(δ) ≥ J B (δ, P 1/2 k (Θ, δ)) such that Ψ(δ)/δ 2 is a nonincreasing function of δ. Then, for some sequence (δ n ) and universal constant c which satisfy √ nδ 2 n ≥ cΨ(δ), we obtain thatP E X h(p Gn (•|X), p G * (•|X)) > δ ≤ c exp(-nδ 2 /c 2 ),for any δ ≥ δ n Proof of Lemma K.2 can be found in [91] . Now, we are ready to provide the proof for convergence rate of density estimation in Theorem J.1 in Appendix K.2.1.It is worth noting that for any t > 0, we haveH B (t, P 1/2 k (Θ, t), ∥ • ∥ 2 ) ≤ H B (t, P k (Ω, t), h).Then, the integral in equation 17 is upper bounded as follows:EQUATIONwhere the second inequality follows from part (ii) of Lemma K.1.As a result, by choosing Ψ(δ) = δ • log(1/δ), we can verify that Ψ(δ)/δ 2 is a non-increasing function of δ. Furthermore, the inequality in equation 18 indicates that Ψ(δ) ≥ J B (δ, P 1/2 k (Θ, δ)). Next, let us consider a sequence (δ n ) defined as δ n := log(n)/n. This sequence can be validated to satisfy the condition √ nδ 2 n ≥ cΨ(δ) for some universal constant c. Therefore, by Lemma K.2, we reach the conclusion of Theorem J.1:P E X [h(p Gn (•|X), p G * (•|X))] > C log(n)/n ≲ n -c ,for some universal constant C depending only on Θ.Part (i). In this part, we will derive the following upper bound for the covering number of metric space (P k (Θ), ∥ • ∥ ∞ ) for any 0 < η < 1/2 given the bounded set Ω:log N (η, P k (Θ), ∥ • ∥ ∞ ) ≲ log(1/η).To start with, we denote Ω := {(a, b, ν) ∈ R d × R × R + : (β, W, a, b, ν) ∈ Ω}. As Θ is a bounded set, the set Ω is also bounded. Therefore, we can find an η-cover of Ω, denoted by Ω η . Additionally, we also define ∆ := {(β, W ) ∈ R × R d : (β, W, a, b, ν) ∈ Ω}, and ∆ η be an η-cover of ∆. Then, it can be validated that|Ω η | ≤ O(η -(d+2)k ), |∆ η | ≤ O(η -(d+1)k ).Next, for each mixing measure G = k i=1 exp(β i )δ (Wi,ai,bi,νi) ∈ G k (Θ), we take into account two other mixing measures. The first measure is G ′ = k i=1 exp(β i )δ (Wi,ai,bi,νi) , where(a i , b i , ν i ) ∈ Ω η is the closest points to (a i , b i , ν i ) in this set for all i ∈ [k]. The second one is G := k i=1 exp(β i )δ (W i,ai,bi,νi) in which (β i , W i ) ∈ ∆ η for any i ∈ [k]. Next, let us define T := {p G ∈ P k (Θ) : (β i , W i ) ∈ ∆ η , (a i , b i , ν i ) ∈ Θ η , ∀i ∈ [k]},then it is obvious that p G ∈ T . Now, we will show that T is an η-cover of metric space (P k (Θ), ∥ • ∥ ∞ ) with a note that it is not necessarily the smallest cover. Indeed, according to the triangle inequality, we haveEQUATIONSince the softmax function is no greater than one, the first term in the right hand side can be upper bounded as follows:EQUATIONSubsequently, we bound the second term ∥p G ′p G ∥ ∞ as follows:EQUATIONIt follows from the results in equation 19, equation 20 and equation 21 that ∥p Gp G ∥ ∞ ≲ η. This result indicates that T is an η-cover of the metric space (P k (Θ), ∥ • ∥ ∞ ). As a consequence, we obtain thatN (η, P k (Θ), ∥ • ∥ ∞ ) ≲ |∆ η | × |Ω η | ≤ O(1/η (2d+3)k ),which leads to the conclusion of this part: logN (η, P k (Θ), ∥ • ∥ ∞ ) ≲ log(1/η).Part (ii). In this part, we provide an upper bound for the bracketing entropy of P k (Θ) under the Hellinger distance h: H B (η, P k (Θ), h) ≲ log(1/η).Since Θ and X are bounded sets, there exist positive constants γ, ℓ, u such that -γ ≤ a ⊤ X + b ≤ γ and ℓ ≤ ν ≤ u. Let us defineB(Y |X) := 1 √ 2πℓ exp -Y 2 8u , for |Y | ≥ 2γ 1 √ 2πℓ , for |Y | < 2γ Then, it can be validated that f (Y |a ⊤ X + b, ν) ≤ B(Y |X) for any (X, Y ) ∈ X × Y.Next, let ζ ≤ η which will be chosen later and {p 1 , . . . , p N } be an ζ-cover of metric space(P k (Θ), ∥ • ∥ ∞ ) with the covering number N := N (ζ, P k (Θ), ∥ • ∥ ∞ ). Additionally, we also consider brackets of the form [Ψ L i (Y |X), Ψ U i (Y |X)] where Ψ L i (Y |X) := max{p i (Y |X) -ζ, 0} Ψ U i (Y |X) := max{p i (Y |X) + ζ, B(Y |X)}. Then, we can check that P k (Θ) ⊆ N i=1 [Ψ L i (Y |X), Ψ U i (Y |X)] and Ψ U i (Y |X) -Ψ L i (Y |X) ≤ min{2ζ, B(Y |X)}. Let S := max{2γ, √ 8u} log(1/ζ), we have for any i ∈ [N ] that ∥Ψ U i -Ψ L i ∥ 1 = |Y |<2γ [Ψ U i (Y |X) -Ψ L i (Y |X)] d(X, Y ) + |Y |≥2γ [Ψ U i (Y |X) -Ψ L i (Y |X)] d(X, Y ) ≤ Sζ + exp - S 2 2u ≤ S ′ ζ,where S ′ is some positive constant. This inequality indicates thatH B (S ′ ζ, P k (Θ), ∥ • ∥ 1 ) ≤ log N (ζ, P k (Θ), ∥ • ∥ ∞ ) ≤ log(1/ζ).By setting ζ = η/S ′ , we obtain that H B (η, P k (Θ), ∥ • ∥ 1 ) ≲ log(1/η). Finally, due to the inequality h 2 ≤ ∥ • ∥ 1 , we reach the conclusion of this part:H B (η, P k (Θ), h) ≲ log(1/η).Hence, the proof is completed.In order to establish the following Total Variation lower bound under the over-specified settings, i.e. when k > k * is unknown:E X [V (p G (•|X), p G * (•|X))] ≳ D 2 (G, G * ), we need to prove two following inequalities:• Inequality A. inf G∈G k (Θ):D2(G,G * )≤ε ′ E X [V (p G (•|X), p G * (•|X))] D 2 (G, G * ) > 0; • Inequality B. inf G∈G k (Θ):D2(G,G * )>ε ′ E X [V (p G (•|X), p G * (•|X))] D 2 (G, G * ) > 0,for some constant ε ′ > 0. As the inequality B can be achieved in the same fashion as in Appendix K.1, we concentrate on showing the inequality A in this proof. For that purpose, it suffices to prove thatEQUATIONAssume that the above claim does not hold true, then there exists a sequence of mixing measuresEQUATION) Since D 2 (G n , G * ) → 0, we deduce that i∈Aj exp(β n i ) → exp(β * j ) and (W n i , a n i , b n i , ν n i ) → (W * j , a * j , b * j , ν * j ) for all i ∈ A j and j ∈ [k * ]. Now, we reuse the three-step framework in Appendix K.1.Firstly, by abuse of notations, let us consider the quantityQ n := k * j=1 exp(-∥W * j -X∥ + β * j ) • [p Gn (Y |X) -p G * (Y |X)].Similar to Step 1 in Appendix K.1, we can express this term asQ n = k * j=1 i∈Aj exp(β n i ) F (Y |X; W n i , a n i , b n i , ν n i ) -F (Y |X; W * j , a * j , b * j , ν * j ) - k * j=1 i∈Aj exp(β n i ) H(Y |X; W n i ) -H(Y |X; W * j ) + k * j=1 i∈Aj exp(β n i ) -exp(β * j ) F (Y |X; W * j , a * j , b * j , ν * j ) -H(Y |X, W * j ) : = A n -B n + E n ,Next, we proceed to decompose A n based on the cardinality of the Voronoi cells as follows:A n = j:|Aj |=1 i∈Aj exp(β n i ) F (Y |X; W n i , a n i , b n i , ν n i ) -F (Y |X; W * j , a * j , b * j , ν * j ) + j:|Aj |>1 i∈Aj exp(β n i ) F (Y |X; W n i , a n i , b n i , ν n i ) -F (Y |X; W * j , a * j , b * j , ν * j ) .By applying the Taylor expansions of order 1 and rj to the first and second terms of A n , respectively, and following the derivation in equation 11, we get thatEQUATION2(rj -|α1|-|α2|) η=0 α3+2α4=η, 0≤α3+α4≤rj -|α1|-|α2| exp(β n i ) 2 α4 α! • (∆W n ij ) α1 (∆a n ij ) α2 × (∆b n ij ) α3 (∆ν n ij ) α4 • X α2 • ∂ |α1| g ∂W α1 (X; W * j ) • ∂ |α2|+η f ∂h |α2|+η 1 (Y |(a * j ) ⊤ X + b * j , ν * j ) + R 4 (X, Y ) where R i (X, Y ) is a Taylor remainder such that R i (X, Y )/D 2 (G n , G * ) → 0 as n → ∞ for i ∈ {3, 4}.Next, we apply the Taylor expansions of order 1 and 2 to the first and second terms of B n , respectively, and following the derivation in equation 12, we get thatB n = j:|Aj |=1 i∈Aj |γ|=1 exp(β n i ) γ! (∆W n ij ) γ • ∂ |γ| g ∂W γ (X; W * j )p Gn (Y |X) + R 5 (X, Y ) j:|Aj |>1 i∈Aj 2 |γ|=1 exp(β n i ) γ! (∆W n ij ) γ • ∂ |γ| g ∂W γ (X; W * j )p Gn (Y |X) + R 6 (X, Y ),where R 5 (X, Y ) and R 6 (X, Y ) are Taylor remainders such that their ratios over D 2 (G n , G * ) approach zero as n → ∞. Subsequently, let us defineS n j,α1,α2,η := i∈Aj α3+2α4=η, 0≤α3+α4≤rj -|α1|-|α2| exp(β n i ) 2 α4 α! • (∆W n ij ) α1 (∆a n ij ) α2 (∆b n ij ) α3 (∆ν n ij ) α4 , T n j,γ := i∈Aj exp(β n i ) γ! (∆W n ij ) γ , for any (α 1 , α 2 , η) ̸ = (0 d , 0 d , 0) and γ ̸ = 0 d . Otherwise, S n j,0 d ,0 d ,0 = T n j,0 d := i∈Aj exp(β n i ) - exp(β * j ).As a consequence, it follows thatEQUATIONEQUATIONStage 2 -Non-vanishing coefficients:In this step, we demonstrate that not all the ratios S n j,α1,α2,η /D 2 (G n , G * ) and T n j,γ /D 2 (G n , G * ) converge to zero as n → ∞. Assume by contrary that all these terms go to zero. Then, by employing arguments for deriving equation 14 and equation 15, we get that1 D 2 (G n , G * ) • k * j=1 i∈Aj exp(β n i ) -exp(β * j ) + j:|Aj |=1 i∈Aj exp(β n i ) ∥∆W n ij ∥ + ∥∆a n ij ∥ + |∆b n ij | + |∆ν n ij | → 0.Taking the summation of j:|Aj |>1 |S n j,α 1 ,α 2 ,η | D2(Gn,G * ) for all (α 1 , α 2 , η) where α 1 ∈ {2e 1 , 2e 2 , . . . , 2e d }, α 2 = 0 d and η = 0, we have1 D 2 (G n , G * ) • j:|Aj |>1 i∈Aj exp(β n i )∥∆W n ij ∥ 2 → 0.Taking the summation ofj:|Aj |>1 |S n j,α 1 ,α 2 ,η | D2(Gn,G * ) for all (α 1 , α 2 , η) where α 1 = 0 d , α 2 ∈ {2e 1 , 2e 2 , . . . , 2e d } and η = 0, we have 1 D 2 (G n , G * ) • j:|Aj |>1 i∈Aj exp(β n i )∥∆a n ij ∥ 2 → 0.Combine the above limit with the formulation of D 2 (G n , G * ) in equation 23, we have that1 D 2 (G n , G * ) • j:|Aj |>1 i∈Aj exp(β n i ) |∆b n ij | rj + |∆ν n ij | rj 2 ̸ → 0.This result implies that we can find some index j′ ∈ [k * ] : |A j ′ | > 1 that satisfies 1 D 2 (G n , G * ) • i∈A j ′ exp(β n i ) |∆b n ij ′ | rj ′ + |∆ν n ij ′ | rj ′ 2 ̸ → 0.For simplicity, we may assume that j ′ = 1. Since S n 1,0 d ,0 d ,η /D 2 (G n , G * ) vanishes as n → ∞ for any 1 ≤ η ≤ rj , we divide this term by the left hand side of the above equation and achieve thatEQUATIONfor any 1 ≤ η ≤ r1 .Subsequently, we define M n := max{|∆b n i1 |, |∆ν n i1 | 1/2 : i ∈ A 1 } and π n := max{exp(β n i ) : i ∈ A 1 }. As a result, the sequence exp(β n i )/π n is bounded, which indicates that we can substitute it with its subsequence that admits a positive limit z 2 5i := lim n→∞ exp(β n i )/π n . Therefore, at least one among the limits z 2 5i equals to one. Furthermore, we also denote(∆b n i1 )/M n → z 3i , (∆ν n i1 )/(2M n ) → z 4i .From the above definition, it follows that at least one among the limits z 3i and z 4i equals to either 1 or -1. By dividing both the numerator and the denominator of the term in equation 25 by π n M η n , we arrive at the following system of polynomial equations:i∈A1 α3+2α4=η, 1≤α3+α4≤r1 z 2 5i z α3 3i z α4 4i α 3 ! α 4 ! = 0,for all 1 ≤ η ≤ r1 . Nevertheless, from the definition of r1 , we know that the above system does not admit any non-trivial solutions, which is a contradiction. Consequently, not all the ratios S n j,α1,α2,η /D 2 (G n , G * ) and T n j,γ /D 2 (G n , G * ) tend to zero as n → ∞. Stage 3 -Fatou's contradiction:Recall that E X [V (p Gn (•|X), p G * (•|X))]/D 2 (G n , G * ) → 0 as n → ∞.Then, by applying the Fatou's lemma, we get 0 = lim n→∞ E X [V (p Gn (•|X), p G * (•|X))] D 2 (G n , G * ) ≥ 1 2 • lim inf n→∞ |p Gn (Y |X) -p G * (Y |X)| D 2 (G n , G * ) d(X, Y ), which implies that |p Gn (Y |X)-p G * (Y |X)|/D 2 (G n , G * ) → 0 as n → ∞ for almost surely (X, Y ).EQUATIONDue to the identifiability of the location-scale Gaussian mixtures [85, 86, 87] , we get that k = k * and softmax(-∥W i -X∥+ β i ) : i ∈ [k] ≡ softmax(-∥W * i -X∥ + β * i ) : i ∈ [k] ,for almost surely X. WLOG, we may assume that softmax(-∥W i -X∥EQUATIONfor almost surely X for any i ∈ [k]. Since the softmax function is invariant to translations, it follows from equation 28 that W i = W * i and β i = β * i + v 0 for some v 0 ∈ R. Notably, from the assumption of the model, we have β k = β * k = 0, which implies that v 0 = 0. As a result, we obtain that β i = β * i for any i ∈ [k * ]. Then, equation 27 can be rewritten asEQUATIONfor almost surely (X, Y ). Next, we denote J 1 , J 2 , . . . , J m as a partition of the index set [k * ], where m ≤ k * , such that exp(β i ) = exp(β i ′ ) for any i, i ′ ∈ J j and j ∈ [m]. On the other hand, when i and i ′ do not belong to the same set J j , we let exp(β i ) ̸ = exp(β i ′ ). Thus, we can reformulate equation 29 as Hence, we reach the conclusion of this lemma.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: computational resources used to run the experiments can be found in Appendix G. Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? Answer: [Yes] Justification: we have reviewed the NeurIPS Code of Ethics and make sure the paper conforms to this. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [Yes] Justification: we have attempted to discuss the impacts of this work in Appendix M. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: our paper has no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: all the assets used in the paper have been properly credited.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.https://github.com/mims-harvard/Raindrop/tree/main/PAMdataThis paper presents research aimed at propelling advancements in the broad domain of machine learning. The implications of our findings are wide-ranging, with potential applications in sectors including healthcare, autonomous driving, and recommendation systems. Based on our current understanding, this research does not warrant an ethics review, and a detailed discussion of the potential societal impacts is not required at the current stage. 
 
 Appendix:  


Decision: This paper presents work on multi-modal data, particularly including time series data.  A mixture of experts framework is developed, with a Laplace gating function.  A variety of fusion/gating strategies is explored as examples of the general MoE / gating approach.

The reviewers appreciated the problem formulation and contribution of the Laplace gating function with its accompanying theoretical analysis.  Concerns were raised over the computation needs, scalability, and effectiveness of the routing/gating mechnanism versus the modality-specific encoders.

After considering the authors' response the reviewers evaluated their positions.  The additional notes on experiments, comparisons, and clarificiations resolved most of the reviewers' concerns.  

Overall, the paper presents interesting work on learning from multi-modal data with a solid contribution in the form of a novel gating function and promising results.  As such, it is recommended for acceptance in NeurIPS.

Reviews: review: paper_summary: paper_summary: This papers proposes a novel MoE arcchitecture for handling and fusing multiple modalities, along with two core contributions:
- a novel router design that can handle missing modalities;
- a laplace gating function that is theoretically proven to ensure better convergence.
strength_weakness: strength_weakness: - 1 The proposed method achieves promising results on the tested benchmark, which are still comparably small scale datasets, e.g, CMU-MOSI is with 2000+ samples. I am wondering if this work can be applicable to training/fine-tuning with larger scale data. 
- 2 Again, for vision task, only CIFAR-10 is used. This cannot fully justify its effectiveness for more general vision or multi-modal tasks;
- 3 The capability in handling missing modalities is promising. The auhor is encouraged to explain why, in some scenarios, mising modalities + the proposed method surposses the variant with full modalities (Fig 4 c).  I am still wondering its generalizability under more scenarios. 
- 4 The visualization and quanlitative analysis of the learned gaing weight would help the reader better understand the method.
questions: questions: - 1 In table 5, Laplace gaiintg shows limited superiority against other variants, any insight for it? 
- 2 In table 4c, why the method achieves better results with less modaltieis? Any further insight or ablations?
limitations: limitations: N
paper_summary: paper_summary: The Paper introduces “FuseMoE”, a novel mixture-of-experts (MoE) framework that can handle multi-modal data even in scenarios with missing elements, sparsity of samples and temporal irregularity.  It proposes an innovative Laplace gating function with theoretical proof to enhance convergence rates and predictive performance across various tasks. The Laplace gating function allows sparser distribution of weights among MoE and encourages more balanced utilization of experts with sharper peaks and heavier tails preventing certain experts to be over dominant. FuseMoE includes modality and irregularity encoder using a discretized multi-time attention (mTAND) module to support generic fleximodal data ingestion with unlimited input modalities instead of other popular pairwise setup that exist in the literature. Paper also outlines various router designs for processing multimodal inputs and discusses the tradeoff between them and demonstrates .

Key:

1. Laplace gating function showed superior performance compared to Softmax gating across multiple tasks
2. Unlike baseline models, FuseMoE demonstrated improved scalability with increasing number of modalities
3. Handles missing modalities with tAND and per modality routers with entropy loss function
4. theoretical analysis of convergence rates for parameter estimation
5. FuseMoe outperformed baseline methods on several fleximodal datasets.
strength_weakness: strength_weakness: 1. Potential over-parameterization when input is small
2. limited discussion of computational needs and scalability to huge datasets
3. Most of the items discussed were already established and proved. Paper seems to be combining existing literature for these datasets.
questions: questions: 1. The paper claims to work with unlimited input modalities. is it really so? how do you anticipate the fuseMoE to handle extremely sparse data across modalities? Any thoughts/estimation on how far we can go here? How does the complexity of setup scale with number of modalities?
2. Is Laplace only tested with CIFAR-10 dataset? Why was it not tested with other bigger datasets?
3. Have you thought about or explored the interpretability side of expert assignement process? 
4. Few other paper discusses and introduces Laplace instead of gaussian for MoE. Can you expand more on how FuseMoE is differnt from these?
    
   - Nguyen, Hien D., and Geoffrey J. McLachlan. "Laplace mixture of linear experts." *Computational Statistics & Data Analysis* 93 (2016): 177-191.
    
   - Wu, Lc., Zhang, Sy. & Li, Ss. Heteroscedastic Laplace mixture of experts regression models and applications. *Appl. Math. J. Chin. Univ.* **36**, 60–69 (2021). https://doi.org/10.1007/s11766-021-3591-2
limitations: limitations: Potential limitation of over-parameterization of networks especially when input data/size is small is adequately discussed in the paper.
paper_summary: paper_summary: This paper proposes an MOE-based model to handle multimodal data fusion. It addresses two challenges: missing modalities and irregularly sampled data trajectories. A Laplace gating function is applied to the MoE Backbone. An entropy regularization loss is proposed to ensure balanced and stable expert utilization. The author validates the method in diverse datasets.
strength_weakness: strength_weakness: Addressing the following weaknesses may improve the paper: 

1. The author should clearly distinguish the proposed method and others’ modules. From the paper, the Laplace gating function is proposed as a new one. I am not sure the author made some contributions to the encoder design, router design, and loss design. The author should make more illustrations about the contributions, not just combine other people’s work together. 

2. The experimental results are not extensive. The author should clearly demonstrate the data modality of the chosen benchmarks. It seems that other modalities, such as text, are not included. The author should explain this. 

3. Please explain more concisely how the gating functions can stabilize the imbalance and sparse multi-modal data. The paper should point this out more concisely and better with experimental results. Too many mathematical proofs and theorems in Sec.3 seem not helpful in illustrating the advantages of the proposed method.
questions: questions: Please see the weakness part to answer the questions.
limitations: limitations: Yes
paper_summary: paper_summary: The paper “FuseMoE: Mixture-of-Experts Transformers for Fleximodal Fusion” introduces a novel framework called FuseMoE, designed to tackle challenges associated with multimodal data in machine learning. This framework addresses key issues such as missing elements, temporal irregularity, and sparsity in data, which are prevalent in fields like healthcare where data is often incomplete and irregularly sampled. The core innovation in FuseMoE is its gating function that integrates various modalities efficiently. The framework enhances predictive performance through: 1. Sparse MoE Layers: Incorporating sparsely gated MoE layers to manage tasks and learn optimal modality partitioning. 2. Laplace Gating Function: A novel gating function theoretically proven to ensure better convergence rates than traditional Softmax functions. 3. FlexiModal Data Handling: Efficiently handling scenarios with missing modalities and irregularly sampled data trajectories. The framework is validated through empirical evaluations on diverse prediction tasks, demonstrating its effectiveness in real-world applications.
strength_weakness: strength_weakness: 1. The comparison works are all very outdated (the latest one UTDE[1] is publicized in 2023, with over half of the compared works predating 2020). Why not compare with more recent works?
2. Fig.1 lacks clear description, as the data processing workflow is not well explained.
3. Gaussian gating and your proposed form are very similar in format, with formula $h(x)=\operatorname{Top} \mathrm{K}\left(-\|W-x\|_{2}^{2}\right)$. Thus, theoretical comparison with original softmax may be unfair as the form is derived from and similar to gaussian gating. Therefore, although theoretical results may be valid, their importance diminishes due to this aspect.

[1] Improving Medical Predictions by Irregular Multimodal Electronic Health Records Modeling
questions: questions: Please refer to the weaknesses.
The biggest question at my side is that the similarity with gaussian gating may lead to some theoretical misunderstanding. Any further explanation on this topic could lead to higher score.
limitations: limitations: The authors adequately addressed the limitations, there is some redundancy in the encoder part that heavily draws inspiration from UTDE's work. From the appendix it can be observed that encoder might have a greater impact on overall performance than what gate function alone contributes to improvement; hence further discussion is needed on whether this gate function must necessarily be tied to methods related to mTAND and imputation for obtaining higher-quality features as inputs for downstream MoE layers from an implementation perspective. Of course, the discussion on this topic is optional.

Comments: title: Kindly Request for Reviewer's Feedback
Summary After Author-Reviewer Discussion
Thank You!
Official Comment from reviewer  6Ay7
comment: Dear Reviewer `sZSN`,

We sincerely appreciate the time you have taken to provide feedback on our work, which has helped us greatly improve its clarity, among other attributes. This is a gentle reminder that the **Author-Reviewer Discussion period ends in just around 12 hours from this comment, i.e., 11:59 pm AoE on August 13.** We are happy to answer any further questions you may have before then, but we will be unable to respond after that time.

If you agree that our responses to your reviews have addressed the questions you listed, we kindly ask that you consider whether raising your score would more accurately reflect your updated evaluation of our paper. Thank you again for your time and thoughtful comments!

Sincerely,

Authors
Dear Area Chairs and Reviewers,

As the Author-Reviewer Discussion period comes to an end, we would like to provide a brief summary of our paper at this point for your reference. We have carefully reviewed all comments and have provided detailed responses to each. These will be comprehensively addressed in the final version of the paper. We extend our sincere gratitude to all reviewers for their insightful feedback and suggestions.

## **Motivation**: 
We introduced FuseMoE, a novel MoE framework specifically designed to enhance multimodal fusion for inputs with arbitrary missingness or irregularity. Multiple reviewers (`P6x6`, `6SUH`, `sZSN`) have **agreed that the problem we are addressing is well-motivated and represents a promising direction for investigation**.

## **Method**: 
FuseMoE incorporates sparsely gated MoE layers in its fusion component, which routes each modality to designated experts specializing in those specific data types while mitigating the influence of missing modalities. **This capability is critical to our method, as it allows for leveraging a vast amount of patient data without relying solely on complete modalities.** As we explained to Reviewer `P6x6`, most other methods that only leverage full modality data result in the exclusion of patient records with incomplete modalities. FuseMoE, by contrast, greatly improves generalization by accessing a broader array of samples.

Additionally, FuseMoE employs various router design mechanisms tailored for multimodal inputs with distinct combinations or interactions, giving users the flexibility to choose based on their specific applications.

## **Theory**: 
FuseMoE integrates a novel Laplace gating function, which is particularly well-suited for multimodal input data, especially in large-scale experiments, compared to the commonly used Softmax gating function. Theoretically, it ensures better convergence rates than Softmax, as discussed in detail with high-level intuitive explanations in our rebuttal to Reviewer `6SUH`. Additionally, we have illustrated the superiority of our gating function over the Gaussian gating function in our response to Reviewer `sZSN`.

## **Experiments**: 
We demonstrate that FuseMoE significantly outperforms existing well-known and representative methodologies within specific categories, in integrating diverse input modality types:
- `Text`: MIMIC-III, MIMIC-IV, CMU-MOSI, CMU-MOSEI
- `Images`: MIMIC-IV, CIFAR-10, ImageNet (new)
- `Time series`: MIMIC-III, MIMIC-IV, PAM
- `ECG waveform`: MIMIC-IV
- `Video frames`: CMU-MOSI, CMU-MOSEI
- `Audio signals`: CMU-MOSI, CMU-MOSEI

with varying missingness and irregular sampling on challenging tasks:

- `Mortality prediction`: MIMIC-III, MIMIC-IV
- `Length-of-stay`: MIMIC-III, MIMIC-IV
- `Phenotyping`: MIMIC-IV
- `Sentiment analysis`: CMU-MOSI, CMU-MOSEI
- `Activity recognition`: PAM
- `Image classification`: CIFAR-10, ImageNet (new)

**The majority of the multimodal datasets we tested are considered large-scale**, including MIMIC-III with 36,212 samples, MIMIC-IV with 73,173 samples, and CMU-MOSEI with 22,777 samples. These datasets cover a sufficient number of benchmarks and research areas, particularly in a context where benchmarks are not abundant.

**We also conducted comprehensive ablation studies and visualizations of FuseMoE** to assess the effectiveness of each building block and enhance the interpretability of our proposed method. These studies revealed that our model architecture consistently outperforms other alternatives and insights into how much each input modality contributes to our tasks of interest.

## **Reproducibility**: 
We have submitted (1) implementations of FuseMoE and its baselines, (2) scripts to process and generate multimodal input data as part of our supplementary material.

Thank you again for your efforts in reviewing this submission. We appreciate all the valuable feedback, which has significantly contributed to improving our manuscript. **We look forward to continuing the discussion with any reviewers who still have questions during the remaining discussion period.**

Sincerely,

Authors
Dear Reviewer `6Ay7`,

Thank you so much for your positive response! Please feel free to reach out if you have any further questions or thoughts. We would be delighted to discuss them with you to improve our paper!

Sincerely,

Authors
I thank the authors for the rebuttal and for taking the time to conduct additional experiments and comparisons.

The rebuttal and the additional clarification address most of my concerns so I am happy to increase my score from 6 to 7.
Cheers!