# Revisiting Incremental Stochastic Majorization-Minimization Algorithms with Applications to Mixture of Experts

Trungkhang Tran, Trungtin Nguyen, Gersende Fort, Tung Doan, Hien Duy Nguyen, Binh T. Nguyen, Florence Forbes, Christopher C. Drovandi

► To cite this version:

Trungkhang Tran, Trungtin Nguyen, Gersende Fort, Tung Doan, Hien Duy Nguyen, et al.. Revisiting Incremental Stochastic Majorization-Minimization Algorithms with Applications to Mixture of Experts. 2026. hal-05497100

HAL Id: hal-05497100

https://hal.science/hal-05497100v1

Preprint submitted on 6 Feb 2026

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Revisiting Incremental Stochastic
Majorization-Minimization Algorithms with
Applications to Mixture of Experts

TrungKhang Tran*,<fnref n="1" />, TrungTin Nguyen*,<fnref n="2" />,<fnref n="3" />, Gersende Fort<fnref n="4" />, Tung Doan<fnref n="5" />,

Hien Duy Nguyen<sup>6,7</sup>, Binh T. Nguyen<sup>8,9</sup>, Florence Forbes<sup>10</sup>, Christopher Drovandi<sup>2,3</sup>

Abstract

Processing high-volume, streaming data is increasingly common in modern statis-
tics and machine learning, where batch-mode algorithms are often impractic-
cal because they require repeated passes over the full dataset. This has moti-
vated incremental stochastic estimation methods, including incremental stochastic
Expectation-Maximization (EM) procedures formulated via stochastic approxima-
tion. In this work, we study an incremental stochastic variant of the Majorization-
Minimization (MM) principle that generalizes incremental stochastic EM as a
special case. Our main methodological contribution is to make this framework
operational for softmax-gated mixture of experts (MoE) models: we construct
tractable majorizer surrogates that yield explicit incremental updates in settings
where no stochastic EM algorithm is available. We complement this development
with consistency guarantees by verifying conditions under which the iterates con-
verge to a stationary point characterized by a vanishing gradient of the objective.
Empirically, for softmax-gated MoE regression, the resulting incremental stochas-
tic MM algorithm consistently outperforms widely used stochastic optimizers,
including stochastic gradient descent, root mean square propagation, adaptive mo-
ment estimation, and second-order clipped stochastic optimization. These results
highlight the practical value of building a principled bridge between incremental
MM theory and modern softmax-gated MoE architectures, given their central role
in contemporary deep neural networks for heterogeneous data modeling and scal-
able conditional computation. Beyond synthetic experiments, we further validate
effectiveness on two real-world datasets, including a bioinformatics study of dent
maize genotypes under drought stress that integrates high-dimensional proteomics
with ecophysiological traits, where incremental stochastic MM yields stable gains
in predictive performance.

<footnote n="1">¹School of Computing, National University of Singapore, Singapore.</footnote>

<footnote n="2">²ARC Centre of Excellence for the Mathematical Analysis of Cellular Systems.</footnote>

<footnote n="3">³School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia.</footnote>

<footnote n="4">⁴Laboratoire d'Analyse et d'Architecture des Systèmes, CNRS, Toulouse, France.</footnote>

<footnote n="5">⁵School of Medicine and Dentistry, Griffith University, Brisbane, Australia.</footnote>

<footnote n="6">⁶Department of Mathematics and Physical Sciences, La Trobe University, Melbourne Australia.</footnote>

<footnote n="7">⁷Institute of Mathematics for Industry, Kyushu University, Fukuoka, Japan.</footnote>

<footnote n="8">⁸Faculty of Mathematics and Computer Science, University of Science, Ho Chi Minh City, Vietnam.</footnote>

<footnote n="9">⁹Vietnam National University, Ho Chi Minh City, Vietnam.</footnote>

<footnote n="10">¹⁰Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000, Grenoble, France.</footnote>

<footnote>*Co-first author.</footnote>
# 1 Introduction

Mixture of Experts (MoE) models, originally introduced by [44, 49], are widely used to capture complex, heterogeneous relationships between inputs and outputs by combining a gating model with multiple expert models, all driven by the covariates. As a prominent example of conditional computation [2, 16], MoE architectures route each input to a small subset of experts, which can substantially increase model capacity while keeping training and inference costs nearly constant by activating only a fraction of parameters per example [97, 27, 21]. Beyond their practical success, MoE models are supported by strong approximation, estimation and model selection theory across mixture models [35, 91, 86, 41, 42, 85, 83, 17], mixture of regression models [26, 43], and more general MoE families [47, 87, 79, 76, 80, 84, 73, 81, 72, 74, 59, 98, 36]. This growing theory also aligns with recent methodological developments on predictive inference for high-dimensional MoE models, which couple $\ell_1$-penalized EM estimation with debiasing to construct prediction sets with conditional coverage guarantees [45]. Beyond uncertainty quantification, these advances also highlight the need for computationally efficient inference routines for modern MoE pipelines. We refer to [101, 69, 75, 82, 7] for broader surveys and applications.

At the same time, modern statistical and machine-learning workflows increasingly operate in streaming or large-scale regimes, where batch-mode methods are often impractical because they require repeated passes over the full dataset. This motivates incremental and mini-batch estimation procedures that update parameters using sequentially arriving data. Many incremental and mini-batch extensions of EM can be interpreted within the classical Stochastic Approximation framework (see, e.g., [6, 55, 25, 31, 100]) and the online-EM line initiated in [12, 11, 29, 32, 78, 54, 51, 20, 50, 88, 89, 15, 33, 30]. In parallel, incremental stochastic and mini-batch adaptations of MM algorithms have primarily been developed using techniques from incremental stochastic learning and convex optimization (see, e.g., [18, 38, 46, 56, 64, 92, 96]), with illustrative examples outlined in [68, 71, 77, 19, 65, 66, 62, 60, 61, 28, 24]. While these online EM and stochastic MM methodologies are effective for many mixture and latent-variable settings, their direct applicability to softmax-gated MoE models is far less clear.

The goal of this work is to establish a principled and practically useful bridge between incremental stochastic MM methodology and modern softmax-gated MoE architectures. Concretely, we revisit the stochastic-approximation formulation of incremental stochastic MM studied in [77] (see also Algorithm 1 in [24]) and develop an incremental stable MM scheme tailored to softmax-gated MoE models. Our analysis imposes regularity assumptions on the surrogate construction together with i.i.d. assumptions on the data stream (see P 1 to P 4). In contrast to EM-based approaches, the incremental stochastic MM framework accommodates general surrogate functions and therefore covers incremental stochastic EM as a special case when a latent-variable exponential-family representation is available [12, 11]. This flexibility is particularly relevant for softmax-gated MoE models and related objectives where a stochastic EM algorithm is not available or is difficult to derive.

**Main contributions and practical implications (informal).** The novelty of the paper lies primarily in the MoE-specific instantiation, which builds a practical bridge between incremental stochastic MM theory and modern softmax-gated MoE training. Concretely:

* We construct explicit exponential-family majorizer surrogates for softmax-gated MoE models and derive concrete incremental stochastic MM updates for both the softmax-gated Gaussian MoE and softmax-gated multinomial logistic MoE settings; see Section 4.

* We establish stability and consistency statements for these MoE instantiations by verifying the required assumptions and invoking stochastic-approximation convergence theory; see Section 3 together with Theorems 2 and 3.

* We explain why existing online EM and incremental stochastic MM variants do not apply off-the-shelf to softmax-gated MoE models, and we make this obstruction explicit by analyzing the failure of standard global regularity conditions (Lipschitz continuity, $L$-smoothness, convexity) for the SGMoE objective; see Section 4.3 and Section 4.3.2.
* We prove consistency under verifiable conditions: assuming stability with probability one, the iterates converge to stationary points (vanishing gradient); see Section 3. Comparable guarantees are typically unavailable for widely used online optimizers for MoE models, including SGD [94], RMSProp [40], Adam [53, 23], and Sophia [63].

* Our incremental stochastic MM algorithms empirically outperform these optimizers on synthetic benchmarks and on real data, including a dent maize drought-stress study combining high-dimensional proteomics with ecophysiological traits; see Section 5.

* We revisit and correct the bound in Lemma 2 of [4], and we leverage the corrected form to build a valid majorizer for the softmax-gated MoE optimization problem; this refinement may be of independent interest for MM constructions in multinomial logistic regression.

Overall, the contribution is not merely an incremental algorithmic variant: it provides a concrete route for bringing incremental stochastic MM methodology to softmax-gated MoE models, together with a clear account of why several seemingly related online EM/MM approaches do not extend to this setting.

**Notation.** Throughout the paper, the set {1, 2, ..., N} is abbreviated as [N] for any positive integer N ∈ N; 1_N and I_N denote the N-dimensional vector of ones and N × N identity matrix, respectively. In accordance with standard conventions, A^T represents the transpose of a given vector or matrix A, and all vectors v = (v_1, ..., v_N)^T ≡ [v_1, ..., v_N]^T are expressed as column vectors. The Loewner order between two matrices is denoted by A ≥ B, implying that the matrix A - B is positive semi-definite (PSD). Additionally, the symbol A ⊗ B refers to the Kronecker product (also known as the matrix direct product) of the matrices A and B. Let ⟨·, ·⟩ denote the Euclidean scalar product. Given any differentiable function f = (f_1, ..., f_M)^T : R^T → R^M, we denote by ∇_θ f = (∂_f / ∂_θ_1, ..., ∂_f / ∂_θ_T)^T the gradient of f when M = 1, and by ∇_θ f^T the T × M matrix whose columns are the gradients, that is, ∇_θ f^T = (∇_θ f_1, ..., ∇_θ f_M). When the variable of differentiation is clear from the context, we may omit the subscript and simply write ∇_θ f instead of ∇_θ f. Accordingly, the symbol ∇_θ f^T is to be understood as either a vector or a matrix, depending on whether the function f is scalar or vector-valued. In the latter case, the usual Jacobian matrix is the transpose of ∇_θ f^T. When f is twice differentiable, we use the notation ∇_θ^2 f as the Hessian matrix which is a T × T matrix with components given by [∇_θ^2 f]_{i,j} = ∂^2 f / ∂_θ_i ∂_θ_j, i, j ∈ [T]. For any r ∈ R^D and any set S ⊂ R^D, we define d(r, S) := inf_{s∈S} ||r - s||, that is, the distance from s to S, where ||·|| denotes the L_2-norm on R^D. Finally, for any d ∈ N and any vector x = (x_1, x_2, ..., x_D)^T ∈ R^D, we define the element-wise d-th power of x by x^d := (x_1^d, x_2^d, ..., x_D^d)^T.

**Paper organization.** The remainder of this paper is organized as follows. In Section 2, we provide an overview of the incremental stochastic MM algorithm, setting the foundation for subsequent discussions. In Section 3, we provide sufficient conditions for the consistency of the incremental stochastic MM algorithm, which are further substantiated through empirical evaluations. Applications to softmax-gated MoE models are presented in Section 4, followed by an experimental study in Section 5 to validate the theoretical findings. Then, we conclude the article and outline directions for future research in Section 6. Finally, all proofs of theoretical results are deferred to the supplementary material in Sections A to C.

# 2 The incremental stochastic MM algorithm

## 2.1 Primitive optimization problem

We begin by examining the following primitive optimization problem in a statistical framework:

$$ \arg\min_{\theta \in \mathbb{T}} \mathbb{E}_{z \sim \mathcal{P}_{\pi}} [f(\theta; z)] = \arg\min_{\theta \in \mathbb{T}} \mathbb{E}_{\pi} [f(\theta; z)]. \quad (1) $$

In this context, the function f: T × Z → R is measurable and represents a family indexed by the parameter θ ∈ T. The set T is defined as a measurable open subset of R^T, and Z ⊆ R^K is a Euclidean subspace endowed with its Borel sigma-algebra. Additionally, z is an Z-valued random variable on the probability space (Ω, A, P), with its probability density function on Z denoted by
$\pi$. The symbols $\mathbb{P}_\pi$ and $\mathbb{E}_\pi$ are used to represent the probability measure and the expected value, respectively, under $\pi$.

In this paper, we focus on the scenario where the expectation $\mathbb{E}_\pi [|f(\theta; z)|]$ is assumed to be finite on $\mathbb{T}$ but cannot be expressed in closed form, and the corresponding optimization problem is addressed using an MM-based algorithm. According to the terminology introduced in [58], the function $g: (\theta, z, \tau) \mapsto g(\theta, z; \tau)$, defined over the domain $\mathbb{T} \times \mathbb{Z} \times \mathbb{T}$, is referred to as a majorizer of $f$ if, for any $\tau \in \mathbb{T}$ and any $(\theta, z) \in \mathbb{T} \times \mathbb{Z}$, the following conditions are satisfied:

$$f(\theta; z) - f(\tau; z) \le g(\theta, z; \tau) - g(\tau, z; \tau). \qquad (2)$$

In this paper, we study the scenario in which the majorizer function $g$ satisfies the following conditions:

**P 1 (Exponential family majorizer surrogate).** The majorizer surrogate $g$ belongs to an exponential family:

$$g(\theta, z; \tau) := -\psi(\theta) + \langle \bar{S}(\tau; z), \phi(\theta) \rangle, \qquad (3)$$

where $\psi: \mathbb{T} \to \mathbb{R}$, $\phi: \mathbb{T} \to \mathbb{R}^D$ and $\bar{S}: \mathbb{T} \times \mathbb{Z} \to \mathbb{R}^D$ are measurable functions. In addition, $\phi$ and $\psi$ are continuously differentiable on $\mathbb{T}$.

**P 2 (Convex set).** There exists a measurable, open, and convex set $\mathbb{S} \subseteq \mathbb{R}^D$ such that

$$\bar{S}(\tau; z) \in \mathbb{S} \quad \text{for all } (\tau, z) \in \mathbb{T} \times \mathbb{Z}.$$

**P 3 (Incrementally independent and identically distributed data with finite expectation).** The expectation $\mathbb{E}_\pi[\bar{S}(\theta; z)]$ exists, lies within $\mathbb{S}$, for any $\theta \in \mathbb{T}$, but it not have a closed-form expression.

Nevertheless, an incremental sequence of independent examples $\{z_n, n \ge 0\}$, distributed under $\pi$, is available. We assume that $\{z_n, n \ge 0\}$ are random variables on $(\Omega, \mathcal{A}, \mathbb{P})$.

$$\mathbb{E}_{\pi}[f(\theta; z)] - \mathbb{E}_{\pi}[f(\tau; z)] \le \psi(\tau) - \psi(\theta) + \langle \mathbb{E}_{\pi}[\bar{S}(\tau; z)], \phi(\theta) - \phi(\tau) \rangle, \quad (4)$$

thereby we establish a valid majorizer function for the variation of the objective function $\theta \mapsto \mathbb{E}_\pi[f(\theta; z)] - \mathbb{E}_\pi[f(\tau; z)]$ as follows: $G(\cdot, \tau) - G(\tau, \tau)$ where $G(\cdot, \tau) := -\psi(\cdot) + \langle \mathbb{E}_\pi[\bar{S}(\tau; z)], \phi(\cdot) \rangle$.

We also assume that the majorizing function $\theta \mapsto g(\theta, z; \tau)$ possesses a unique minimizer.

**P 4 (Unique global minimum).** For any $s \in \mathbb{S}$, there exists a unique root in $\mathbb{T}$ to the function $\theta \mapsto -\nabla\psi(\theta) + \nabla\phi(\theta)^\top s$. This root corresponds to the unique global minimum of the function $\theta \mapsto h(s; \theta) := -\psi(\theta) + \langle s, \phi(\theta) \rangle$ over $\mathbb{T}$. The root is denoted by $\bar{\theta}(s)$, where

$$\bar{\theta}(s) := \arg\min_{\theta \in \mathbb{T}} [-\psi(\theta) + \langle s, \phi(\theta) \rangle], \quad -\nabla\psi(\bar{\theta}(s)) + \nabla\phi(\bar{\theta}(s))^\top s = 0. \qquad (5)$$

A regularity assumption is also assumed on the objective function.

**P 5.** The function $\theta \mapsto \mathbb{E}_\pi[f(\theta; z)]$ is continuously differentiable on $\mathbb{T}$.

## 2.2 Incremental MM algorithm

Given that the expectation in Equation (1) may lack a closed form, but infinitely large datasets $\{z_n, n \ge 0\}$ are available (as stated in P 3), we consider the Incremental (Online) MM algorithm introduced in [77] (see also [24], which is more general as it does not assume a specific form of the oracle of $\mathbb{E}_\pi[\bar{S}(\theta_n; z)]$) and is described in Algorithm 1.

This algorithm defines a sequence $\{s_n, n \ge 0\}$, where the update mechanism, detailed in Equation (6), is a Stochastic Approximation iteration. The approach involves constructing a sequence of majorizer functions by defining their parameter $s_n$ within the functional space spanned by $-\psi, \phi_1, \dots, \phi_D$.
**Algorithm 1 Incremental MM algorithm**

**Require:** Incremental data points {$z_n : n \ge 0$}, an initial value $s_0 \in \mathbb{S}$ and positive step sizes {$\gamma_{n+1}, n \ge 0$} in (0, 1).

**Ensure:** A $\mathbb{T}$-valued sequence {$\theta_n, n \ge 1$} and a $\mathbb{S}$-valued sequence {$s_n, n \ge 0$}.

1: **for** $n = 0, 1, \dots, \mathbf{do}$

2: Compute

$$s_{n+1} = s_n + \gamma_{n+1} \{ \bar{S}(\bar{\theta}(s_n); z_{n+1}) - s_n \}, \quad (6)$$

$$\theta_{n+1} = \bar{\theta}(s_{n+1}) := \arg\min_{\theta \in \mathbb{T}} [-\psi(\theta) + \langle s_{n+1}, \phi(\theta) \rangle]. \quad (7)$$

3: **end for**

Under P 2, it is readily seen that since $s_0 \in \mathbb{S}$ and $\gamma_n \in (0, 1)$

$$\mathbb{P}(\forall n \ge 0, s_n \in \mathbb{S}) = 1; \quad (8)$$

with probability one, the algorithm generates an $\mathbb{S}$-valued sequence.

## 3 Convergence properties of the incremental stochastic MM algorithm

In this section, we investigate the convergence properties of Algorithm 1. Throughout this section, it is assumed that the MM framework described by Equation (2) and P 1 to P 5 holds.

### 3.1 Limiting points of Algorithm 1 and stationary points of the original optimization problem

Equation (6) in Algorithm 1 is a Stochastic Approximation algorithm [94]. Such iterative schemes are designed to solve root-finding (fixed-point) problems, that is, to find $s^0$ such that $E_\pi [\bar{S}(\bar{\theta}(s^0); z)] = s^0$. Thus, Algorithm 1 is formulated to approximate the intractable expectation at $\bar{\theta}(s^0)$, where $s^0$ adheres to this fixed-point condition.

Proposition 1 was initially proposed and proved as Lemma 1 in [77]. For the sake of completeness and to maintain consistency with the notations used in our paper, we have also included the proof of Proposition 1 in Section A.1. Proposition 1 elucidates the connection between the limiting points of Equation (6) and the optimization problem defined in Equation (1). Notably, it demonstrates that any limiting value $s^0$ corresponds to a stationary point $\theta^0 := \bar{\theta}(s^0)$ of the objective function $E_\pi [f(\theta; z)]$, indicating that $\theta^0$ is a root of its derivative. The proof leverages the methodology presented in Proposition 1 of [12], which has also been used in Lemma 2 of [22], Proposition 1 of [29], Proposition 2.1 of [31], and the contemporaneous Proposition 1 of [24], to establish convergence rates for stochastic EM-type algorithms, including the incremental stochastic EM algorithm studied here. We commence by defining the following notations:

$$\eta(s) := E_\pi [\bar{S}(\bar{\theta}(s); z)] - s, \quad F := \{s \in \mathbb{S} : \eta(s) = 0\}. \quad (9)$$

**Proposition 1.** Let $\mathbb{L}$ denote the set of stationary points of this function, defined as

$$\mathbb{L} := \{\theta \in \mathbb{T} : \nabla_{\theta} E_{\pi}[f(\theta; z)] = 0\}.$$

If $s^0 \in \mathbb{F}$, then $\bar{\theta}(s^0) \in \mathbb{L}$. Conversely, if $\theta^0 \in \mathbb{L}$, it follows that $s^0 := E_\pi [\bar{S}(\theta^0; z)] \in \mathbb{F}$.

**Additional regularity assumptions for Lyapunov function.** Leveraging the results from [22] on the asymptotic convergence of stochastic approximation algorithms, along with the additional regularity conditions specified in Assumptions 1 to 3 for $\psi$, $\phi$, and $\bar{\theta}$, Proposition 2 (whose proof is detailed in Section A.2) establishes that the algorithm described in Equation (6) admits a continuously differentiable Lyapunov function $V$, defined on $\mathbb{S}$ as

$$s \mapsto V(s) := E_\pi [f(\bar{\theta}(s); z)],$$

which satisfies $\langle \nabla V(s), \eta(s) \rangle \le 0$, with strict inequality holding outside the set $\mathbb{F}$ (refer to [12, Prop. 2]). Furthermore, the additional regularity assumptions detailed below must also be satisfied, in addition to P 1 to 5.
**Assumption 1** (Twice continuously differentiable exponential family surrogate). The parameter space $\mathcal{T}$ is defined as a convex and open subset of $\mathbb{R}^T$, with the functions $\psi$ and $\phi$ in Equation (3) assumed to be twice continuously differentiable on $\mathcal{T}$.

**Assumption 2** (Continuously differentiable global minimum function). The global minimizer function $s \mapsto \theta(s)$, defined in Equation (5), is continuously differentiable on $\mathcal{S}$.

**Proposition 2.** Under Assumption 1 and Assumption 2, the following holds:

(a) The function $V$ is continuously differentiable on $\mathcal{S}$.

(b) For any $s \in \mathcal{S}$, we have $\langle \nabla V(s), \eta(s) \rangle \le 0$, and $\{s \in \mathcal{S} : \langle \nabla V(s), \eta(s) \rangle = 0\} \supseteq \mathbb{F}$.

## 3.2 Consistency of the incremental stochastic MM algorithm

We now establish the almost-sure convergence of the sequence $\{s_n : n \ge 0\}$ stated in Theorem 1. The detailed proof is provided in Section A.3, where we verify the sufficient conditions for convergence of a general stochastic approximation algorithm as given in [1]. Moreover, the sequence $\{\theta(s_n), n \ge 0\}$ is also proven to converge almost surely to the set $\mathbb{L}$, which corresponds to the stationary points of the objective function $\theta \mapsto \mathbb{E}_{\pi}[f(\theta; z)]$.

To this end, we require the following assumptions on moment conditions on the statistics $\bar{S}$, on the model at hand, on the stability of the algorithm, and on the step sizes:

**Assumption 3** (Second finite moment). For all compact subsets $\mathbb{K} \subset \mathcal{S}$,

$$ \sup_{s \in \mathbb{K}} \mathbb{E}_{\pi} [\| \bar{S} (\bar{\theta}(s); z) \|^2] < \infty. $$

**Assumption 4.** The function $s \mapsto \eta(s)$ is continuous on $\mathcal{S}$. In addition, $\{s \in \mathcal{S} : \langle \nabla V(s), \eta(s) \rangle = 0\} \subset \mathbb{F}$.

Combined with Proposition 2, Assumption 4 implies that the two sets are equal and are closed, since the set $\mathbb{F}$ collects the zeros of a continuous function. It is established in the proof of Proposition 2 (see Section A.2) that

$$ \langle \nabla V(s), \eta(s) \rangle = - \left( \{ (\nabla \phi)(\bar{\theta}(s)) \}^{\top} \eta(s) \right)^{\top} A(s) \left( \{ (\nabla \phi)(\bar{\theta}(s)) \}^{\top} \eta(s) \right), $$

where $A(s)$ is a positive definite matrix under P 4. Consequently, $\langle \nabla V(s), \eta(s) \rangle = 0$ iff $\{(\nabla \phi)(\bar{\theta}(s))\}^{\top} \eta(s) = 0$, so that if $\{(\nabla \phi)(\bar{\theta}(s))\}$ has full row rank, then $\eta(s) = 0$. Thus, this rank condition is a sufficient condition for the two sets $\{s : \langle \nabla V(s), \eta(s) \rangle = 0\}$ and $\mathbb{F}$ to be equal.

**Assumption 5.** The closure of $V(\mathbb{F})$ has an empty interior.

**Assumption 6** (Sublevel sets of the Lyapunov function). There exist $M_0 < M_1$ such that $\mathbb{F} \subset \{s \in \mathcal{S} : V(s) < M_0\}$ and $\{s \in \mathcal{S} : V(s) \le M_1\}$ is compact.

**Assumption 7** (Stability assumption). With probability one, there exists a compact set $\mathbb{K}$ of $\mathcal{S}$ such that $\mathbb{K} \cap \mathbb{F} \ne \emptyset$ and $s_n \in \mathbb{K}$ for all $n \ge 0$.

**Assumption 7** states that, with probability one, $\limsup_n \|s_n\| < \infty$ and $\liminf_n d(s_n, \mathcal{S}^c) > 0$. It imposes a stability assumption, which is not straightforward. In general, the stability of the algorithm can be guaranteed by truncating the algorithm updates, either within a fixed set, as outlined in [55, Chapter 2], or within a progressively expanding sequence of sets, as demonstrated in [1]. To maintain conciseness in the exposition, we refrain from explicitly performing these constructions, which also impact the convergence analysis, and leave them as directions for future work.

