Paper deep dive
REALM: Regime-Switching, Explainable, and Activation-Induced Linear Models
Xiaoran Cheng, Sen Na, Jia Li
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 9/29/2026, 12:45:59 PM
Summary
The paper introduces REALM, a Regime-Switching, Explainable, and Activation-Induced Linear Model framework designed to create interpretable surrogates for Deep Neural Networks (DNNs). REALM distills a deep teacher network into a Wide Shallow Student Network (WSSN), clusters the WSSN's binarized activation patterns to define stable regimes, and fits local linear models within each regime. An explanatory gate (multiclass logistic regression) is trained to interpret regime assignments, providing dual interpretability at the feature and regime levels while maintaining competitive predictive performance.
Entities (8)
Relation Signals (6)
Deep Neural Network → hasstructure → Piecewise-Affine
confidence 95% · Deep ReLU networks are piecewise-affine mappings that partition the input space into cells
REALM → includescomponent → Explanatory Gate
confidence 95% · we train a multiclass logistic regression, the explanatory gate, to reproduce the regime assignments.
WSSN → isdistilledfrom → Deep Neural Network
confidence 95% · we first distill a deep teacher into a wide, shallow student network (WSSN)
REALM → uses → Wide Shallow Student Network
confidence 95% · we first distill a deep teacher into a wide, shallow student network (WSSN)... REALM then clusters the WSSN’s binarized hidden activation patterns
REALM → achieves → Interpretability
confidence 90% · This structure motivates fitting a local linear model within each cell to preserve predictive accuracy while improving interpretability.
Explanatory Gate → interprets → Regime Assignment
confidence 90% · To make regime assignment interpretable, we train a multiclass logistic regression, the explanatory gate, to reproduce the regime assignments.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Deep ReLU networks are piecewise-affine mappings that partition the input space into cells, each characterized by a distinct activation pattern. This structure motivates fitting a local linear model within each cell to preserve predictive accuracy while improving interpretability. The challenge is to identify regimes that are stable, data-adaptive, and easy to explain. We propose REALM, a mixture of linear models whose regimes are induced by neural activation patterns. Because the number of activation cells in a deep neural network (DNN) can grow rapidly with depth, we first distill a deep teacher into a wide, shallow student network (WSSN), then binarize and cluster its hidden-layer activations to define the regimes and fit a linear model within each regime. Since the regimes are discovered from internal structure, the router does not carry the predictive burden. To make regime assignment interpretable, we train a multiclass logistic regression, the explanatory gate, to reproduce the regime assignments. The two-level structure is interpretable at both stages in terms of raw tabular or learned convolutional features: the gate identifies features that determine regime assignments, while the linear models identify features that drive predictions within each regime. We analyze an idealized setting that illustrates a trade-off between partition complexity and stability: as the number of regimes grows, finer partitions can improve approximation but may reduce regime-assignment stability. Experiments on tabular and image datasets show that REALM achieves competitive predictive performance relative to other DNN-guided mixture surrogates and inherently interpretable models while producing stable regime-level explanations.
Tags
Links
- Source: https://arxiv.org/abs/2609.32141v1
- Canonical: https://arxiv.org/abs/2609.32141v1
Trouble viewing inline? Open PDF directly →
Full Text
91,187 characters extracted from source content.
Expand or collapse full text
REALM: Regime-Switching, Explainable, and Activation-Induced Linear Models Xiaoran Cheng Sen Na Jia Li †thanks: X. Cheng and J. Li are with the Department of Statistics, The Pennsylvania State University, University Park, PA 16802 USA (e-mail: xjc5161@psu.edu; jol2@psu.edu).†thanks: S. Na is with the H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA 30332 USA (e-mail: senna@gatech.edu). Abstract Deep ReLU networks are piecewise-affine mappings that partition the input space into cells, each characterized by a distinct activation pattern. This structure motivates fitting a local linear model within each cell to preserve predictive accuracy while improving interpretability. The challenge is to identify regimes that are stable, data-adaptive, and easy to explain. We propose REALM, a mixture of linear models whose regimes are induced by neural activation patterns. Because the number of activation cells in a deep neural network (DNN) can grow rapidly with depth, we first distill a deep teacher into a wide, shallow student network (WSSN), then binarize and cluster its hidden-layer activations to define the regimes and fit a linear model within each regime. Since the regimes are discovered from internal structure, the router does not carry the predictive burden. To make regime assignment interpretable, we train a multiclass logistic regression, the explanatory gate, to reproduce the regime assignments. The two-level structure is interpretable at both stages in terms of raw tabular or learned convolutional features: the gate identifies features that determine regime assignments, while the linear models identify features that drive predictions within each regime. We analyze an idealized setting that illustrates a trade-off between partition complexity and stability: as the number of regimes grows, finer partitions can improve approximation but may reduce regime-assignment stability. Experiments on tabular and image datasets show that REALM achieves competitive predictive performance relative to other DNN-guided mixture surrogates and inherently interpretable models while producing stable regime-level explanations. Index Terms: Adaptive modeling, interpretable models, knowledge distillation, mixture of experts, ReLU networks. I Introduction DNNs are now standard tools in vision, language, and scientific computing, but their black-box nature remains a major barrier in settings that demand auditing, debugging, or accountability [1, 2, 3, 4]. A common strategy is to replace a trained DNN with an interpretable surrogate, for example by distilling the DNN into a soft decision tree [5]. However, a single global surrogate often struggles to capture the complex decision logic of a highly nonlinear model without sacrificing predictive accuracy. A more expressive alternative is to use a collection of simple models, each responsible for a distinct regime of the input space. This divide-and-conquer view motivates adaptive surrogate modeling: a complex model can be approximated by partitioning the input space into regimes and fitting a simple local model within each regime. The central design question then becomes how to construct the partition: what criterion should define the regimes so that the surrogate remains accurate while both the regime partitioning and the local models remain interpretable? Mixture-of-experts (MoE) models [6, 7] provide a natural framework for adaptive surrogate modeling. An MoE uses a router to assign each input to one or a few simple local experts. However, in classical MoE models, the router and experts are jointly optimized for predictive accuracy, so the resulting partition is not explicitly designed for interpretability [6, 7]. Subsequent variants have introduced explicit input partitioning [8], Bayesian non-parametric mixtures for explanation [9], boosted mixtures of experts [10], and sparsely gated layers for conditional computation [11]. These methods expand the flexibility and scalability of MoE models, but they still do not provide a fully interpretable prediction path for a trained DNN. The central challenge is that the router can absorb substantial predictive logic. Because the router and experts are trained jointly, the resulting decision logic can become complex [12, 13], making it difficult to disentangle whether a prediction is driven by the router, the experts, or their interaction. Large sparse MoE models can also suffer from router-induced training instability [14]. Constraining either component to be inherently interpretable [15, 16, 17] often comes at the cost of predictive accuracy. Interpretability may further deteriorate under diffuse routing, where a single input activates many experts with similar weights, obscuring feature attribution. Moreover, when the partition is learned directly from the input, it may define a new predictive structure of its own rather than reflect the internal regime structure of the trained DNN [6, 7, 16, 11, 15]. In this case, the router itself carries part of the predictive burden. These limitations suggest that an interpretable DNN surrogate should not learn the partition solely through the inputs by optimizing prediction accuracy jointly with the experts. Instead, the regimes should align with the DNN’s own internal geometric structure. ReLU networks provide a natural signal for this purpose. A trained ReLU network partitions the input space into polyhedral cells, within each of which the network reduces to an affine map [18, 19, 20, 21, 22, 23]. Thus, its activation patterns provide a natural basis for defining regimes that support locally linear modeling. Since ReLU is widely used in modern architectures [24], this structure can be exploited to construct interpretable surrogates. However, directly using the exact partition induced by activation patterns is impractical. The number of activation cells grows rapidly with network depth [22, 23], and a trained network may generate nearly instance-specific activation patterns [19, 25]. This limits interpretability because the resulting explanations become highly localized and difficult to compare across instances. Several families of interpretable surrogates address parts of this challenge. IME [16] constructs a mixture in which both the experts and the assignment module are interpretable, while the route–interpret–repeat framework [17] iteratively decomposes a black-box model into interpretable models, each covering the samples it can explain. However, in both approaches, the partition is still learned through routing on the inputs, so the resulting regimes are determined by the routing objective rather than derived from the internal structure of the trained network. More closely tied to ReLU network geometry, ReLU Region Reasoning (Re3) [25] computes the exact per-sample affine map associated with each activation pattern, while MASALA [26] fits a surrogate over an adaptively chosen neighborhood around each instance. Although both methods exploit local structure, neither coarsens the exact activation cells into a shared set of regimes. This suggests a natural direction for improving the interpretability of geometry-based surrogates: aggregate the exact activation cells into a small number of shared, interpretable regimes. The methods most closely related to ours are DNN-guided mixture surrogates, particularly MLM (Mixture of Linear Models Co-supervised by DNNs) [27] and SEE-Net (Synced Explanation-Enhanced N) [28]. Both methods cluster continuous hidden-layer representations of the DNN layer by layer, combine the resulting per-layer clusters through Cartesian products, and fit local linear models within the resulting regimes. Because these methods cluster hidden-layer values rather than activation patterns, inputs that share the same activation pattern may still be assigned to different regimes. SEE-Net further transforms the DNN into a pre-interpretable architecture, which can incur additional training cost. In contrast, our proposed REALM uses the trained DNN only as a teacher for distillation and data augmentation, and directly clusters the binarized activation patterns of a WSSN. Consequently, its regimes are formed as unions of the WSSN’s activation cells. Overview of REALM We propose Regime-Switching, Explainable, and Activation-Induced Linear Models (REALM). Given a trained DNN teacher fTf_T, REALM first distills fTf_T into a WSSN fSf_S to match the teacher’s soft predictions. This step is motivated by the fact that wide and shallow networks can faithfully mimic deep networks under knowledge distillation [29, 30, 31]. As a ReLU network, the WSSN also partitions its input space into polyhedral activation cells within which it is affine. Its single hidden layer induces far fewer cells than the deep teacher [22, 23, 32], and its activation patterns are lower-dimensional than the teacher’s while remaining aligned with the teacher network. The reduced dimension stabilizes regime discovery since distance-based clustering becomes unreliable as the dimension grows [33, 34]. REALM then clusters the WSSN’s binarized hidden activation patterns into regimes and fits a regularized local linear expert within each regime. The final prediction is a weighted average of the expert outputs, with weights given by the posterior from the clustering step. To make interpretation more faithful, we employ an inverse temperature that sharpens the posterior toward hard routing, so that each prediction relies on only a few experts. To make regime assignment easy to interpret, we additionally train a regularized multiclass logistic regression, called the explanatory gate, that predicts the regime assignment label directly from the input representation. In the basic REALM formulation used for the main results in Section V, the explanatory gate serves only to explain regime assignments rather than to determine routing. In Section V-C, we also evaluate a variant, REALM-EG, that uses the explanatory gate itself for routing and achieves performance comparable to that obtained by directly clustering the WSSN’s binarized activation patterns. REALM achieves competitive predictive performance relative to SEE-Net across six datasets and to MLM and common inherently interpretable models on tabular data, while producing stable regime-specific explanations through its linear experts and explanatory gate. Main Contributions • Interpretable partitioning from ReLU networks’ internal structure. We propose REALM, an interpretable MoE framework that distills a trained DNN into a WSSN and clusters its low-dimensional, teacher-aligned ReLU activation patterns into a small set of regimes, yielding dual interpretability with respect to raw tabular or learned convolutional features: 1) local feature attribution within each regime; 2) an explanatory gate that attributes each regime assignment to these features. • An idealized analysis of the complexity–stability trade-off. In an idealized setting, we analyze when regime assignments stay stable under input perturbations, and, assuming sparse within-regime coefficients, when LASSO recovers the correct support. This analysis reveals a trade-off between partition complexity and stability: finer partitions can improve approximation but may reduce regime-assignment stability. • Practical pipeline with cross-modality validation. We benchmark REALM against SEE-Net on both tabular and image tasks, and against MLM and other interpretable models on tabular tasks, demonstrating competitive accuracy with stable, regime-specific interpretations. Ablations show that 1) the router does not carry the predictive burden; 2) the explanatory gate can serve as the routing signal with little performance loss; 3) the results for different K values are mostly consistent with the complexity–stability trade-off. The remainder of this paper is organized as follows. Section I reviews the piecewise-affine structure of ReLU networks. Section I introduces REALM and its learning algorithm. Section IV illustrates a trade-off between partition complexity and stability by analyzing an idealized setting. Section V reports experiments and ablation studies, and Section VI concludes the paper. I Preliminaries: Piecewise-Affine Structure of ReLU DNNs Although a ReLU network is globally nonlinear, once the active or inactive status of every hidden unit is fixed, each activation pattern indexes a cell of the input space on which the network reduces to a single affine map [18, 21, 19, 25]. Thus, instead of learning an arbitrary partition through joint routing–prediction optimization, REALM derives interpretable partitioning using activation patterns to discover the latent regimes and fits simple local models on the resulting regimes. Consider a ReLU network with input x∈=ℝpx =R^p, hidden layers ℓ=1,…,L =1,…,L, and output dimension O: h(0)(x) h^(0)(x) =x, =x, a(ℓ)(x) a^( )(x) =W(ℓ)h(ℓ−1)(x)+b(ℓ), =W^( )h^( -1)(x)+b^( ), h(ℓ)(x) h^( )(x) =ρ(a(ℓ)(x)), =ρ (a^( )(x) ), where ρ(t)=max0,tρ(t)= \0,t\. The output layer is f(x)=W(L+1)h(L)(x)+b(L+1).f(x)=W^(L+1)h^(L)(x)+b^(L+1). Here n0=pn_0=p denotes the input dimension, nℓn_ is the width of hidden layer ℓ , and nL+1=On_L+1=O is the output dimension. Thus, h(ℓ−1)(x)∈ℝnℓ−1h^( -1)(x) ^n_ -1, a(ℓ)(x),h(ℓ)(x),b(ℓ)∈ℝnℓa^( )(x),h^( )(x),b^( ) ^n_ , W(ℓ)∈ℝnℓ×nℓ−1W^( ) ^n_ × n_ -1, and f(x),b(L+1)∈ℝOf(x),b^(L+1) ^O, W(L+1)∈ℝO×nLW^(L+1) ^O× n_L. Each hidden unit is either on (ai(ℓ)(x)>0a^( )_i(x)>0) or off (ai(ℓ)(x)≤0a^( )_i(x)≤ 0). To track these states across the network, we use the activation pattern [23, 21]. Definition I.1 (Activation Pattern [21]). For each hidden layer, define the activation mask s(ℓ)(x)=a(ℓ)(x)>0,D(ℓ)(x)=diag(s(ℓ)(x)), s^( )(x)=1 \a^( )(x)>0 \,\;D^( )(x)=diag (s^( )(x) ), (1) where s(ℓ)(x)∈0,1nℓs^( )(x)∈\0,1\^n_ records which units in layer ℓ are active and D(ℓ)(x)∈ℝnℓ×nℓD^( )(x) ^n_ × n_ is the corresponding diagonal selector matrix. Each ReLU layer can then be written as h(ℓ)(x)=D(ℓ)(x)a(ℓ)(x)=D(ℓ)(x)(W(ℓ)h(ℓ−1)(x)+b(ℓ)), splith^( )(x)&=D^( )(x)\,a^( )(x)\\ &=D^( )(x) (W^( )h^( -1)(x)+b^( ) ), split so the only source of nonlinearity is the variation of the masks D(ℓ)(x)D^( )(x). Once they are fixed, every layer becomes an affine map [18, 21, 19, 25]. Let s(x)=(s(1)(x),…,s(L)(x))∈0,1∑ℓ=1Lnℓs(x)= (s^(1)(x),…,s^(L)(x) )∈\0,1\ _ =1^Ln_ denote the full activation pattern, and let exp≜s(x):x∈S_exp \s(x):x \ be the set of distinct patterns realized on X [21, 25]. Since forward propagation assigns every input a unique pattern s(x)s(x), each pattern s selects an activation cell ℛs=x∈:s(x)=sR_s=\x :s(x)=s\, an intersection of finitely many half-spaces and hence a convex polyhedral region. Since distinct patterns give disjoint cells, they partition the input space such that =⨆s∈expℛsX= _s _expR_s [21, 22, 25]. However, this exact partition ℛsR_s can be too fine-grained to support interpretable regime assignment, since the number of cells can grow rapidly with network depth [22]. An upper bound on the number of cells is (nmaxpL)O(n_ ^pL), where nmax=max1≤ℓ≤Lnℓn_ = _1≤ ≤ Ln_ denotes the maximum hidden-layer width [23]. This fine partition may produce many tiny or even singleton cells [22, 21, 19, 25]. Practical local surrogate modeling therefore requires collapsing these cells into a small number of stable, interpretable regimes, each a union of activation cells on which the network is approximately affine. This motivates the activation-induced clustering of REALM developed in Section I. I REALM: Regime-Switching, Explainable, and Activation-Induced Linear Models REALM has three stages: 1) distilling a DNN teacher fTf_T into a WSSN fSf_S, which has teacher-aligned but lower-dimensional activation patterns; 2) discovering K latent regimes by clustering the WSSN’s activation patterns; 3) fitting a mixture of local experts together with a separate explanatory gate, making both feature attributions and regime assignments interpretable. A schematic illustration of REALM is provided in Fig. 1. Fig. 1: Overview of the REALM pipeline. Stage 1 distills a deep teacher fTf_T into the WSSN fSf_S. Stage 2 clusters the WSSN’s binarized activation patterns into K regimes, giving the routing posterior πk(x) _k(x). Stage 3 fits a mixture of local linear experts under πk(x) _k(x) to produce the prediction f^(x) f(x), and trains the explanatory gate q^(x) q(x) to explain the assignments. Stage 1: Teacher-to-Student Distillation Let u(x)u(x) denote the input to the student WSSN: u(x)=x,for tabular data x∈ℝp,φ(x),for image data x∈ℝH×W×C.u(x)= casesx,&for tabular data x ^p,\\ (x),&for image data x ^H× W× C. cases Here φ(x)∈ℝp (x) ^p denotes the output of the WSSN’s convolutional feature extractor. Throughout, we write u(x)∈ℝpu(x) ^p. The extractor φ is trained jointly with the WSSN during distillation and is frozen afterwards. Therefore, regime discovery, the local experts, and the explanatory gate all operate on the fixed representation u(x)u(x). Note that for image data, local experts and the explanatory gate are linear in u(x)u(x) rather than in the raw pixels. We then pass u(x)u(x) through the one-hidden-layer ReLU head, with pre-activation a(x)=W1u(x)+b1a(x)=W_1u(x)+b_1 and hidden activation h(x)=ρ(a(x))h(x)=ρ(a(x)), giving fS(x)=W2h(x)+b2.f_S(x)=W_2h(x)+b_2. We detail the construction of φ(x) (x) for each image dataset in the supplementary material. For both tabular and image data, the WSSN has a single hidden layer (excluding the feature-extraction layers for image data), so it induces only a manageable number of activation cells. For classification, fSf_S is trained by knowledge distillation from fTf_T: ℒKDcls _KD^cls =αT2ℒKL(softmax(fT(x)T)∥softmax(fS(x)T)) =α T^2\,L_KL (softmax ( f_T(x)T )\, \|\,softmax ( f_S(x)T ) ) (2) +(1−α)ℒCE(fS(x),y), +(1-α)\,L_CE (f_S(x),y ), where ℒKLL_KL is the Kullback–Leibler divergence between the WSSN and the teacher distributions, and ℒCEL_CE is the cross-entropy against the label y. The temperature T>0T>0 controls the softness of the softmax: T>1T>1 smooths the output distribution and facilitates transfer of dark knowledge, while T→0+T→ 0^+ sharpens it toward a one-hot vector. The weight α∈[0,1]α∈[0,1] balances distillation and label-matching. For regression, we apply the mean-squared error: ℒKDreg=α‖fS(x)−fT(x)‖22+(1−α)‖fS(x)−y‖22.L_KD^reg=α\,\|f_S(x)-f_T(x)\|_2^2+(1-α)\,\|f_S(x)-y\|_2^2. Stage 2: Soft Clustering for Latent Regime Discovery We freeze the WSSN after distillation, so its activation pattern s(x)s(x) is fixed as defined in Definition I.1. Therefore, the WSSN is an affine function of its input u(x)u(x) within each exact activation cell. We then discover K regimes by soft clustering with learnable centroids μ1,…,μK∈ℝnhid _1,…, _K ^n_hid, where nhidn_hid is the hidden width of the WSSN, and the weighted distance: Dk(x)=1nhid∑j=1nhid[1+(ω−1)sj(x)](sj(x)−μk,j)2,D_k(x)= 1n_hid _j=1^n_hid [1+(ω-1)s_j(x) ] (s_j(x)- _k,j )^2, where ω>1ω>1 upweights mismatches on active units (weight ω when sj(x)=1s_j(x)=1 and weight 11 otherwise). The WSSN-induced posterior uses an inverse temperature βclus>0 _clus>0: πk(x)=exp−βclusDk(x)∑j=1Kexp−βclusDj(x). _k(x)= \- _clusD_k(x)\ _j=1^K \- _clusD_j(x)\. (3) A small βclus _clus produces flatter regime assignments, whereas a large βclus _clus concentrates the assignments on a few regimes. Near-one-hot assignment enhances interpretability, since each prediction only draws from a few regimes and its feature attribution localizes to a small set of experts. To avoid imposing near-one-hot assignments at the initialization of the clustering algorithm, we progressively increase the inverse temperature according to βclus(t+1)=minηβclus(t),βmax, _clus^(t+1)= \η\, _clus^(t),\, _ \, (4) where t indexes clustering epochs, η>1η>1 is the sharpening factor, and βmax _ is the upper bound on the inverse temperature. For a fixed K, the centroids are learned by minimizing ℒclus=1n∑i=1n∑k=1Kπk(xi)Dk(xi),L_clus= 1n _i=1^n _k=1^K _k(x_i)D_k(x_i), where centroids are initialized from K randomly chosen training activation patterns. We update the centroids by differentiating through both Dk(xi)D_k(x_i) and πk(xi) _k(x_i). Define the regime assignment as z(x)=argmaxkπk(x)z(x)= argmax_k _k(x). Thus, each regime k=x:z(x)=kC_k=\x:z(x)=k\ is a union of activation cells on which the WSSN is piecewise-affine. Since clustering merges cells with similar activation patterns, the WSSN remains approximately affine within each regime, as supported empirically by the strong performance in Section V. Because the WSSN is frozen, clustering merely summarizes how the distilled WSSN organizes the input space rather than absorbing the predictive burden into the router. This is shown in the ablation study in Section V-C. Stage 3: Fixed-Router Mixture of Regularized Local Experts Local Experts. With the posterior πk(x) _k(x) fixed, we train a mixture of local linear experts. Each regime-k expert is: ek(x)=softmax(Wk⊤u(x)+bk),classification,wk⊤u(x)+bk,regression.e_k(x)= casessoftmax(W_k u(x)+b_k),&classification,\\ w_k u(x)+b_k,®ression. cases (5) For O-class classification, Wk∈ℝp×OW_k ^p× O and bk∈ℝOb_k ^O, whereas for regression wk∈ℝpw_k ^p and bk∈ℝb_k . The final prediction aggregates the experts under the posterior πk(x) _k(x): f^(x)=∑k=1Kπk(x)ek(x). f(x)= _k=1^K _k(x)\,e_k(x). (6) To optimize the expert parameters, let ℒyL_y be the predictive loss of the fixed-router mixture, with ℒy=1n∑i=1nℒCE(f^(xi),yi),classification,1n∑i=1nℒMSE(f^(xi),yi),regression.L_y= cases 1n _i=1^nL_CE ( f(x_i),y_i ),&classification,\\ 1n _i=1^nL_MSE ( f(x_i),y_i ),®ression. cases We estimate the expert parameters through minθyℒy(θy)+λyΩy(θy), _ _y _y( _y)+ _y _y( _y), (7) where θy _y collects the local linear experts’ parameters in Eq. (5) across all K regimes, θy=(wk,bk)k=1K,regression,(Wk,bk)k=1K,classification, _y= cases\(w_k,b_k)\_k=1^K,®ression,\\ \(W_k,b_k)\_k=1^K,&classification, cases and Ωy _y penalizes only the coefficients, Ωy(θy)=∑k=1K∑j=1p|wk,j|,regression,∑k=1K∑j=1p∑o=1O|(Wk)jo|,classification. _y( _y)= cases _k=1^K _j=1^p|w_k,j|,®ression,\\[6.0pt] _k=1^K _j=1^p _o=1^O|(W_k)_jo|,&classification. cases Explanatory Gate. The posterior πk(x) _k(x) routes inputs but does not explain why an input is assigned to a regime. To make the assignment rule interpretable, we train a separate multiclass logistic regression to predict the hard regime label zi=argmaxkπk(xi)z_i= argmax_k _k(x_i): q^(x)=softmax(Γu(x)+c), q(x)=softmax( u(x)+c), (8) where Γ=[γ1,…,γK]⊤∈ℝK×p =[ _1,…, _K] ^K× p, γk∈ℝp _k ^p, and c∈ℝKc ^K are trainable gate parameters. We collect the gate parameters as θg=Γ,c _g=\ ,c\. Since a common shift of all γk _k leaves q q unchanged, coefficients are read as contrasts between regimes. The explanatory gate is trained using the cross-entropy loss ℒg(θg)=1n∑i=1nℒCE(q^(xi),zi).L_g( _g)= 1n _i=1^nL_CE ( q(x_i),z_i ). We estimate the gate parameters through minθgℒg(θg)+λg‖Γ‖1, _ _g _g( _g)+ _g\| \|_1, (9) where ‖Γ‖1=∑k=1K∑j=1p|Γkj|\| \|_1= _k=1^K _j=1^p| _kj| denotes the entrywise ℓ1 _1 norm. Since both the experts and the explanatory gate are linear models, REALM yields dual interpretability: 1) the expert coefficients explain how each prediction is formed; 2) the gate coefficients explain why a sample is assigned to its regime. Teacher-Assisted Data Augmentation. To stabilize expert estimation in low-support regimes, we augment small regimes with teacher-labeled synthetic samples, giving the DNN a second role as a source of supervision beyond distillation. Let ℐk=i:z(xi)=kI_k=\i:z(x_i)=k\ denote the set of training samples assigned to regime k, and let nk=|ℐk|n_k=|I_k| be the number of original samples in regime k. If nk<τn_k<τ for a predefined minimum regime size τ, we generate nkaug=τ−nkn_k^aug=τ-n_k synthetic samples. We split the variables into a continuous block contJ_cont and a categorical (one-hot) block catJ_cat: xk,jaug[cont] x^aug_k,j[J_cont] =μ¯k+ϵσk⊙ξj,ξj∼(0,I|cont|), = μ_k+ε\, _k _j, _j (0,I_|J_cont|), (10) xk,jaug[cat] x^aug_k,j[J_cat] =xij[cat],ij∼Unif(ℐk), =x_i_j[J_cat], i_j (I_k), where μ¯k,σk∈ℝ|cont| μ_k, _k ^|J_cont| are the within-regime mean and standard deviation of the continuous variables, and ϵ>0ε>0 controls the variance magnitude. Continuous variables are thus simulated around the regime mean, while categorical variables are copied from a uniformly sampled real member of the regime to preserve valid encodings. Each synthetic sample is labeled by the teacher, yk,jaug=argmaxcfT(xk,jaug)c,classification,fT(xk,jaug),regression,y^aug_k,j= cases argmax_cf_T(x^aug_k,j)_c,&classification,\\[2.84526pt] f_T(x^aug_k,j),®ression, cases and appended to the training set. The experts and explanatory gate are then estimated through the objectives in Eqs. (7) and (9), respectively, with the augmented samples routed by the same posterior πk _k. Samples generated for regime k may be assigned to another regime. Regimes with nk=0n_k=0 are skipped. We augment only the tabular datasets because every regime exceeds the minimum size τ on the image datasets. REALM-EG. In REALM, predictions are weighted by the WSSN-induced posterior πk(x) _k(x) in Eq. (3). The explanatory gate q^(x) q(x) in Eq. (8) explains these assignments post hoc but never routes inputs in the main results in Section V. Promoting the gate to the router yields a variant we call REALM-EG (EG for explanatory gate). REALM-EG is then a two-layer linear model on the representation u(x)u(x): the gate reassigns every sample, zig=argmaxkq^k(xi),ℐkg=i:zig=k,z_i^g= argmax_k q_k(x_i), _k^g=\i:z_i^g=k\, and predictions are mixed by the gate posterior, f^EG(x)=∑k=1Kq^k(x)ekg(x), f_EG(x)= _k=1^K q_k(x)\,e_k^g(x), which is Eq. (6) with q^k q_k in place of πk _k. The local experts ekge_k^g are then refitted on all training samples by solving the same objective in Eq. (7) with f^(xi) f(x_i) replaced by f^EG(xi) f_EG(x_i). The hard labels zigz_i^g thus define the regimes ℐkgI_k^g of REALM-EG, whereas the refit uses the same soft mixture as REALM, so q^k(x) q_k(x) weights the experts both during training and at test time. Both the assignment rule and the within-regime predictions of REALM-EG are then directly readable from linear coefficients. The ablation study in Section V-C shows that REALM-EG performs competitively with REALM. IV Theoretical Analysis of Structural Stability and Support Recovery We analyze the joint reliability of regime assignment and sparse support recovery in REALM under an idealized setting where the data follow a latent piecewise sparse linear model whose regimes are separated by a linear winner-take-all regime selector. Here x represents u(x)u(x), the winner-take-all rule is a linear abstraction of regime assignment, and λ controls within-regime sparsity. The setting is idealized in three ways: 1) the regimes are treated as known; 2) the model is exactly linear within each of them, whereas REALM discovers regimes by clustering activation patterns and is only approximately linear inside each regime; 3) the true coefficients are sparse, whereas the fitted experts need not be. We allow the test-time regime selector to be learned, with error bounded by ηn _n, while retaining true regime membership for expert training. A recurring theme of the analysis is to make explicit how the number of regimes K affects both stability and overall support recovery. Note that this analysis also assumes a regression model with additive Gaussian noise. IV-A Model Setup Let x∈=ℝpx =R^p denote the input signal and suppose the data are generated from a latent piecewise sparse linear model: y=x⊤βk∗+ε,y=x _k^*^*+ , where k∗=g∗(x)k^*=g^*(x) denotes the true regime, βk∗∈ℝp _k^* ^p are sparse vectors with support Sk=supp(βk∗),|Sk|=sk,S_k=supp( _k^*), |S_k|=s_k, and ε∼(0,σ2) (0,σ^2). The true regime assignment is determined by a linear regime selector g∗g^* with gating vectors vk∗k=1K⊂ℝp\v_k^*\_k=1^K ^p, specified in Assumptions IV.1 and IV.3. For simplicity, the idealized selector omits intercepts. Define the pairwise margin: m(x)=minj≠k∗(vk∗−vj∗)⊤x.m(x)= _j≠ k^*(v_k^*^*-v_j^*) x. (11) We analyze stability at inputs with m(x)>0m(x)>0. This margin quantifies how much the winning regime’s gating score exceeds that of its closest competitor. IV-B Stability of Regime Assignment We study how reliably the true regime selector preserves its decision under additive perturbations of a fixed input signal. Throughout this subsection we treat x as a deterministic signal and take all probabilities with respect to the perturbation δ alone. For the probabilistic bounds, we assume K≥2K≥ 2 and σg>0 _g>0. Assumption IV.1 (Linear Gating, Deterministic Signal). The true regime selector is the linear winner-take-all rule g∗(x)=argmaxk∈1,…,Kvk∗⊤x,g^*(x)= argmax_k∈\1,…,K\v_k^* x, (12) with gating vectors vk∗k=1K\v_k^*\_k=1^K fixed. Assumption IV.2 (Additive Isotropic Perturbation). We observe x+δx+δ where δ∼(0,σg2Ip)δ (0, _g^2I_p) is independent of x. Proposition IV.1 (Deterministic Stability (fixed x)). Suppose Assumption IV.1 holds. Fix x with k∗=g∗(x)k^*=g^*(x) and margin m(x)>0m(x)>0. For any additive perturbation δ∈ℝpδ ^p, if ‖δ‖2<m(x)maxj≠k∗‖vk∗−vj∗‖2,\|δ\|_2< m(x) _j≠ k^*\|v_k^*^*-v_j^*\|_2, then g∗(x+δ)=g∗(x)g^*(x+δ)=g^*(x). Proof. For every j≠k∗j≠ k^*, the Cauchy–Schwarz inequality gives (vk∗−vj∗)⊤(x+δ) (v_k^*^*-v_j^*) (x+δ) ≥m(x)−‖vk∗−vj∗‖2‖δ‖2>0. ≥ m(x)-\|v_k^*^*-v_j^*\|_2\,\|δ\|_2>0. Hence g∗(x+δ)=k∗g^*(x+δ)=k^*. ∎ We now strengthen this deterministic result to a probabilistic stability bound under Gaussian perturbations. To make the dependence on the number of regimes K explicit, we additionally assume bounded gating energy. Assumption IV.3 (Bounded Gating Energy). ‖vk∗‖2≤Bv\|v_k^*\|_2≤ B_v for all k∈1,…,Kk∈\1,…,K\. Fix x and let k∗=g∗(x)k^*=g^*(x). For each j≠k∗j≠ k^*, define the gating-vector difference aj≜vj∗−vk∗a_j v_j^*-v_k^*^*. Define the Gaussian-max statistic (random in δ) as Mk∗(δ)≜maxj≠k∗aj⊤δM_k^*(δ) _j≠ k^*a_j δ. A misassignment occurs only if some competitor’s perturbed gain exceeds the margin: g∗(x+δ)≠k∗⊆Mk∗(δ)≥m(x).\g^*(x+δ)≠ k^*\ \M_k^*(δ)≥ m(x) \. The following result bounds this misassignment probability via Gaussian-max concentration [3]. Proposition IV.2 (Probabilistic Stability Bound via Gaussian-Max Concentration). Under Assumptions IV.1–IV.3, the misassignment probability (with respect to δ) satisfies ℙδ(g∗(x+δ)≠k∗∣x)≤exp(−(m(x)−δMk∗(δ))+22σg2maxj≠k∗‖aj‖22),P_δ\! (g^*(x+δ)≠ k^* x )≤ \! (- (m(x)-E_δM_k^*(δ) )_+^22\, _g^2 _j≠ k^*\|a_j\|_2^2 ), (13) where (u)+≜maxu,0(u)_+ \u,0\ and δE_δ denotes expectation with respect to δ∼(0,σg2Ip)δ (0, _g^2I_p). Moreover, using δMk∗(δ)≤σgmaxj≠k∗‖aj‖22log(K−1),E_δM_k^*(δ)≤ _g _j≠ k^*\|a_j\|_2 2 (K-1), and ‖aj‖2≤2Bv\|a_j\|_2≤ 2B_v, we obtain the explicit bound ℙδ(g∗(x+δ)≠k∗∣x)≤exp(−12(m(x)2σgBv−2log(K−1))+2).P_δ\! (g^*(x+δ)≠ k^* x )\\ ≤ \! (- 12 ( m(x)2 _gB_v- 2 (K\!-\!1) )_+^2 ). (14) For any prescribed margin threshold mmin>0m_ >0, Eq. (14) holds uniformly for all x∈x satisfying m(x)≥mminm(x)≥ m_ , with m(x)m(x) replaced by mminm_ , since the bound is decreasing in m(x)m(x). Proof. The proof is given in the supplementary material. ∎ Three factors compete in the explicit bound in Eq. (14): a larger margin m(x)m(x) promotes stability, a higher noise level σg _g erodes it, and the regime-complexity term 2log(K−1) 2 (K-1) imposes an additional cost that grows only slowly with K. Remark IV.1 (When the bound is nonvacuous). The explicit bound in Proposition IV.2 yields a non-trivial guarantee (i.e., strictly less than 11) only when m(x)>2σgBv2log(K−1).m(x)>2 _gB_v 2 (K-1). This condition requires the margin to exceed an upper bound on the expected maximum Gaussian fluctuation. Otherwise, Eq. (14) becomes vacuous. IV-C Joint Regime and Support Recovery with Learning Error We study support recovery with a learned regime selector and oracle-trained experts. This isolates regime uncertainty and does not account for the effect of assignment errors on expert training. All probabilities in this support-recovery analysis are conditional on the training inputs. Let nD_n denote the training data, independent of the test perturbation δ. Let g g be a learned regime selector constructed from nD_n, and define k^=g^(x+δ) k= g(x+δ) at test time. For each regime k, define the oracle index set ℐk∗=i:g∗(xi)=k,nk=|ℐk∗|.I_k^*=\i:g^*(x_i)=k\, n_k=|I_k^*|. Let β^k β_k be the LASSO estimator fitted using the samples indexed by ℐk∗I_k^*, and let S^k=supp(β^k) S_k=supp( β_k). These oracle index sets are unavailable to the practical algorithm because g∗g^* is unknown. We assume that observations in regime k follow y=x⊤βk∗+ε,ε∼(0,σ2),y=x _k^*+ , (0,σ^2), independently across samples, with σ>0σ>0. To bring classical LASSO support-recovery guarantees to bear within each regime, we impose some standard conditions on the within-regime design and signal: a restricted eigenvalue condition (Assumption 0.1), an irrepresentable condition with incoherence parameter γk∈(0,1] _k∈(0,1] (Assumption 0.2), a beta-min condition (Assumption 0.3), and the penalty choice λ=Aλσlogp/nk∗λ=A_λσ p/n_k^* with Aλ≥4/γk∗A_λ≥ 4/ _k^*, where γk∗ _k^* is the incoherence parameter of the true regime k∗k^*. Under these conditions, the LASSO result for given training inputs [1] guarantees support recovery within regime k∗k^* with probability at least 1−c1exp(−c2nk∗λ2/σ2)1-c_1 (-c_2n_k^*λ^2/σ^2) for some constants c1,c2>0c_1,c_2>0. We state the corresponding result as Theorem 0.1 in the supplementary material. Since the regime selector is itself estimated, misassignment of the regime is a third source of error, alongside LASSO failure within the correct regime and the perturbation-driven switching of Section IV-B. We assume that the learned selector g g satisfies supx∈ℙn(g^(x)≠g∗(x))≤ηn _x P_D_n ( g(x)≠ g^*(x) )≤ _n for a sequence ηn→0 _n→ 0 (Assumption 0.4). All of these conditions are stated formally in Section IV of the supplementary material. Combining the three sources yields an end-to-end guarantee for oracle-trained experts. For fixed x as in Section IV-B, under Assumptions IV.1–IV.3 and 0.1–0.4 and the stated penalty choice, the probability that the support recovered in the predicted regime differs from the true support satisfies ℙ(S^k^≠Sk∗) ( S_ k≠ S_k^* ) ≤ηn⏟regime learning+c1exp(−c2nk∗λ2/σ2)⏟support recovery ≤ _n_regime learning+ c_1 \! (-c_2n_k^*λ^2/σ^2 )_support recovery (15) +exp(−12(m(x)2σgBv−2log(K−1))+2)⏟perturbation stability. + \! (- 12 ( m(x)2 _gB_v- 2 (K-1) )_+^2 )_perturbation stability. The formal statement and proof are given in the supplementary material as Proposition 0.1. Section IV-D explains the implications of this bound for the choice of K. IV-D Design Implications and Structural Trade-Offs Regime Learning Accuracy. The term ηn _n captures the generalization error of the learned regime selector g g. Assumption 0.4 abstracts the regime-identification performance of g g through this single quantity. Since a full analysis of the regime-learning procedure lies beyond our scope, we instead take ηn _n as given and combine it additively with the switching probability and the within-regime recovery error. This yields an end-to-end reliability bound. A smaller ηn _n tightens the overall support-recovery bound. Perturbation Stability. The perturbation-driven instability term is controlled by the margin m(x)m(x) relative to the noise-complexity radius 2σgBv2log(K−1)2 _gB_v 2 (K-1). This exposes a trade-off: increasing the number of regimes K raises approximation capacity, but it also enlarges the noise-complexity radius at rate logK K, thereby weakening the stability guarantee. This comparison fixes m(x)m(x), σg _g, and BvB_v. In REALM, changing K may also alter margins and improve approximation, but these effects are not captured by the fixed-regime analysis. Per-Regime Sample Complexity. Recall from Section IV-C that λ≍σlogp/nk∗λ σ p/n_k^*. With balanced regimes (nk∗≈n/Kn_k^*≈ n/K) and all else fixed, a larger K tightens the beta-min requirement: minj∈Sk∗|βk∗,j∗|≳σKlogpn. _j∈ S_k^*| _k^*,j^*| σ K pn. Consequently, a larger K raises the signal strength required by the support-recovery guarantee. Complexity–Stability Trade-Off. Overall, these effects make K a single knob that trades representational capacity against statistical reliability. A larger K allows a finer partition and can reduce approximation bias when the regression function is smooth, but it also shrinks the per-regime sample size and can increase switching instability, potentially increasing estimation variability and making the support-recovery conditions harder to satisfy. Under balanced sampling, we use the sample-size heuristic nK≫smaxlogp, nK s_ p, (16) where smax=max1≤k≤K|Sk|s_ = _1≤ k≤ K|S_k|. Together with reliable regime learning and adequate margins, this heuristic motivates a moderate number of well-separated regimes to balance approximation accuracy and explanation reliability. We examine this qualitative trade-off in the ablation on K in Section V-C. V Experiments We evaluate REALM on six datasets, comparing it with SEE-Net [28], MLM [27], other inherently interpretable models, and standard MLP/CNN models. MLM and the inherently interpretable models are evaluated only on tabular datasets. SEE-Net (Pre-intp) and SEE-Net (Intp) denote SEE-Net’s pre-interpretable N, which predicts through the guidance N, and its final interpretable N, respectively. MLM-CELL is the cell-level mixture obtained from layer-wise clustering, whereas MLM-EPIC merges cells into larger Explainable Prediction-induced Input Clusters (EPIC). SEE-Net (Intp) and MLM-EPIC are the interpretable models intended for deployment in practice. The six datasets comprise three tabular datasets (Bike Sharing, California Housing, Covertype) and three image datasets (MNIST, Flowers, CIFAR-10). We detail the full dataset descriptions, implementation and training details, and additional per-dataset interpretability analyses in the supplementary material. In the main experiments, REALM routes inputs using the WSSN-induced posterior πk(x) _k(x) rather than the explanatory gate, so that the experiments examine the usefulness of the WSSN activation patterns themselves. Section V-C examines REALM-EG, the derived two-layer linear model in which the explanatory gate takes over routing. V-A Synthetic Data Study To illustrate the regime-recovery behavior of REALM, we study a synthetic dataset with known ground-truth regimes. We draw n=2,000n=2,000 points uniformly over the square [−1,1]2[-1,1]^2 with independent covariates X1,X2X_1,X_2, and let the response Y follow a piecewise-linear trapezoidal-bowl surface with an affine interior and additive Gaussian noise. Let r=max(|X1|,|X2|)r= (|X_1|,|X_2|). Then Y=β0+β1X1+β2X2+εif r≤t,X2−t+εif r>t,X2≥|X1|,−X2−t+εif r>t,X2≤−|X1|,X1−t+εif r>t,X1>|X2|,−X1−t+εif r>t,X1<−|X2|,Y= cases _0+ _1X_1+ _2X_2+ &if r≤ t,\\ X_2-t+ &if r>t,\ X_2≥|X_1|,\\ -X_2-t+ &if r>t,\ X_2≤-|X_1|,\\ X_1-t+ &if r>t,\ X_1>|X_2|,\\ -X_1-t+ &if r>t,\ X_1<-|X_2|, cases where t=0.5t=0.5, ε∼(0,0.022) (0,0.02^2) is the Gaussian noise, and the center regime uses β0=0.5 _0=0.5, β1=0.5 _1=0.5, β2=−0.5 _2=-0.5. For this experiment we fix the final number of regimes at K=5K=5 for REALM, MLM, and SEE-Net. Since each regime is linear in the inputs, the optimal decomposition follows the edges of the trapezoidal bowl. Regime boundaries are marked by jumps or changes in slope. Fig. 2 compares MLM, SEE-Net, and REALM on the same synthetic draw. MLM in Fig. 2 misses almost the entire right and middle regions, and SEE-Net in Fig. 2 cuts across the true boundaries and recovers the regimes only partially, whereas REALM in Fig. 2 aligns most closely with the true geometry. (a) MLM (b) SEE-Net (c) REALM Fig. 2: Predicted regime assignments for MLM, SEE-Net, and REALM. Ground-truth boundaries in red, predicted boundaries as black outlines. V-B Real-World Data Study Prediction Performance. Tables I and I report the predictive performance of REALM against inherently interpretable models, SEE-Net, MLM, and MLP/CNN models, where applicable. We also report Gate AUC, which measures how accurately the explanatory gate predicts the WSSN-induced regime assignments. Table I shows that REALM is the strongest interpretable model on all tabular datasets, outperforming SEE-Net (Intp), MLM-CELL, MLM-EPIC, and the inherently interpretable models. The MLM baselines exhibit substantial variation on Bike Sharing, suggesting instability under our experimental setup. Strong black-box models such as MLP and XGBoost still perform better than REALM, but REALM remains simple to read because the routing posterior is near one-hot, so each prediction can be explained mostly by a single linear expert, as Fig. 10– in the supplementary material indicates. In Table I, REALM nearly matches the teacher CNN on MNIST. On Flowers, the distilled WSSN achieves a higher AUC than the teacher, illustrating a potential predictive benefit of the distillation stage. REALM retains part of this gain while providing regime-level explanations. Overall, REALM improves over SEE-Net (Intp) on all three image datasets. This shows that REALM’s performance advantage carries over from tabular features to high-dimensional images. More importantly, these image datasets demonstrate that the interpretability benefits come at little cost in accuracy, although a gap to the strongest black-box models remains on the tabular datasets. REALM is also stable across runs since the standard deviations are small and only larger on harder datasets such as California Housing and Flowers. This consistency indicates that the distillation, clustering, and fixed-router fitting pipeline yields robust final predictions. The gate-recovery results in both Tables I and I provide a complementary interpretability view. The WSSN-induced regime labels are recovered almost perfectly by the linear gate on Bike Sharing, Covertype, MNIST, and CIFAR-10, suggesting that the switching logic is both effective and explainable. The lower Gate AUC on California Housing reflects the fact that a larger, more fragmented partition is harder to explain with a simple linear assignment rule, as this dataset uses the largest partition, K=150K=150, among all datasets. Among the image datasets, Flowers has the lowest Gate AUC, suggesting that its regime assignments are harder to recover with a linear gate. Nevertheless, its Gate AUC of 0.974 remains high, indicating that the assignments are still largely recoverable from the learned representation. Table I: Prediction performance on the three tabular datasets, reported as mean (std) over five runs. Gate AUC measures how accurately the linear gate recovers the discovered regime assignments. We fit the same explanatory gate to the regimes of SEE-Net and the MLM variants, so the column is comparable across all regime-switching methods. All AUC values are macro-averaged one-vs-rest over classes or regimes. Empty cell means not applicable. Bold black text marks the best result for performance column (excluding SEE-Net (Pre-intp), WSSN, and Gate AUC); bold dark blue, the best interpretable model. Bike Sharing Cal. Housing Covertype Method Test RMSE ↓ Test Gate AUC ↑ Test RMSE ↓ Test Gate AUC ↑ Test AUC ↑ Test Gate AUC ↑ Logistic Reg. 0.937 (0.000) SVM 0.925 (0.000) LASSO Reg. 203.485 (0.000) 0.830 (0.008) GAM 241.067 (0.000) 1.050 (0.009) 0.943 (0.000) Random Forest 94.510 (0.456) 0.536 (0.007) 0.956 (0.001) SVR 127.276 (0.000) 0.568 (0.006) XGBoost 89.355 (0.236) 0.461 (0.007) 0.986 (0.000) SEE-Net (Pre-intp) 269.171 (0.215) 0.995 (0.005) 0.663 (0.008) 0.943 (0.021) 0.988 (0.002) 0.948 (0.013) SEE-Net (Intp) 262.994 (2.442) 0.997 (0.003) 3.893 (3.934) 1.000 (0.000) 0.924 (0.008) 1.000 (0.000) MLM-CELL 875.494 (379.166) 0.943 (0.013) 0.651 (0.022) 0.880 (0.010) 0.816 (0.014) 0.994 (0.002) MLM-EPIC 823.153 (369.246) 0.739 (0.092) 0.682 (0.130) 0.879 (0.014) Teacher MLP 89.498 (3.812) 0.515 (0.006) 0.995 (0.000) WSSN 90.367 (5.497) 0.525 (0.006) 0.991 (0.000) REALM 142.273 (4.241) 0.988 (0.008) 0.608 (0.022) 0.738 (0.113) 0.963 (0.001) 1.000 (0.000) Table I: Prediction performance on the three image datasets, reported as mean (std) over five runs. Inherently interpretable models and MLM-CELL/MLM-EPIC are omitted, as they are not designed for image data. Conventions follow Table I. MNIST Flowers CIFAR-10 Method Test AUC ↑ Test Gate AUC ↑ Test AUC ↑ Test Gate AUC ↑ Test AUC ↑ Test Gate AUC ↑ SEE-Net (Pre-intp) 0.999 (0.000) 0.979 (0.008) 0.875 (0.007) 0.804 (0.137) 0.920 (0.005) 0.779 (0.015) SEE-Net (Intp) 0.985 (0.002) 1.000 (0.000) 0.841 (0.053) 0.999 (0.001) 0.773 (0.029) 0.994 (0.013) Teacher CNN 1.000 (0.000) 0.850 (0.018) 0.961 (0.001) WSSN 1.000 (0.000) 0.893 (0.015) 0.927 (0.001) REALM 0.999 (0.000) 0.999 (0.000) 0.862 (0.019) 0.974 (0.007) 0.908 (0.002) 0.993 (0.000) Regime Recovery and Interpretability Analysis. Beyond prediction, REALM provides interpretable regime assignments and regime-specific prediction rules, explaining both why an instance is routed to a regime and how its prediction is formed. We examine this dual interpretability on the image datasets from three perspectives: 1) whether routing is effective without collapsing to a single regime; 2) what the explanatory gate and local expert coefficients reveal about regime assignment and prediction; 3) whether the resulting regimes and explanations are human-readable. The corresponding diagnostics for the tabular datasets are reported in the supplementary material. (a) MNIST (b) CIFAR-10 (c) Flowers Fig. 3: Routing sharpness and regime size balance on MNIST (a), CIFAR-10 (b), and Flowers (c). Top row: distribution of the maximum routing posterior maxkπk(x) _k _k(x) over training samples. Bottom row: distribution of regime sizes, i.e., the number of training samples assigned to each non-empty regime. Bars are means over five independent runs, and error bars are ±1.96± 1.96 standard errors across runs. Fig. 4: Single-run explanatory gate coefficients on MNIST (a), CIFAR-10 (b), and Flowers (c). Rows are the 10 largest regimes by training-set size. Columns show feature units of u(x)u(x), labeled g_j. For MNIST, a feature unit is one spatial location in a convolutional channel; for CIFAR-10 and Flowers, it is a max-pooled channel. Color represents the gate coefficient of the labeled feature unit for that regime. Fig. 3 shows that routing is effective on all three datasets. On MNIST, the maximum routing posterior concentrates near one and the regime sizes are balanced across all the regimes (Fig. 3), indicating sharp, noncollapsed routing. On CIFAR-10, routing is less sharp than on MNIST, which is consistent with the greater visual complexity, while the regime sizes remain balanced with few oversized or undersized regimes (Fig. 3). On Flowers, routing is likewise less sharp than on MNIST, while the regime sizes stay relatively balanced (Fig. 3). Fig. 4 illustrates regime-assignment rules for the image datasets. For MNIST, g_613 has a coefficient of +0.67+0.67 for regime 1 against −0.70-0.70 and −0.68-0.68 for regimes 2 and 3, so stronger activation of this feature-map location favors regime 1 over regimes 2 and 3. For CIFAR-10, g_31 has a coefficient of +0.73+0.73 for regime 1 against −0.40-0.40 for regime 14, so stronger activation of this channel favors regime 1 over regime 14. For Flowers, some pooled features receive coefficients of opposite signs across regimes, reflecting different regime-assignment rules. (a) Regime 1 (b) Regime 2 (c) Regime 6 (d) Regime 7 (e) Regime 8 (f) Regime 0 (g) Regime 3 (h) Regime 6 (i) Regime 7 (j) Regime 9 (k) Regime 0 (l) Regime 1 (m) Regime 6 (n) Regime 7 (o) Regime 9 Fig. 5: Training images grouped by regime. (a)–(e) MNIST digits 1 and 3 in five selected regimes, revealing variations in writing style among same-digit images. (f)–(j) CIFAR-10 dogs and trucks in five selected regimes, revealing variations in color, shape, and style among same-class images. (k)–(o) Dandelions and sunflowers from Flowers in five selected regimes, revealing variations in visual tone, framing, and bloom stage among same-class images. Within each regime, each class is shown as a 3×33× 3 grid of randomly sampled training images. Fig. 5 shows that the learned regimes are human-readable. For MNIST in Fig. 5–, each regime isolates a coherent writing style. For example, digit 1 has a more formal style in regime 1, a strongly slanted version in regime 2, and thin upright strokes in regime 8. Digit 3 also shows distinct variations across regimes. For CIFAR-10 in Fig. 5–, the partition captures interpretable structure even at 32×3232× 32 resolution. Dogs separate largely by pose and framing: regime 0 collects head-only close-ups facing the camera; regime 7 gathers side-facing views; regimes 3 and 6 contain full-body shots. Trucks are instead organized by viewpoint, with regimes 3, 6, and 7 each grouping trucks facing a consistent direction. Regime 9 is color-driven across both classes, with predominantly white dogs on blue backgrounds and blue-toned trucks. For Flowers in Fig. 5–, regimes 6 and 9 share a warm visual tone and differ mainly in dandelions, with white seed heads in regime 6 and yellow blooms in regime 9; regime 1 gathers sunflowers on lighter backgrounds, and regime 0 leans toward distant, wide-field shots. The visual differences between the displayed regimes are less pronounced for Flowers than for MNIST and CIFAR-10. Overall, REALM organizes images by writing style, pose, framing, and color, suggesting that it partitions the input space along human-readable features. Fig. 6: CIFAR-10 local expert contribution maps. The first column shows the original image and its index (“idx”). The next four columns show the filter groups with the largest total absolute weight×activation. Each group sums the expert-weighted pre-pooling feature maps across the three RGB input channels. The final two columns sum contributions over these four groups (Top-4) and all channels, respectively. All maps are computed for the ground-truth class using the expert with the highest routing posterior. Larger positive contributions appear darker in the per-group maps and brighter in the aggregate maps. After analyzing regime discovery, we now turn to regime-specific prediction logic. Fig. 6 visualizes the feature-map locations that are most influential for each prediction. In row 1, f60 and f12 trace the car’s silhouette while f64 focuses on the windshield. The Top-4 sum concentrates on the silhouette. In row 2, f101 outlines the silhouette, f85 responds to the upper body and f64 to the lower body, while f51 instead focuses on the background. The Top-4 sum emphasizes the lower body together with the full silhouette. Row 5 is the most consistent example in which every displayed filter group, together with the Top-4 and All-channel sums, concentrates on the shadowed part of the car. For the frog in row 6, f74, f1, and f44 spread over the whole animal while f5 focuses on the background, and the Top-4 sum covers the body and traces its outline. Notably, several filter groups recur across images: f64 on both cars (rows 1–2), f44 and f5 on both frogs (rows 3 and 6), and f60 on a car (row 1) and the deer (row 4), which suggests that these filter groups capture reusable visual patterns. Across rows, the Top-4 and All-channel sums remain similar, indicating that predictive contributions are concentrated in a few dominant filter groups. V-C Ablation Studies Does the Router Carry the Predictive Burden? A key design claim of REALM is that the routing stage discovers a meaningful partition while the predictive burden is carried mostly by the local linear experts rather than by the router itself. To verify this claim, we replace the local linear experts with a majority-vote baseline that fits no model within each regime. For classification, the majority vote replaces each local expert prediction ek(x)e_k(x) in Eq. (6) with the empirical class-frequency vector dkd_k of regime k. Each component dk,cd_k,c is the proportion of samples in class c among all training samples in regime k. The final prediction aggregates them as in Eq. (6), giving the weighted vote f^mv(x)=argmaxc∑k=1Kπk(x)dk,c f_mv(x)= argmax_c _k=1^K _k(x)\,d_k,c. For AUC evaluation, we use the mixed class-frequency vector ∑k=1Kπk(x)dk _k=1^K _k(x)d_k before taking the argmax. For regression, we use y¯k y_k, the mean response for training samples in regime k, so the prediction becomes f^mv(x)=∑k=1Kπk(x)y¯k f_mv(x)= _k=1^K _k(x)\, y_k. Since linear modeling is skipped in all regimes, this experiment reveals the standalone prediction capability of the routing stage. Fig. 7: Effect of replacing the local linear experts with majority-vote predictions. Blue bars give the performance degradation (percentage increase in RMSE for regression, or percentage drop in AUC for classification). Orange bars show the increase in the cross-run standard deviation of the test metric, expressed as a percentage of the full model’s test metric. Negative orange bars indicate a decrease in standard deviation. Positive values mean the ablation is worse than the full REALM. Fig. 7 summarizes the results across all six datasets. Replacing the local experts with majority-vote predictions reduces predictive performance on every dataset and increases the cross-run standard deviation of the test metric on four of the six datasets, most noticeably on MNIST and Bike Sharing. This indicates that the discovered regimes alone are insufficient for strong predictive performance. These results show that the linear experts recover predictive performance that within-regime constants cannot, suggesting that the learned regimes successfully decompose a nonlinear prediction task into subproblems that can be effectively modeled by linear models. Overall, these findings indicate that REALM identifies a meaningful two-level structure that achieves both interpretability and strong predictive accuracy. Routing with the Explanatory Gate. Section I introduced REALM-EG, the two-layer linear model in which the explanatory gate is used as the router instead of as a post hoc explanation of the router. Since the gate already recovers the discovered regimes with high Gate AUC, as shown in Tables I and I, we examine whether REALM-EG predicts as well as REALM. Fig. 8: Effect of replacing the router with the explanatory gate and refitting the local linear experts on all training samples under the gate posterior q^k(x) q_k(x). Bars and sign convention as in Fig. 7. As shown in Fig. 8, REALM-EG performs comparably to REALM. The largest performance degradation is 1.6% on Flowers, and REALM-EG slightly improves on Bike Sharing and California Housing. The differences on Covertype, MNIST, and CIFAR-10 are negligible. These results indicate that REALM-EG achieves predictive performance comparable to REALM, demonstrating that a two-layer linear model can achieve strong predictive performance while remaining inherently interpretable. We retain WSSN-based routing in the main results to evaluate the predictive usefulness of activation-induced regimes directly. Effect of the Number of Regimes K. Section IV-D motivates a trade-off between representational capacity and statistical reliability as K varies. We therefore examine how REALM performs with varying K on the shared grid K∈5, 10, 50, 100, 150, 200K∈\5,\,10,\,50,\,100,\,150,\,200\ for each dataset. We additionally include results for CIFAR-10 and Covertype at K=20K=20 and K=30K=30, respectively, which are the values of K used for these datasets in the main experiments. Fig. 9 illustrates the results for each dataset, with regression and classification tasks separated into two plots. Eq. (16) provides a viewpoint for interpreting Fig. 9. Since the fitted experts retain nearly all features (Table IV of the supplementary material), we take smax=ps_ =p in Eq. (16) and report the ratio n/(Kplogp)n/(K\,p p) as a heuristic diagnostic of sample size relative to model complexity as K varies. Larger values suggest a more favorable balance between the available samples and the complexity of the local models. For Covertype and California Housing, the ratio remains above 1.01.0 over the entire grid, decreasing from 345.3345.3 to 8.68.6 on Covertype (n=371,847n=371,847, p=54p=54) and from 158.8158.8 to 4.04.0 on California Housing (n=13,209n=13,209, p=8p=8). Performance improves monotonically on both datasets. As K=200K=200 yields little gain over K=150K=150 on California Housing, we chose K=150K=150 for better interpretability, and likewise K=30K=30 for Covertype. Bike Sharing (n=11,122n=11,122, p=19p=19) shows that this ratio alone does not identify the best K. Its RMSE is minimized at K=10K=10 (mean 134.707, std 14.026), where the ratio is still 19.919.9, and increases thereafter. The ratio decreases to approximately 1.01.0 at K=200K=200. Nevertheless, Bike Sharing clearly shows a U-shaped curve, with the best performance at an intermediate K, consistent with the trade-off between representational capacity and statistical reliability. For the image datasets, the representation is high-dimensional, and the ratio is below 1.01.0 over most of the grid: the ratio drops below 1.01.0 at K=20K=20 on CIFAR-10 (n=45,000n=45,000, p=384p=384), and at K=5K=5 it is already 0.540.54 on MNIST (n=48,000n=48,000, p=2,304p=2,304) and 0.220.22 on Flowers (n=2,569n=2,569, p=384p=384). Fig. 9 shows that CIFAR-10’s AUC decreases steadily, with a sharper decline beyond K=50K=50, while Flowers shows the largest AUC decrease among the three image datasets. MNIST maintains a near-perfect AUC, illustrating that this ratio alone does not determine predictive performance. Except for Covertype and California Housing, K was matched to the applicable baseline settings. Fig. 9: Effect of the number of regimes K on test performance. Left: regression RMSE (Bike Sharing on the left axis, California Housing on the right). Right: classification AUC. VI Conclusion We introduced REALM, regime-switching, explainable, and activation-induced linear models that derive latent regimes from the piecewise-affine activation geometry of a distilled wide, shallow ReLU student network (WSSN) and fit a mixture of local linear experts in the discovered regimes. REALM provides dual interpretability: 1) regime-specific feature attributions explaining each prediction; 2) an explanatory gate explaining the regime-switching logic. By analyzing an idealized setting in which the true regimes are known and the model is exactly linear within each, we characterized how regime assignment responds to input perturbations and when within-regime supports are recoverable, revealing a complexity–stability trade-off in the number of regimes K, consistent with the performance curves in the ablation study in Section V-C. Experiments on tabular and image datasets showed that REALM attains competitive predictive accuracy relative to prior DNN-guided mixture surrogates and inherently interpretable models, with local explanations illustrated qualitatively through gate coefficient heat maps and regime-level analysis. Ablations further showed that the router does not carry the predictive burden and that the explanatory gate can itself serve as the routing signal with little loss in accuracy. References [1] C. Rudin, “Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead,” Nature Machine Intelligence, vol. 1, no. 5, p. 206–215, May 2019. [2] F. Mumuni and A. Mumuni, “Explainable artificial intelligence (XAI): from inherent explainability to large language models,” arXiv:2501.09967, 2025. [3] R. Caruana, Y. Lou, J. Gehrke, P. Koch, M. Sturm, and N. Elhadad, “Intelligible models for HealthCare: Predicting pneumonia risk and hospital 30-day readmission,” in ACM SIGKDD, Aug. 2015, p. 1721–1730. [4] E. Tabassi, “Artificial intelligence risk management framework (AI RMF 1.0),” National Institute of Standards and Technology, Gaithersburg, MD, USA, NIST Trustworthy and Responsible AI NIST AI 100-1, Jan. 2023. [5] N. Frosst and G. Hinton, “Distilling a neural network into a soft decision tree,” in Workshop on Comprehensibility and Explanation in AI and ML (CEX), 2017. [6] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton, “Adaptive mixtures of local experts,” Neural Computation, vol. 3, no. 1, p. 79–87, Feb. 1991. [7] M. I. Jordan and R. A. Jacobs, “Hierarchical mixtures of experts and the EM algorithm,” Neural Computation, vol. 6, no. 2, p. 181–214, Mar. 1994. [8] B. Tang, M. Heywood, and M. Shepherd, “Input partitioning to mixture of experts,” in International Joint Conference on Neural Networks, 2002, p. 227–232. [9] W. Guo, S. Huang, Y. Tao, X. Xing, and L. Lin, “Explaining deep learning models—a Bayesian non-parametric approach,” in Advances in Neural Information Processing Systems, vol. 31, 2018. [10] S. E. Yuksel, J. N. Wilson, and P. D. Gader, “Twenty years of mixture of experts,” IEEE Transactions on Neural Networks and Learning Systems, vol. 23, no. 8, p. 1177–1193, Aug. 2012. [11] 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. [12] B. Badjie, J. Cecílio, and A. Casimiro, “Decoupling mixture-of-experts routing from gradient noise: A framework for structured specialization and soft generalization toward robust and efficient inference,” Expert Systems with Applications, vol. 309, p. 131128, May 2026. [13] B. Schwab and J. C. Schneider, “Interpretable mixture-of-experts via soft decision tree routing and additive networks,” SSRN preprint 5810123, Nov. 2025. [14] B. Zoph, I. Bello, S. Kumar, N. Du, Y. Huang, J. Dean, N. Shazeer, and W. Fedus, “ST-MoE: Designing stable and transferable sparse expert models,” arXiv:2202.08906, 2022. [15] M. Vasić, A. Petrović, K. Wang, M. Nikolić, R. Singh, and S. Khurshid, “MoËT: Mixture of expert trees and its application to verifiable reinforcement learning,” Neural Networks, vol. 151, p. 34–47, Jul. 2022. [16] A. A. Ismail, S. Ö. Arik, J. Yoon, A. Taly, S. Feizi, and T. Pfister, “Interpretable mixture of experts,” Transactions on Machine Learning Research, 2023. [17] S. Ghosh, K. Yu, F. Arabshahi, and K. Batmanghelich, “Dividing and conquering a BlackBox to a mixture of interpretable models: Route, interpret, repeat,” in International Conference on Machine Learning, 2023, p. 11 360–11 397. [18] W.-L. Hwang and A. Heinecke, “Un-rectifying non-linear networks for signal representation,” IEEE Transactions on Signal Processing, vol. 68, p. 196–210, 2020. [19] N. Picchiotti and M. Gori, “Clustering-based interpretation of deep ReLU network,” in AIxIA 2021 – Advances in Artificial Intelligence, ser. Lecture Notes in Computer Science, vol. 13196. Springer International Publishing, 2022, p. 403–412. [20] J. He, L. Li, J. Xu, and C. Zheng, “ReLU deep neural networks and linear finite elements,” Journal of Computational Mathematics, vol. 38, no. 3, p. 502–527, Jun. 2020. [21] A. Sudjianto, W. Knauth, R. Singh, Z. Yang, and A. Zhang, “Unwrapping the black box of deep ReLU networks: interpretability, diagnostics, and simplification,” arXiv:2011.04041, 2020. [22] G. Montúfar, R. Pascanu, K. Cho, and Y. Bengio, “On the number of linear regions of deep neural networks,” in Advances in Neural Information Processing Systems, vol. 27, 2014. [23] M. Raghu, B. Poole, J. Kleinberg, S. Ganguli, and J. Sohl-Dickstein, “On the expressive power of deep neural networks,” in International Conference on Machine Learning, 2017, p. 2847–2854. [24] S. R. Dubey, S. K. Singh, and B. B. Chaudhuri, “Activation functions in deep learning: A comprehensive survey and benchmark,” Neurocomputing, vol. 503, p. 92–108, Sep. 2022. [25] A. Barua, M. U. Ahmed, and S. Begum, “Mechanistic interpretability of ReLU neural networks through piecewise-affine mapping,” Machine Learning, vol. 115, p. 17, Jan. 2026. [26] S. Anwar, N. Griffiths, A. Bhalerao, and T. Popham, “MASALA: Model-agnostic surrogate explanations by locality adaptation,” in KDD Workshop on Human-Interpretable AI, 2024. [27] B. Seo, L. Lin, and J. Li, “Mixture of linear models co-supervised by deep neural networks,” Journal of Computational and Graphical Statistics, vol. 31, no. 4, p. 1303–1317, Sep. 2022. [28] B. Seo and J. Li, “Explainable machine learning by SEE-Net: closing the gap between interpretable models and DNNs,” Scientific Reports, vol. 14, no. 1, p. 26302, Nov. 2024. [29] G. Hinton, O. Vinyals, and J. Dean, “Distilling the knowledge in a neural network,” arXiv:1503.02531, 2015. [30] C. Buciluă, R. Caruana, and A. Niculescu-Mizil, “Model compression,” in ACM SIGKDD, Aug. 2006, p. 535–541. [31] L. J. Ba and R. Caruana, “Do deep nets really need to be deep?” in Advances in Neural Information Processing Systems, vol. 27, 2014. [32] T. Serra, C. Tjandraatmadja, and S. Ramalingam, “Bounding and counting linear regions of deep neural networks,” in International Conference on Machine Learning, 2018, p. 4558–4566. [33] C. C. Aggarwal, A. Hinneburg, and D. A. Keim, “On the surprising behavior of distance metrics in high dimensional space,” in Database Theory – ICDT 2001, ser. Lecture Notes in Computer Science, vol. 1973. Springer Berlin Heidelberg, 2001, p. 420–434. [34] M. Verleysen and D. François, “The curse of dimensionality in data mining and time series prediction,” in Computational Intelligence and Bioinspired Systems (IWANN 2005), ser. Lecture Notes in Computer Science, vol. 3512. Springer Berlin Heidelberg, 2005, p. 758–770. [35] R. Vershynin, High-Dimensional Probability: An Introduction with Applications in Data Science. Cambridge University Press, Sep. 2018. [36] M. J. Wainwright, “Sharp thresholds for high-dimensional and noisy sparsity recovery using ℓ1 _1-constrained quadratic programming (LASSO),” IEEE Transactions on Information Theory, vol. 55, no. 5, p. 2183–2202, May 2009. Supplementary Material for REALM: Regime-Switching, Explainable, and Activation-Induced Linear Models Xiaoran Cheng, Sen Na, and Jia Li, I Datasets Bike Sharing. This dataset records hourly bike rentals in the Capital Bikeshare system during 2011 and 2012. The target is the total number of rented bikes. The inputs are temporal and weather covariates, including year, month, hour, holiday indicator, weekday, working-day indicator, temperature, apparent temperature, humidity, wind speed, season, and weather condition. Hour, weekday, and month are encoded as sine and cosine pairs; season and weather condition are dummy encoded with the first level dropped, yielding 1919 input features. California Housing. Derived from the 1990 U.S. Census, this dataset contains census block groups in California, with median house value as the target. The 8 covariates are median income, house age, average numbers of rooms and bedrooms, population, average occupancy, latitude, and longitude. Covertype. This dataset is used to classify forest cover types in the Roosevelt National Forest of northern Colorado. Each instance represents a 30m×30m30\,m× 30\,m patch of forested land. The inputs comprise 10 continuous variables (elevation, aspect, slope, horizontal and vertical distances to hydrology, horizontal distance to roadways, hillshade at 9 a.m., noon, and 3 p.m., and horizontal distance to fire points), 4 binary wilderness-area indicators, and 40 binary soil-type indicators, yielding 5454 features. MNIST. This handwritten digit classification dataset contains 28×2828× 28 grayscale images of digits 00 through 99. Flowers. This dataset contains 180×180180× 180 color images (after resizing) from five flower categories: daisy, dandelion, roses, sunflowers, and tulips. CIFAR-10. This dataset contains 32×3232× 32 color images from ten object classes: airplane, automobile, bird, cat, deer, dog, frog, horse, ship, and truck. I Implementation Details Across all datasets, we first train a DNN teacher and then follow the three-stage pipeline in Section I of the main paper: 1) we distill the teacher into the WSSN with distillation weight α=0.7α=0.7 and, for classification, temperature T=2T=2 in Eq. (2); 2) we freeze the WSSN, fixing both s(x)s(x) and u(x)u(x), and optimize only the K centroids to discover regimes; 3) with the routing posterior πk(x) _k(x) fixed, we fit the local experts and explanatory gate using Eqs. (7) and (9), respectively. Table I lists the dataset-specific hyperparameters. Table I: Dataset-specific hyperparameters. Shared hyperparameters are given in the text. The minimum regime size τ is used only for tabular data augmentation. Bike Sharing Cal. Housing Covertype MNIST Flowers CIFAR-10 Train / validation / test 64/16/20% 64/16/20% 64/16/20% 48k/12k/10k 2,569/550/551 45k/5k/10k Batch size 256 256 512 300 32 200 Learning rate (teacher) 10−310^-3 10−310^-3 10−310^-3 10−310^-3 10−310^-3 10−310^-3 Learning rate (Stage 1) 10−310^-3 10−310^-3 3×10−43× 10^-4 10−310^-3 10−310^-3 10−310^-3 Learning rate (Stages 2–3) 10−310^-3 10−310^-3 10−310^-3 10−310^-3 10−310^-3 10−310^-3 Representation dimension p 19 8 54 2,304 384 384 WSSN width nhidn_hid 256 256 256 512 512 512 Regimes K 50 150 30 10 10 20 Active-unit weight ω 10.0 2.0 2.0 2.0 2.0 2.0 Sharpening factor η 1.02 1.02 1.02 1.1 1.1 1.1 Maximum inverse temperature βmax _ 150 150 150 100 300 300 Minimum regime size τ 50 100 50 Optimization. All training stages use Adam for at most 200200 epochs with an early-stopping patience of 1010 epochs. Learning rates are listed in Table I. The expert and gate penalty weights are selected from 0,10−3,10−2,10−1\0,10^-3,10^-2,10^-1\ using the validation set. For final predictions, we use the checkpoint with the highest validation AUC for classification or the lowest validation RMSE for regression. For all datasets, the inverse temperature starts at βclus=1 _clus=1 and follows Eq. (4). Augmentation. For tabular data, we generate τ−nkτ-n_k teacher-labeled synthetic samples for each regime with 0<nk<τ0<n_k<τ using Eq. (10) with ϵ=0.05ε=0.05. Image data are not augmented because sufficient training samples are already available in each regime. Architectures. For tabular data, the teacher is a ReLU MLP with two hidden layers of widths (256,64)(256,64). The WSSN has one ReLU hidden layer and operates directly on the original features, so u(x)=xu(x)=x. For MNIST, the teacher has two 5×55× 5 convolutional layers with 44 and 88 channels, respectively, each followed by average pooling, and then a fully connected layer of width 120120 and a final classification layer. The WSSN feature extractor uses a single 5×55× 5 convolution with 44 channels. Note that, for MNIST, the 4×24×244× 24× 24 feature maps are flattened without pooling. For Flowers and CIFAR-10, the teacher is a LeNet-style CNN with two 3×33× 3 convolutional layers, each with 3232 channels. The WSSN feature extractor uses a depthwise 6×66× 6 convolution with 128128 filters per input channel, followed by adaptive max pooling to one value per channel. For image datasets, the WSSN’s ReLU head, local experts, and explanatory gate all operate on the extracted representation u(x)u(x). Data Splits and Preprocessing. For tabular datasets, we randomly split the data into training, validation, and test sets in each run, using class-stratified sampling for classification tasks. For MNIST and CIFAR-10, we use the official test sets and randomly split the training sets into training and validation subsets. For Flowers, we randomly divide the 3,6703,670 images into training, validation, and test sets and resize them to 180×180180× 180. All image inputs are normalized separately for each channel. I Additional Interpretability Analyses We provide tabular routing and coefficient diagnostics, followed by additional local expert contribution maps for MNIST and Flowers. Routing and Regime Sizes. Fig. 10 shows sharp routing across all three datasets, with maximum posteriors concentrated near one and no single dominant regime. California Housing has many small regimes and a long right tail in regime size, indicating a more fragmented partition. (a) Covertype (b) Bike Sharing (c) California Housing Fig. 10: Routing sharpness and regime sizes on Covertype (a), Bike Sharing (b), and California Housing (c). Conventions follow Fig. 3 of the main paper. (a) Covertype (b) Bike Sharing (c) California Housing Fig. 11: Single-run local expert coefficients (top) and explanatory gate coefficients (bottom) on Covertype (a), Bike Sharing (b), and California Housing (c). Covertype experts show class-0 logit coefficients. Conventions follow Fig. 4 of the main paper. Local Expert Coefficients. The top row of Fig. 11 shows that the same feature can have coefficients of opposite signs across regimes. These sign changes reflect distinct local prediction rules. For example, on Covertype, Soil_Type29 has a coefficient of −7.77-7.77 in regime 12 and +7.76+7.76 in regime 15. On Bike Sharing, hr_cos has a coefficient of +1.90+1.90 in regime 36 and −1.36-1.36 in regime 33. On California Housing, the coefficient of HouseAge changes from +0.74+0.74 in regime 16 to −0.33-0.33 in regime 77. As shown in Table IV, local experts retain most or all features in the largest regime, indicating limited feature sparsity. Explanatory Gate Coefficients. The bottom row of Fig. 11 shows regime-assignment rules. On Covertype, the coefficient of Elevation is positive in regime 29 and negative in regime 26, so higher Elevation favors regime 29 over regime 26 in the explanatory gate. On Bike Sharing and California Housing, many cells are near white, indicating greater feature sparsity in the explanatory gate. These patterns agree with Table IV. In the largest regime, the gate retains an average of 5454 out of 5454 features on Covertype, 11.911.9 out of 1919 on Bike Sharing, and 1.41.4 out of 88 on California Housing. Table IV: Feature-selection statistics for the main results in Tables I and I, reported as mean (std) over five runs. Selected features counts features with absolute coefficient ≥10−3≥ 10^-3 in the largest regime by mean training-set size, separately for the explanatory gate and the local expert. Total features is the representation dimension p. Selected features Dataset Explanatory gate Local expert Total features Bike Sharing 11.9 (5.4) 17.3 (1.6) 19 Cal. Housing 1.4 (2.6) 8.0 (0.0) 8 Covertype 54.0 (0.0) 54.0 (0.0) 54 MNIST 2,286.6 (9.8) 2,290.1 (5.6) 2,304 Flowers 237.3 (127.1) 359.9 (49.6) 384 CIFAR-10 381.9 (0.9) 383.3 (0.8) 384 Local Expert Contribution Maps. Fig. 12 shows how the WSSN’s four convolutional channels capture different parts of MNIST digits. Since the feature extractor has only four channels, the Top-4 and All-channel sums coincide. For the “3”, channel 2 emphasizes the bends, channels 1 and 3 follow the horizontal strokes, and channel 0 highlights the left side. For the “7”, channel 1 highlights the top bar, channel 0 follows the vertical stroke, and channel 2 emphasizes the inner corner. For Flowers in Fig. 12, all four selected filter groups emphasize the radial structure of the dandelion seed head in row 1, a pattern also visible in both aggregate maps. For the daisy in row 3, the selected filter groups trace the white petals and yellow disc. (a) MNIST (b) Flowers Fig. 12: Local expert contribution maps on MNIST (a) and Flowers (b). For MNIST, ch j denotes convolutional channel j. Other conventions follow Fig. 6. IV Assumptions and Auxiliary Results for Section IV-C of the Main Paper We continue to use the notation and idealized regression setting of Section IV-C of the main paper. Below, we state the assumptions and the LASSO support-recovery bound used in that analysis. All probabilities in this support-recovery analysis are conditional on the training inputs. Assumption 0.1 (Restricted Eigenvalue Condition). Consider regime k with true support Sk=supp(βk∗)S_k=supp( _k^*). Let Δ∈ℝp ^p denote a generic vector, and in particular the estimation error Δ=β^k−βk∗ = β_k- _k^*. For any index set A⊆1,…,pA \1,…,p\, let ΔA _A denote the restriction of Δ to coordinates in A. We assume that Xk∈ℝnk×pX_k ^n_k× p satisfies the following condition: there exists a constant κk>0 _k>0 such that 1nk‖XkΔ‖22≥κk‖ΔSk‖22 1n_k\|X_k \|_2^2\;≥\; _k\| _S_k\|_2^2 for all vectors Δ satisfying the cone constraint ‖ΔSkc‖1≤3‖ΔSk‖1.\| _S_k^c\|_1≤ 3\| _S_k\|_1. We also assume maxj‖Xk,j‖2/nk≤1 _j\|X_k,j\|_2/ n_k≤ 1, where Xk,jX_k,j is column j of XkX_k. Assumption 0.2 (Irrepresentable Condition). Let Σk=1nkXk⊤Xk _k= 1n_kX_k X_k denote the within-regime Gram matrix. There exists γk∈(0,1] _k∈(0,1] such that ‖Σk,SkcSkΣk,SkSk−1‖∞≤ 1−γk, \| _k,S_k^cS_k\, _k,S_kS_k^-1 \|_∞\;≤\;1- _k, where, for a matrix A, ‖A‖∞=max∑ji|Aij|\|A\|_∞= _i _j|A_ij| is the ℓ∞ _∞ operator norm (maximum absolute row sum). Assumption 0.3 (Beta-Min Condition). For each regime k, minj∈Sk|βk,j∗|≥Ckλ _j∈ S_k| _k,j^*|≥ C_kλ for Ck>‖Σk,SkSk−1‖∞+4/κkC_k>\| _k,S_kS_k^-1\|_∞+4/ _k, where λ is the LASSO regularization hyperparameter. These conditions are standard in high-dimensional sparse regression. The restricted eigenvalue condition rules out degenerate design directions, the irrepresentable condition limits dependence between on-support and off-support variables, and the beta-min condition ensures that non-zero coefficients are large enough to be detected. Under Gaussian noise, these assumptions and the penalty choice below yield the support-recovery bound in Theorem 0.1 (see, e.g., [1, 2]). Since the regime selector is estimated, we additionally quantify how accurately it recovers the regime assignments. Assumption 0.4 (Regime Learning Accuracy). Let g g be a regime selector learned from nD_n. For the sequence of training inputs considered here, there exists a sequence ηn→0 _n→ 0 such that supx∈ℙn(g^(x)≠g∗(x))≤ηn. _x P_D_n ( g(x)≠ g^*(x) )≤ _n. Theorem 0.1 (LASSO support recovery; adapted from Wainwright [1]). Consider the linear model within the true regime k∗k^*: y=Xβ∗+ε,ε∼(0,σ2In),y=Xβ^*+ , (0,σ^2I_n), where n=nk∗n=n_k^*, β∗=βk∗β^*= _k^*^*, σ>0σ>0, and X=Xk∗∈ℝn×pX=X_k^* ^n× p is fixed. The LASSO estimator minimizes ‖y−Xb‖22/(2n)+λ‖b‖1\|y-Xb\|_2^2/(2n)+λ\|b\|_1 over b∈ℝpb ^p. Suppose Assumptions 0.1–0.3 hold for this regime, with λ=Aλσlogp/nλ=A_λσ p/n and Aλ≥4/γk∗A_λ≥ 4/ _k^*. Then there exist constants c1,c2>0c_1,c_2>0, depending only on the design constants, such that ℙ(S^k∗≠Sk∗)≤c1e−c2nk∗λ2/σ2.P ( S_k^*≠ S_k^* )≤ c_1e^-c_2n_k^*λ^2/σ^2. (17) Here S^k∗ S_k^* is the support estimated from the nk∗n_k^* oracle training samples in regime k∗k^*. The bound follows from Theorem 1 of [1] applied to y/σy/σ with penalty λ/σλ/σ. Proposition 0.1 (Oracle End-to-End Support-Recovery Bound). Fix a deterministic signal x with true regime k∗=g∗(x)k^*=g^*(x) and margin m(x)>0m(x)>0, and let k^=g^(x+δ) k= g(x+δ) denote the predicted regime. Each expert is fitted on the oracle index set ℐk∗I_k^* defined in Section IV-C of the main paper. The training data nD_n are independent of the test perturbation δ. Under Assumptions IV.1–IV.3 of the main paper, Assumptions 0.1–0.4 above, and the stated penalty choice, Proposition IV.2 and Theorem 0.1 give ℙ(S^k^≠Sk∗) ( S_ k≠ S_k^* ) ≤ηn⏟regime learning+c1e−c2nk∗λ2/σ2⏟support recovery ≤ _n_regime learning+ c_1e^-c_2n_k^*λ^2/σ^2_support recovery (18) +exp(−12(m(x)2σgBv−2log(K−1))+2)⏟perturbation stability. + \! (- 12 ( m(x)2 _gB_v- 2 (K-1) )_+^2 )_perturbation stability. Here c1,c2>0c_1,c_2>0 are the constants in Theorem 0.1. The probability is over training noise, learner randomness, and the independent perturbation δ, conditional on the training inputs and with x fixed. V Proofs V-A Proof of Proposition IV.2 Since k∗=g∗(x)k^*=g^*(x) and m(x)=minj≠k∗(vk∗−vj∗)⊤x,m(x)= _j≠ k^*(v_k^*^*-v_j^*) x, misassignment implies that some competitor j≠k∗j≠ k^* satisfies vj∗⊤(x+δ)≥vk∗⊤(x+δ).v_j^* (x+δ)≥ v_k^*^* (x+δ). With aj=vj∗−vk∗a_j=v_j^*-v_k^*^*, this implies aj⊤δ≥m(x),a_j δ≥ m(x), so g∗(x+δ)≠k∗⊆Mk∗(δ)≥m(x).\g^*(x+δ)≠ k^*\ \M_k^*(δ)≥ m(x)\. The function δ↦Mk∗(δ)δ M_k^*(δ) is L-Lipschitz with L=maxj≠k∗‖aj‖2,L= _j≠ k^*\|a_j\|_2, since each map δ↦aj⊤δ a_j δ is linear and ‖aj‖2\|a_j\|_2-Lipschitz. By the Borell–Tsirelson–Ibragimov–Sudakov inequality for Lipschitz functions of Gaussian vectors [3], ℙδ(Mk∗(δ)−δMk∗(δ)≥t)≤exp(−t22σg2L2),P_δ\! (M_k^*(δ)-E_δM_k^*(δ)≥ t )≤ \! (- t^22 _g^2L^2 ), for all t≥0t≥ 0. Set t=m(x)−δMk∗(δ),t=m(x)-E_δM_k^*(δ), and apply the concentration inequality for t≥0t≥ 0 to obtain Eq. (13). For t<0t<0, the positive-part term is zero and the right-hand side of Eq. (13) equals 11, so the bound follows from the trivial probability bound. The stated bound on δMk∗(δ)E_δM_k^*(δ) follows from a standard maximal inequality for sub-Gaussian random variables [3]. Finally, ‖aj‖2=‖vj∗−vk∗‖2≤‖vj∗‖2+‖vk∗‖2≤2Bv.\|a_j\|_2=\|v_j^*-v_k^*^*\|_2≤\|v_j^*\|_2+\|v_k^*^*\|_2≤ 2B_v. Substituting these two bounds into Eq. (13) yields Eq. (14). V-B Proof of Proposition 0.1 We decompose the event S^k^≠Sk∗\ S_ k≠ S_k^*\ (the estimated support in the predicted regime differs from the true support in the true regime) into three error events: S^k^≠Sk∗ \ S_ k≠ S_k^*\ ⊆g^(x+δ)≠g∗(x+δ)⏟(i) learned regime selector errors \ g(x+δ)≠ g^*(x+δ)\_(i) learned regime selector errors ∪g∗(x+δ)≠g∗(x)⏟(i) true regime selector switches ∪\; \g^*(x+δ)≠ g^*(x)\_(i) true regime selector switches ∪S^k∗≠Sk∗,k^=k∗⏟(i) LASSO fails in the correct regime. ∪\; \ S_k^*≠ S_k^*,\, k=k^*\_(i) LASSO fails in the correct regime. We bound the three probabilities as follows: • Assumption 0.4 bounds the learned selector’s error by ηn _n uniformly over inputs. Independence of the training randomness and δ lets us apply this conditional bound at x+δx+δ. • Proposition IV.2 bounds the probability that the true regime switches under δ. • The third event is contained in S^k∗≠Sk∗\ S_k^*≠ S_k^*\. For the oracle-trained expert in regime k∗k^*, Theorem 0.1 bounds its probability by c1exp(−c2nk∗λ2/σ2)c_1 (-c_2n_k^*λ^2/σ^2). The union bound yields Eq. (18). References [1] M. J. Wainwright, “Sharp thresholds for high-dimensional and noisy sparsity recovery using ℓ1 _1-constrained quadratic programming (LASSO),” IEEE Transactions on Information Theory, vol. 55, no. 5, p. 2183–2202, May 2009. [2] P. Bühlmann and S. van de Geer, Statistics for High-Dimensional Data: Methods, Theory and Applications. Springer Berlin Heidelberg, 2011. [3] R. Vershynin, High-Dimensional Probability: An Introduction with Applications in Data Science. Cambridge University Press, Sep. 2018.