**Assumption 8.** The step size sequence is a monotone nonincreasing sequence, $\gamma_n \in (0, 1)$, $\sum_n \gamma_n = +\infty$ and $\sum_n \gamma_n^2 < \infty$.

**Assumption 8** represents a conventional requirement for stochastic approximation methods employing vanishing step sizes, as highlighted in [55]. This condition is satisfied, for instance, by selecting the step size $\gamma_n = \gamma_0 n^{-\alpha}$, with $\alpha$ belonging to the interval $(\frac{1}{2}, 1]$. The supplementary conditions that $\gamma_n \in (0, 1)$ and that $s_0 \in \mathcal{S}$, together with P 2, ensure that the entire sequence $\{s_n\}_{n \ge 0}$ remains in $\mathcal{S}$ for all $n$.
**Theorem 1 (Consistency).** Let $\mathbb{F}$ and $\mathbb{L}$ be given by Equation (9) and Proposition 1. Under Assumption 1 to Assumption 8, with probability 1,

$$ \limsup_{n \to \infty} d(\mathbf{s}_n, \mathbb{K} \cap \mathbb{F}) = 0 \quad \text{and} \quad \limsup_{n \to \infty} d(\boldsymbol{\theta}_n, \mathbb{L}) = 0. \qquad (10) $$

## 4 An application to the softmax-gated MoE models

To demonstrate the effectiveness of the proposed approach, we consider a regression model that, as outlined in Section 2.1, represents a scenario where the expectation $\mathbb{E}_{\pi}[f(\boldsymbol{\theta}; \mathbf{z})]$ in Equation (1) can not be explicitly computed in closed form. The associated optimization problem is therefore tackled using the MM-framework-based Algorithm 1 in Section 2.2.

**Softmax-gated MoE models for heterogeneous data.** We refer to the outputs $\mathbf{y} \in \mathbb{Y}$, as the target or response variables and the inputs $\mathbf{x} \in \mathbb{X} \subset \mathbb{R}^P$, $P \in \mathbb{N}$, as the fixed explanatory or predictor variables. We consider an incremental sequence $\{(\mathbf{x}_n, \mathbf{y}_n) : n \ge 1\}$ of i.i.d. copies of $(\mathbf{x}, \mathbf{y}) =: \mathbf{z}$, which follows a distribution characterized by the true but unknown probability density function $\pi(\cdot | \mathbf{x} = \mathbf{x})$. Here, the corresponding observed values are denoted by $(\mathbf{x}, \mathbf{y}) =: \mathbf{z}$. Motivated by universal approximation theorems for MoE models, the function $\pi$ can be estimated using a softmax-gated MoE model, defined as:

$$ s_{\theta}(\mathbf{y}|\mathbf{x}) = \sum_{k=1}^{K} g_k(\mathbf{w}(\mathbf{x}))e_k(\mathbf{y}; \mathbf{v}(\mathbf{x})). \qquad (11) $$

The softmax gating networks are defined as follows:

$$ g_k(\mathbf{w}(\mathbf{x})) = \frac{\exp(w_k(\mathbf{x}))}{\sum_{l=1}^{K} \exp(w_l(\mathbf{x}))}, \quad \forall k \in [K]. $$

Here, $\mathbf{w}(\mathbf{x}) = (w_1(\mathbf{x}), \dots, w_K(\mathbf{x}))$ represents the weight functions, and the expert networks $e_k(\mathbf{y}; \mathbf{v}(\mathbf{x}))$, for all $k \in [K]$, are specified later in Sections 4.1 and 4.2 for continuous and discrete output variables respectively, and $K$ is the number of mixture of expert components. In Equation (11), $\boldsymbol{\theta}$ collects all finite-dimensional parameters that specify the gating and expert networks, i.e., the weight functions $\{w_k(\cdot)\}_{k \in [K]}$ and the expert predictors $\{\mathbf{v}(\cdot)\}$ (and hence $\{e_k(\cdot; \mathbf{v}(\cdot))\}$) are assumed to belong to a prescribed parametric family indexed by $\boldsymbol{\theta}$.

**From batch MM to incremental stochastic MM for softmax-gated MoE models.** Since the true conditional density $\pi(\cdot | \mathbf{x} = \mathbf{x})$ is typically unknown, we formulate parameter learning for the softmax-gated MoE model in Equation (11) as an expected-risk minimization problem of the form Equation (1). We consider

$$ \boldsymbol{\theta}^0 \in \arg \min_{\boldsymbol{\theta} \in \mathbb{T}} \mathbb{E}[f(\boldsymbol{\theta}; \mathbf{z})], \quad f(\boldsymbol{\theta}; \mathbf{z}) := \ell(\mathbf{y}, s_{\theta}(\cdot | \mathbf{x})), \qquad (12) $$

where $s_{\theta}(\mathbf{y} | \mathbf{x}) = \sum_{k=1}^{K} g_k(\mathbf{w}(\mathbf{x})) e_k(\mathbf{y}; \mathbf{v}(\mathbf{x}))$ is defined in Equation (11), and $\ell(\cdot, \cdot)$ is a chosen loss (e.g., a proper scoring rule). Here, we consider the negative log predictive density, i.e.,

$$ \ell(\mathbf{y}, s_{\theta}(\cdot | \mathbf{x})) = -\log s_{\theta}(\mathbf{y} | \mathbf{x}). \qquad (13) $$

The expectation in Equation (12) is taken with respect to the joint law of $\mathbf{z} = (\mathbf{x}, \mathbf{y})$: equivalently, $\mathbb{E}[f(\boldsymbol{\theta}; \mathbf{z})] = \mathbb{E}_{\mathbf{x}}[\mathbb{E}_{\mathbf{y}|\mathbf{x}}[-\log s_{\theta}(\mathbf{y} | \mathbf{x}) | \mathbf{x}]]$, so $\boldsymbol{\theta}^0$ is a single global parameter (it does not depend on $\mathbf{x}$). In general, the expectation $\mathbb{E}[f(\boldsymbol{\theta}; \mathbf{z})]$ is not available in closed form. A direct batch strategy would replace $\mathbb{E}[f(\boldsymbol{\theta}; \mathbf{z})]$ by an empirical average based on $\{\mathbf{z}_n\}_{1 \le n \le N}$ and then apply an MM procedure to this finite-sum objective. However, such a batch approach implicitly assumes the availability of the entire dataset $\{\mathbf{z}_n\}_{1 \le n \le N}$ for each $N$, which is infeasible in practice when observations are accessed as a stream and when memory constraints make storing the full sequence impossible for large $N$. Accordingly, we work under P 3, where an incremental stream of i.i.d. observations $\{\mathbf{z}_n = (\mathbf{x}_n, \mathbf{y}_n)\}_{n \ge 0}$ is available. We therefore apply the incremental (online) MM scheme in Algorithm 1. This algorithm maintains an auxiliary statistic sequence $\{s_n\}_{n \ge 0}$ updated by the stochastic approximation recursion Equation (6), which constructs data-driven majorizer surrogates within the functional class spanned by $-\psi, \phi_1, \dots, \phi_D$ (see P 1). The parameter iterate is then obtained via the deterministic map $\boldsymbol{\theta}_{n+1} = \tilde{\boldsymbol{\theta}}(s_{n+1})$ in Equation (7).
## 4.1 Softmax-gated Gaussian MoE (SGMoE) models for continuous output data

In this case, the response variable $y$ is continuous, i.e., $y \in \mathbb{Y} \subset \mathbb{R}^Q$, $Q \in \mathbb{N}$, and the Gaussian expert, $e_k(y; v(x)) := \mathcal{N}(y, \mu_k(x), \Sigma_k)$, is parameterized by the mean function $\mu_k(x)$ and the covariance matrix $\Sigma_k$, i.e.,

$$e_k(y; v(x)) := \mathcal{N}(y; \mu_k(x), \Sigma_k) = \frac{\exp\left(-\frac{1}{2}(y - \mu_k(x))^\top \Sigma_k^{-1} (y - \mu_k(x))\right)}{\sqrt{(2\pi)^Q |\Sigma_k|}}$$

For simplicity, we assume that $\Sigma_k$ is diagonal, with diagonal entries $\sigma_{k,q}^2$ for $q \in [Q]$. We leave the extension to a general covariance matrix for future work.

### 4.1.1 Motivation for polynomial regression

To address the heterogeneous regression problem, certain authors have employed softmax-gated Gaussian MoE models under specific simplifying assumptions. Notably, [13, 14] explored Gaussian MoE models for multiple regression with univariate output variables, where both the weights and means are modeled as linear functions of the input variables. However, this linear assumption restricts the capacity of MoE models. Indeed, in the context of convolutional neural networks, [16] empirically demonstrated that while mixtures of linear experts outperform single-expert models, they fall significantly short when compared to mixtures of non-linear experts. Motivated by these findings, we aim to enhance the flexibility of MoE models by incorporating nonlinearities. Specifically, we define the gating weights and expert means as polynomials of the input variables. Regarding the convergence properties of such polynomial-based MoE models, we refer to [70], which provides insights into the optimal convergence rates of MoE models where each expert employs a polynomial regression framework. Motivated by regression with polynomial features, both the gating weights $w_k$ and the expert means $\mu_k$ are modeled as polynomial functions of the input variables $x$, as follows: given any $\omega_{k,d} = (\omega_{k,d,p})_{p \in [P]} \in \mathbb{R}^P$,

$$w_k(x) := w_k(x; \omega_k) = \sum_{d=0}^{D_W} \omega_{k,d}^\top x^d = \sum_{d=0}^{D_W} \left( \sum_{p=1}^{P} \omega_{k,d,p} x_p^d \right),$$

$$\mu_k(x) := \mu_k(x; \mathbf{\Upsilon}_k) = \sum_{d=0}^{D_V} \mathbf{\Upsilon}_{k,d} x^d, \text{ with } \mathbf{\Upsilon}_{k,d} \in \mathbb{R}^{Q \times P}.$$

Then, let $\omega = (\omega_k)_{k \in [K]}$, $\omega_k = (\omega_{k,d,p})_{d \in \{0, \dots, D_W\}, p \in [P]}$ and $\mathbf{\Upsilon}_k = (\mathbf{\Upsilon}_{k,d})_{d \in \{0, 1, \dots, D_V\}}$ be the tuples of unknown coefficients with the maximum degrees $D_W$ and $D_V$ of polynomials for the weight and mean functions, respectively. Subsequently, the unknown parameters of the model are denoted as follows: $\theta = (\omega_k, \mathbf{\Upsilon}_k, \Sigma_k)_{k \in [K]}$.

At this point, the global regression function is defined as the conditional expectation of the mixed probability density function in Equation (11), that is

$$m(x) = \mathbb{E}(y|x) = \sum_{k=1}^{K} g_k(w(x)) \mu_k(x), \quad (14)$$

### 4.1.2 Identifiability of softmax-gated Gaussian MoE models

It is worth noting that if we translate $\omega_{k,d}$ to $\omega_{k,d} + t$, $t \in \mathbb{R}^P$, the values of the softmax gating function remain unchanged. This implies that the gating parameters $\omega_{k,d}$ are identifiable only up to a translation. This issue has also been discussed and addressed in the context of establishing the convergence rate of parameter estimation via MLE for softmax-gated Gaussian MoE models in [73]. To alleviate this problem, based on [39, 48], we parameterize the gating parameters via the constraints, without loss of generality, $\{\omega_{K,d,p}\}_{d \in \{0, \dots, D_W\}, p \in [P]} = 0$ such that

$$g_K(x; \omega) := 1 - \sum_{k=1}^{K-1} g_k(x, \omega), \quad g_k(x; \omega) := \frac{\exp(w_k(x; \omega_k))}{1 + \sum_{l=1}^{K-1} \exp(w_l(x; \omega_l))}, \quad \forall k \in [K-1].$$
### 4.1.3 Incremental stochastic MM algorithm for softmax-gated Gaussian MoE models

We begin by introducing the following notations:

* $s_{n,4,k:q}$ is the segment of the vector $s_{n,4}$ spanning from position $\{[(k-1)Q+q]-1\} P(D_V+1)+1$ to $[(k-1)Q+q]P(D_V+1)$. Similarly, $s_{n,5,k:q}$ denotes the segment of the vector $s_{n,5}$ ranging from position $\{[(k-1)Q+q]-1\} [P(D_V+1)]^2+1$ to $[(k-1)Q+q][P(D_V+1)]^2$.

* $s_{n,3,kq}$ is the $kq$-th component of the vector $s_{n,3}$, while $s_{n,6,kq}$ denotes the $kq$-th component of the vector $s_{n,6}$.

* $\mathrm{mat}(s_{n+1,2})$ is the process of reshaping the vector $s_{n+1,2}$ into a matrix with dimensions $(K-1)(D_W+1)P \times (K-1)(D_W+1)P$. Similarly, $\mathrm{mat}(s_{n+1,5,k})$ denotes the reshaping of the vector $s_{n+1,5,k:q}$ into a square matrix of dimensions $P(D_V+1) \times P(D_V+1)$.

* The posterior probability that the data point $(\mathbf{x}_n, \mathbf{y}_n)$ belongs to the $k$-th expert is

$$ \tau_{n,k} = \frac{g_k(\mathbf{x}_n; \omega_n) \mathcal{N}(\mathbf{y}_n; \boldsymbol{\mu}_{n,k}(\mathbf{x}_n; \mathbf{\Upsilon}_{n,k}), \boldsymbol{\Sigma}_{n,k})}{\sum_{l=1}^{K} g_l(\mathbf{x}_n; \omega_n) \mathcal{N}(\mathbf{y}_n; \boldsymbol{\mu}_{n,l}(\mathbf{x}_n; \mathbf{\Upsilon}_{n,l}), \boldsymbol{\Sigma}_{n,l})}, \quad (15) $$

where, $\theta_n = (\omega_{n,k}, \mathbf{\Upsilon}_{n,k}, \boldsymbol{\Sigma}_{n,k})_{k \in [K]}$ denotes the parameter vector at the $n$-th iteration (time point).

* $\xi_n = [\tau_{n,1}x_{n,1}^0, \tau_{n,1}x_{n,1}^1, \dots, \tau_{n,1}x_{n,1}^{D_W}, \dots, \tau_{n,K-1}x_{n,1}^{D_W}, \dots, \tau_{n,K-1}x_{n,P}^{D_W}]^\top$.

* $\tau_{n,:} = [\tau_{n,1}, \tau_{n,2}, \dots, \tau_{n,K}]^\top$, $\mathbf{r}_n = \mathrm{vec}(\{\mathbf{x}_n^0, \mathbf{x}_n^1, \dots, \mathbf{x}_n^{D_V}\})$.

* $\mathbf{\Upsilon}_k = [\mathbf{\Upsilon}_{k,0}, \dots, \mathbf{\Upsilon}_{k,D_V}]$, $\mathbf{\Upsilon}_{k,q,:}$ is the vector containing the $q$-th row of $\mathbf{\Upsilon}_k$ for $q \in [Q]$.

* Given $\hat{\mathbf{x}}_n = [x_{n,1}^0, \dots, x_{n,1}^{D_W}, \dots, x_{n,P}^0, \dots, x_{n,P}^{D_W}]^\top$, we define the following matrix:

$$ \mathbf{B}_{n,K} = \left( \frac{3}{4} \mathbf{I}_{K-1} - \frac{\mathbf{1}_{K-1} \mathbf{1}_{K-1}^\top}{2(K-1)} \right) \otimes \hat{\mathbf{x}}_n \hat{\mathbf{x}}_n^\top + \varepsilon^* \mathbf{I}_{(K-1)P(D_W+1)}, \quad (16) $$

where $\varepsilon^*$ is a positive real number chosen sufficiently small.

With $\mathbf{B}_{n,K}$ defined in Equation (16), it can be proven that any eigenvalue of $\mathbf{B}_{n,K}$ is bounded in the interval $(\frac{1}{M_0}, M_0)$ for some positive real number $M_0$. Then we can define $\mathcal{S}$ as follows:

$$ \mathcal{S} = \{ (\mathbf{s}_1, \mathrm{vec}(\mathbf{S}_2), \mathbf{s}_3, \mathbf{s}_4, \mathbf{s}_5, \mathbf{s}_6) : \begin{array}{l} \mathbf{s}_1 \in \mathbb{R}^{(K-1)P(D_W+1)}; \\ \mathbf{S}_2 \in \mathbb{R}^{(K-1)P(D_W+1) \times (K-1)P(D_W+1)}; \\ \mathbf{S}_2 \succeq \mathbf{0}, \lambda \in \left(\frac{1}{M_0}, M_0\right) \text{ for all eigenvalues } \lambda \text{ of } \mathbf{S}_2; \\ \mathbf{s}_3 \in \mathbb{R}^{KQ}; \quad \mathbf{s}_4 \in \mathbb{R}^{KPQ(D_V+1)}; \quad \mathbf{s}_6 \in \mathbb{R}_+^{KQ}; \\ \mathbf{s}_5 \in \mathbb{R}^{KQ(P(D_V+1))^2}, \mathrm{mat}(\mathbf{s}_{5,k:q}) \succeq \mathbf{0} \text{ for all } k \in [K], q \in [Q] \end{array} \}. \quad (17) $$

Through essential inequalities for constructing surrogate functions as outlined in Section C, and following the procedure in Section 2.2, we derive the exponential family majorizer surrogate function for Equation (13). This is formally presented in Proposition 3 and its proof is provided in Section A.4.1.

**Proposition 3.** The construction of the exponential family majorizer surrogate function for the softmax-gated Gaussian MoE models is as follows:

$$ \begin{aligned} -\log s_{\theta}(\mathbf{y}_n | \mathbf{x}_n) \le & \sum_{k=1}^{K} \tau_{n,k} \log(\tau_{n,k}) + \mathbf{g}(\omega_n; \mathbf{x}_n) - \omega_n^\top \nabla \mathbf{g}(\omega_n; \mathbf{x}_n) - \sum_{k=1}^{K} \tau_{n,k} w_k(\mathbf{x}_n) \\ & + \left\{ \omega^\top \nabla \mathbf{g}(\omega_n; \mathbf{x}_n) + \frac{1}{2}(\omega - \omega_n^\top) \mathbf{B}_{n,K} (\omega - \omega_n) \right\} \\ & - \sum_{k=1}^{K} \tau_{n,k} \left[ \log(\mathcal{N}(\mathbf{y}_n, \boldsymbol{\mu}_k(\mathbf{x}_n), \boldsymbol{\Sigma}_k)) \right] =: -\log[g(\theta, \mathbf{x}_n, \mathbf{y}_n; \theta_n)]. \quad (18) \end{aligned} $$
Here $\underline{g}(\omega_n; \mathbf{x}_n) = \log \left( 1 + \sum_{k=1}^{K-1} \exp(w_{n,k}(\mathbf{x}_n)) \right)$ where $w_{n,k}(\mathbf{x}_n) = \sum_{d=0}^{D_W} \omega_{n,k,d}^{\top} \mathbf{x}_n^d$.

We adapt the incremental stochastic MM algorithm Algorithm 1 for the softmax-gated Gaussian MoE models as in Algorithm 2 as follows.

**Algorithm 2** Incremental stochastic MM algorithm for softmax-gated Gaussian MoE models

**Require:** Incremental data points $\{(\mathbf{x}_n, \mathbf{y}_n) : n \ge 0\}$, positive step sizes $\{\gamma_{n+1}, n \ge 1\}$ in $(0, 1)$ and an initial value $s_0 = (s_{0,1}, \dots, s_{0,6}) \in \mathbb{S}$, where $\mathbb{S}$ is defined in Equation (17), such that:

$$ \begin{aligned} \text{mat}(s_{0,2}) &\succeq 0, \quad s_{0,3} = \mathbf{1}_{KQ}, \quad \hat{\mathbf{s}} \in \mathbb{R}^{P(D_V+1)}, \quad s_{0,4} = \mathbf{1}_{KQ} \otimes \hat{\mathbf{s}}, \\ s_{0,5} &= \mathbf{1}_{KQ} \otimes \text{vec}(\hat{\mathbf{s}}\mathbf{s}^{\top} + \mathbf{I}_{P(D_V+1)}), \quad s_{0,6} > 0. \end{aligned} $$

**Ensure:** A $\mathbb{T}$-valued sequence $\{\theta_n, n \ge 0\}$ and a $\mathbb{S}$-valued sequence $\{s_n, n \ge 0\}$.

1: **for** $n = 0, 1, \dots, \text{do}$

2: Compute

$$ \begin{aligned} s_{n+1,1} &= s_{n,1} + \gamma_{n+1} (-\xi_{n+1} + \nabla g(\omega_n; \mathbf{x}_n) - \mathbf{B}_{n+1,K} \omega_n - s_{n,1}), \\ s_{n+1,2} &= s_{n,2} + \gamma_{n+1} \left[ \frac{1}{2} \text{vec}(\mathbf{B}_{n+1,K}) - s_{n,2} \right], \\ s_{n+1,3} &= s_{n,3} + \gamma_{n+1} (\tau_{n+1,:} \otimes \mathbf{y}_{n+1}^2 - s_{n,3}), \\ s_{n+1,4} &= s_{n,4} + \gamma_{n+1} (-2\tau_{n+1,:} \otimes (\mathbf{y}_{n+1} \otimes \mathbf{r}_{n+1}) - s_{n,4}), \\ s_{n+1,5} &= s_{n,5} + \gamma_{n+1} (\tau_{n+1,:} \otimes (\text{vec}(\mathbf{r}_{n+1}\mathbf{r}_{n+1}^{\top}) \otimes \mathbf{1}_Q) - s_{n,5}), \\ s_{n+1,6} &= s_{n,6} + \gamma_{n+1} (\tau_{n+1,:} \otimes \mathbf{1}_Q - s_{n,6}), \\ \omega_{n+1} &= -(2\text{mat}(s_{n+1,2}))^{-1} s_{n+1,1}, \\ \mathbf{\Upsilon}_{n+1,k,q,:} &= -(2\text{mat}(s_{n+1,5,k,q}))^{-1} s_{n+1,4,k,q}, \\ \sigma_{n+1,k,q}^2 &= \frac{s_{n+1,3,kq} + s_{n+1,4,k,q}^{\top} \mathbf{\Upsilon}_{k,q,:}^{\top} + s_{n+1,5,k,q}^{\top} \text{vec}(\mathbf{\Upsilon}_{k,q,:}^{\top} \mathbf{\Upsilon}_{k,q,:})}{s_{n+1,6,kq}}, \\ \theta_{n+1} &= (\omega_{n+1,k,p}, \mathbf{\Upsilon}_{n+1,k,q,:}, \sigma_{n+1,k,q}^2)_{k \in [K]; q \in [Q]; p \in [P]} \end{aligned} $$

3: **end for**

A proof of Theorem 2, which provides the theoretical guarantees for Algorithm 2, is given in Section A.4. Meanwhile, a detailed derivation of the update procedure in Algorithm 2 is presented in Section A.5.

**Theorem 2** (Stability and consistency for Algorithm 2). Consider the softmax-gated Gaussian MoE setting described in Section 4.1, and assume that the covariates and responses are bounded. With the surrogate construction given in Equation (18) and the choice of $\mathbf{B}_{n,K}$ in Equation (16), the model-specific assumptions P 1 to 5 and Assumptions 1 to 3 and 7 are verified (see Section A.4.2). The remaining assumptions required to invoke the general stochastic-approximation convergence results are treated as standing assumptions for this MoE instantiation, namely Assumptions 4 to 6 (Lyapunov-level regularity, stability, and compact sublevel sets for almost-sure convergence). In addition, the step-size condition Assumption 8 can be enforced by design, for example by taking $\gamma_n = \gamma_0 n^{-\alpha}$ with $\alpha \in (\frac{1}{2}, 1]$ and $\gamma_0 \in (0, 1)$. Under the verified assumptions above, together with these standing assumptions, the conclusions of Propositions 1 and 2 and Theorem 1 apply to Algorithm 2.

## 4.2 Softmax-gated multinomial logistic MoE models for discrete output data

In this situation, the output variable $y \in \mathbb{Y} := [M], M \in \mathbb{N}$, is a discrete response variable. Furthermore, the multinomial (multiclass) logistic regression expert network is defined via the
conditional probability function as follows:

$$e_k(y = m; v(x; v)) := e_k(v_{y_n}(x_n; v)) = \frac{\exp(v_{m,k}(x))}{\sum_{l=1}^{M} \exp(v_{l,k}(x))}, \forall k \in [K], \forall m \in [M].$$

Here the expert function $v = (v_m)_{m \in [M]} = (v_{m,k})_{m \in [M], k \in [K]}$ is modeled as polynomial function of the input variable $x$ as follows:

$$v_{m,k}(\boldsymbol{x}) = \sum_{d=0}^{D_V} \boldsymbol{v}_{m,k,d}^{\top} \boldsymbol{x}^{d} = \sum_{d=0}^{D_V} \left( \sum_{p=1}^{P} v_{m,k,d,p} \boldsymbol{x}_p^d \right) \quad \text{for any } \boldsymbol{v}_{m,k,d} = (v_{m,k,d,p})_{p \in [P]} \in \mathbb{R}^P.$$

Then, let $\omega = (\omega_k)_{k \in [K]}$ where $\omega_k = (\omega_{k,d,p})_{d \in \{0, \dots, D_W\}, p \in [P]}$ and $v = (v_m)_{m \in [M]} = (v_{m,k})_{m \in [M], k \in [K]}$ where $v_{m,k} = (v_{m,k,d,p})_{d \in \{0, 1, \dots, D_V\}, p \in [P]}$ be the tuples of unknown coefficients with the maximum degrees $D_W$ and $D_V$ of polynomials for the softmax-gated and expert functions, respectively. Then, we define

$$e_k(y = m; v(x; v)) := e_k(v_{y_n}(x_n; v)), \forall k \in [K].$$

Finally, the unknown parameters of the softmax-gated multinomial logistic MoE models are denoted as follows: $\theta = (\omega_k, v_{m,k})_{m \in [M], k \in [K]}$.

### 4.2.1 Identifiability of softmax-gated multinomial logistic MoE models

Apart from addressing the identifiability of the parameters from in the softmax-gated Gaussian MoE models in Section 4.1.2, we also propose an additional constraint on the parameters, which is also considered in the study of the convergence rate of parameter estimation for softmax-gated multinomial logistic MoE models in [72]. Specially, for the multinomial logistic regression expert network, we impose the following constraint: without loss of generality, $\{v_{M,k,d,p}\}_{d \in \{0, \dots, D_W\}, p \in [P]} = \mathbf{0}$ for all $k \in [K]$ such that

$$e_k(v_M(\boldsymbol{x})) = 1 - \sum_{m=1}^{M-1} e_k(v_m(\boldsymbol{x})), \quad e_k(v_m(\boldsymbol{x})) = \frac{\exp(v_{m,k}(\boldsymbol{x}))}{1 + \sum_{l=1}^{M-1} \exp(v_{l,k}(\boldsymbol{x}))}, \forall m \in [M-1].$$

### 4.2.2 Incremental MM algorithm for softmax-gated multinomial logistic MoE models

Let $\mathbb{1}(z, l)$ be an indicator that maps $[M] \times [M] \to \{0, 1\}$ where $\mathbb{1}(z, l) = 1$ if $z = l$ and 0, otherwise. We can choose that function as follows: $\mathbb{1}(z, l) = [\prod_{q=1}^{M}(z - q)] / [(z - l)(z - 1)(M - z)!( -1)^{M-z}]$. Then, given $r_{d,l,p,k} = \mathbb{1}(y, l)\tau_k x_p^d$, we denote

$$\mathbf{r}_n = [r_{n,0,1,1,1}, \dots, r_{n,D_V,1,1,1}, \dots, r_{n,D_V,M,1,1}, \dots, r_{n,D_V,M,P,1}, \dots, r_{n,D_V,M,P,K}]^{\top}.$$

Given any $t \in \mathbb{N}$, $n \in \mathbb{N}$, and any matrix $A \in \mathbb{R}^{nt \times tn}$, we define $bdiag_t(A) \in \mathbb{R}^{nt \times t}$ as follows:

$$bdiag_t(\mathbf{A})_{[it+1:it+t], [1:t]} = \mathbf{A}_{[it+1:it+t], [it+1:it+t]}, \quad \forall i \in \{0, \dots, n-1\}.$$

For any $K \in \mathbb{N}$ and $M \in \mathbb{N}$, let $B_{n,K}$, $B_{n,M}$ and $\tau_{n,\cdot}$ to have the same definition as the previous Section 4.1.3. Now we go into the main part of finding the surrogate function of Equation (13) via Proposition 4, which is proved in Section A.6.1.

**Proposition 4.** The construction of the exponential family majorizer surrogate construction for softmax-gated multinomial logistic MoE models is as follows:

$$\begin{aligned} - \log s_{\theta}(y_n | \mathbf{x}_n) & \le C_{n,k} - (\boldsymbol{\omega}^{\top} \boldsymbol{\xi}_n) - (\boldsymbol{v}^{\top} \mathbf{r}_n) \\ & \quad + \left\{ \underline{g}(\omega_n; \mathbf{x}_n) + (\omega - \omega_n)^{\top} \nabla \underline{g}(\omega_n; \mathbf{x}_n) + \frac{1}{2}(\omega - \omega_n)^{\top} \mathbf{B}_{n,K}(\omega - \omega_n) \right\} \\ & \quad + \sum_{k=1}^{K} \tau_{n,k} \left\{ \underline{e}(\mathbf{c}_{n,k}; \mathbf{x}_n) + (\mathbf{c}_k - \mathbf{c}_{n,k})^{\top} \nabla \underline{e}(\mathbf{c}_{n,k}; \mathbf{x}_n) \right. \\ & \quad \left. + \frac{1}{2}(\mathbf{c}_k - \mathbf{c}_{n,k})^{\top} \mathbf{B}_{n,M}(\mathbf{c}_k - \mathbf{c}_{n,k}) \right\} =: - \log[g(\boldsymbol{\theta}, \mathbf{x}_n, y_n; \boldsymbol{\theta}_n)], \quad (19) \end{aligned}$$
where:

$$c_{n,k} = [\nu_{n,1,k,0,1}, \dots, \nu_{n,1,k,D_V,1}, \dots, \nu_{n,M,k,D_V,1}, \dots, \nu_{n,M,k,D_V,P}]^\top,$$

$$\underline{e}(c_{n,k}; \mathbf{x}_n) = \log \left( 1 + \sum_{l=1}^{M-1} \exp(v_{n,l,k}(\mathbf{x}_n)) \right), \quad \underline{g}(\omega_n; \mathbf{x}_n) = \log \left( 1 + \sum_{k=1}^{K-1} \exp(w_{n,k}(\mathbf{x}_n)) \right),$$

$$\tau_{n,k} = \frac{g_k(\mathbf{x}_n; \omega_n) e_k(v_{y_n}(\mathbf{x}_n; v_{y_n}))}{\sum_{l=1}^{K} g_l(\mathbf{x}_n; \omega_n) e_l(v_{y_n}(\mathbf{x}_n; v_{y_n}))},$$

**Algorithm 3 Incremental MM algorithm for softmax-gated multinomial logistic MoE models**

**Require:** Incremental data points {($\mathbf{x}_n, y_n$) : $n \ge 0$}, positive step sizes {$\gamma_{n+1}, n \ge 1$} in (0, 1), and an initial value $s_0 = (s_{0,1}, \dots, s_{0,6}) \in \mathbb{S}$ such that:

$$\mathrm{mat}(s_{0,2}) \succeq 0, \quad \mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{0,4})) \succeq 0, \quad \text{where } \iota = \mathrm{MP}(D_V + 1).$$

**Ensure:** A T-valued sequence {$\theta_n, n \ge 0$} and a S-valued sequence {$s_n, n \ge 0$}:

1. **for** $n = 0, 1, \dots, \mathbf{do}$

2. Compute

$$s_{n+1,1} = s_{n,1} + \gamma_{n+1} (-s_{n+1} + \nabla \underline{g}(w_n; \mathbf{x}_n) - B_{n+1,K} w_n - s_{n,1}),$$

$$s_{n+1,2} = s_{n,2} + \gamma_{n+1} \left[ \frac{1}{2} \mathrm{vec}(B_{n+1,K}) - s_{n,2} \right],$$

$$s_{n+1,3} = s_{n,3} + \gamma_{n+1} (-\mathbf{r}_{n+1} + \mathrm{vec}([\tau_{n+1,k} \nabla \underline{e}(v_n; \mathbf{x}_n) - \tau_{n+1,k} B_{n+1,M} v_n]_{k=1,\dots,K}) - s_{n,3}),$$

$$s_{n+1,4} = s_{n,4} + \gamma_{n+1} \left[ \frac{1}{2} \mathrm{vec}(\tau_{n+1,:} \otimes B_{n+1,M}) - s_{n,4} \right],$$

$$\mathbf{w}_{n+1} = - (\mathrm{mat}(s_{n+1,2}) + \mathrm{mat}(s_{n+1,2})^\top)^{-1} s_{n+1,1},$$

$$\mathbf{v}_{n+1} = - [\mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n+1,4})) + \mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n+1,4}))^\top]^{-1} s_{n+1,3},$$

$$\theta_{n+1} = (\mathbf{w}_{n+1,k}, \mathbf{v}_{n+1,k})_{k \in [K]}.$$

3. **end for**

We provide the complete derivation of the stochastic recursion underlying Algorithm 3 in Section A.7. The resulting sequence of iterates fits naturally within the general incremental MM framework. Consequently, the convergence properties of the algorithm follow from standard stochastic-approximation arguments, after verifying the model-dependent assumptions. This leads to the stability and consistency result stated below, whose proof is given in Section A.6.2.

**Theorem 3 (Stability and consistency for Algorithm 3).** Consider the softmax-gated multinomial logistic MoE model introduced in Section 4.2, and assume that both the covariates and the responses are bounded. With the surrogate construction specified in Equation (19) and the parameter-update rules defined in Algorithm 3, the MoE-specific assumptions P 1 to 5 and Assumptions 1 to 3 and 7 are verified (see Section A.6.2). The remaining hypotheses required to apply the general stochastic-approximation convergence theory are treated as standing assumptions for this model instantiation, namely Assumptions 4 to 6 (Lyapunov regularity, stability, and compactness of sublevel sets). Moreover, the step-size condition Assumption 8 can be enforced by construction, for instance by choosing $\gamma_n = \gamma_0 n^{-\alpha}$ with $\alpha \in (\frac{1}{2}, 1]$ and $\gamma_0 \in (0, 1)$. Under the verified assumptions together with these standing conditions, all conclusions of Propositions 1 and 2 and Theorem 1 hold for Algorithm 3.

## 4.3 Why have other variants of incremental stochastic MM and EM algorithms not been successfully implemented in practice for the softmax-gated MoE models?

In the literature, numerous incremental stochastic EM and MM algorithms have been proposed to optimize complex latent-variable models. However, most of these methods are not directly applicable to the softmax-gated MoE model in Equation (11). A common limitation is that many online EM variants are formulated under an exponential-family assumption on the target model,
whereas Equation (11) does not generally belong to an exponential family. For example, [10] present an online EM procedure assuming that the complete-data model is in an exponential family (see Assumption 1(a) therein). While this is stated as a structural assumption on the model itself, a careful reading of their analysis shows that the key ingredient exploited in the convergence argument is the availability of a surrogate/objective representation with an exponential-family-like structure, namely a decomposition in terms of sufficient statistics and a parameter update given by a deterministic maximization map. In our setting, we do not impose an exponential-family form on the softmax-gated MoE model Equation (11); instead, we require this structure only for the majorizer surrogate (cf. P 1 to 5), which yields a strictly weaker requirement and allows a more general treatment. Subsequent extensions, such as [9] (Hidden Markov Models) and [8] (state-space models), build on the same framework and retain the exponential-family condition, and thus they likewise do not directly cover Equation (11). Other EM-based adaptations exhibit similar constraints: the Fast iEM algorithm of [51] introduces a bias-correction step but still leverages exponential-family structure; the SAEM approach in [54] employs an ergodic Markov sampling kernel yet preserves the same requirement; and the Monte Carlo E-step strategy of [95], which samples from distributions defined by the current parameter iterate, again relies on exponential-family assumptions.

MM-based schemes often fail for various reasons, primarily due to their reliance on overly restrictive regularity conditions. For instance, the stochastic MM algorithm proposed in [67] requires strict convexity and Lipschitz continuity conditions that are typically violated by Equation (11). The incremental stochastic MM algorithm proposed in [28] relies on expectation bounds stated in Assumptions A2(a) and A3(a) of their paper (discussed in Condition 1), which are generally not satisfied by the surrogate functions defined in Equation (18) and Equation (19). These surrogate functions correspond to the original objective function in Equation (11) for the softmax-gated Gaussian and multinomial logistic MoE models for continuous and discrete output data, respectively. In [68], a variety of surrogate functions are proposed; however, aside from the quadratic majorizer (also adopted in our work), all other surrogates require properties such as L-smoothness, Lipschitz continuity, or convexity. Lastly, the large-scale incremental stochastic MM approach in [52] employs the MISSO method, which is based on the MISO method proposed by [68] and updates surrogate functions only partially at each iteration. However, they require that the surrogate function satisfies Assumption H2 (see Condition 2), and even though this condition is true if $\hat{e}(\theta; \theta_n)$ is L-smooth. Nonetheless, formulating such a surrogate for Equation (11) is still extremely challenging.

In summary, there is no existing incremental stochastic variant of the EM or MM algorithm that can be directly applied to SGMoE models, according to our knowledge. The failure of SGMoE to satisfy strict conditions such as Lipschitz continuity, L-smoothness, or convexity is demonstrated in Section 4.3.2. While these conditions can be rigorously shown to not hold in general, Condition 1 and Condition 2 remain unproven. This is not due to their validity but rather the inherent difficulty of establishing the non-existence of a suitable surrogate function. At best, we can only suggest that constructing such a surrogate appears to be intractable in practice.

### 4.3.1 Restrictive regularity conditions in existing incremental stochastic MM frameworks

We now examine the regularity assumptions imposed by several recent variants of incremental stochastic MM algorithms and explain why these conditions are not satisfied by the SGMoE objective considered in this work.

**Condition 1** (Assumptions A2(a) and A3(a) in [28]). Let $g(\theta, z, u)$ denote a surrogate function for the objective, where $u$ belongs to an auxiliary index set $\mathcal{U}$. Assumption A2(a) in [28] requires

$$ \mathbb{E}_{\pi} [ |g(\theta, z, u)| ] < \infty, \quad \forall (\theta, u) \in \mathbb{T} \times \mathcal{U}, $$

while Assumption A3(a) further assumes the existence of some $p > 1$ such that

$$ \mathbb{E}_{\pi} \left[ \sup_{(\theta, u) \in \mathbb{T} \times \mathcal{U}} |g(\theta, z, u)|^p \right] < \infty. $$

**Condition 2** (Assumption H2 in [52]). Consider objective functions $\mathcal{L}_n(\theta)$ indexed by $n \in [N]$, and let $\hat{\mathcal{L}}_n(\theta; \bar{\theta}_n)$ be a majorizing surrogate of $\mathcal{L}_n$ at $\bar{\theta}_n \in \mathbb{T}$. Define

$$ \hat{e}(\theta; \bar{\theta}_n) := \hat{\mathcal{L}}_n(\theta; \bar{\theta}_n) - \mathcal{L}_n(\theta). $$

Assumption H2 in [52] postulates the existence of a constant $L > 0$ such that

$$ \|\nabla \hat{e}(\theta; \bar{\theta}_n)\|^2 \le 2L \hat{e}(\theta; \bar{\theta}_n), \quad \forall \theta \in \mathbb{T}. $$
**Remark 1.** The surrogate function constructed in Equation (18) (see Proposition 3) does not satisfy the moment conditions in Condition 1. Indeed, if the covariance matrices $\Sigma_k$ are allowed to grow unbounded while all other components of $\theta$ are fixed, the Gaussian log-likelihood term

$$-\sum_{k=1}^{K} \tau_{n,k} \log(\mathcal{N}(\mathbf{y}_n; \boldsymbol{\mu}_k(\mathbf{x}_n), \boldsymbol{\Sigma}_k))$$

diverges due to the $\log|\boldsymbol{\Sigma}_k|$ contribution. Consequently, $\sup_{(\theta, u)} |g(\theta, z, u)|^p$ is not integrable for any $p > 1$, violating both Assumptions A2(a) and A3(a) of [28].

Assumption H2 of [52] is likewise not applicable in our setting. That framework requires the construction of surrogates directly majorizing the objective $\mathcal{L}_n$, which in turn relies on strong global regularity properties of $f(\theta; z)$ such as Lipschitz continuity of the gradient, L-smoothness, or convexity. As we show below, these properties do not generally hold for the SGMoE objective.

### 4.3.2 Failure of Lipschitz continuity, smoothness, and convexity for SGMoE

Many existing incremental stochastic MM algorithms rely on strong global regularity assumptions on the objective function, such as Lipschitz continuity, L-smoothness, or convexity. We now show that these properties fail, in general, for the SGMoE negative log predictive density

$$f(\theta; z) = -\log s_{\theta}(\mathbf{y} | \mathbf{x}),$$

even when covariates and responses are bounded.

Let $\boldsymbol{\mu}(\mathbf{x})$ be a polynomial function of $\mathbf{x}$. Then there exists a parameter matrix $\mathbf{\Upsilon}_0$ such that $\boldsymbol{\mu}_0(\mathbf{x}) = \mathbf{y}$. Define a sequence of parameter vectors

$$\theta_n := (\omega, \mathbf{\Upsilon}_0, \boldsymbol{\Sigma}_n), \quad \theta^* := (\omega, \mathbf{\Upsilon}_0, \boldsymbol{\Sigma}^*),$$

where $\boldsymbol{\Sigma}_n \succeq 0$ with $\boldsymbol{\Sigma}_n \downarrow 0$ and $\boldsymbol{\Sigma}^* \succeq 0$ is fixed.

**Failure of global Lipschitz continuity.** For Gaussian experts, when $\boldsymbol{\mu}_0(\mathbf{x}) = \mathbf{y}$, the SGMoE density satisfies

$$s_{\theta_n}(\mathbf{y} | \mathbf{x}) \geq g_k(\mathbf{w}(\mathbf{x})) \mathcal{N}(\mathbf{y}; \mathbf{y}, \boldsymbol{\Sigma}_n)$$

for some $k \in [K]$, and therefore

$$f(\theta_n; z) = -\log s_{\theta_n}(\mathbf{y} | \mathbf{x}) \geq -\log \mathcal{N}(\mathbf{y}; \mathbf{y}, \boldsymbol{\Sigma}_n) = \frac{1}{2} \log |\boldsymbol{\Sigma}_n| + C \xrightarrow{n \to \infty} +\infty.$$

At the same time, the parameter difference $\|\theta_n - \theta^*\|$ remains bounded for any norm on $\mathbb{R}^T$. Hence,

$$\frac{|f(\theta_n; z) - f(\theta^*; z)|}{\|\theta_n - \theta^*\|} \xrightarrow{n \to \infty} +\infty,$$

which rules out global Lipschitz continuity of $f(\theta; z)$.

**Failure of L-smoothness.** The gradient of the Gaussian negative log-likelihood with respect to $\boldsymbol{\Sigma}$ contains terms of the form $\boldsymbol{\Sigma}^{-1}$ and $\boldsymbol{\Sigma}^{-2}$. Consequently,

$$\|\nabla_{\theta} f(\theta_n; z)\| \xrightarrow{n \to \infty} +\infty \quad \text{as} \quad \boldsymbol{\Sigma}_n \downarrow 0,$$

showing that no global Lipschitz constant on the gradient can exist. Thus, $f$ is not L-smooth.

**Failure of convexity.** Let $\boldsymbol{\mu}_1(\mathbf{x}) \neq \mathbf{y}$ and $\boldsymbol{\mu}_2(\mathbf{x}) \neq \mathbf{y}$ satisfy

and define

$$\theta_i := (\omega, \mathbf{\Upsilon}_i, \boldsymbol{\Sigma}_n), \quad i = 1, 2,$$

where $\mathbf{\Upsilon}_i$ corresponds to $\boldsymbol{\mu}_i$. As $\boldsymbol{\Sigma}_n \downarrow 0$, both $f(\theta_1; z)$ and $f(\theta_2; z)$ remain finite, since $\boldsymbol{\mu}_i(\mathbf{x}) \neq \mathbf{y}$. However,

$$f(\frac{\theta_1+\theta_2}{2}; z) = -\log \mathcal{N}(\mathbf{y}; \mathbf{y}, \boldsymbol{\Sigma}_n) \xrightarrow{n \to \infty} +\infty.$$

This violates Jensen's inequality and shows that the SGMoE objective is not convex, and hence not strongly convex, in general.

These results demonstrate that the SGMoE negative log predictive density does not satisfy global Lipschitz continuity, L-smoothness, or convexity. As a consequence, several existing incremental stochastic MM frameworks that rely on these restrictive regularity conditions are not directly applicable to SGMoE models.
# 5 Numerical experiments

In this section, we conduct comprehensive experiments on both synthetic and real-world datasets to evaluate the effectiveness and robustness of the proposed method. All experiments were conducted using Python 3.12.2 on a standard Unix-based system.

## 5.1 Experimental settings

**Assessment of initialization schemes.** We consider two initialization strategies:

(i) Perturbed ground truth initialization: The initial parameters are generated by adding Gaussian noise to the true data-generating parameters. This strategy mostly serves to test the local stability of the algorithm.

(ii) Data-driven initialization: This strategy involves partitioning the first batch of incoming data into clusters using a standard K-means clustering procedure and then compute an initial set of parameters. For real-world data, this is the realistic initialization approach.

**Assessment of regression performance.** Furthermore, we benchmark the regression performance of our algorithm against five baseline stochastic models using a standard train-test evaluation framework. With the regression function in Equation (14), there are two approaches for assessment. Firstly, for strictly synthetic data, we compute the regression values using the estimated parameters and compare it to the ones computed with true parameters, this is called estimation error. Secondly, for both synthetic and real-world data, we measure the “fit” of the regression line on the testing samples, similarly to how regression models are traditionally bench-marked, this is called prediction error. For both approaches, we utilize the following point-wise error metrics: mean squared error (MSE), mean absolute percentage error (MAPE), and normalized root mean squared error (NRMSE). It is also worth noting that, throughout this section, the maximum log-likelihood at each iteration of the incremental stochastic MM algorithm is computed by substituting the parameter estimates from that iteration into Equation (13).

**Comparative methods.** In this simulation study on MoE models, we compare our proposed method against several widely-used stochastic online or incremental stochastic optimization algorithms, including stochastic approximation or SGD [94], RMSProp [40], Adam [53], and Sophia [63].

## 5.2 Experiments on synthetic data

We now present experiments conducted on synthetic data. Specifically, we consider a dataset with two generating Gaussian distributions, which correspond to two clusters. Each cluster consists of $N = 1,000$ samples. We set the model configurations as follows: $D_W = 1$, $D_V = 1$, $P = 2$, and $Q = 1$. The true parameters, denoted by $\omega^0$, $v^0$, and $\sigma^0$, are specified as

$$ \omega^0 = \begin{bmatrix} [0 \quad 8] \\ [0 \quad 0] \\ [0 \quad 0] \end{bmatrix}, \quad v^0 = \begin{bmatrix} \begin{bmatrix} [0 \quad 0] \\ [-2.5 \quad 0] \\ [2.5 \quad 0] \end{bmatrix} \end{bmatrix}^T, \quad \sigma^0 = [[1], [1]]. $$

A visualization of the simulated dataset is provided in Figure 1. The dataset is divided into training and testing subsets, with 80% of the data allocated to training. It is important to note that under the perturbed ground truth initialization setting, the focus is to analyze the convergence behavior of the MLE, and thus only the training data is used. In contrast, the data-driven initialization scheme utilize both the training and testing data to evaluate generalization performance.

### 5.2.1 Evaluation on synthetic data with perturbed ground truth initialization

We first evaluate the performance of our algorithm under the perturbed ground truth initialization setting. More precisely, the noise variables are independently drawn from a Gaussian distribution with mean 0 and variance 1, and then scaled by a factor of 0.005. These small perturbations are added to the true parameters to generate the initial values. For each series, the initial latent state $s_{i,0}$ is computed as the average of the corresponding $S(z)$ values over the first 85 observations of $z$, using the initialized parameters.

In Figure 2a, we present the progression of the maximum log-likelihood value over $N = 1600$ iterations. Furthermore, we apply Polyak averaging, a widely adopted technique used to reduce
<figure><img src="image_2.png" /><figcaption>Figure 1: Visualization of the synthetic dataset generated from the true parameters and estimated clusters and regression functions by our proposed model with perturbed ground truth initialization in softmax-gated MoE model.</figcaption></figure>

variance in optimization trajectories and to dampen the influence of initially volatile parameter estimates. Starting from a designated iteration $N_0$ (we choose $N_0 = 100$), the averaged sequence is initialized as $\theta_{N_0}^P = 0$. For all $n \ge N_0$, the updates follow the rule:

$$ \theta_{n+1}^P = \theta_n^P + \alpha_{n-N_0+1} (\theta_n - \theta_n^P), \quad \text{where } \alpha_n \text{ is typically set to } \frac{1}{n}. $$

In the figure, the true MLE trajectory is plotted in black, and the MLE trajectory with Polyak averaging is shown in red. It can be observed that both trajectories converge, and the red line, which incorporates Polyak averaging, exhibits greater overall stability. Additionally, we plot the trajectories of selected parameter dimensions, including $\omega_{1,0,1}^0$, $\nu_{0,1,0,0}^0$, and $\sigma_{0,0}^0$, as shown in Figure 2b. These trajectories clearly demonstrate that our optimization method effectively drives the parameters to converge to their true values. Furthermore, in Figure 1b, we visualize the predicted mean values from our model compared to the true mean, providing further insight into the predictive accuracy of the estimated parameters.

### 5.2.2 Evaluation on synthetic data with data-driven initialization

To further demonstrate the robustness of our algorithm, we employ an alternative initialization strategy based on K-means clustering. Specifically, K-means is executed using randomly selected centroids, repeated up to ten times, and the configuration with the lowest distortion is selected as the initialization. The initial latent state values for each series, denoted by $s_{i,0}$, are computed in the same manner as previously described. The progression of the maximum log-likelihood function over iterations is presented in Figure 2c, while the convergence of selected parameter values is depicted in Figure 2d. Despite the inferior starting point provided by K-means initialization compared to the ideal setting, our algorithm exhibits strong robustness, successfully guiding the parameters toward their true values.

Additionally, Figure 3 illustrates the trajectory of the distance between the estimated and true log-likelihood functions over iterations. This plot clearly demonstrates that our proposed method achieves the fastest convergence rate among all compared baselines. To mitigate the impact of numerical instability observed during early iterations, we report the maximum log-likelihood values starting from iteration 400. The results, shown in Table 1, indicate that our proposed method consistently outperforms all baseline approaches across all metrics. Second, each algorithm is trained solely on the training data and evaluated on the test set, with prediction accuracy again measured against the ground truth using the same metrics. Details of the train/test data split are provided at the beginning of this section. The corresponding results are summarized in Table 2.
<figure><img src="image_3.png" /><figcaption>(a) Loss function with perturbed ground truth initialization. (b) Parameter estimation with perturbed ground truth initialization.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) Loss function with data-driven initialization. (d) Parameter estimation with data-driven initialization.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Distance to the true maximum log-likelihood of all algorithms for each iteration.</figcaption></figure>

<table><caption>Table 1: Estimation errors with data-driven initialization</caption><thead><tr><th>Model</th><th>MSE↓</th><th>MAPE↓</th><th>NMRSE↓</th></tr></thead><tbody><tr><td>Inc. MM</td><td><strong>0.014</strong></td><td><strong>0.228</strong></td><td><strong>0.009</strong></td></tr><tr><td>SGD</td><td>0.415</td><td>2.230</td><td>0.050</td></tr><tr><td>Adam</td><td>0.824</td><td>1.272</td><td>0.061</td></tr><tr><td>Adam W</td><td>0.192</td><td>1.816</td><td>0.035</td></tr><tr><td>RMSProp</td><td>0.322</td><td>0.545</td><td>0.042</td></tr><tr><td>Sophia</td><td>0.272</td><td>1.633</td><td>0.042</td></tr></tbody></table>

<table><caption>Table 2: Prediction errors with data-driven initialization</caption><thead><tr><th>Model</th><th>MSE↓</th><th>MAPE↓</th><th>NMRSE↓</th></tr></thead><tbody><tr><td>Inc. MM</td><td><strong>0.950</strong></td><td>0.480</td><td><strong>0.059</strong></td></tr><tr><td>SGD</td><td>1.517</td><td>0.332</td><td>0.075</td></tr><tr><td>Adam</td><td>1.789</td><td>0.422</td><td>0.082</td></tr><tr><td>Adam W</td><td>1.139</td><td><strong>0.219</strong></td><td>0.065</td></tr><tr><td>RMSProp</td><td>1.505</td><td>0.478</td><td>0.075</td></tr><tr><td>Sophia</td><td>1.256</td><td>0.257</td><td>0.069</td></tr></tbody></table>
## 5.3 Experiments on higher-dimensional synthetic data

In the previous section, we demonstrated that our proposed methods perform exceptionally well on low-dimensional synthetic data. However, real-world datasets are often higher-dimensional, and the performance observed in low-dimensional settings may not necessarily generalize. To evaluate the robustness of our approach under more realistic conditions, we conduct experiments on simulated data with significantly higher dimensionality. As we show in this section, even in these more challenging scenarios, our model consistently outperforms baseline methods, underscoring its robustness and effectiveness in high-dimensional settings.

Similar to the previous section, we set the number of data points and the number of clusters to 1,000 and 2, respectively. The parameters $D_W$ and $D_V$ are set to 1, while $P$ and $Q$ are chosen to be 10 and 1, respectively. Additionally, the true parameters, denoted by $\omega^0$, $\nu^0$, and $\sigma^0$, are initialized using Algorithm 4.

**Algorithm 4 Initialization of $\omega^0$, $\nu^0$, and $\sigma^0$ for Higher-Dimensional Simulation**

1: **Input:** Number of experts $K$, dimensions $D_V$, $D_W$, parameters $P$, $Q$, seed $s_e$.

2: Initialize tensors:

• $\nu^0 \in \mathbb{R}^{K \times (D_V+1) \times Q \times P} \leftarrow 0$.

• $\sigma^0 \in \mathbb{R}^{K \times Q \times Q} \leftarrow 0$.

• $\omega^0 \in \mathbb{R}^{P \times K \times (D_W+1)} \leftarrow 0$.

3: **Step 1: Initialization of variance components**

4: **for** $k=1$ **to** $K$ **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;**for** $q=1$ **to** $Q$ **do**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set random seed $s_e + kQ + q$.

7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Sample $\sigma^0[k, q, q] \sim \text{Uniform}(0.5, 1.5)$.

8: &nbsp;&nbsp;&nbsp;&nbsp;**end for**

9: **end for**

10: **Step 2: Initialization of expert parameters**

11: **for** $(k, d, q, p)$ over their respective ranges $([K] \times \{\{0\} \cup [D_V]\} \times [Q] \times [P])$ **do**

12: &nbsp;&nbsp;&nbsp;&nbsp;Set random seed: $s_e + s_e + k(D_V + 1)QP + dQP + qP + p$.

13: &nbsp;&nbsp;&nbsp;&nbsp;**if** $p=1$ **then**

14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Sample $\nu^0[k, d, q, p] \sim \text{Uniform}(1.5, 5)$.

15: &nbsp;&nbsp;&nbsp;&nbsp;**end if**

16: &nbsp;&nbsp;&nbsp;&nbsp;**if** $k > 1$ and $d > 0$ **then**

17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\nu^0[k, d, q, p] \leftarrow -\nu^0[k, d, q, p]$.

18: &nbsp;&nbsp;&nbsp;&nbsp;**end if**

19: **end for**

20: **Step 3: Initialization of gating parameters**

21: **for** $(p, k, d)$ over their respective ranges $([P] \times [K-1] \times \{\{0\} \cup [D_V]\})$ **do**

22: &nbsp;&nbsp;&nbsp;&nbsp;**if** $p=1, k=1, d=1$ **then**

23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\omega^0[p, k, d] \leftarrow 8$.

24: &nbsp;&nbsp;&nbsp;&nbsp;**else**

25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\omega^0[p, k, d] \leftarrow 0$.

26: **end if**

27: **end for**

### 5.3.1 Evaluation on higher dimensional synthetic data with perturbed ground truth initialization

Following the same procedure as in Section 5.2.1, we begin our evaluation under the setting using perturbed ground truth initialization. Specifically, the initialization is performed by perturbing the ground-truth parameters with Gaussian noise, and the latent series $s_i$ is initialized accordingly. The progression of the maximum log-likelihood values across epochs is illustrated in Figure 4. In this plot, the trajectory corresponding to the true log-likelihood is depicted in black, while the log-likelihood obtained via Polyak averaging is shown in red. The figure clearly demonstrates that, as the number of iterations increases, the proposed Incremental MM algorithm effectively maximizes the log-likelihood function. Moreover, the inclusion of Polyak averaging enhances the stability of the learning process
and leads to consistently superior performance in terms of log-likelihood values. As summarized in Table 3, the error values across all parameters remain low under all metrics, thereby confirming the ability of our algorithm to recover the true model parameters with high precision.

<table><caption>Table 3: The distance between each parameter and its ground truth using different metrics.</caption><thead><tr><th>Parameters</th><th>MSE</th><th>MAPE</th><th>NMRSE</th></tr></thead><tbody><tr><td>Omega</td><td>0,0051</td><td>0,4750</td><td>0,0087</td></tr><tr><td>Upsilon</td><td>0,1849</td><td>1,3460</td><td>0,0917</td></tr><tr><td>Sigma</td><td>0,0015</td><td>0,0183</td><td>0,0435</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 4: Maximum log-likelihood function with and without applying the Polyak averaging method over each iteration for higher dimensional synthetic data.</figcaption></figure>

### 5.3.2 Evaluation on higher dimensional synthetic data with data-driven initialization

In this section, we evaluate the performance of our method on high-dimensional synthetic data using data-driven initialization. This initialization strategy is generally more adaptable and better aligned with real-world scenarios, thus serving as a robust framework for assessing the resilience of clustering algorithms. To ensure consistency and comparability, we adopt the same configuration as described in Section 5.2.2, using the synthetic dataset introduced at the beginning of this section. As in Section 5.2.2, we utilize two complementary evaluation strategies. The first involves comparing the predicted means against the true mean trajectory of the data, while the second employs a train-test split to assess the generalization performance. The results of the first evaluation are summarized in Table 4, and the results of the train-test evaluation are presented in Table 5. These findings further demonstrate the robustness and superior accuracy of our proposed method in high-dimensional settings.

<table><caption>Table 4: Estimation errors with data-driven initialization for high dimensional simulated data</caption><thead><tr><th>Model</th><th>MSE↓</th><th>MAPE↓</th><th>NMRSE↓</th></tr></thead><tbody><tr><td>Inc. MM</td><td>2,809</td><td>0,197</td><td>0,081</td></tr><tr><td>SGD</td><td>3,278</td><td>0,243</td><td>0,073</td></tr><tr><td>Adam</td><td>3,402</td><td>0,245</td><td>0,074</td></tr><tr><td>Adam W</td><td>2,904</td><td>0,233</td><td>0,080</td></tr><tr><td>RMSProp</td><td>2,840</td><td>0,226</td><td>0,079</td></tr><tr><td>Sophia</td><td>2,912</td><td>0,256</td><td>0,077</td></tr></tbody></table>

<table><caption>Table 5: Prediction errors with data-driven initialization for high dimensional simulation data</caption><thead><tr><th>Model</th><th>MSE↓</th><th>MAPE↓</th><th>NMRSE↓</th></tr></thead><tbody><tr><td>Inc. MM</td><td>4,397</td><td>0,293</td><td>0,098</td></tr><tr><td>SGD</td><td>4,607</td><td>0,304</td><td>0,099</td></tr><tr><td>Adam</td><td>4,725</td><td>1,148</td><td>0,087</td></tr><tr><td>Adam W</td><td>4,434</td><td>0,439</td><td>0,090</td></tr><tr><td>RMSProp</td><td>4,478</td><td>0,359</td><td>0,106</td></tr><tr><td>Sophia</td><td>4,674</td><td>0,314</td><td>0,105</td></tr></tbody></table>
## 5.4 Experiments on real-world data

To further assess the robustness and practical applicability of our proposed model, we conducted empirical evaluations on two real-world datasets.

### 5.4.1 Dent maize genotypes dataset

The first dataset used in our study aims to elucidate the genetic and molecular mechanisms underlying drought-responsive traits in maize [90, 3]. It comprises 254 dent maize genotypes that capture the genetic diversity of the species. These genotypes were grown under two distinct watering conditions and evaluated for seven ecophysiological traits, although some measurements were incomplete [90].

In parallel, proteomic profiling was performed on leaf samples collected from the same plants, yielding quantification for 2,055 proteins. Of these, 973 proteins were measured using continuous data, while the remaining 1,082 were recorded as count data [3]. For our analysis, we focused on a subset of the dataset, specifically two ecophysiological traits (late leaf area and water use together with the 973 continuously quantified proteins, all collected under water deficit conditions).

After removing samples with missing values, the final dataset consisted of 233 maize genotypes (i.e., *N* = 233), each with measurements for the two selected traits (*Q* = 2) and 973 protein expression levels (*P* = 973). To mitigate numerical issues associated with high latent dimensions *P* and *Q*, we applied Lasso regression with feature selection (from scikit-learn), see also [34, 99], to identify the 11 most relevant features from *x*. These 11 features were used as inputs, and the first dimension of *y* was selected as the output for the experiment.

The dataset was split into a training set with 185 samples and a test set with 47 samples. We set both $D_V$ and $D_W$ to 1. Because the true underlying parameters are unknown in this real-world dataset, we adopted K-means for initialization for all the baseline. Each method was trained on the training set and evaluated on the test set using MSE, MAPE, and NRMSE as performance metrics. To further ensure robustness, we conducted 5-fold cross-validation, reporting the average performance over the folds along with the corresponding standard deviations to reflect result stability.

Learning rates for all solvers were finely tuned to ensure optimal performance. For initialization of $\omega$, we considered two strategies: (1) setting all entries to zero, and (2) applying logistic regression between the inputs $x$ and the labels generated via K-means. The results are summarized in Table 6, where each cell presents the mean metric value, with the standard deviation shown in parentheses. Incremental MM refers to the approach with zero-initialized $\omega$, while Incremental MM* denotes initialization via logistic regression. Overall, Incremental MM and Incremental MM* consistently outperform the baseline methods, even on real-world data, and maintains comparable standard deviations, highlighting both accuracy and stability.

### 5.4.2 Communities and Crime dataset

The second dataset we consider is the Communities and Crime dataset, available at [93]. This dataset integrates socio-economic data from the 1990 U.S. Census, law enforcement data from the 1990 Law Enforcement Management and Administrative Statistics survey, and crime statistics from the 1995 FBI Uniform Crime Reporting program. It comprises 1,994 instances, each corresponding to a distinct U.S. community, and includes 127 attributes capturing a range of factors such as demographics, economic indicators, housing conditions, and characteristics of local law enforcement agencies. The primary target variable is the per capita rate of violent crime, computed based on population size and the aggregated number of reported incidents of murder, rape, robbery, and assault.

Similar to the previous experiment, we employ Lasso regression with feature selection from scikit-learn, see also [34, 99], to identify the ten most informative attributes. These ten features are used as input variables $x$, alongside the output variable $y$, which corresponds to the per capita rate of violent crimes. The dataset is subsequently divided into a training set comprising 1600 samples and a test set containing 397 samples. We set $D_W = 1$ by default and evaluate the performance of Incremental MM under two configurations: one with $D_V = 1$ and another with $D_V = 2$. As the true underlying parameters are unknown in this context, we employ K-means initialization. To assess robustness, we apply five-fold cross-validation, reporting the mean and standard deviation of each metric across the folds.
To ensure fair comparisons, the learning rates of all solvers are carefully fine-tuned. Consistent with our previous experimental settings, we examine two strategies for initializing $\omega$: (1) setting all entries to zero, and (2) a “warm-start” approach where we cluster the first data batch and utilize the assignment labels as targets to solve a logistic regression problem to get the initial gating parameters. The final evaluation results are presented in Table 7.

<table><caption>Table 6: Prediction errors for the Dent maize genotypes dataset with data-driven initialization. The best result for each metric is shown in **bold**, and the second-best is **underlined**. *Incremental MM* denotes the case where $\omega$ is initialized using logistic regression.</caption><thead><tr><th>Model</th><th>MSE↓</th><th>MAPE↓</th><th>NMRSE↓</th></tr></thead><tbody><tr><td>Incremental MM</td><td><u>0.1190</u> (0.021)</td><td><u>0.1073</u> (0.013)</td><td><u>0.1421</u> (0.016)</td></tr><tr><td>Incremental MM*</td><td><b>0.1156</b> (0.017)</td><td><b>0.1065</b> (0.013)</td><td><b>0.1397</b> (0.017)</td></tr><tr><td>SGD</td><td>0.1332 (0.030)</td><td>0.1147 (0.018)</td><td>0.1514 (0.021)</td></tr><tr><td>Adam</td><td>0.1329 (0.026)</td><td>0.1164 (0.015)</td><td>0.1522 (0.024)</td></tr><tr><td>Adam W</td><td>0.1327 (0.026)</td><td>0.1163 (0.015)</td><td>0.1521 (0.024)</td></tr><tr><td>RMSProp</td><td>0.1324 (0.028)</td><td>0.1168 (0.017)</td><td>0.1514 (0.023)</td></tr><tr><td>Sophia</td><td>0.1332 (0.030)</td><td>0.1147 (0.018)</td><td>0.1514 (0.021)</td></tr></tbody></table>

<table><caption>Table 7: Prediction errors for the Communities and Crime dataset with data-driven initialization. The best result for each metric is highlighted in **bold**, while the second-best is **underlined**. Incremental MM$^i$ denotes the Incremental MM algorithm with $D_V$ set to $i$. Incremental MM* indicates the case where $\omega$ is initialized using logistic regression.</caption><thead><tr><th>Model</th><th>MSE↓</th><th>MAPE↓</th><th>NMRSE↓</th></tr></thead><tbody><tr><td>Incremental MM<sup>2</sup></td><td>0.0321 (0.004)</td><td><u>0.9537</u> (0.116)</td><td><b>0.1307</b> (0.016)</td></tr><tr><td>Incremental MM<sup>1</sup></td><td><b>0.0283</b> (0.003)</td><td>1.0440 (1.083)</td><td>0.1461 (0.012)</td></tr><tr><td>Incremental MM*</td><td><u>0.0316</u> (0.010)</td><td><b>0.9500</b> (0.229)</td><td>0.1594 (0.010)</td></tr><tr><td>SGD</td><td>0.0462 (0.012)</td><td>1.8288 (0.231)</td><td>0.1591 (0.010)</td></tr><tr><td>Adam</td><td>0.0314 (0.005)</td><td>2.2066 (0.914)</td><td><u>0.1413</u> (0.007)</td></tr><tr><td>Adam W</td><td>0.0315 (0.004)</td><td>1.5257 (0.507)</td><td>0.1417 (0.006)</td></tr><tr><td>RMSProp</td><td>0.0322 (0.005)</td><td>2.126 (1.184)</td><td>0.1415 (0.006)</td></tr><tr><td>Sophia</td><td>0.0309 (0.005)</td><td>1.9419 (1.059)</td><td>0.1430 (0.007)</td></tr></tbody></table>

# 6 Conclusion and future investigation

This work introduces a novel incremental stochastic MM algorithm that generalizes the incremental stochastic EM algorithm and provides a robust framework for handling high-volume and incremental data. By relaxing the stringent assumptions typically required by EM-based methods, such as convexity and explicit latent variable representations, our approach broadens the applicability of incremental stochastic estimation to a wider class of models. Notably, we demonstrate the effectiveness of the proposed method in fitting softmax-gated MoE models, a setting in which the incremental stochastic EM algorithm fails, achieving superior empirical performance over widely-used stochastic optimization algorithms. Looking forward, several promising directions for future research emerge. First, our incremental stochastic MM framework naturally accommodates extensions beyond Gaussian experts, allowing for more flexible modeling with non-Gaussian components. This generalization opens the door to integrating deep neural network architectures within each expert, enabling powerful hybrid models that combine statistical rigor with the representational capacity of deep learning. Second, while the current formulation processes one data point per iteration, extending the algorithm to support mini-batch updates could improve both stability and scalability, particularly in high-throughput or distributed environments. Developing a theoretically grounded mini-batch MM algorithm with provable convergence guarantees represents an important step toward scalable incremental stochastic learning. Third, addressing the challenges posed by high-dimensional datasets requires incorporating dimensionality reduction techniques or sparsity-inducing regularization within the MM framework to mitigate overfitting and enhance computational efficiency. Lastly, extending the algorithm to
handle mixed-type data, including both quantitative and qualitative covariates and parameters, would broaden its applicability to diverse real world problems. This could be achieved by designing expert components capable of modeling discrete structures or employing data transformation techniques compatible with the MM principle. Collectively, these avenues underscore the potential of the proposed incremental stochastic MM framework as a versatile and foundational tool in modern statistical and machine learning pipelines.

## Acknowledgments

TrungTin Nguyen, Hien Duy Nguyen, Florence Forbes, and Gersende Fort acknowledge funding from the Australian Research Council grant DP230100905, and from Inria Project WOMBAT. Christopher Drovandi and TrungTin Nguyen acknowledge funding from the Australian Research Council Centre of Excellence for the Mathematical Analysis of Cellular Systems (CE230100001). Christopher Drovandi was also supported by an Australian Research Council Future Fellowship (FT210100260).

# Supplementary Materials for
## “Revisiting Incremental Stochastic Majorization–Minimization
## Algorithms with Applications to Mixture of Experts”

In this supplementary material, we present the proofs of the main results in Section A and the technical proofs and results in Section B and Section C, respectively.

## A Proofs of main results

We start with a preliminary lemma

**Lemma 1.** Under Equation (2), P 1 and P 5, it holds,

$$1. \text{ for all } \tau \in \mathbb{T} \text{ and } z \in \mathbb{Z}, \nabla_{\theta} f(\cdot; z)|_{\theta=\tau} + \nabla\psi(\tau) - \nabla\phi(\tau)^{\top}\bar{S}(\tau; z) = 0.$$

$$2. \text{ for all } \tau \in \mathbb{T}, \nabla_{\theta} \mathbb{E}_{\pi} [f(\cdot; z)]|_{\theta=\tau} + \nabla\psi(\tau) - \nabla\phi(\tau)^{\top} \mathbb{E}_{\pi} [\bar{S}(\tau; z)] = 0.$$

*Proof.* Let $z \in \mathbb{Z}$, and $\tau \in \mathbb{T}$. Applying Equation (2) and Equation (3), we have

$$f(\cdot; z) - f(\tau; z) - \psi(\tau) + \psi(\cdot) + [\phi(\tau) - \phi(\cdot)]^{\top} \bar{S}(\tau; z) \le 0 \quad (20)$$

and this function is equal to zero at $\tau$. This implies that the derivative at $\tau$ is zero, thus yielding the first property.

We assumed that the objective function is finite on $\mathbb{T}$ (see Section 2.1), i.e. for every $\theta \in \mathbb{T}$, $\mathbb{E}_{\pi} [|f(\theta; z)|] < \infty$. Take the expectation w.r.t. $z$ in Equation (20); observing again that the function

$$\mathbb{E}_{\pi} [f(\cdot; z)] - \mathbb{E}_{\pi} [f(\tau; z)] - \psi(\tau) + \psi(\cdot) + [\phi(\tau) - \phi(\cdot)]^{\top} \mathbb{E}_{\pi} [\bar{S}(\tau; z)] \le 0$$

is non-positive on $\mathbb{T}$ and equal to zero at $\tau$, we conclude that its derivative is zero at $\tau$, which yields the second property. $\square$

### A.1 Proof of Proposition 1

Let $s^0 \in \mathbb{F}$. Apply Lemma 1 with $\tau \leftarrow \bar{\theta}(s^0)$. It then follows that

$$\nabla_{\theta} \mathbb{E}_{\pi} [f(\cdot; z)]|_{\theta=\bar{\theta}(s^0)} + \nabla\psi(\bar{\theta}(s^0)) - \nabla\phi(\bar{\theta}(s^0))^{\top} s^0 = 0,$$

which implies $\bar{\theta}(s^0) \in \mathbb{L}$ by P 4. Conversely, let $\theta^0 \in \mathbb{L}$; then by Lemma 1, we have

$$\nabla\psi(\theta^0) - \nabla\phi(\theta^0)^{\top}\mathbb{E}_{\pi}[\bar{S}(\theta^0; z)] = 0,$$

which, by P 3 and P 4, implies that $\theta^0 = \bar{\theta}(\mathbb{E}_{\pi}[\bar{S}(\theta^0; z)])$. By definition of $s^0$, this yields $\theta^0 = \bar{\theta}(s^0)$, from which it follows that $s^0 = \mathbb{E}_{\pi}[\bar{S}(\bar{\theta}(s^0); z)];$ i.e., $s^0 \in \mathbb{F}$.
## A.2 Proof of Proposition 2

Recall that

$$V(s) := \mathbb{E}_{\pi}[f(\bar{\theta}(s); z)], \quad \eta(s) := \mathbb{E}_{\pi}[\bar{S}(\bar{\theta}(s); z)] - s,$$

and that, by definition of $\bar{\theta}(s)$,

$$\bar{\theta}(s) := \arg\min_{\theta \in \mathbb{T}} h(s; \theta), \quad h(s; \theta) := -\psi(\theta) + \langle s, \phi(\theta) \rangle, \quad (21)$$

so that the first-order optimality condition reads

$$-\nabla\psi(\bar{\theta}(s)) + \nabla\phi(\bar{\theta}(s))^{\top}s = 0, \quad s \in \mathbb{S}. \quad (22)$$

(a) Differentiability of V. By Assumption 2 and P 5, V is continuously differentiable on $\mathbb{S}$.

(b) Sign of $\langle \nabla V(s), \eta(s) \rangle$ Using Lemma 1 with $\tau = \bar{\theta}(s)$ gives

$$\nabla_{\theta} \mathbb{E}_{\pi}[f(\bar{\theta}(s); z)] = -(\nabla\psi)(\bar{\theta}(s)) + (\nabla\phi)(\bar{\theta}(s))^{\top} \mathbb{E}_{\pi}[\bar{S}(\bar{\theta}(s); z)]. \quad (23)$$

Therefore, by Equation (22) and Equation (23),

$$\nabla_{\theta} \mathbb{E}_{\pi}[f(\cdot; z)]|_{\theta=\bar{\theta}(s)} = \{(\nabla\phi)(\bar{\theta}(s))\}^{\top} \eta(s). \quad (24)$$

Applying the chain rule to $V(s) = \mathbb{E}_{\pi}[f(\bar{\theta}(s); z)]$ and using Equation (23) and Equation (24) yields

$$\nabla V(s) = \nabla \bar{\theta}(s)^{\top} \nabla_{\theta} \mathbb{E}_{\pi}[f(\cdot; z)]|_{\theta=\bar{\theta}(s)} = \nabla \bar{\theta}(s)^{\top} \{(\nabla \phi)(\bar{\theta}(s))\}^{\top} \eta(s). \quad (25)$$

Next we differentiate the optimality condition in Equation (22) with respect to $s$. Writing $h(s; \theta)$ as in Equation (21), P 4 and Assumption 1 imply that, for each fixed $s \in \mathbb{S}$, the function $\theta \mapsto h(s; \theta)$ is twice continuously differentiable and has a unique minimizer at $\theta(s)$, with positive definite Hessian $H(s) := \nabla^2_{\theta}h(s; \theta)|_{\theta=\bar{\theta}(s)}$. Differentiating $s \mapsto \nabla_{\theta}h(s; \theta)|_{\theta=\bar{\theta}(s)} = 0$ with respect to $s$ gives

$$\nabla^2_{\theta}h(s; \theta)|_{\theta=\bar{\theta}(s)} \nabla \bar{\theta}(s) + \nabla_s[\nabla_{\theta}h(s; \theta)]|_{\theta=\bar{\theta}(s)} = 0.$$

Using $\nabla_{\theta}h(s; \theta) = -\nabla\psi(\theta) + \nabla\phi(\theta)^{\top}s$ and the fact that $\nabla_s s^{\top}$ is the identity matrix, we obtain

$$\nabla_s[\nabla_{\theta}h(s; \theta)]|_{\theta=\bar{\theta}(s)} = \nabla_{\theta}\phi(\bar{\theta}(s))^{\top}.$$

Hence

$$\nabla \bar{\theta}(s)^{\top} = -\nabla \phi(\bar{\theta}(s)) H(s)^{-1}. \quad (26)$$

Plugging Equation (26) into Equation (25) yields

$$\begin{aligned} \langle \nabla V(s), \eta(s) \rangle &= \eta(s)^{\top} \nabla V(s) \\ &= -\eta(s)^{\top} (\nabla \phi)(\bar{\theta}(s)) H(s)^{-1} \{(\nabla \phi)(\bar{\theta}(s))\}^{\top} \eta(s) \\ &= -\left(\{(\nabla \phi)(\bar{\theta}(s))\}^{\top} \eta(s)\right)^{\top} H(s)^{-1} \left(\{(\nabla \phi)(\bar{\theta}(s))\}^{\top} \eta(s)\right). \end{aligned} \quad (27)$$

By P 4, $H(s)$ is positive definite, so its inverse is also positive definite. Consequently, $\langle \nabla V(s), \eta(s) \rangle \le 0$ for all $s \in \mathbb{S}$, as claimed.

If $s \in \mathbb{F}$, then $\eta(s) = 0$ by definition, and $\langle \nabla V(s), \eta(s) \rangle = 0$ follows immediately.

## A.3 Proof of Theorem 1

**Definitions.** Define the filtration $\{\mathcal{F}_n, n \ge 0\}$ by

$$\mathcal{F}_n := \sigma(s_0, z_1, \dots, z_n).$$

Then $s_n \in \mathcal{F}_n$ for all $n \ge 0$. Define the martingale increment

$$\xi_{n+1} := \bar{S}(\bar{\theta}(s_n); z_{n+1}) - \mathbb{E}[\bar{S}(\bar{\theta}(s_n); z_{n+1}) | \mathcal{F}_n].$$

Note that for any positive measurable function $\tau$, $\mathbb{E}[\tau(s_n, z_{n+1}) | \mathcal{F}_n] = \mathbb{E}_{\pi}[\tau(s, z)]|_{s=s_n}$.
Our algorithm Equation (6) can be rewritten as:

$$s_{n+1} = s_n + \gamma_{n+1} \{\eta(s_n) + \xi_{n+1}\}. \quad (28)$$

The proof consists in checking the assumptions of [1, Theorem 2.3], which provide a general result for the almost-sure convergence of a stable Stochastic Approximation algorithm.

**Checking the assumption A1 in [1].** It is verified under Assumption 4, Assumption 5, Assumption 6 and by Proposition 2.

**A stable algorithm w.p.1.** By Assumption 7, with probability one, the sequence {$s_n, n \ge 0$} remains in a compact set $\mathbb{K}$ of $\mathbb{S}$ ($\mathbb{K}$ is path-dependent). Hence, the stability condition is verified. Furthermore, $\mathbb{K} \cap \mathbb{F} \neq \emptyset$.

**The conditions on the step size sequence.** By Assumption 8, the step size sequence is monotone nonincreasing, $\lim_k \gamma_k = 0$ and $\sum_k \gamma_k = +\infty$. For any $\lambda_0 > 0$, there exists $K$ such that $\gamma_K \le \lambda_0$; since we are interested in the convergence when $k \to \infty$, we can assume without loss of generality that the algorithm is studied for all $k \ge K$.

**Show that with probability one,** $\limsup_n \sup_{k \ge n} \|\sum_{j=n}^k \gamma_j \xi_j\| = 0$. We prove the equivalent property: with probability one, $\lim_n M_n$ exists (and is finite), where $M_n := \sum_{j=1}^n \gamma_j \xi_j$.

Let {$\mathbb{K}_q, q \ge 1$} be a sequence of compact sets of $\mathbb{S}$ such that $\mathbb{K}_q \subseteq \mathbb{K}_{q+1}$ and $\cup_{q \in \mathbb{N}} \mathbb{K}_q = \mathbb{S}$; for example, set $\mathbb{K}_q := \{s \in \mathbb{S} : \|s\| \le q, d(s, \mathbb{S}^c) \ge \frac{1}{q}\}$. By Assumption 7,

$$\mathbb{P}\left(\bigcup_{q \in \mathbb{N}} \{\forall n \ge 0, s_n \in \mathbb{K}_q\}\right) = 1. \quad (29)$$

We write

$$\begin{aligned} \mathbb{P}\left(\lim_n M_n \text{ exists}\right) &= \mathbb{P}\left(\left\{\lim_n M_n \text{ exists}\right\} \cap \bigcup_{q \in \mathbb{N}} \{\forall n \ge 0, s_n \in \mathbb{K}_q\}\right) \\ &= \lim_q \uparrow \mathbb{P}\left(\left\{\lim_n M_n \text{ exists}\right\} \cap \{\forall n \ge 0, s_n \in \mathbb{K}_q\}\right) \\ &= \lim_q \uparrow \mathbb{P}\left(\left\{\lim_n M_{n,q} \text{ exists}\right\} \cap \{\forall n \ge 0, s_n \in \mathbb{K}_q\}\right) \end{aligned} \quad (30)$$

where $M_{n,q} := \sum_{j=1}^n \gamma_j \xi_j \mathbf{1}_{s_{j-1} \in \mathbb{K}_q}$. For fixed $q \ge 0$, the random variables {$M_{n,q}, n \ge 1$} is a $\mathcal{F}_n$-adapted martingale sequence, upon noting that by definition of $\xi_n$, $\mathbb{E}[\xi_n \mathbf{1}_{s_{n-1} \in \mathbb{K}_q} | \mathcal{F}_{n-1}] = 0$. By Assumption 3, it is a square-integrable martingale:

$$\mathbb{E}[\|M_{n,q}\|^2] = \sum_{j=1}^{n} \gamma_j^2 \mathbb{E}[\|\xi_j\|^2 \mathbf{1}_{s_{j-1} \in \mathbb{K}_q}] \le \sup_{s \in \mathbb{K}_q} \mathbb{E}_{\pi} [\|\bar{\mathcal{S}}(\bar{\theta}(s); z)\|^2] \sum_{j=1}^{n} \gamma_j^2 < \infty.$$

Therefore, by [37, Theorem 2.15]

$$\{\lim_n M_{n,q} \text{ exists}\} \supset \{\sum_{j \ge 1} \gamma_j^2 \mathbb{E}[\|\xi_j\|^2 \mathbf{1}_{s_{j-1} \in \mathbb{K}_q} | \mathcal{F}_{j-1}] < \infty\}. \quad (31)$$

Using again Assumption 3 and since $\sum_j \gamma_j^2 < \infty$ (see Assumption 8), we have

$$\mathbb{P}\left(\sum_{j \ge 1} \gamma_j^2 \mathbb{E}[\|\xi_j\|^2 \mathbf{1}_{s_{j-1} \in \mathbb{K}_q} | \mathcal{F}_{j-1}] < \infty\right) = 1. \quad (32)$$

By Equation (30), Equation (31), Equation (32) and Equation (29), we have

$$\mathbb{P}\left(\lim_n M_n \text{ exists}\right) \ge \lim_q \uparrow \mathbb{P}(\forall n, s_n \in \mathbb{K}_q) = 1.$$

The proof is concluded.
**Convergence of the sequence** {$s_n, n \ge 0$}. By [1, Theorem 2.3], with probability one, $\limsup_n d(s_n, \{\langle \nabla V(s), \eta(s) \rangle = 0\} \cap \mathbb{K}) = 0$. By Assumption 4, this yields $\limsup_n d(s_n, \mathbb{F} \cap \mathbb{K}) = 0$.

**Convergence of the sequence** {$\theta_n, n \ge 0$}. We prove an inclusion of sets: $\{\limsup_n d(s_n, \mathbb{F}) = 0\} \subseteq \{\limsup_n d(\theta_n, \mathbb{L}) = 0\}$; since the set on the LHS is almost-sure, this will imply that the set on the RHS is almost-sure. Let $\omega \in \{\omega : \limsup_n d(s_n(\omega), \mathbb{F}) = 0\}$. Set $\theta_n(\omega) := \theta(s_n(\omega))$. Let $\{\theta_{n_q}(\omega), q \ge 0\}$ be a converging subsequence with limiting value $\theta_\star$; we have to prove that $\theta_\star \in \mathbb{L}$.

From the subsequence {$s_{n_q}(\omega), q \ge 0$}, we extract a converging subsequence {$s_{n_{\rho(q)}}(\omega), q \ge 0$} with limit $s_\star$; by Assumption 7, such a subsequence exists and $s_\star \in \mathbb{S}$. Since $\omega \in \{\limsup_n d(s_n, \mathbb{F}) = 0\}$ and $\mathbb{F}$ is closed (as a consequence of Assumption 4), then $s_\star \in \mathbb{F}$. By Assumption 2, $\theta$ is continuous so that $\theta_\star = \lim_q \bar{\theta}(s_{n_q}(\omega)) = \lim_q \bar{\theta}(s_{n_{\rho(q)}}(\omega)) = \bar{\theta}(s_\star)$. Therefore, by Proposition 1, $\theta_\star \in \mathbb{L}$. This concludes the proof.

## A.4 Proof of Theorem 2

### A.4.1 Proof of exponential family majorizer surrogate via Equation (18) in Proposition 3

Denote:

$$ \underline{g}(\omega; \mathbf{x}) = \log \left( 1 + \sum_{k=1}^{K-1} \exp(w_k(\mathbf{x})) \right) \quad \text{where} \quad w_k(\mathbf{x}) = \sum_{d=0}^{D_W} \omega_{k,d}^{\top} \mathbf{x}^d. $$

We have $\underline{g}$ is $C^2$, then Taylor's theorem with Lagrange remainder implies that there exists some $\alpha \in (0, 1)$ such that

$$ \underline{g}(\omega; \mathbf{x}) - \underline{g}(\omega_n; \mathbf{x}) = (\omega - \omega_n)^{\top} \nabla \underline{g}(\omega_n; \mathbf{x}) + \frac{1}{2} (\omega - \omega_n)^{\top} \nabla^2 \underline{g}(\omega_n + \alpha(\omega - \omega_n); \mathbf{x}_n) (\omega - \omega_n). \quad (33) $$

**Calculate the gradient of g.** For some $k \in [K-1]$, $p \in [P]$ and $d \in \{0, 1, \dots, D_W\}$, the gradient can be calculated as follows:

$$ \nabla_{\omega_{k,d,p}} \underline{g} = \frac{x_p^d \exp \left( \sum_{d=0}^{D_W} \omega_{k,d}^{\top} \mathbf{x}^d \right)}{1 + \sum_{l=1}^{K-1} \exp \left( \sum_{d=0}^{D_W} \omega_{l,d}^{\top} \mathbf{x}^d \right)} = x_p^d g_k(\mathbf{w}(\mathbf{x})). $$

Hence, the gradient of $\underline{g}$ is given as:

$$ \nabla \underline{g} = [x_1^0 g_1(\omega(\mathbf{x})), \dots, x_1^{D_W} g_1(\omega(\mathbf{x})), \dots, x_P^0 g_{K-1}(\omega(\mathbf{x})), \dots, x_P^{D_W} g_{K-1}(\omega(\mathbf{x}))]^{\top} = \hat{g} \otimes \hat{\mathbf{x}}, $$

where $\hat{\mathbf{x}} = [x_1^0, \dots, x_1^{D_W}, \dots, x_P^0, \dots, x_P^{D_W}]^\top$, $\hat{g} = [g_1(\mathbf{w}(\mathbf{x})), \dots, g_{K-1}(\mathbf{w}(\mathbf{x}))]^\top$ and $\otimes$ is the Kronecker product.
**Calculate for the Hessian the gradient of g.** For some $k_1, k_2 \in [K-1]$, $d_1, d_2 \in \{0, 1, \dots, D_W\}$ and $p_1, p_2 \in [P]$, the Hessian can be calculated as follow:

$$ \nabla_{\omega_{k_1, d_1, p_1} \omega_{k_2, d_2, p_2}} \underline{g} = \frac{-x_{p_1}^{d_1} \exp(w_{k_1}(\mathbf{x})) x_{p_2}^{d_2} \exp(w_{k_2}(\mathbf{x}))}{\left[1 + \sum_{l=1}^{K-1} \exp(w_l(\mathbf{x}))\right]^2} = -x_{p_1}^{d_1} x_{p_2}^{d_2} \mathbf{g}_{k_1}(\mathbf{w}(\mathbf{x})) \mathbf{g}_{k_2}(\mathbf{w}(\mathbf{x})) \text{ for } k_1 \neq k_2, $$

$$ \begin{aligned} \nabla_{\omega_{k_1, d_1, p_1} \omega_{k_1, d_2, p_2}} \underline{g} &= \frac{u_{p_1}^{d_1} x_{p_2}^{d_2} \exp(w_{k_1}(\mathbf{x})) \left( \sum_{l=1}^{K} \exp(w_l(\mathbf{x})) - \exp(w_{k_1}(\mathbf{x})) \right)}{\left[1 + \sum_{l=1}^{K-1} \exp(w_l(\mathbf{x}))\right]^2} \\ &= u_{p_1}^{d_1} x_{p_2}^{d_2} \mathbf{g}_{k_1}(\mathbf{w}(\mathbf{x})) (1 - \mathbf{g}_{k_1}(\mathbf{w}(\mathbf{x}))). \end{aligned} $$

Hence, the Hessian is given as:

$$ \nabla^2 \underline{g} = (\mathrm{diag}(\hat{g}) - \hat{g} \hat{g}^\top) \otimes \hat{x} \hat{x}^\top. $$

**Lemma 2** (See, e.g., [4]). If $A \succeq B$ then for symmetric, nonnegative definite $C$:

$$ A \otimes C \succeq B \otimes C. $$

**Proposition 5.** Let $\pi_k$ be be real values in $(0, 1)$ for $k \in [K+1]$ and $\sum_{k=1}^{K+1} \pi_k = 1$. Denote $\hat{\pi} = [\pi_1, \pi_2, \dots, \pi_K]^\top$. We have:

$$ \mathrm{diag}(\hat{\pi}) - \hat{\pi} \hat{\pi}^\top \le \left( \frac{3}{4} \mathbf{I}_{K-1} - \frac{\mathbf{1}_{K-1} \mathbf{1}_{K-1}^\top}{2(K-1)} \right). $$

The proofs of Lemma 2 are provided in Section B.1, while the proofs of Proposition 5 are detailed in Section B.2.

**Proof of exponential family majorizer surrogate via Equation (18).** Applying Lemma 2 and Proposition 5 to Equation (33) we receive the majorization:

$$ \underline{g}(\omega; x_n) \le \underline{g}(\omega_n; x_n) + (\omega - \omega_n)^\top \nabla \underline{g}(\omega_n; x_n) + \frac{1}{2} (\omega - \omega_n)^\top \mathbf{B}_{n,K} (\omega - \omega_n). $$

Hence, given $\omega := (\omega_{k,d,p} |_{p \in [P], k \in [K-1], d \in \{0, \dots, D_W\}})$, the overall majorization function is defined:

$$ \begin{aligned} - \log s_{\theta}(\mathbf{y}_n | \mathbf{x}_n) &\le \sum_{k=1}^{K} \tau_{n,k} \log(\tau_{n,k}) - \sum_{k=1}^{K} \tau_{n,k} w_k(\mathbf{x}_n) - \sum_{k=1}^{K} \tau_{n,k} [\log(\mathcal{N}(\mathbf{y}_n; \boldsymbol{\mu}_k(\mathbf{x}_n), \boldsymbol{\Sigma}_k))] \\ &\quad + \left\{ \underline{g}(\omega_n; \mathbf{x}_n) + (\omega - \omega_n)^\top \nabla \underline{g}(\omega_n; \mathbf{x}_n) + \frac{1}{2} (\omega - \omega_n)^\top \mathbf{B}_{n,K} (\omega - \omega_n) \right\}, \\ &= \sum_{k=1}^{K} \tau_{n,k} \log(\tau_{n,k}) + \underline{g}(\omega_n; \mathbf{x}_n) - \omega_n^\top \nabla \underline{g}(\omega_n; \mathbf{x}_n) - \sum_{k=1}^{K} \tau_{n,k} w_k(\mathbf{x}_n) \\ &\quad + \left\{ \omega^\top \nabla \underline{g}(\omega_n; \mathbf{x}_n) + \frac{1}{2} (\omega - \omega_n^\top) \mathbf{B}_{n,K} (\omega - \omega_n) \right\} \\ &\quad - \sum_{k=1}^{K} \tau_{n,k} [\log(\mathcal{N}(\mathbf{y}_n; \boldsymbol{\mu}_k(\mathbf{x}_n), \boldsymbol{\Sigma}_k))] =: - \log[g(\boldsymbol{\theta}, \mathbf{x}_n, \mathbf{y}_n; \boldsymbol{\theta}_n)]. \end{aligned} $$

### A.4.2 Verification of assumptions in Theorem 2

To verify the model-specific assumptions P 1 to 5 and Assumptions 1 to 3, we will assume that $x$ and $y$ are bounded by $u^\#$ and $y^\#$, respectively.
**Verification of P 1.** For the problem of softmax-gated MoE, we consider the majorizer surrogate $g$ to have the following form that belongs to the exponential family

$$ g(\theta, z; \tau) := -\psi(\theta) + \langle \bar{S}(\tau; z), \phi(\theta) \rangle, \quad (34) $$

where,

$$ \phi(\theta) := \begin{bmatrix} \omega \\ \operatorname{vec}(\omega\omega^\top) \\ \kappa \\ \zeta \\ \Delta \\ \Sigma \end{bmatrix}, \quad \bar{S}(\tau; z) = \begin{bmatrix} -\xi + \nabla f(\tau) - B\tau \\ \frac{1}{2}\operatorname{vec}(B) \\ \tau \otimes y^2 \\ -2\tau \otimes (y \otimes r) \\ \tau \otimes (\operatorname{vec}(rr^\top) \otimes 1_Q) \\ \tau \otimes 1_Q \end{bmatrix}, \quad \psi(\theta) = C. $$

Here, $C$ represents the constants that are independent of $\omega$ and $\kappa, \zeta, \Delta, \bar{\Sigma}$ and $\tau$ are defined as follows:

$$
\begin{align*}
\bar{\Sigma} &= [\log \sigma_{1,1}^2, \dots, \log \sigma_{1,Q}^2, \log \sigma_{2,1}^2, \dots, \log \sigma_{K,Q}^2]^\top, \quad \kappa = \left[ \frac{1}{\sigma_{1,1}^2}, \dots, \frac{1}{\sigma_{1,Q}^2}, \frac{1}{\sigma_{2,1}^2}, \dots, \frac{1}{\sigma_{K,Q}^2} \right]^\top, \\
\zeta &= \operatorname{vec} \left( \left[ \frac{\mathbf{\Upsilon}_{1,1}^\top}{\sigma_{1,1}^2}, \dots, \frac{\mathbf{\Upsilon}_{1,Q}^\top}{\sigma_{1,Q}^2}, \frac{\mathbf{\Upsilon}_{2,1}^\top}{\sigma_{2,1}^2}, \dots, \frac{\mathbf{\Upsilon}_{K,Q}^\top}{\sigma_{K,Q}^2} \right] \right), \quad \tau = [\tau_1, \tau_2, \dots, \tau_K]^\top, \\
\Delta &= \operatorname{vec} \left( \frac{\operatorname{vec}(\mathbf{\Upsilon}_{1,1}^\top \mathbf{\Upsilon}_{1,1})}{\sigma_{1,1}^2}, \dots, \frac{\operatorname{vec}(\mathbf{\Upsilon}_{1,Q}^\top \mathbf{\Upsilon}_{1,Q})}{\sigma_{1,Q}^2}, \frac{\operatorname{vec}(\mathbf{\Upsilon}_{2,1}^\top \mathbf{\Upsilon}_{2,1})}{\sigma_{2,1}^2}, \dots, \frac{\operatorname{vec}(\mathbf{\Upsilon}_{K,Q}^\top \mathbf{\Upsilon}_{K,Q})}{\sigma_{K,q}^2} \right).
\end{align*}
$$

Furthermore, since $\phi, \psi$ and $\bar{S}$ are composed of elementary functions applied component-wise to their parameters, they are trivially continuously differentiable and measurable due to the smoothness and measurability properties of elementary functions.

**Verification of P 2.** For $B_{n,K}$ defined in Equation (16), it can be proven that for any eigenvalue $\lambda_{B_{n,K}}$ of $B_{n,K}$, $\lambda_{B_{n,K}}$ is bounded in the interval $(\frac{1}{M_0}, M_0)$ for some positive real number $M_0$.

Additionally, since $\lambda_{B_{n,K}}^{\max} = \frac{1}{\lambda_{B_{n,K}}^{\min}} < \frac{1}{M_0}$ and $\lambda_{B_{n,K}}^{\min} = \frac{1}{\lambda_{B_{n,K}}^{\max}} > \frac{1}{M_0}$. Thus, for all eigenvalues $\lambda_{B_{n,K}}^{-1}$ of $B_{n,K}^{-1}$, $\lambda_{B_{n,K}}^{-1} \in (\frac{1}{M_0}, M_0)$.

Knowing this value $M_0$, we can define $\mathcal{S}$ as follows:

$$
\begin{align}
\mathcal{S} ={}& \{(s_1, \operatorname{vec}(S_2), s_3, s_4, s_5, s_6) : s_1 \in \mathbb{R}^{(K-1)P(D_w+1)}; \nonumber \\
& \qquad S_2 \in \mathbb{R}^{(K-1)P(D_w+1) \times (K-1)P(D_w+1)}, S_2 \succeq 0, \lambda \in (\frac{1}{M_0}, M_0) \text{ for all eigenvalues } \lambda \text{ of } S_2; \nonumber \\
& \qquad s_3 \in \mathbb{R}^{KQ}; \quad s_4 \in \mathbb{R}^{KPQ(D_v+1)}; \quad s_6 \in \mathbb{R}_+^{KQ}; \nonumber \\
& \qquad s_5 \in \mathbb{R}^{KQ(P(D_v+1))^2}, \operatorname{mat}(s_{5,k:q}) \succeq 0 \text{ for all } k \in [K], q \in [Q]\}. \tag{35}
\end{align}
$$

It is trivial that $\mathcal{S}$ denoted above satisfies the requirements of P 2.

**Verification of P 3.** $\bar{S}(\theta, z)$ is composed of elementary functions applied component-wise to $x$ and $y$, hence its expectation exist and lies in $\mathcal{S}$. Furthermore, since $x$ and $y$ are bounded, it can be shown that $\tau, B, r, \xi$ and $\nabla g$ are all bounded. Thus, the expectation $\mathbb{E}_\pi[\bar{S}(\theta; z)]$ exists, lies within $\mathcal{S}$, for any $\theta \in \mathbb{T}$.

**Verification of P 4.** This assumption is proven in Section A.5.3 where we have shown that $\omega, \mathbf{\Upsilon}, \sigma^2$ of the next iteration are the unique solution to the minimization problem of the current iteration. Indeed, since the objective is strictly convex in each block, the minimizer is unique and measurable.

**Verification of P 5.** Since $f(\theta; z)$ is continuously differentiable in $\theta$ for all $z$ and dominated by an integrable envelope (due to bounded covariates and responses), differentiation under the expectation is valid, hence $\theta \mapsto \mathbb{E}[f(\theta; z)]$ is $C^1$ on $\mathbb{T}$.
**Verification of Assumption 1.** Since we have $\omega \in \mathbb{R}^{P \times K \times (D_w+1)}$, $\mathbf{\Upsilon} \in \mathbb{R}^{K \times (D_v+1) \times Q \times P}$ and $\sigma^2 \in \mathbb{R}_+^{K \times Q \times Q}$, thus $\mathbb{T} = \mathbb{R}^{P \times K \times (D_w+1)} \times \mathbb{R}^{K \times (D_v+1) \times Q \times P} \times \mathbb{R}_+^{K \times Q \times Q}$ is a convex and open set. Furthermore, similar to the verification of P 1, we can show that $\psi$ and $\phi$ is twice continuously differentiable throughout $\mathbb{T}$.

**Verification of Assumption 2.** From Section A.5.3, we have the function $\bar{\theta}(s)$ is defined as:

$$ \bar{\theta}(s) = \left[ - \frac{(2 \mathrm{mat}(s_2))^{-1} s_1}{s_{3,kq} + s_{4,k:q}^{\top} \mathbf{\Upsilon}_{k,q:}^{\top} + s_{5,k:q}^{\top} \mathrm{vec}(\mathbf{\Upsilon}_{k,q:}^{\top} \mathbf{\Upsilon}_{k,q:})} \right]^{-1} s_1, \quad s = (s_1, s_2, s_3, s_4, s_5, s_6). $$

It is trivial to show that $\bar{\theta}(s)$ is continuously differentiable for $s_1, s_3, s_4, s_6$. For $s_2, s_5$, we have $v \to \mathrm{mat}(v) + \mathrm{mat}(v)^{\top}$ is an affine transformation, thus it is continuously differentiable. Furthermore, from Section A.5.3, we know that $2 \mathrm{mat}(s_2)$ and $\mathrm{mat}(s_5) + \mathrm{mat}(s_5)^{\top}$ are positive definite, thus their inverts are non-singular, thus their inverts are continuously differentiable. This concludes our proof.

**Verification of Assumption 3.** By applying Lemma 3 in addition with the fact that $\kappa, \zeta, \Delta, \bar{\Sigma}$ (defined in Section A.5.2) are all bounded. We can conclude that for all compact set $\mathbb{K} \subset \mathbb{S}$, $\mathcal{S}(\theta(s), x)$ is bounded by a number $\beta$. Hence,

$$ \sup_{s \in \mathbb{K}} (\mathbb{E}_{\pi} [ |\bar{\mathcal{S}}(\theta(s); x)|^2 ] ) < \beta^2 < \infty. $$

**Lemma 3.** Let $\mathbb{K}$ be a compact subset of $\mathbb{S}$. Then for all $s_1, s_2 \in \mathbb{K}$, we have $(2 \mathrm{mat}(s_2))^{-1} s_1$ is bounded.

*Proof of Lemma 3.* Since $s_2 \in \mathbb{K} \subset \mathbb{S}$, we have $\mathrm{mat}(s_2)$ eigenvalues are bounded in $(\frac{1}{M_0}, M_0)$. This leads to $\mathrm{mat}(s_2)$ being bounded. Furthermore, we have $s_1 \in \mathbb{K}$ and $\mathbb{K}$ is compact, so $s_1$ is also bounded. Combining these two reasons we conclude the proof. $\square$

**Verification of Assumption 7.** First, notice that since $x$ and $y$ are bounded, $\kappa, \zeta, \Delta, \bar{\Sigma}$ are also bounded. Thus by applying Lemma 4, we can show that $s_{i,n}$ is bounded for all $i \in \{2, \dots, 6\}$. Additionally, since $\kappa, \zeta, \Delta, \bar{\Sigma}$ are bounded, due to how we choose $\mathbb{S}$, for all $i \in \{2, \dots, 6\}$, there exist a compact set $\mathcal{J}_i \subset \mathbb{S}$ that satisfy the condition of Lemma 5 for sequence $s_i$. Hence, by applying Lemma 5, we can show that

$$ \liminf d(s_{i,n}, \mathbb{S}^c) > 0, \quad \forall i \in \{2, \dots, 6\}. $$

Since, $x$ and $y$ are taken from compact sets, we can apply Lemma 5 for $s_i$, $\forall i \in \{2, \dots, 6\}$, and thus having

$$ \liminf d(s_{i,n}, \mathbb{S}^c) > 0, \quad \forall i \in \{2, \dots, 6\}. $$

**Lemma 4.** Given a sequence of vectors $(v_n)_{n=0}^{+\infty}$ satisfying

$$ v_{n+1} = v_n + \gamma_n(f(n) - v_n), $$

where $0 < \gamma_n < 1$. If $f$ is bounded (element-wise), then, $v_n$ is bounded (element-wise).

**Lemma 5.** Let $\mathbb{H}$ be a convex, open set. Given a sequence of vectors $(v_n)_{n=0}^{+\infty}$ in $\mathbb{H}$ satisfying

$$ v_{n+1} = v_n + \gamma_n(f(n) - v_n), $$

where $0 < \gamma_n < 1$. Then, if there exist a compact set $\mathcal{J} \subset \mathbb{H}$ such that $f(n) \in \mathcal{J}$ $\forall n$, we have

$$ \liminf \{d(v_n, \mathbb{S}^c)\} > 0. $$

*Proof of Lemma 4.* Let say $f$ is bounded (element-wise) by $\alpha$. Denote $\alpha_0 = \max(\alpha, t_0)$ (element-wise). We will show by induction that $v_n$ is upper bounded by $\alpha_0$ (element-wise). We have $t_0 \le \alpha_0$. Let say for some $n$, we have $v_n \le \alpha_0$ (element-wise), then:

$$ v_{n+1} = (1 - \gamma_{n+1})v_n + \gamma_{n+1}f(n) \le (1 - \gamma_{n+1})\alpha_0 + \gamma_{n+1}\alpha_0 = \alpha_0. $$

Hence, by induction, we have shown that $v_n$ is upper bounded. Doing similarly for the lower bound, we complete the proof. $\square$
**Proof of Lemma 5.** Let $\hat{\mathcal{J}}$ be the smallest compact, convex set that contains $\mathcal{J} \cup \{t_0\}$ (i.e $\hat{\mathcal{J}}$ is the convex hull of $\mathcal{J} \cup \{t_0\}$). We know that $t_0 \in \mathbb{H}$ and $\mathcal{J} \subseteq \mathbb{H}$, hence $\mathcal{J} \cup \{t_0\} \subseteq \mathbb{H}$. Since $\mathbb{H}$ is convex and open, we have $\hat{\mathcal{J}} \subset \mathbb{H}$. We will show by induction that $v_n \in \hat{\mathcal{J}}$ for all $n$. We have $t_0 \in \hat{\mathcal{J}}$. Assume that $v_n \in \hat{\mathcal{J}}$, we have:

$$v_{n+1} = (1 - \gamma_n)v_n + \gamma_n f(n)$$

is a convex combination of $v_n$ and $f(n)$. But since both $v_n$ and $f(n)$ are elements of $\hat{\mathcal{J}}$, $v_{n+1} \in \hat{\mathcal{J}}$. Hence, by induction we have $v_n \in \hat{\mathcal{J}} \quad \forall n$. Lastly, since $\hat{\mathcal{J}}$ is a compact subset of the open set $\mathbb{S}$, there exist $\varepsilon > 0$ such that for all $a \in \hat{\mathcal{J}}$ and $b \in \mathbb{S}^c$,

$$d(a, b) > \varepsilon.$$

Combining this with the previous proven facts, we yield

$$\liminf\{d(v_n, \mathbb{S}^c)\} > \varepsilon > 0.$$

This completes the proof.

<figure><img src="image_7.png" /></figure>

**Verification of Assumption 8.** By taking $\gamma_n = \gamma_0 n^{-\alpha}$ with $\alpha \in (\frac{1}{2}, 1]$ and $\gamma_0 \in (0, 1)$, all requirements of Assumption 8 are met.

**Remark 2.** Although formal proofs for Assumption 7 (in the case of $s_{n,1}$) and Assumptions 4 and 5 remain elusive, we adopt these as working assumptions for the MoE model. Nonetheless, our experiments on both synthetic and real-world datasets consistently demonstrate that our proposed method achieves convergence reliably, with a faster convergence rate and higher accuracy compared to the baseline models.

## A.5 Derivation of the Algorithm 2

### A.5.1 Construction of the incremental stochastic MM algorithm for $O_{gate}$

$$O_{gate}(\omega, \theta_n) = -\sum_{k=1}^{K} \tau_k w_k(x) + \left\{ \omega^\top \nabla \underline{g}(\omega_n; x_n) + \frac{1}{2} (\omega - \omega_n)^\top B_{n,K} (\omega - \omega_n) \right\}.$$

Let us denote vector:

$$\xi = [\tau_1 x_1^0, \tau_1 x_1^1, \dots, \tau_1 x_1^{D_W}, \dots, \tau_{K-1} x_1^{D_W}, \dots, \tau_{K-1} x_P^{D_W}]^\top.$$

Hence:

$$O_{gate}(\omega, \theta_n) = (-\xi^\top + \nabla \underline{g}(\omega_n; x_n)^\top - \omega_n^\top B) \omega + \frac{1}{2} \omega^\top B \omega + C, \\ B = \left( \frac{3}{4} I_{K-1} - \frac{1}{2(K-1)} \mathbf{1}_{K-1} \right) \otimes \hat{x}\hat{x}^\top,$$

where $C$ represents the constants that are independent of $\omega$.

Here the incremental stochastic MM algorithm can be written as:

$$\phi(\omega) := \left[ \frac{\omega}{\mathrm{vec}(\omega\omega^\top)} \right], \quad \bar{S}(\tau; z) = \left[ \frac{-\xi + \nabla \underline{g}(\tau) - B\tau}{\frac{1}{2} \mathrm{vec}(B)} \right], \quad \psi(\omega) = C.$$

### A.5.2 Construction of the incremental stochastic MM algorithm for $O_{expert}$

Recall that the expert function to be minimized is given by:

$$O_{expert}(\mathbf{\Upsilon}, \sigma^2, \theta) = \sum_{k=1}^{K} \left[ \tau_k (\mathbf{y} - v_k(\mathbf{x}))^\top \Sigma_k^{-1} (\mathbf{y} - v_k(\mathbf{x})) + \tau_k \log|\Sigma_k| \right]^\top \\ = \sum_{k=1}^{K} \left[ \tau_k (\mathbf{y} - \mathbf{\Upsilon}_k \mathbf{r})^\top \Sigma_k^{-1} (\mathbf{y} - \mathbf{\Upsilon}_k \mathbf{r}) + \tau_k \log|\Sigma_k| \right]^\top.$$
Where $\mathbf{\Upsilon}_k = [\mathbf{\Upsilon}_{k,0}, \dots, \mathbf{\Upsilon}_{k,D_V}]$ is a $Q \times P(D_V + 1)$ matrix and $\mathbf{r} = \text{vec}([\mathbf{x}^0, \dots, \mathbf{x}^{D_V}])$ is a $P(D_V + 1)$-vector. We further denote $\mathbf{\Upsilon}_{k,q}$, to be the vector containing the q-th row of $\mathbf{\Upsilon}_k$ and $\sigma_{k,q}^2$ to be the values on the diagonal of $\mathbf{\Sigma}_k$ for $q \in [Q]$.

Here the incremental stochastic MM algorithm can be written as:

$$
\begin{align*}
\phi(v, \sigma) &:= \begin{bmatrix} \kappa \\ \zeta \\ \frac{\Delta}{\bar{\Sigma}} \end{bmatrix}, \quad \bar{S}(\tau; x) = \begin{bmatrix} \tau \otimes y^2 \\ -2\tau \otimes (y \otimes r) \\ \tau \otimes (\text{vec}(rr^\top) \otimes 1_Q) \\ \tau \otimes 1_Q \end{bmatrix}, \quad \text{where} \\
\bar{\Sigma} &= [\log \sigma_{1,1}^2, \dots, \log \sigma_{1,Q}^2, \log \sigma_{2,1}^2, \dots, \log \sigma_{K,Q}^2]^\top, \quad \kappa = \begin{bmatrix} \frac{1}{\sigma_{1,1}^2}, \dots, \frac{1}{\sigma_{1,Q}^2}, \frac{1}{\sigma_{2,1}^2}, \dots, \frac{1}{\sigma_{K,Q}^2} \end{bmatrix}^\top, \\
\zeta &= \text{vec} \left( \begin{bmatrix} \mathbf{\Upsilon}_{1,1}^\top, \dots, \mathbf{\Upsilon}_{1,Q}^\top, \mathbf{\Upsilon}_{2,1}^\top, \dots, \mathbf{\Upsilon}_{K,Q}^\top \end{bmatrix} \right), \quad \tau = [\tau_1, \tau_2, \dots, \tau_K]^\top, \\
\Delta &= \text{vec} \left( \frac{\text{vec}(\mathbf{\Upsilon}_{1,1}^\top \mathbf{\Upsilon}_{1,1})}{\sigma_{1,1}^2}, \dots, \frac{\text{vec}(\mathbf{\Upsilon}_{1,Q}^\top \mathbf{\Upsilon}_{1,Q})}{\sigma_{1,Q}^2}, \frac{\text{vec}(\mathbf{\Upsilon}_{2,1}^\top \mathbf{\Upsilon}_{2,1})}{\sigma_{2,1}^2}, \dots, \frac{\text{vec}(\mathbf{\Upsilon}_{K,Q}^\top \mathbf{\Upsilon}_{k,q;})}{\sigma_{k,q}^2} \right).
\end{align*}
$$

Results from Sections A.5.1 and A.5.2 imply the desired updates of Equation (6) presented in Algorithm 2.

### A.5.3 Parameters update

**Lemma 6.** Consider the series $\mathbf{u}_m$ defined as follows:

$$ \mathbf{u}_0 = \mathbf{A} \in \mathbb{R}^{m \times m}, \mathbf{u}_{m+1} = \mathbf{u}_m + \Lambda(\mathbf{B} - \mathbf{u}_m), \forall m \ge 1, \text{ here } \mathbf{B} \in \mathbb{R}^{m \times m}. $$

Then, $\mathbf{u}_m \succeq \mathbf{0}$ for all $m$ if the following conditions are satisfied: $\mathbf{u}_0 \succeq \mathbf{0}$, $\mathbf{B} \succeq \mathbf{0}$ and $1 > \Lambda > 0$.

The proofs for Lemma 6 are shown in Section B.3.

**Gating update:**

At the $n+1$ iteration, our gating parameters are calculated as follows:

$$ \omega_{n+1} = \arg\min_{\omega} [\omega^{\top} s_{n+1,1} + \text{vec}(\omega\omega^{\top})^{\top} s_{n+1,2}] . $$

The first and second derivatives of $r_1(\omega) = \omega^\top s_{n+1,1} + \text{vec}(\omega\omega^\top)^\top s_{n+1,2}$ are:

$$ \nabla r_1(\omega) = s_{n+1,1} + (\text{mat}(s_{n+1,2}) + \text{mat}(s_{n+1,2})^\top)\omega, \quad (36) $$

$$ \nabla^2 r_1(\omega) = \text{mat}(s_{n+1,2}) + \text{mat}(s_{n+1,2})^\top. \quad (37) $$

Since

$$ s_{n+1,2} = s_{n,2} + \gamma_{n+1} \left[ \frac{1}{2} \text{vec} \mathbf{B}_{n,K} - s_{n,2} \right], $$

hence:

$$ \text{mat}(s_{n+1,2}) = \text{mat}(s_{n,2}) + \gamma_{n+1} \left[ \frac{1}{2} \mathbf{B}_{n,K} - \text{mat}(s_{n,2}) \right]. $$

Now by applying Lemma 6 (knowing $\text{mat}(s_{0,2}) \succeq \mathbf{0}$), we have $\text{mat}(s_{n,2})$ is positive semi-definite for all $n$. Combine this with Equation (37), we have $\omega_{n+1}$ is the solution of Equation (36). Moreover, since $s_{0,2}$ and $\mathbf{B}_{n,K}$ are PSD matrices for all integer $n$, by induction, we have $\text{mat}(s_{n,2})$ is a PSD matrix for all $n$.

Hence,

$$ \omega_{n+1} = -(2 \text{mat}(s_{n+1,2}))^{-1} s_{n+1,1}, $$

where $\text{mat}(s_{n+1,2})$ represents the reshaping of $s_{n+1,2}$ into a $K(D_W + 1)P \times K(D_W + 1)P$ square matrix.
**Expert mean functions update:**

At the $n + 1$ iteration, our expert mean functions parameters are calculated as followed:

$$ \mathbf{\Upsilon}_{n+1,k,q,\cdot} = \arg\min_{\mathbf{\Upsilon}_{k,q,\cdot}} \left[ \mathbf{s}_{n+1,4,k;q}^{\top} \frac{\mathbf{\Upsilon}_{k,q,\cdot}^{\top}}{\sigma_{k,q}^{2}} + \mathbf{s}_{n+1,5,k;q}^{\top} \frac{\mathrm{vec}(\mathbf{\Upsilon}_{k,q,\cdot}^{\top}) \mathbf{\Upsilon}_{k,q,\cdot}}{\sigma_{k,q}^{2}} \right]. $$

The first and second derivatives of

$$ r_2(\mathbf{\Upsilon}_{k,q,\cdot}) = \mathbf{s}_{n+1,4,k;q}^{\top} \frac{\mathbf{\Upsilon}_{k,q,\cdot}^{\top}}{\sigma_{k,q}^{2}} + \mathbf{s}_{n+1,5,k;q}^{\top} \frac{\mathrm{vec}(\mathbf{\Upsilon}_{k,q,\cdot}^{\top}) \mathbf{\Upsilon}_{k,q,\cdot}}{\sigma_{k,q}^{2}} $$

are:

$$ \nabla r_2(\mathbf{\Upsilon}_{k,q,\cdot}) = \frac{\mathbf{s}_{n+1,4,k;q}}{\sigma_{k,q}^2} + \frac{\left( \mathrm{mat}(\mathbf{s}_{n+1,5,k;q}) + \mathrm{mat}(\mathbf{s}_{n+1,5,k;q})^{\top} \right)}{\sigma_{k,q}^2} \mathbf{\Upsilon}_{k,q,\cdot}^{\top}, \quad \text{and} \quad (38) $$

$$ \nabla^2 r_2(\mathbf{\Upsilon}_{k,q,\cdot}) = \frac{\left( \mathrm{mat}(\mathbf{s}_{n+1,5,k;q}) + \mathrm{mat}(\mathbf{s}_{n+1,5,k;q})^{\top} \right)}{\sigma_{k,q}^2}. \quad (39) $$

Since

$$ \mathbf{s}_{n+1,5,k;q} = \mathbf{s}_{n,5,k;q} + \gamma_{n+1} (\mathrm{vec}(\mathbf{r}\mathbf{r}^{\top}) - \mathbf{s}_{n,5,k;q}), $$

hence:

$$ \mathrm{mat}(\mathbf{s}_{n+1,5,k;q}) = \mathrm{mat}(\mathbf{s}_{n,5,k;q}) + \gamma_{n+1} (\mathbf{r}\mathbf{r}^{\top} - \mathrm{mat}(\mathbf{s}_{n,5,k;q})). $$

Now by applying Lemma 6, knowing that:

$$ \mathrm{mat}(\mathbf{s}_{0,5,k;q}) = \hat{\mathbf{s}}\hat{\mathbf{s}}^{\top} + \mathbf{I}_{D_V+1} \succeq \mathbf{0}, $$

we have $\mathrm{mat}(\mathbf{s}_{n,5,k;q})$ is positive definite for all $n$ and $k \in [K]$, $q \in [Q]$. Combine this with Equation (39), we have the solution of Equation (38) is the global minimum of $r_2(\mathbf{\Upsilon}_{k,q,\cdot})$. Additionally, it can be shown by induction that $\mathrm{mat}(\mathbf{s}_{n,5,k;q})$ is always a symmetric matrix.

Hence:

$$ \mathbf{\Upsilon}_{n+1,k,q,\cdot} = -(2\mathrm{mat}(\mathbf{s}_{n+1,5,k;q}))^{-1} \mathbf{s}_{n+1,4,k;q}, $$

where $\mathrm{mat}(\mathbf{s}_{n+1,5,k;q})$ is $\mathbf{s}_{n+1,5,k;q}$ reshape into $P(D_V + 1) \times P(D_V + 1)$ square matrix.

**Expert variance update:**

At the $n + 1$ iteration, our expert variance parameters are calculated as followed:

$$ \sigma_{n+1,k,q}^2 = \arg\min_{\sigma_{k,q}^2} \left[ s_{n+1,3,kq} \frac{1}{\sigma_{k,q}^2} + \mathbf{s}_{n+1,4,k;q}^{\top} \frac{\mathbf{\Upsilon}_{k,q,\cdot}^{\top}}{\sigma_{k,q}^2} + \mathbf{s}_{n+1,5,k;q}^{\top} \frac{\mathrm{vec}(\mathbf{\Upsilon}_{k,q,\cdot}^{\top}) \mathbf{\Upsilon}_{k,q,\cdot}}{\sigma_{k,q}^2} + s_{n+1,6,kq} \log(\sigma_{k,q}^2) \right]. $$

The first and second derivatives of:

$$ r_3(\sigma_{k,q}^2) = s_{n+1,3,kq} \frac{1}{\sigma_{k,q}^2} + \mathbf{s}_{n+1,4,k;q}^{\top} \frac{\mathbf{\Upsilon}_{k,q,\cdot}^{\top}}{\sigma_{k,q}^2} + \mathbf{s}_{n+1,5,k;q}^{\top} \frac{\mathrm{vec}(\mathbf{\Upsilon}_{k,q,\cdot}^{\top}) \mathbf{\Upsilon}_{k,q,\cdot}}{\sigma_{k,q}^2} + s_{n+1,6,kq} \log(\sigma_{k,q}^2) $$

are:

$$ \nabla r_3(\sigma_{k,q}^2) = - \frac{s_{n+1,3,kq} + \mathbf{s}_{n+1,4,k;q}^{\top} \mathbf{\Upsilon}_{k,q,\cdot}^{\top} + \mathbf{s}_{n+1,5,k;q}^{\top} \mathrm{vec}(\mathbf{\Upsilon}_{k,q,\cdot}^{\top})}{\sigma_k^4} + \frac{s_{n+1,6,kq}}{\sigma_{k,q}^2}, \quad (40) $$

$$ \nabla^2 r_3(\sigma_{k,q}^2) = 2 \frac{s_{n+1,3,kq} + \mathbf{s}_{n+1,4,k;q}^{\top} \mathbf{\Upsilon}_{k,q,\cdot}^{\top} + \mathbf{s}_{n+1,5,k;q}^{\top} \mathrm{vec}(\mathbf{\Upsilon}_{k,q,\cdot}^{\top})}{\sigma_{k,q}^6} - \frac{s_{n+1,6,kq}}{\sigma_{k,q}^4}. $$
We shall prove that the solution to Equation (40) is a global minimum of:

$$s_{n+1,3,kq} \frac{1}{\sigma_{k,q}^2} + s_{n+1,4,k;q}^{\top} \frac{\mathbf{\Upsilon}_{k,q,:}^{\top}}{\sigma_{k,q}^2} + s_{n+1,5,k;q}^{\top} \frac{\mathrm{vec}(\mathbf{\Upsilon}_{k,q,:}^{\top} \mathbf{\Upsilon}_{k,q,:})}{\sigma_{k,q}^2} + s_{n+1,6,kq} \log(\sigma_{k,q}^2).$$

Denote $\sigma_{n+1,k,q}^2 = \frac{s_{n+1,3,kq} + s_{n+1,4,k;q}^{\top} \mathbf{\Upsilon}_{k,q,:}^{\top} + s_{n+1,5,k;q}^{\top} \mathrm{vec}(\mathbf{\Upsilon}_{k,q,:}^{\top} \mathbf{\Upsilon}_{k,q,:})}{s_{n+1,6,kq}}$ be such a solution, we have:

$$\begin{aligned}
\sigma_{n+1,k,q}^6 \nabla^2 r_3(\sigma_{n+1,k,q}^2) &= s_{n+1,3,kq} + s_{n+1,4,k;q}^{\top} \mathbf{\Upsilon}_{k,q,:}^{\top} + s_{n+1,5,k;q}^{\top} \mathrm{vec}(\mathbf{\Upsilon}_{k,q,:}^{\top} \mathbf{\Upsilon}_{k,q,:}) \\
&= s_{n,3,kq} + \gamma_{n+1} (\tau_k y_q^2 - s_{n,3,kq}) + [s_{n,4,k;q} + \gamma_{n+1} ((-2y_q \tau_k \mathbf{r}) - s_{n,4,k;q})]^{\top} \mathbf{\Upsilon}_{k,q,:}^{\top} \\
&\quad + [s_{n,5,k;q} + \gamma_{n+1} (\tau_k (\mathrm{vec}(\mathbf{r}\mathbf{r}^{\top})) - s_{n,5,k;q})]^{\top} \mathrm{vec}(\mathbf{\Upsilon}_{k,q,:}^{\top} \mathbf{\Upsilon}_{k,q,:}) \\
&= (1 - \gamma_{n+1}) [s_{n,3,kq} + s_{n,4,k;q}^{\top} \mathbf{\Upsilon}_{k,q,:}^{\top} + s_{n,5,k;q}^{\top} \mathrm{vec}(\mathbf{\Upsilon}_{k,q,:}^{\top} \mathbf{\Upsilon}_{k,q,:})] + \gamma_{n+1} \tau_k (y_q - \mathbf{\Upsilon}_{k,q,:} \mathbf{r})^2.
\end{aligned}$$

Now notice that for all value of $\mathbf{\Upsilon}_{k,q,:}$ we have:

$$s_{0,3,kq} + s_{0,4,kq}^{\top} \mathbf{\Upsilon}_{k,q,:}^{\top} + s_{0,5,k;q}^{\top} \mathrm{vec}(\mathbf{\Upsilon}_{k,q,:}^{\top} \mathbf{\Upsilon}_{k,q,:}) = \left(\frac{1}{2} + \mathbf{\Upsilon}_{k,q,:} \hat{s}\right)^2 + \mathbf{\Upsilon}_{k,q,:}^{\top} \mathbf{\Upsilon}_{k,q,:} + \frac{3}{4} > 0.$$

Hence using induction we can prove that $\sigma_{n+1,k,q}^6 \nabla^2 r_3(\sigma_{n+1,k,q}^2) > 0$ for all $n, k, q, \mathbf{\Upsilon}_{k,q,:}$. Which directly leading to $\sigma_{n+1,k,q}^2$ being a local minimum. Now since $\sigma_{n+1,k,q}^2$ is the only solution of $\nabla r_3(\sigma_{n+1,k,q}^2)$, to prove that it is the global minimum, we only need to show that

$$\lim_{\sigma_{k,q}^2 \to 0^+} r_3(\sigma_{k,q}^2) \ge r_3(\sigma_{n+1,k,q}^2) \quad \text{and} \quad \lim_{\sigma_{k,q}^2 \to +\infty} r_3(\sigma_{k,q}^2) \ge r_3(\sigma_{n+1,k,q}^2).$$

Denote $U = s_{n+1,3,kq} + s_{n+1,4,k;q}^{\top} \mathbf{\Upsilon}_{k,q,:}^{\top} + s_{n+1,5,k;q}^{\top} \mathrm{vec}(\mathbf{\Upsilon}_{k,q,:}^{\top} \mathbf{\Upsilon}_{k,q,:}) > 0$, we have:

$$\lim_{\sigma_{k,q}^2 \to +\infty} r_3(\sigma_{k,q}^2) = \lim_{\sigma_{k,q}^2 \to +\infty} \frac{U}{\sigma_{k,q}^2} + \lim_{\sigma_{k,q}^2 \to +\infty} s_{n+1,6,kq} \log(\sigma_{k,q}^2) = 0 + \infty = +\infty > r_3(\sigma_{n+1,k,q}^2),$$

$$\lim_{\sigma_{k,q}^2 \to 0^+} r_3(\sigma_{k,q}^2) = \lim_{\sigma_{k,q}^2 \to +\infty} r_3(\frac{1}{\sigma_{k,q}^2}) = \lim_{\sigma_{k,q}^2 \to +\infty} (U\sigma_{k,q}^2 - s_{n+1,6,kq} \log(\sigma_{k,q}^2)) = +\infty > r_3(\sigma_{n+1,k,q}^2).$$

Hence this has proven that $\sigma_{n+1,k,q}^2$ is the global minimum, thus:

$$\sigma_{n+1,k,q}^2 = \frac{s_{n+1,3,kq} + s_{n+1,4,k;q}^{\top} \mathbf{\Upsilon}_{k,q,:}^{\top} + s_{n+1,5,k;q}^{\top} \mathrm{vec}(\mathbf{\Upsilon}_{k,q,:}^{\top} \mathbf{\Upsilon}_{k,q,:})}{s_{n+1,6,kq}}.$$

### Notes on the updates of $\mathbf{\Upsilon}$ and $\Sigma$:

**Lemma 7.** Consider the function

$$f(u, v) = \frac{g(v)}{u} + c \log(u).$$

Where $v$ is a vector, $u$ is a positive real number, $c$ is a positive constant and $g(v) : \mathbb{R}^n \to \mathbb{R}^+$ be a function that has exactly one global minimum. Denote $v_0 = \arg\min_{v} g(v)$ be the global minimum of

$g(v)$ and $u_0 = \arg\min_{u} \left( \frac{g(v_0)}{u} + \log(u) \right)$. We have $f(u_0, v_0)$ is the global minimum of $f$ over all $u, v$.

**Proof of Lemma 7.** For all $v$, we have:

$$\nabla_u f(u, v) = \frac{1}{u^2} (-g(v) + cu), \quad \nabla_u^2 f(u, v) = \frac{1}{u^3} (2g(v) - cu).$$
Hence $u = \frac{g(v)}{c}$ is the local minimum of $f(u, v)$ for all $v$. But since $\nabla_u f(u, v)$ has exactly one solution for $u$ and both $\lim_{u \to +\infty} f(u, v)$, $\lim_{u \to +0} f(u, v)$ is $+\infty$. We can conclude that $u = \frac{g(v)}{c}$ is the global minimum of $f(u, v)$ for all given $v$. Assume that $(u_1, v_1)$ is the minimum of $f(u, v)$, we have

$$f(u_1, v_1) \ge f\left(\frac{g(v_1)}{c}, v_1\right) = c \log\left(\frac{g(v_1)}{c}\right) \ge c \log\left(\frac{g(v_0)}{c}\right) = f\left(\frac{g(v_0)}{c}, v_0\right).$$

Which lead to $g(v_1) = g(v_0)$. But since $v_0$ is the only global minimum of $g(v)$, we must have $v_1 = v_0$. Thus, leading to $u_1$ also is equal to $u_0$. $\square$

**Remark 3.** In our problem, $v = \Upsilon_{k,q;:}^{\top}$, $u = \sigma_{k,q}^2$ and $g(v) = \sigma_{k,q}^2 r_2(\Upsilon_{k,q;:}^{\top})$. It has been proven in the expert mean functions update section that $\sigma_{k,q}^2 r_2(\Upsilon_{k,q;:}^{\top})$ has exactly unique global minimum (since it is strictly convex). Hence, Lemma 7 shows that our way of updating $\Upsilon$ and then $\Sigma$ is the same as updating the expert function all at once.

## A.6 Proof of Theorem 3

### A.6.1 Proof of surrogate function construction in Proposition 4

By applying Lemma 8 for the function $f = -\log$ and:

$$ \begin{aligned} \hat{c}_k &= [1, 1, \dots, 1]^{\top}, \\ a &= [g_k(x_n; \omega) e_k(v_{y_n}(x_n; v))]_{k \in [K]}^{\top}, \\ b &= [g_k(x_n; \omega_n) e_k(v_{y_n}(x_n; v_n))]_{k \in [K]}^{\top}. \end{aligned} $$

Here we $\hat{c}$ is used as $c$ in Lemma 8 to not be mistaken for the vectorization of $v$ in Proposition 4. We obtain:

$$-\log \left[ \sum_{k=1}^{K} g_k(x_n; \omega) e_k(v_{y_n}(x_n; v)) \right] \le -\sum_{k=1}^{K} \tau_{n,k} \log \left[ \frac{1}{\tau_{n,k}} g_k(x_n; \omega) e_k(v_{y_n}(x_n; v)) \right],$$

$$\text{where } \tau_{n,k} = \frac{g_k(x_n; \omega_n) e_k(v_{y_n}(x_n; v_n))}{\sum_{l=1}^{K} g_l(x_n; \omega_n) e_l(v_{y_n}(x_n; v_n))}.$$

Hence:

$$ \begin{aligned} -\log s_{\theta}(y_n | x_n) &\le -\sum_{k=1}^{K} \tau_{n,k} \log \left[ \frac{1}{\tau_{n,k}} g_k(x_n; \omega) e_k(v_{y_n}(x_n; v)) \right] \\ &= \sum_{k=1}^{K} \tau_{n,k} \log(\tau_{n,k}) - \sum_{k=1}^{K} \tau_{n,k} [\log(g_k(x_n; \omega)) + \log(e_k(v_{y_n}(x_n; v)))] \end{aligned} $$

Similar to Section A.4.1, by applying Lemma 2 and Proposition 5, we obtain the first part of the majorization function as follows:

$$ \begin{aligned} -\log s_{\theta}(y_n | x_n) &\le \sum_{k=1}^{K} \tau_{n,k} \log(\tau_{n,k}) + g(\omega_n; x_n) - \omega_n^{\top} \nabla g(\omega_n; x_n) - \sum_{k=1}^{K} \tau_{n,k} w_k(x_n) \\ &\quad + \left\{ \omega^{\top} \nabla g(\omega_n; x_n) + \frac{1}{2}(\omega - \omega_n^{\top}) B_{K,n} (\omega - \omega_n) \right\} \\ &\quad - \sum_{k=1}^{K} \tau_{n,k} \log(e_k(v_{y_n}(x_n; v))). \end{aligned} \quad (41) $$
Hence, using the same idea of the majorization construction for $-\sum_{k=1}^{K} \tau_{n,k} \log(e_k(v_{y_n}(x_n; v)))$, our overall majorization function takes the form:

$$ 
\begin{aligned}
& - \log s_{\theta}(y_n | x_n) \\
\le & \sum_{k=1}^{K} \tau_{n,k} \log(\tau_{n,k}) - \sum_{k=1}^{K} \tau_{n,k} w_k(x_n) - \sum_{k=1}^{K} \tau_{n,k} \left( \sum_{l=1}^{M} 1(y_n, l) v_{l,k}(x_n) \right) \\
& + \left\{ \underline{g}(\omega_n; x_n) + (\omega - \omega_n)^{\top} \nabla \underline{g}(\omega_n; x_n) + \frac{1}{2}(\omega - \omega_n)^{\top} \mathbf{B}_{n,K}(\omega - \omega_n) \right\} \\
& + \sum_{k=1}^{K} \tau_{n,k} \left\{ \underline{e}(c_{n,k}; x_n) + (c_{n,k} - c_{n,k})^{\top} \nabla \underline{e}(c_{n,k}; x_n) + \frac{1}{2}(c_{n,k} - c_{n,k})^{\top} \mathbf{B}_{n,M}(c_{n,k} - c_{n,k}) \right\} \\
= & \sum_{k=1}^{K} \tau_{n,k} \log(\tau_{n,k}) - (\omega^{\top} s_n) - (v^{\top} r_n) \\
& + \left\{ \underline{g}(\omega_n; x_n) + (\omega - \omega_n)^{\top} \nabla \underline{g}(\omega_n; x_n) + \frac{1}{2}(\omega - \omega_n)^{\top} \mathbf{B}_{n,K}(\omega - \omega_n) \right\} \\
& + \sum_{k=1}^{K} \tau_{n,k} \left\{ \underline{e}(c_{n,k}; x_n) + (c_{n,k} - c_{n,k})^{\top} \nabla \underline{e}(c_{n,k}; x_n) + \frac{1}{2}(c_{n,k} - c_{n,k})^{\top} \mathbf{B}_{n,M}(c_{n,k} - c_{n,k}) \right\}.
\end{aligned}
$$

Here $v$ is ordered as $(v_{m,k,d,p}|_{k \in [K], p \in [P], m \in [M-1], d \in \{0, \dots, D_W\}})$ and

$$ \mathbf{c}_{n,k} = [v_{n,1,k,0,1}, \dots, v_{n,1,k,D_V,1}, \dots, v_{n,M,k,D_V,1}, \dots, v_{n,M,k,D_V,P}]^{\top}, $$

$$ \underline{e}(c_{n,k}; x_n) = \log \left( 1 + \sum_{l=1}^{M-1} \exp(v_{n,l,k}(x_n)) \right), \quad \underline{g}(\omega_n; x_n) = \log \left( 1 + \sum_{k=1}^{K-1} \exp(w_{n,k}(x_n)) \right), $$

$$ \tau_{n,k} = \frac{\underline{g}_k(x_n; \omega_n) e_k(v_{y_n}(x_n; v_n))}{\sum_{l=1}^{K} \underline{g}_l(x_n; \omega_n) e_l(v_{y_n}(x_n; v_n))}. $$

### A.6.2 Verification of P 1 to 5, and Assumptions 1, 2 and 8 for the softmax-gated multinomial logistic MoE models

We will assume that $x$ and $y$ are bounded by $x^\#$ and $y^\#$, respectively.

**Verification of P 1.** For the problem of softmax-gated MoE, we consider the majorizer surrogate $g$ to have the following form that belongs to the exponential family

$$ g(\theta, z; \tau) := -\psi(\theta) + \langle \bar{S}(\tau; z), \phi(\theta) \rangle, \quad (42) $$

where $\psi(\omega) = C_{n,k}$, and

$$ \phi(\omega) := \begin{bmatrix} \omega \\ \mathrm{vec}(\omega\omega^{\top}) \\ v \\ \mathrm{vec}((\mathrm{diag}_{MP(D_V+1)}(vv^{\top}))) \end{bmatrix}, \quad \bar{S}(\tau; z) = \begin{bmatrix} -s + \nabla \underline{g}(\tau) - \mathbf{B}_{K-1}\tau \\ \frac{1}{2}\mathrm{vec}(\mathbf{B}_{K-1}) \\ -\mathbf{r} + \mathrm{vec}([\tau_k \nabla \underline{e}(\tau) - \tau_k \mathbf{B}_{M-1}\tau]_{k=1,\dots,K}) \\ \frac{1}{2}\mathrm{vec}(\tau \otimes \mathbf{B}_{M-1}) \end{bmatrix}. $$

where $\tau$, $B$ are defined Section 4.1.3 and $v, \omega, r, s$ are defined in Section 4.2.2. Furthermore, since $\phi, \psi$ and $\bar{S}$ are composed of elementary functions applied component-wise to their parameters, they are trivially continuously differentiable and measurable due to the smoothness and measurability properties of elementary functions.

**Verification of P 2.** We can define $\mathcal{S}$ as follows:

$$ \mathcal{S} = \{ (s_1, \mathrm{vec}(S_2), s_3, \mathrm{vec}(S_4)) : s_1 \in \mathbb{R}^{(K-1)P(D_W+1)}; S_2 \in \mathbb{R}^{(K-1)P(D_W+1) \times (K-1)P(D_W+1)}, \\ S_2 \succeq 0; s_3 \in \mathbb{R}^{(D_V+1)MPK}; S_4 \in \mathbb{R}^{(D_V+1)MPK \times (D_V+1)MPK}, S_4 \succeq 0 \}. $$

It is trivial that $\mathcal{S}$ denoted above satisfies the requirements of P 2.

**Verification of P 3.** $\bar{S}(\theta, z)$ is composed of elementary functions applied component-wise to $x$ and $y$, hence its expectation exist and lies in $\mathcal{S}$. Furthermore, since $x$ and $y$ are bounded, it can be shown that $\tau, B, r, s, \nabla g, \nabla e$ are all bounded. Thus, $\mathbb{E}[\bar{S}(\theta, z)]$ is finite.
**Verification of P 4.** This assumption is proven in Section A.7.2, where we have shown that $v, \omega$ of the next iteration are the unique solution to the minimization problem of the current iteration.

**Verification of P 5.** Since $f(\theta; z)$ is continuously differentiable in $\theta$ for all $z$ and dominated by an integrable envelope (due to bounded covariates and responses), differentiation under the expectation is valid, hence $\theta \mapsto E[f(\theta; z)]$ is $C^1$ on $\mathbb{T}$.

**Verification of Assumption 1.** Since we have $\omega \in \mathbb{R}^{P \times K \times (D_W+1)}$ and $v \in \mathbb{R}^{(D_V+1) \times M \times P \times K}$, thus $\mathbb{T} = \mathbb{R}^{P \times K \times (D_W+1)} \times \mathbb{R}^{(D_V+1) \times M \times P \times K}$ is a convex and open set.

Furthermore, similar to the verification of P 1, we can show that $\psi$ and $\phi$ is twice continuously differentiable throughout $\mathbb{T}$.

**Verification of Assumption 2.** From Section A.7.2, given $s = (s_1, s_2, s_3, s_4)$, we have the function $\bar{\theta}(s)$ is defined as:

$$ \bar{\theta}(s) = \begin{bmatrix} - (\mathrm{mat}(s_{2,n+1}) + \mathrm{mat}(s_{2,n+1})^{\top})^{-1} s_{1,n+1} \\ - [\mathrm{bdiag}_{\ell}^{-1}(\mathrm{mat}_{\ell}(s_{4,n+1})) + \mathrm{bdiag}_{\ell}^{-1}(\mathrm{mat}_{\ell}(s_{4,n+1}))^{\top}]^{-1} s_{3,n+1} \end{bmatrix}. $$

It is trivial to show that $\bar{\theta}(s)$ is continuously differentiable for $s_1, s_3$.

We have $A \to \mathrm{mat}(A) + \mathrm{mat}(A)^{\top}$ and $A \to \mathrm{bdiag}_{\ell}^{-1}(\mathrm{mat}_{\ell}(A)) + \mathrm{bdiag}_{\ell}^{-1}(\mathrm{mat}_{\ell}(A))^{\top}$ are affine transformations, thus they are continuously differentiable. Furthermore, from Section A.7.2, we know that $\mathrm{mat}(s_2) + \mathrm{mat}(s_2)^{\top}$ and $\mathrm{bdiag}_{\ell}^{-1}(\mathrm{mat}_{\ell}(s_4)) + \mathrm{bdiag}_{\ell}^{-1}(\mathrm{mat}_{\ell}(s_4))^{\top}$ are positive definite, thus their inverts are non-singular, thus their inverts are continuously differentiable. This concludes our proof.

**Verification of Assumption 8.** By simply choosing $\gamma_n = \gamma_0 n^{-\alpha}$ with $\alpha \in (\frac{1}{2}, 1]$ and $\gamma_0 \in (0, 1)$, all requirements of Assumption 8 will be met.

**Remark 4.** Assumption 3 cannot be verified due to the ambiguity in the definition of the bdiag function. Similarly to softmax-gated Gaussian MoE, Assumption 7 for scenarios $s_1$ and $s_3$, as well as Assumptions 4 and 5, also remain unproven.

## A.7 Derivation of Algorithm 3

### A.7.1 Construction of the incremental MM algorithm

From here the majorization from the previous section, the incremental MM algorithm can be constructed as:

$$ \phi(\omega) := \begin{bmatrix} \omega \\ \mathrm{vec}(\omega\omega^{\top}) \\ v \\ \mathrm{vec}(\mathrm{bdiag}_{MP(D_V+1)}(vv^{\top})) \end{bmatrix}, \quad \bar{S}(\tau; z) = \begin{bmatrix} -s + \nabla g(\tau) - B_{K-1}\tau \\ -r + \mathrm{vec}([\tau_k \nabla e(\tau) - \tau_k B_{M-1}\tau]_{k=1,...,K}) \\ -\frac{1}{2} \mathrm{vec}(\tau \otimes B_{M-1}) \end{bmatrix}. $$

and $\psi(\omega) = C_{n,k}$, which is the independent constant with respect to the parameter $\omega$. Here, recall that we have $\tau = [\tau_1, ..., \tau_k]^\top$ and for any matrix $A \in \mathbb{R}^{n \times n}$ (n is a positive integer) we have $\mathrm{bdiag}_t(A) \in \mathbb{R}^{n \times t}$ defined as:

$$ \mathrm{bdiag}_t(A)_{[it+1:it+t];[1:t]} = A_{[it+1:it+t];[it+1:it+t]}, \quad \forall i \in \{0, \dots, n-1\}. $$

### A.7.2 Parameters update

**Gating parameters update:**

At the $n+1$ iteration, our gating parameters are calculated as follows:

$$ \omega_{n+1} = \arg\min_{\omega} [\omega^{\top} s_{n+1,1} + \mathrm{vec}(\omega\omega^{\top})^{\top} s_{n+1,2}]. $$

The first and second derivatives of $r_1(\omega) = \omega^{\top} s_{n+1,1} + \mathrm{vec}(\omega\omega^{\top})^{\top} s_{n+1,2}$ are:

$$ \nabla r_1(\omega) = s_{n+1,1} + (\mathrm{mat}(s_{n+1,2}) + \mathrm{mat}(s_{n+1,2})^{\top})\omega, \quad (43) $$

$$ \nabla^2 r_1(\omega) = \mathrm{mat}(s_{n+1,2}) + \mathrm{mat}(s_{n+1,2})^{\top}. \quad (44) $$
Since

$$s_{n+1,2} = s_{n,2} + \gamma_{n+1} \left[ \frac{1}{2} \mathrm{vec} \, B_{n,K} - s_{n,2} \right],$$

hence:

$$\mathrm{mat}(s_{n+1,2}) = \mathrm{mat}(s_{n,2}) + \gamma_{n+1} \left[ \frac{1}{2} B_{n,K} - \mathrm{mat}(s_{n,2}) \right].$$

Now by applying Lemma 6 (knowing $\mathrm{mat}(s_{0,2}) \succeq \mathbf{0}$), we have $\mathrm{mat}(s_{n,2})$ is positive definite for all $n$. Combine this with Equation (44), we have $\omega_{n+1}$ is the solution of Equation (43).

Hence,

$$\omega_{n+1} = - (\mathrm{mat}(s_{n+1,2}) + \mathrm{mat}(s_{n+1,2})^{\top})^{-1} s_{n+1,1},$$

where $\mathrm{mat}(s_{n+1,2})$ represents the reshaping of $s_{n+1,2}$ into a $K(D_W + 1)P \times K(D_W + 1)P$ square matrix.

### Expert parameters update:

Let us denote $\iota = MP(D_V + 1)$ and for any positive integer $k$, let us define the two functions $\mathrm{mat}_k : \mathbb{R}^{ks} \to \mathbb{R}^{s \times k}$ and $\mathrm{bdiag}_k^{-1} : \mathbb{R}^{g \times k} \to \mathbb{R}^{g \times gk}$ satisfying:

$$\mathrm{vec}(\mathrm{mat}_k(A)) = A \quad \forall A \in \mathbb{R}^{ks},$$

$$\mathrm{bdiag}_k(\mathrm{bdiag}_k^{-1}(A)) = A \quad \forall A \in \mathbb{R}^{ks \times k}.$$

At the $n+1$ iteration, our expert parameters are calculated as followed:

$$v_{n+1} = \arg\min_{v} \left[ v^{\top} s_{n+1,3} + \mathrm{vec} \left( \mathrm{bdiag}_{\iota}(v v^{\top}) \right)^{\top} s_{n+1,4} \right].$$

The first and second derivatives of $r_1(v) = v^{\top}s_{n+1,3} + \mathrm{vec} \left( \mathrm{bdiag}_{\iota}(v v^{\top}) \right)^{\top}s_{n+1,4}$ are:

$$\nabla r_1(v) = s_{n+1,3} + [\mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n+1,4})) + \mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n+1,4}))^{\top}] v, \quad (45)$$

$$\nabla^2 r_1(v) = \mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n+1,4})) + \mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n+1,4}))^{\top}. \quad (46)$$

Since

$$s_{n+1,4} = s_{n,4} + \gamma_{n+1} \left[ \frac{1}{2} \mathrm{vec} \left[ \tau \otimes B_{n,M} \right] - s_{n,4} \right],$$

hence:

$$\mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n+1,4})) =$$

$$\mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n,4})) + \gamma_{n+1} \left[ \frac{1}{2} \mathrm{bdiag}_{\iota}^{-1} [\tau \otimes B_{n,M}] - \mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n,4})) \right].$$

Now by applying Lemma 6 (knowing $\mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{0,4})) \succeq \mathbf{0}$), we have $\mathrm{mat}(s_{n,4})$ is positive definite for all $n$. Combine this with Equation (46), we have $v_{n+1}$ is the solution of Equation (45).

Hence,

$$v_{n+1} = - [\mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n+1,4})) + \mathrm{bdiag}_{\iota}^{-1}(\mathrm{mat}_{\iota}(s_{n+1,4}))^{\top}]^{-1} s_{n+1,3}.$$

# B Technical proofs

## B.1 Proof of Lemma 2

Using the spectral decomposition we can write:

$$C = M_1^{\top} D_1 M_1, \quad B - A = M_2^{\top} D_2 M_2.$$

Hence:

$$C \otimes (B - A) = (M_1^{\top} D_1 M_1) \otimes (M_2^{\top} D_2 M_2) = (M_1 \otimes M_2)^{\top} (D_1 \otimes D_2) (M_1 \otimes M_2),$$

which lead to $C \otimes (B - A)$ are nonnegative definite, thus proven the lemma.
## B.2 Proof of Proposition 5

Denote $\epsilon = \frac{1}{\pi_1 + \pi_2 + \dots + \pi_N}$. Hence $\epsilon\hat{p}$ is a stochastic vector which according to Lemma 11:

$$ \epsilon \mathrm{diag}(\hat{\pi}) - \epsilon^2 \hat{\pi} \hat{\pi}^{\top} \leq \frac{1}{2} \left( \mathbf{I}_K - \frac{\mathbf{1}\mathbf{1}^{\top}}{K} \right). $$

Hence:

$$ \begin{aligned} \mathrm{diag}(\hat{\pi}) - \hat{\pi}\hat{\pi}^{\top} &\leq \frac{1}{2\epsilon^2} \left( \mathbf{I}_K - \frac{\mathbf{1}\mathbf{1}^{\top}}{K} \right) + \left( 1 - \frac{1}{\epsilon} \right) \mathrm{diag}(\hat{\pi}) \\ &\leq \frac{1}{2} \left( \mathbf{I}_K - \frac{\mathbf{1}\mathbf{1}^{\top}}{K} \right) + \pi_{N+1} \max_{n \in [N]} \pi_n \mathbf{I}_K \\ &\leq \frac{1}{2} \left( \mathbf{I}_K - \frac{\mathbf{1}\mathbf{1}^{\top}}{K} \right) + \frac{1}{4} \mathbf{I}_K \\ &\leq \frac{1}{2} \left( \frac{3}{2} \mathbf{I}_K - \frac{\mathbf{1}\mathbf{1}^{\top}}{K} \right). \end{aligned} $$

## B.3 Proof of Lemma 6

Let say that all the conditions are met, we shall prove that $u_m \succeq 0$ for all $m$ using induction.

Since $u_1 \succeq 0$, we assume that $u_i \succeq 0$ for all $i \leq m$. Our target is to prove that $u_{m+1} \succeq 0$.

This is true since the term $u_m(1 - \Lambda) \succeq 0$ and the term $\Lambda B \succeq 0$, hence their sum which is $u_{m+1}$ is also positive definite.

# C Technical results

Surrogate function construction plays a pivotal role in optimization techniques such as MM algorithms. By leveraging mathematical inequalities, it is possible to approximate or bound complex objective functions with simpler, more tractable alternatives. Below are several fundamental inequalities commonly employed in the development of surrogate functions:

**Lemma 8** (See, e.g., page 191 in [57]). Given three vectors $\mathbf{a}, \mathbf{b}, \mathbf{c} \in \mathbb{R}^D$ such that all components of $\mathbf{a}, \mathbf{b}, \mathbf{c}$ are positive and a convex function $f$. The inequality below holds:

$$ f(\mathbf{c}^{\top} \mathbf{a}) \leq \sum_{d=1}^{D} \frac{c_d y_d}{\mathbf{c}^{\top} \mathbf{b}} f\left(\frac{\mathbf{c}^{\top} \mathbf{b}}{y_d} x_d\right), \text{ where the equality holds when } \mathbf{a} = \mathbf{b}. $$

**Lemma 9** (See, e.g., page 191 in [57]). Given four vectors $\mathbf{a}, \mathbf{b}, \mathbf{c}, \alpha \in \mathbb{R}^D$, and a convex function $f$, such that all $\alpha_d \geq 0$, $\sum_{d=1}^{D} \alpha_d = 1$, and $\alpha_d > 0$ whenever $c_d \neq 0$ for all $d \in [D]$. The inequality below holds:

$$ f(\mathbf{c}^{\top} \mathbf{a}) \leq \sum_{d=1}^{D} \alpha_d f\left(\frac{c_d}{\alpha_d}(x_d - y_d)\right) + \mathbf{c}^{\top} \mathbf{b}, \text{ where the equality holds when } \mathbf{a} = \mathbf{b}. $$

**Lemma 10** (See, e.g., page 206 in [57]). Given 2 positive reals $u, y$. The inequality below holds:

$$ \log(u) \leq \frac{u}{y} + \log(y) - 1, \text{ where the equality holds when } u = y. $$

**Lemma 11** (Theorem 5.3 in [5]). Let $\pi_N$ be a stochastic vector of dimension $N$. We have:

$$ \mathrm{diag}(\pi_N) - \pi_N \pi_N^{\top} \leq \frac{1}{2} \left( \mathbf{I}_K - \frac{\mathbf{1}_N \mathbf{1}_N^{\top}}{N} \right). $$
References

[1] C. Andrieu, E. Moulines, and P. Priouret. Stability of Stochastic Approximation under Verifiable Conditions. *SIAM Journal on Control and Optimization*, 44(1):283–312, 2005. (Cited on pages 6, 24, and 25.)

[2] Y. Bengio. Deep Learning of Representations: Looking Forward. In *Statistical Language and Speech Processing*, pages 1–37, Berlin, Heidelberg, 2013. (Cited on page 2.)

[3] M. Blein-Nicolas, S. S. Negro, T. Balliau, C. Welcker, L. Cabrera-Bosquet, S. D. Nicolas, A. Charcosset, and M. Zivy. A systems genetics approach reveals environment-dependent associations between snps, protein coexpression, and drought-related traits in maize. *Genome Research*, 30(11):1593–1604, 2020. (Cited on page 20.)

[4] D. Böhning. Multinomial logistic regression algorithm. *Annals of the institute of Statistical Mathematics*, 44(1):197–200, 1992. (Cited on pages 3 and 26.)

[5] D. Böhning and B. G. Lindsay. Monotonicity of quadratic-approximation algorithms. *Annals of the Institute of Statistical Mathematics*, 40(4):641–663, 1988. (Cited on page 37.)

[6] V. S. Borkar. *Stochastic approximation: a dynamical systems viewpoint*, volume 9. Springer, 2008. (Cited on page 2.)

[7] W. Cai, J. Jiang, F. Wang, J. Tang, S. Kim, and J. Huang. A Survey on Mixture of Experts in Large Language Models. *IEEE Transactions on Knowledge and Data Engineering*, 37(7):3896–3915, 2025. (Cited on page 2.)

[8] O. Cappé. Online sequential monte carlo em algorithm. In *2009 IEEE/SP 15th Workshop on Statistical Signal Processing*, pages 37–40. IEEE, 2009. (Cited on page 13.)

[9] O. Cappé. Online em algorithm for hidden markov models. *Journal of Computational and Graphical Statistics*, 20(3):728–749, 2011. (Cited on page 13.)

[10] O. Cappé and E. Moulines. On-line expectation-maximization algorithm for latent data models. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 71(3):593–613, 2009. (Cited on page 13.)

[11] O. Cappé. Online EM Algorithm for Hidden Markov Models. *Journal of Computational and Graphical Statistics*, 20(3):728–749, Jan. 2011. (Cited on page 2.)

[12] O. Cappé and E. Moulines. On-Line Expectation-Maximization Algorithm for latent Data Models. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 71(3):593–613, Feb. 2009. (Cited on pages 2 and 5.)

[13] F. Chamroukhi and B. T. Huynh. Regularized maximum-likelihood estimation of mixture-of-experts for regression and clustering. In *2018 International Joint Conference on Neural Networks (IJCNN)*, pages 1–8, 2018. (Cited on page 8.)

[14] F. Chamroukhi and B. T. Huynh. Regularized maximum likelihood estimation and feature selection in mixtures-of-experts models. *Journal de la Société Française de Statistique*, 160(1):57–85, 2019. (Cited on page 8.)

[15] J. Chen, J. Zhu, Y. W. Teh, and T. Zhang. Stochastic Expectation Maximization with Variance Reduction. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, *Advances in Neural Information Processing Systems*, volume 31. Curran Associates, Inc., 2018. (Cited on page 2.)

[16] Z. Chen, Y. Deng, Y. Wu, Q. Gu, and Y. Li. Towards Understanding the Mixture-of-Experts Layer in Deep Learning. In *Advances in Neural Information Processing Systems*, 2022. (Cited on pages 2 and 8.)

[17] M. C. Chong, H. D. Nguyen, and T. Nguyen. Risk Bounds for Mixture Density Estimation on Compact Domains via the h-Lifted Kullback-Leibler Divergence. *Transactions on Machine Learning Research*, 2024. (Cited on page 2.)

[18] E. Chouzenoux and J.-B. Fest. Sabrina: A stochastic subspace majorization-minimization algorithm. *Journal of Optimization Theory and Applications*, 195:919–952, 2022. (Cited on page 2.)

[19] E. Chouzenoux and J.-C. Pesquet. A Stochastic Majorize-Minimize Subspace Algorithm for Online Penalized Least Squares Estimation. *IEEE Transactions on Signal Processing*, 65(18):4770–4783, 2017. (Cited on page 2.)
[20] S. L. Corff and G. Fort. Online Expectation Maximization based algorithms for inference in Hidden Markov Models. *Electronic Journal of Statistics*, 7(none):763–792, 2013. (Cited on page 2.)

[21] D. Dai, C. Deng, C. Zhao, R. Xu, H. Gao, D. Chen, J. Li, W. Zeng, X. Yu, Y. Wu, Z. Xie, Y. Li, P. Huang, F. Luo, C. Ruan, Z. Sui, and W. Liang. DeepSeekMoE: Towards Ultimate Expert Specialization in Mixture-of-Experts Language Models. In L.-W. Ku, A. Martins, and V. Srikumar, editors, *Proceedings of the 62nd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)*, pages 1280–1297, Bangkok, Thailand, Aug. 2024. Association for Computational Linguistics. (Cited on page 2.)

[22] B. Delyon, M. Lavielle, and E. Moulines. Convergence of a stochastic approximation version of the EM algorithm. *Annals of statistics*, pages 94–128, 1999. (Cited on page 5.)

[23] S. Dereich, A. Jentzen, and S. Kassing. ODE approximation for the Adam algorithm: General and overparametrized setting. *arXiv preprint arXiv:2511.04622*, 2025. (Cited on page 3.)

[24] A. Dieuleveut, G. Fort, M. Hegazy, and H.-T. Wai. Federated Majorize-Minimization: Beyond Parameter Aggregation. *arXiv preprint arXiv:2507.17534*, 2025. (Cited on pages 2, 4, and 5.)

[25] A. Dieuleveut, G. Fort, E. Moulines, and H.-T. Wai. Stochastic Approximation Beyond Gradient for Signal Processing and Machine Learning. *IEEE Transactions on Signal Processing*, 71:3117–3148, 2023. (Cited on page 2.)

[26] D. Do, L. Do, and X. Nguyen. Strong identifiability and parameter learning in regression with heterogeneous response. *arXiv preprint arXiv:2212.04091*, 2022. (Cited on page 2.)

[27] W. Fedus, B. Zoph, and N. Shazeer. Switch Transformers: Scaling to Trillion Parameter Models with Simple and Efficient Sparsity. *Journal of Machine Learning Research*, 23(120):1–39, 2022. (Cited on page 2.)

[28] G. Fort, F. Forbes, and H. D. Nguyen. Sequential Sample Average Majorization–Minimization. *Statistics and Computing*, 2025. Accepted for publication. (Cited on pages 2, 13, and 14.)

[29] G. Fort, P. Gach, and E. Moulines. Fast incremental expectation maximization for finite-sum optimization: nonasymptotic convergence. *Statistics and Computing*, 31(4):48, June 2021. (Cited on pages 2 and 5.)

[30] G. Fort and E. Moulines. The Perturbed Prox-Preconditioned Spider Algorithm: Non-Asymptotic Convergence Bounds. In *2021 IEEE Statistical Signal Processing Workshop (SSP)*, pages 96–100, 2021. (Cited on page 2.)

[31] G. Fort and E. Moulines. Stochastic variable metric proximal gradient with variance reduction for non-convex composite optimization. *Statistics and Computing*, 33(3), Apr. 2023. (Cited on pages 2 and 5.)

[32] G. Fort, E. Moulines, and H.-T. Wai. A Stochastic Path Integral Differential Estimator Expectation Maximization Algorithm. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, *Advances in Neural Information Processing Systems*, volume 33, pages 16972–16982. Curran Associates, Inc., 2020. (Cited on page 2.)

[33] G. Fort, E. Moulines, and H.-T. Wai. Geom-Spider-EM: Faster Variance Reduced Stochastic Expectation Maximization for Nonconvex Finite-Sum Optimization. In *ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, pages 3135–3139, 2021. (Cited on page 2.)

[34] J. H. Friedman, T. Hastie, and R. Tibshirani. Regularization Paths for Generalized Linear Models via Coordinate Descent. *Journal of Statistical Software*, 33(1):1–22, Feb. 2010. Section: Articles. (Cited on page 20.)

[35] C. R. Genovese and L. Wasserman. Rates of convergence for the Gaussian mixture sieve. *The Annals of Statistics*, 28(4):1105–1127, aug 2000. (Cited on page 2.)

[36] D. T. Hai, T. N. Mai, T. Nguyen, N. Ho, B. T. Nguyen, and C. Drovandi. Dendrograms of Mixing Measures for Softmax-Gated Gaussian Mixture of Experts: Consistency without Model Sweeps. In *Proceedings of The 29th International Conference on Artificial Intelligence and Statistics*, Proceedings of Machine Learning Research. PMLR, 2026. (Cited on page 2.)

[37] P. Hall and C. C. Heyde. *Martingale limit theory and its applications* / P. Hall, C.C. Heyde. Probability and mathematical statistics. Academic Press, New York;, 1980. (Cited on page 24.)
[38] E. Hazan. *Introduction to online convex optimization*, volume 2. Now Publishers, Inc., 2016. (Cited on page 2.)

[39] C. Hennig. Identifiability of Models for Clusterwise Linear Regression. *Journal of Classification*, 17(2):273–296, 2000. (Cited on page 8.)

[40] G. Hinton, N. Srivastava, and K. Swersky. Neural networks for machine learning: lecture 6a overview of mini-batch gradient descent. *University of Toronto, Technical Report*, 14(8):2, 2012. (Cited on pages 3 and 15.)

[41] N. Ho and X. Nguyen. Convergence rates of parameter estimation for some weakly identifiable finite mixtures. *The Annals of Statistics*, 44(6):2726–2755, 2016. (Cited on page 2.)

[42] N. Ho and X. Nguyen. On strong identifiability and convergence rates of parameter estimation in finite mixtures. *Electronic Journal of Statistics*, 10(1):271–307, 2016. (Cited on page 2.)

[43] N. Ho, C.-Y. Yang, and M. I. Jordan. Convergence Rates for Gaussian Mixtures of Experts. *Journal of Machine Learning Research*, 2022. (Cited on page 2.)

[44] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton. Adaptive mixtures of local experts. *Neural computation*, 3(1):79–87, 1991. (Cited on page 2.)

[45] A. Javanmard, S. Shao, and J. Bien. Prediction sets for high-dimensional mixture of experts models. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, page qkae117, Jan. 2025. (Cited on page 2.)

[46] H. Jiang and J. Xu. The stochastic proximal distance algorithm. *Statistics and Computing*, 34:210, 2024. (Cited on page 2.)

[47] W. Jiang and M. A. Tanner. Hierarchical mixtures-of-experts for exponential family regression models: approximation and maximum likelihood estimation. *Annals of Statistics*, pages 987–1011, 1999. (Cited on page 2.)

[48] W. Jiang and M. A. Tanner. On the identifiability of mixtures-of-experts. *Neural Networks*, 12(9):1253–1258, 1999. (Cited on page 8.)

[49] M. I. Jordan and R. A. Jacobs. Hierarchical mixtures of experts and the EM algorithm. *Neural computation*, 6(2):181–214, 1994. (Cited on page 2.)

[50] B. Karimi, B. Miasojedow, E. Moulines, and H.-T. Wai. Non-asymptotic Analysis of Biased Stochastic Approximation Scheme. In A. Beygelzimer and D. Hsu, editors, *Proceedings of the Thirty-Second Conference on Learning Theory*, volume 99 of *Proceedings of Machine Learning Research*, pages 1944–1974. PMLR, June 2019. (Cited on page 2.)

[51] B. Karimi, H.-T. Wai, E. Moulines, and M. Lavielle. On the Global Convergence of (Fast) Incremental Expectation Maximization Methods. In *Advances in Neural Information Processing Systems*, volume 32, 2019. (Cited on pages 2 and 13.)

[52] B. Karimi, H.-T. Wai, E. Moulines, and P. Li. Minimization by Incremental Stochastic Surrogate Optimization for Large Scale Nonconvex Problems. In S. Dasgupta and N. Haghtalab, editors, *Proceedings of The 33rd International Conference on Algorithmic Learning Theory*, volume 167 of *Proceedings of Machine Learning Research*, pages 606–637. PMLR, Apr. 2022. (Cited on pages 13 and 14.)

[53] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In *The 3rd International Conference for Learning Representations*, 2015. (Cited on pages 3 and 15.)

[54] E. Kuhn, C. Matias, and T. Rebafka. Properties of the stochastic approximation EM algorithm with mini-batch sampling. *Statistics and Computing*, 30(6):1725–1739, Nov. 2020. (Cited on pages 2 and 13.)

[55] H. Kushner and G. G. Yin. *Stochastic Approximation and Recursive Algorithms and Applications*, volume 35. Springer Science & Business Media, 2003. (Cited on pages 2 and 6.)

[56] G. Lan. *First-order and stochastic optimization methods for machine learning*, volume 1. Springer, 2020. (Cited on page 2.)

[57] K. Lange. *Numerical analysis for statisticians*, volume 1. Springer, 2010. (Cited on page 37.)

[58] K. Lange. *MM Optimization Algorithms*. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2016. (Cited on page 4.)
[59] M. Le, A. Nguyen, H. Nguyen, T. Nguyen, T. Pham, L. Van Ngo, and N. Ho. Mixture of Experts Meets Prompt-Based Continual Learning. In *Advances in Neural Information Processing Systems*, 2024. arXiv preprint arXiv:2405.14124. (Cited on page 2.)

[60] H. A. Le Thi and V. T. Ho. Online Learning Based on Online DCA and Application to Online Classification. *Neural Computation*, 32(4):759-793, Apr. 2020. (Cited on page 2.)

[61] H. A. Le Thi and V. T. Ho. DCA for online prediction with expert advice. *Neural Computing and Applications*, 33(15):9521-9544, Aug. 2021. (Cited on page 2.)

[62] H. A. Le Thi, H. P. H. Luu, and T. P. Dinh. Online Stochastic DCA With Applications to Principal Component Analysis. *IEEE Transactions on Neural Networks and Learning Systems*, 35(5):7035-7047, 2024. (Cited on page 2.)

[63] H. Liu, Z. Li, D. L. W. Hall, P. Liang, and T. Ma. Sophia: A Scalable Stochastic Second-order Optimizer for Language Model Pre-training. In *The Twelfth International Conference on Learning Representations*, 2024. (Cited on pages 3 and 15.)

[64] D. Lupu and I. Necoara. Convergence analysis of stochastic higher-order majorization-minimization algorithms. *Optimization Methods and Software*, 39:384-413, 2024. (Cited on page 2.)

[65] H. Lyu. Stochastic regularized majorization-minimization with weakly convex and multi-convex surrogates. *Journal of Machine Learning Research*, 25:1-83, 2024. (Cited on page 2.)

[66] H. Lyu, C. Strohmeier, and D. Needell. Online Nonnegative CP-dictionary Learning for Markovian Data. *Journal of Machine Learning Research*, 23(148):1-50, 2022. (Cited on page 2.)

[67] J. Mairal. Stochastic majorization-minimization algorithms for large-scale optimization. *Advances in Neural Information Processing Systems*, 26, 2013. (Cited on page 13.)

[68] J. Mairal. Incremental Majorization-Minimization Optimization with Application to Large-Scale Machine Learning. *SIAM Journal on Optimization*, 25(2):829-855, 2015. (Cited on pages 2 and 13.)

[69] S. Masoudnia and R. Ebrahimpour. Mixture of experts: a literature survey. *Artificial Intelligence Review*, 42(2):275-293, 2014. (Cited on page 2.)

[70] E. F. Mendes and W. Jiang. On convergence rates of mixtures of polynomial experts. *Neural computation*, 24(11):3025-3051, 2012. (Cited on page 8.)

[71] A. Mensch, J. Mairal, B. Thirion, and G. Varoquaux. Stochastic Subsampling for Factorizing Huge Matrices. *IEEE Transactions on Signal Processing*, 66(1):113-128, 2018. (Cited on page 2.)

[72] H. Nguyen, P. Akbarian, T. Nguyen, and N. Ho. A General Theory for Softmax Gating Multinomial Logistic Mixture of Experts. In *Proceedings of The 41st International Conference on Machine Learning*, 2024. (Cited on pages 2 and 11.)

[73] H. Nguyen, T. Nguyen, and N. Ho. Demystifying Softmax Gating in Gaussian Mixture of Experts. In *Advances in Neural Information Processing Systems*, Dec. 2023. (Cited on pages 2 and 8.)

[74] H. Nguyen, T. Nguyen, K. Nguyen, and N. Ho. Towards Convergence Rates for Parameter Estimation in Gaussian-gated Mixture of Experts. In *Proceedings of The 27th International Conference on Artificial Intelligence and Statistics*, volume 238, pages 2683-2691, May 2024. (Cited on page 2.)

[75] H. D. Nguyen and F. Chamroukhi. Practical and theoretical aspects of mixture-of-experts modeling: An overview. *Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery*, 8(4):e1246, 2018. (Cited on page 2.)

[76] H. D. Nguyen, F. Chamroukhi, and F. Forbes. Approximation results regarding the multiple-output Gaussian gated mixture of linear experts model. *Neurocomputing*, 366:208-214, 2019. (Cited on page 2.)

[77] H. D. Nguyen, F. Forbes, G. Fort, and O. Cappé. An Online Minorization-Maximization Algorithm. In *Classification and Data Science in the Digital Age - 17th Conference of the International Federation of Classification Societies, IFCS 2022, Proceedings*, pages 263-271. Springer, 2023. (Cited on pages 2, 4, and 5.)
[78] H. D. Nguyen, F. Forbes, and G. J. McLachlan. Mini-batch learning of exponential family finite mixture models. *Statistics and Computing*, 30(4):731–748, July 2020. (Cited on page 2.)

[79] H. D. Nguyen, L. R. Lloyd-Jones, and G. J. McLachlan. A universal approximation theorem for mixture-of-experts models. *Neural computation*, 28(12):2585–2593, 2016. (Cited on page 2.)

[80] H. D. Nguyen, T. Nguyen, F. Chamroukhi, and G. J. McLachlan. Approximations of conditional probability density functions in Lebesgue spaces via mixture of experts models. *Journal of Statistical Distributions and Applications*, 8(1):13, Aug. 2021. (Cited on page 2.)

[81] H. D. Nguyen, T. Nguyen, and F. Forbes. Bayesian Likelihood Free Inference using Mixtures of Experts. In *2024 International Joint Conference on Neural Networks (IJCNN)*, pages 1–8, 2024. (Cited on page 2.)

[82] T. Nguyen. *Model Selection and Approximation in High-dimensional Mixtures of Experts Models: from Theory to Practice*. PhD Thesis, Normandie Université, Dec. 2021. (Cited on page 2.)

[83] T. Nguyen, F. Chamroukhi, H. D. Nguyen, and G. J. McLachlan. Approximation of probability density functions via location-scale finite mixtures in Lebesgue spaces. *Communications in Statistics - Theory and Methods*, pages 1–12, May 2022. (Cited on page 2.)

[84] T. Nguyen, H. D. Nguyen, F. Chamroukhi, and F. Forbes. A non-asymptotic approach for model selection via penalization in high-dimensional mixture of experts models. *Electronic Journal of Statistics*, 16(2):4742–4822, 2022. (Cited on page 2.)

[85] T. Nguyen, H. D. Nguyen, F. Chamroukhi, and G. J. McLachlan. Approximation by finite mixtures of continuous density functions that vanish at infinity. *Cogent Mathematics & Statistics*, 7(1):1750861, 2020. (Cited on page 2.)

[86] X. Nguyen. Convergence of latent mixing measures in finite and infinite mixture models. *The Annals of Statistics*, 41(1):370–400, 2013. (Cited on page 2.)

[87] A. Norets. Approximation of conditional densities by smooth mixtures of regressions. *The Annals of Statistics*, 38(3):1733–1766, 2010. (Cited on page 2.)

[88] G. Oudoumanessah, T. Coudert, C. Lartizien, M. Dojat, T. Christen, and F. Forbes. Scalable magnetic resonance fingerprinting: Incremental inference of high dimensional elliptical mixtures from large data volumes, 2024. (Cited on page 2.)

[89] G. Oudoumanessah, T. Coudert, L. Meyer, A. Delphin, M. Dojat, C. Lartizien, and F. Forbes. Cluster globally, reduce locally: Scalable efficient dictionary compression for magnetic resonance fingerprinting. In *IEEE International Symposium on Biological Imaging (ISBI)*, 2025. (Cited on page 2.)

[90] S. A. Prado, L. Cabrera-Bosquet, A. Grau, A. Coupel-Ledru, E. J. Millet, C. Welcker, and F. Tardieu. Phenomics allows identification of genomic regions affecting maize stomatal conductance with conditional effects of water deficit and evaporative demand. *Plant, Cell & Environment*, 41(2):314–326, 2018. (Cited on page 20.)

[91] A. Rakhlin, D. Panchenko, and S. Mukherjee. Risk bounds for mixture density estimation. *ESAIM: PS*, 9:220–229, 2005. (Cited on page 2.)

[92] M. Razaviyayn, M. Sanjabi, and Z.-Q. Luo. A Stochastic Successive Minimization Method for Nonsmooth Nonconvex Optimization with Applications to Transceiver Design in Wireless Communication Networks. *Mathematical Programming*, 157(2):515–545, June 2016. (Cited on page 2.)

[93] M. Redmond. Communities and Crime. UCI Machine Learning Repository, 2002. DOI: https://doi.org/10.24432/C53W3X. (Cited on page 20.)

[94] H. Robbins and S. Monro. A Stochastic Approximation Method. *The Annals of Mathematical Statistics*, 22(3):400–407, 1951. Publisher: Institute of Mathematical Statistics. (Cited on pages 3, 5, and 15.)

[95] D. Rohde and O. Cappé. Online maximum-likelihood estimation for latent factor models. In *2011 IEEE Statistical Signal Processing Workshop (SSP)*, pages 565–568. IEEE, 2011. (Cited on page 13.)

[96] S. Shalev-Shwartz. Online Learning and Online Convex Optimization. *Foundations and Trends® in Machine Learning*, 4(2):107–194, 2012. (Cited on page 2.)
[97] N. Shazeer, A. Mirhoseini, K. Maziarz, A. Davis, Q. Le, G. Hinton, and J. Dean. Outrageously Large Neural Networks: The Sparsely-Gated Mixture-of-Experts Layer. In *International Conference on Learning Representations*, 2017. (Cited on page 2.)

[98] T. Thai, T. Nguyen, D. Do, N. Ho, and C. Drovandi. Model Selection for Gaussian-gated Gaussian Mixture of Experts Using Dendrograms of Mixing Measures. *arXiv preprint arXiv:2505.13052*, 2025. (Cited on page 2.)

[99] R. Tibshirani. Regression shrinkage and selection via the Lasso. *Journal of the Royal Statistical Society: Series B (Methodological)*, 58(1):267–288, 1996. Publisher: Wiley Online Library. (Cited on page 20.)

[100] Z. Wang, K. Ji, Y. Zhou, Y. Liang, and V. Tarokh. SpiderBoost and Momentum: Faster Variance Reduction Algorithms. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d. Alché-Buc, E. Fox, and R. Garnett, editors, *Advances in Neural Information Processing Systems*, volume 32. Curran Associates, Inc., 2019. (Cited on page 2.)

[101] S. E. Yuksel, J. N. Wilson, and P. D. Gader. Twenty Years of Mixture of Experts. *IEEE Transactions on Neural Networks and Learning Systems*, 23(8):1177–1193, 2012. (Cited on page 2.)