Paper deep dive
Beyond the Mean: Distribution-Aware Loss Functions for Bimodal Regression
Abolfazl Mohammadi-Seif, Carlos Soares, Rita P. Ribeiro, Ricardo Baeza-Yates
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 94%
Last extracted: 3/26/2026, 1:30:26 AM
Summary
The paper introduces a family of distribution-aware loss functions for bimodal regression, addressing the 'mean-collapse' behavior of standard regression models (like MSE) when dealing with bimodal error distributions. By integrating normalized RMSE with Wasserstein and Cramér distances, the proposed framework enables deep regression models to recover bimodal distributions without the optimization instability associated with Mixture Density Networks (MDNs). The approach is validated across four experimental stages, demonstrating improved fidelity and robustness in aleatoric uncertainty estimation.
Entities (5)
Relation Signals (3)
Mean Squared Error → causes → mean-collapse
confidence 95% · standard regression models trained with Mean Squared Error (MSE) fail to capture this bimodal reality... leading to mean-collapse behavior
Mixture Density Networks → exhibits → optimization instability
confidence 95% · they suffer from severe optimization instability
Wasserstein loss → reduces → Jensen-Shannon Divergence
confidence 90% · reducing Jensen-Shannon Divergence by 45% on complex bimodal datasets
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Despite the strong predictive performance achieved by machine learning models across many application domains, assessing their trustworthiness through reliable estimates of predictive confidence remains a critical challenge. This issue arises in scenarios where the likelihood of error inferred from learned representations follows a bimodal distribution, resulting from the coexistence of confident and ambiguous predictions. Standard regression approaches often struggle to adequately express this predictive uncertainty, as they implicitly assume unimodal Gaussian noise, leading to mean-collapse behavior in such settings. Although Mixture Density Networks (MDNs) can represent different distributions, they suffer from severe optimization instability. We propose a family of distribution-aware loss functions integrating normalized RMSE with Wasserstein and Cramér distances. When applied to standard deep regression models, our approach recovers bimodal distributions without the volatility of mixture models. Validated across four experimental stages, our results show that the proposed Wasserstein loss establishes a new Pareto efficiency frontier: matching the stability of standard regression losses like MSE in unimodal tasks while reducing Jensen-Shannon Divergence by 45% on complex bimodal datasets. Our framework strictly dominates MDNs in both fidelity and robustness, offering a reliable tool for aleatoric uncertainty estimation in trustworthy AI systems.
Tags
Links
- Source: https://arxiv.org/abs/2603.22328v1
- Canonical: https://arxiv.org/abs/2603.22328v1
Trouble viewing inline? Open PDF directly →
Full Text
63,686 characters extracted from source content.
Expand or collapse full text
Beyond the Mean: Distribution-Aware Loss Functions for Bimodal Regression Abolfazl Mohammadi-Seif 1* , Carlos Soares 2 , Rita P. Ribeiro 3,4 , Ricardo Baeza-Yates 1 1* Department of Engineering, Universitat Pompeu Fabra, Roc Boronat St., Barcelona, 08018, Barcelona, Spain. 2 Department of Engineering, University of Porto, Pra ̧ca de Gomes Teixeira, Porto, 4099-002, Portugal. 3 Faculty of Sciences, University of Porto, 4169-007 Porto, Portugal. 4 INESC TEC, 4200-465 Porto, Portugal. *Corresponding author(s). E-mail(s): abolfazl.mohammadiseif@upf.edu; Contributing authors: csoares@fe.up.pt; rpribeiro@fc.up.pt; ricardo.baeza@upf.edu; Abstract Despite the strong predictive performance achieved by machine learning models across many application domains, assessing their trustworthiness through reliable estimates of predictive confidence remains a critical challenge. This issue arises in scenarios where the likelihood of error inferred from learned representations follows a bimodal distribution, resulting from the coexistence of confident and ambiguous predictions. Standard regression approaches often struggle to ade- quately express this predictive uncertainty, as they implicitly assume unimodal Gaussian noise, leading to mean-collapse behavior in such settings. Although Mixture Density Networks (MDNs) can represent different distributions, they suffer from severe optimization instability. We propose a family of distribution- aware loss functions integrating normalized RMSE with Wasserstein and Cram ́er distances. When applied to standard deep regression models, our approach recov- ers bimodal distributions without the volatility of mixture models. Validated across four experimental stages, our results show that the proposed Wasserstein loss establishes a new Pareto efficiency frontier: matching the stability of stan- dard regression losses like MSE in unimodal tasks while reducing Jensen-Shannon Divergence by 45% on complex bimodal datasets. Our framework strictly domi- nates MDNs in both fidelity and robustness, offering a reliable tool for aleatoric uncertainty estimation in trustworthy AI systems. 1 arXiv:2603.22328v1 [cs.LG] 20 Mar 2026 Keywords: Bimodal Regression, Distributional Loss, Trustworthy AI 1 Introduction Deep learning models, such as EfficientNet, have achieved remarkable success in image classification. However, a model’s confidence is often not calibrated to its actual cor- rectness. To build trustworthy systems, we often need a secondary complexity predictor that estimates the likelihood of error for a given input based on its embeddings. Stan- dard approaches typically treat this complexity as a continuous scalar value to be regressed, implicitly assuming a unimodal Gaussian distribution of difficulty. In practice, this assumption is not well supported. As we show in our analysis, the distribution of prediction errors is not Gaussian but exhibits a clear bimodal shape, reflecting a split between confidently correct and systematically hard samples. This mismatch suggests that regression-based complexity predictors may be inherently mis-specified, motivating the need for a more faithful modeling approach. To better understand whether this assumption holds in practice, we conducted a preliminary investigation using standard image classification benchmarks. We trained an EfficientNet-V2-M model [1] on the CIFAR-100 and Food-101 datasets using stan- dard cross-entropy loss. Post-training, we computed the prediction error for each sample x i with true class label c i . We denote this prediction error by y i , defined as the complement of the model’s predicted probability assigned to the ground truth class: y i = 1− ˆ p(c i | x i )(1) Analysis of the distribution of these y i values reveals a distinct bimodal shape rather than a Gaussian distribution (Fig. 1a). The modes correspond to two regimes: easy, where the classifier is correct and confident ( ˆ p ≈ 1,e ≈ 0), and hard, where the classifier misclassifies or has low confidence ( ˆ p≪ 1,e > 0). This empirical observation forms the primary motivation for our work. As illus- trated in Figure 1 (b), standard regression models trained with Mean Squared Error (MSE) fail to capture this bimodal reality. Because MSE minimizes the squared dif- ference from the mean, the model’s predictions collapse to a single mode in the center, the sparse region between the two actual clusters. Consequently, the model predicts a state of average difficulty that rarely exists in practice, rather than identifying the distinct easy or hard nature of the samples. To bridge this gap, we propose a distribution-aware loss framework. By integrat- ing normalized RMSE with statistical distance metrics, we treat the regression target not as a single point, but as a continuous probability measure. This formulation com- pels the model to align its predictive distribution with the true target distribution, theoretically enabling the recovery of bimodal distribution without requiring complex architectural changes. We validate this framework through a comprehensive four-stage experimental protocol to rigorously test its stability and fidelity. We formulate the following research questions to guide our investigation. 2 (a)(b) Fig. 1: Bimodal distribution of targets (a). Standard regression models trained with MSE fail to capture the underlying bimodal distribution (b). RQ1 (Fidelity): Can distance-based loss functions recover bimodal error distribu- tions where standard point-estimate losses suffer from mode collapse? RQ2 (Stability): How do these losses behave during the transition from simple to complex noise distributions? Can they avoid the optimization instability often observed in mixture models? RQ3 (Robustness): Are the proposed losses robust when applied to real-world tabular tasks with unknown target distributions? RQ4 (Application): In the motivating context of instance difficulty assessment, can the framework accurately distinguish between easy and hard samples? The remainder of this paper is organized as follows. Section 2 reviews related work. Section 3 details the proposed loss framework and its mathematical for- mulation. Section 4 presents the four-stage experimental evaluation, ranging from synthetic benchmarks to high-dimensional computer vision tasks, while Section 5 presents the results. Section 6 discusses the implications of the findings regarding the stability-fidelity trade-off. Finally, Section 7 concludes the study and outlines future directions. 2 Related Work This research lies at the intersection of three domains: probabilistic modeling of aleatoric uncertainty, regression techniques for imbalanced or non-parametric distri- butions, and the motivating application of complexity assessment in deep learning. 2.1 Aleatoric Uncertainty and Mixture Models To generalize beyond scalar point estimates, probabilistic regression methods explicitly model predictive uncertainty. Heteroscedastic Regression [2] addresses this by mod- eling Aleatoric Uncertainty, the noise inherent in the data itself. Since our goal is to 3 estimate the likelihood of error for a specific input, an intrinsic property of the sam- ple relative to the model, our work is grounded in the principles of heteroscedastic regression. As formalized in deep learning by Kendall and Gal [3], these models output both a mean μ(x) and a variance σ 2 (x), optimized via the Gaussian Negative Log-Likelihood (NLL). This allows the network to attenuate the loss for ambiguous samples by pre- dicting high variance, effectively learning to express lack of confidence rather than forcing a mean prediction. To theoretically approximate any conditional probability density function, includ- ing multimodal ones, Mixture Density Networks (MDNs) [4] are often employed. By outputting the parameters of a Gaussian Mixture Model, MDNs can capture multi- ple modes explicitly. We note that while implicit density estimation via Energy-Based Models (EBMs) [5, 6] offers an alternative class of state-of-the-art performance, they typically require computationally expensive sampling (e.g., MCMC) during inference. However, MDNs are notoriously unstable in practice. Makansi et al. [7] highlight that MDNs often suffer from mode collapse and training instabilities. They proposed an Evolving Winner-Takes-All (EWTA) loss to prevent hypothesis collapse during training. Despite these challenges, we compare against the standard MDN to evaluate the trade-off between the theoretical flexibility of mixtures and the stability of our proposed approach. 2.2 Deep Imbalanced and Non-Parametric Regression Our specific focus on bimodal error distributions relates closely to the broader field of Deep Imbalanced Regression (DIR). While standard regression assumes uniform coverage of the target space, real-world data is often highly skewed. Yang et al. [8] formally define DIR as learning from data with continuous targets where the distribu- tion is non-uniform, emphasizing that accuracy in rare and extreme values is crucial for generalization. In the context of deep learning, recent work has attempted to mitigate this imbal- ance through density estimation and statistical adjustment. Yang et al. [8] proposed Label Distribution Smoothing (LDS), which uses kernel density estimation to smooth empirical label frequencies, allowing the model to weight samples based on their effective density. Similarly, Ren et al. [9] introduced Balanced MSE, which statisti- cally adjusts the training objective to compensate for the shift between the training distribution and the test metric. Other deep learning approaches focus on sample generation or explicit weighting. Tian et al. [10] proposed DIRVAE, a generative framework combining Variational Autoencoders with LSTMs to synthesize missing minority samples. Steininger et al. [11] introduced DenseWeight, which scales gradients inversely to sample density using Kernel Density Estimation (KDE). While effective, these methods primarily focus on re-weighting or augmenting data to fix the imbalance, rather than modeling the distributional shape of the prediction error itself. Outside of deep learning, statistical modeling has explicitly addressed bimodality through parametric families. For instance, Vasconcelos et al. [12] proposed the Odd Log-Logistic Exponential Gaussian (OLLExGa) for agricultural data. By extending 4 the normal distribution with additional shape parameters, this method achieves the flexibility required to fit two distinct modes. However, applying such parametric mod- els to deep neural networks is challenging, as the latent error distribution of a network is unknown a priori and may not conform to specific rigid functional forms. To address distribution shape without enforcing strict parametric assumptions, Quantile Regression [13] offers a robust alternative. Rather than estimating a single mean, this approach minimizes the Pinball Loss to learn conditional quantiles. For- mally, for a given quantile level τ ∈ (0, 1), the model learns a value q such that the probability of the target being less than q is τ . By estimating multiple quantiles simul- taneously, the model can approximate the bounds and median of complex, multi-modal distributions without assuming normality [14]. Similarly, Parzen [15] established foundational work on consistent non-parametric density estimation using kernel functions. In the context of deep learning, Frogner et al. [16] introduced a loss function based on the Wasserstein distance (Earth Mover’s Distance) for multi-label learning, demonstrating that penalizing predictions based on semantic distance encourages smoothness and structural alignment. Our work adapts these distributional distances to regression, using them to explicitly align the predicted error distribution with the bimodal ground truth. 2.3 Application: Complexity and Difficulty Assessment Quantifying the difficulty of a training sample is central to paradigms such as Cur- riculum Learning [17], which posits that training on samples in increasing order of complexity improves convergence and generalization. Hacohen and Weinshall [18] fur- ther formalized this by decomposing curriculum learning into difficulty scoring and pacing functions, showing that ideal curricula modify the optimization landscape to be steeper around the global minimum. While traditional curricula often rely on loss magnitude, recent approaches define difficulty more intrinsically. For instance, Toneva et al. [19] characterize difficulty via forgetting events, while Baldock et al. [20] use prediction depth to rank samples from easy to hard. Crucially, however, these methods effectively flatten difficulty into a con- tinuous scalar index. We argue that this assumes a smooth spectrum of hardness that obscures the underlying structure of the problem: difficulty in deep classification is often structurally bimodal (distinctly easy vs. distinctly confusing). Standard regres- sion losses like MSE fail in this regime because they are minimized by the arithmetic mean, causing predictions to collapse into the sparse region between the two modes rather than identifying the distinct clusters. While methods such as Deep Ensembles [21] and MC-Dropout [22] provide robust uncertainty quantification to address this, they typically require multiple forward passes or extensive memory overhead. Our work focuses specifically on single-pass complexity assessment, where the goal is to imbue a standard regression model with distributional awareness without the computational cost of ensembling. 5 3 Methodology To address the limitations of existing methods, specifically the tendency of standard regression to collapse to the mean and the optimization instability of mixture models, we propose a family of distribution-aware loss functions. These losses combine point- wise accuracy with distributional alignment and range matching. Unlike parametric approaches (e.g., OLLExGa [12]) that assume rigid distributions, or Mixture Density Networks (MDN) that suffer from stability issues, our approach directly minimizes empirical divergences without parametric assumptions. Our framework is modular, consisting of three components: a pointwise error metric, a distributional divergence metric, and a range alignment term. 3.1 Component 1: Pointwise Fidelity The first component represents the typical regression objective: minimizing the distance between the prediction and the target value. While our framework accom- modates any standard metric (e.g., Mean Absolute Error, Huber Loss), we utilize the Root Mean Squared Error (RMSE) to penalize large deviations. For numerical stability, we include a small constant ε (e.g., 10 −8 ) inside the rad- ical to ensure gradients remain bounded when the error approaches zero. Since D represents a distance-like quantity between predictions and targets, we define the raw pointwise metric as: D RMSE = v u u t 1 n n X i=1 (y i − ˆ y i ) 2 ! + ε(2) 3.2 Component 2: Distributional Fidelity Standard pointwise losses (like Eq. 2) treat samples independently. Consequently, a model can achieve low MSE by predicting the conditional mean for every sample, even if the underlying problem structure is bimodal. This collapse to the mean ignores the aggregate shape of the error distribution, which is the primary signal of interest in complexity assessment. To recover this structure, we augment the loss with a divergence term D(F ˆy ,F y ) that penalizes differences between the distribution of predictions F ˆy and the distribu- tion of targets F y within a batch. A common choice for matching distributions is the Kullback-Leibler (KL) Diver- gence. However, KL-Divergence requires the estimation of continuous probability densities from discrete batches. In regression settings, this necessitates Kernel Density Estimation (KDE) or histogram binning, which can be sensitive to hyperparameter choice. Furthermore, KL-Divergence is theoretically undefined or infinite when the support of the predicted and target distributions do not overlap, a common occurrence during the early stages of training. To avoid these issues, we select Statistical Distance Metrics over density-based divergences. Unlike KL, metrics such as the Wasserstein and Cram ́er distances operate directly on the geometry of the CDFs or quantiles [23]. This provides consistent, 6 non-vanishing gradients even when the distributions are disjoint. We investigate two specific metrics: 1. Cram ́er Distance: Defined as the integrated squared difference between the Cumulative Distribution Functions (CDFs): D Cram ́er = Z ∞ −∞ (F ˆy (x)− F y (x)) 2 dx(3) 2. Wasserstein Distance: Also known as the Earth Mover’s Distance, defined for 1D distributions as the L 1 distance between CDFs [24]: D Wasserstein = Z ∞ −∞ |F ˆy (x)− F y (x)|dx(4) These metrics offer complementary strengths. As noted in [25], the squared term in the Cram ́er distance makes it highly sensitive to higher-order moment differences, which is beneficial for separating distinct modes. Conversely, the Wasserstein distance provides a linear gradient proportional to the transport cost, excelling at detecting smooth shifts in probability mass. 3.3 Component 3: Range Alignment In bimodal regression, models often shrink the predicted range to minimize variance. To explicitly counter this, we introduce a range alignment penalty: D Range =| ˆ y max − ˆ y min − (y max − y min )|(5) 3.4 Scale Invariance and Final Loss To ensure balanced optimization between error types, we apply a non-linear normal- ization to map all raw metrics D into a bounded interval [0, 1]. The normalization function N is given by: N (D) = 1− 1 1 + D (6) We note that minimizing the normalized RMSE component (N (D RMSE )) remains equivalent to minimizing the standard MSE. Since N (x) and f (x) = √ x + ε are both strictly monotonically increasing functions for x≥ 0, their composition preserves the global minimizer of the standard squared error. The final loss function is the weighted sum of the normalized components: L =N (D RMSE ) + α·N (D Distribution ) + β·N (D Range )(7) where D Distribution is either D Cram ́er or D Wasserstein . The hyperparameter α ≥ 0 controls the distributional alignment, and β ≥ 0 controls the range constraint. 7 4 Experimental Setup 4.1 Datasets To rigorously validate our hypothesis, we designed a hierarchical experimental protocol using 17 distinct datasets divided into four specific stages, as summarized in Table 1. This progression allows us to evaluate the loss functions in increasingly complex envi- ronments, ranging from controlled synthetic distributions to high-dimensional visual feature spaces. Table 1: Summary of the 17 datasets used across the four-stage experimental protocol. IDDatasetTypeTarget Modality Stage I: Synthetic Proof of Concept 1Inverse SquareSynthetic TabularConditional Bimodal 2Two PathSynthetic TabularConditional Bimodal Stage I: Controlled Separation (Unimodal→ Bimodal) 3AirfoilReal-World TabularUnimodal* 4Bike Sharing (Hour)Real-World TabularUnimodal* 5Concrete StrengthReal-World TabularUnimodal* Stage I: Natural Bimodality 6Houses (OpenML)Real-World TabularBimodal 7Protein StructureReal-World TabularBimodal 8Energy Efficiency (Heating Load)Real-World TabularBimodal 9Energy Efficiency (Cooling Load)Real-World TabularBimodal 10Bike Sharing (Day)Real-World TabularTrimodal Stage IV: Image Complexity Assessment 11CIFAR-10Image Error PredictionUnimodal 12Fashion-MNISTImage Error PredictionUnimodal 13SVHNImage Error PredictionUnimodal 14CIFAR-100Image Error PredictionBimodal 15Oxford Flowers 102Image Error PredictionBimodal 16Food-101Image Error PredictionBimodal 17Caltech-256Image Error PredictionBimodal * Datasets subjected to separation injection parameter s ∈ [0, 1] to transition from unimodal to bimodal distributions. 4.2 Synthetic and Real-World Data (Stages I–I) We begin by evaluating the loss functions on tabular benchmarks categorized by their source and complexity structure: Stage I: Synthetic Proof of Concept. We utilize two synthetic datasets, Inverse Square and Two Path (Rows 1–2), explicitly constructed to exhibit strong conditional bimodality. To ensure reproducibility, we define their generation processes as follows: 8 • Inverse Square: Targets are sampled uniformly y ∼ U [−3, 3]. The informative fea- ture is generated as x = y 2 +ε, with Gaussian noise ε∼N (0, 0.5). The model must learn the inverse mapping x7→ y ≈± √ x, representing two symmetric modes. • Two Path: Data is generated from the noisy parametric equations of a circle with radius r = 5: x = r cos(φ) + ε x and y = r sin(φ) + ε y , where φ ∼ U [0, 2π] and ε∼N (0, 0.3). For any input x, the conditional density P (y|x) contains two modes corresponding to the upper and lower arcs. Both datasets utilize 2-dimensional feature vectors consisting of the generated signal x and an independent Gaussian noise channel to simulate irrelevant features. Stage I: Controlled Separation. To study the transition from unimodal to bimodal distributions, we employ three standard regression datasets: Airfoil [26], Bike Sharing (Hour) [27], and Concrete Strength [28] (Rows 3–5). Since these datasets are originally unimodal, we introduce a controllable separation parameter S ∈ [0, 1] to artificially induce bimodality. The transformation mechanism operates by first identifying the natural latent clus- ters in the target variable y using K-Means (K = 2). Considering c mid as the midpoint between the cluster centroids, we apply a divergent linear transformation: y ′ = ( y· (1− S) + y min · Sif y ≤ c mid y· (1− S) + y max · S if y > c mid (8) where S = 0 retains the original unimodal distribution, and S → 1 forces the clusters toward the domain boundaries (y min ,y max ), creating a distinct bimodal structure with a clean separation gap. Stage I: Natural Bimodality. We include five real-world datasets (Rows 6– 10) that exhibit naturally occurring bimodality without modification. These datasets consist of California Housing [29, 30], Protein Structure, Energy Efficiency [31] and Bike Sharing (Day) [27] where unobserved latent variables often create distinct regimes in the target variable, making them ideal benchmarks for distribution-aware losses. 4.3 Image Data Generation (Stage IV) A critical contribution of this work is evaluating complexity assessment on high- dimensional visual data. For the image complexity assessment, we utilize seven established computer vision datasets: CIFAR-10 [32], Fashion-MNIST [33], Street View House Numbers (SVHN) [34], CIFAR-100 [32], Oxford Flowers 102 [35], Food-101 [36], and Caltech-256 [37]. Unlike tabular data where targets are intrinsic, the regression targets for image datasets must be derived from a classification task. We utilized a two-step procedure for the image benchmarks (Rows 11–17): Target Generation (Classification). We trained an EfficientNet-V2-M classifier [1] on the training split of each dataset. The model was optimized using Adam with a weight decay of 10 −4 . After convergence, we utilized the unseen test split to calculate the scalar error y i = 1− P (c i |x i ). This scalar serves as the ground truth complexity (y) for the subsequent regression task. We utilize this bounded metric [0, 1] rather than Cross-Entropy (− logP ) to ensure numerical stability. 9 Feature Extraction. To isolate the performance of the loss function from the capacity of the feature extractor, we froze the backbone of an EfficientNet-B0 pre- trained on ImageNet. We replaced the classification head with an identity mapping to extract a 1280-dimensional embedding vector for every image. Thus, for the regression experiments, the model receives a fixed 1280-dimensional feature vector (x) and must predict the scalar classification error (y) derived from the test samples in previous step. 4.4 Data Preprocessing To ensure training stability and comparability across datasets with varying physical units and magnitudes, we applied standard scaling (z-score normalization) to both the input features and the regression targets. This transformation standardizes the data by subtracting the arithmetic mean and dividing by the standard deviation. We note that, to prevent data leakage, these statistics are derived exclusively from the training split and subsequently applied to the validation and test sets. This prepro- cessing strategy was applied uniformly across all dataset types. For tabular datasets, all numerical input features were standardized. For the image error datasets, the same standardization was applied to the 1280-dimensional feature embeddings before they were fed into the MLP. Finally, all regression targets (e.g., prediction error, physical loads) were standardized to zero mean and unit variance. 4.5 Proposed Variants and Ablation To isolate the individual contributions of the distance metrics and the range constraint, we evaluate four distinct configurations of our framework. This decomposition allows us to address RQ2 (Stability) by comparing the gradient behavior of Cram ́er vs. Wasserstein, and RQ1 (Fidelity) by ablating the range term. For brevity, we abbreviate Wasserstein as Wasser in the reported results: • Cram ́er-Simple / Wasser-Simple (β = 0): These variants focus solely on the trade-off between pointwise error and distributional matching. By setting β = 0, we evaluate the capability of the raw distance metrics to recover bimodal structures without explicit boundary constraints. • Cram ́er-Range / Wasser-Range (β = α/2): These variants incorporate the range alignment term. We fix β = α/2 as a heuristic stability constraint. This ensures that the distributional fidelity term (α) remains the primary optimization objective, while the range term acts as a soft regularizer to prevent support collapse. 4.6 Network Architecture To ensure that performance variations are strictly attributable to the proposed loss functions rather than architectural discrepancies, we employ a unified experimental framework across all datasets. For all 17 datasets, ranging from 2D synthetic data to 1280D image embeddings, we utilize an identical Multi-Layer Perceptron (MLP) architecture. This acts as a consistent probe to evaluate the loss functions. 10 The architecture consists of an input layer adapting to the feature dimension (1 for synthetic, 1280 for images), followed by three hidden layers with 512, 256, and 128 units, respectively, and a final bottleneck layer of 64 units before the single scalar output. To mitigate overfitting, we employ Batch Normalization and GELU activation after each hidden layer, along with Dropout rates of 0.3 for the initial layers and 0.2 for deeper layers. All models are implemented in PyTorch [38] and trained using the Adam optimizer [39] with a learning rate of 10 −3 and a batch size of 64 for 50 epochs. For the MDN baseline comparison, we employ a network with a shared hidden rep- resentation comprising two fully connected layers with Tanh activations. This shared trunk branches into three separate heads to predict the mixing coefficients (π), means (μ), and standard deviations (σ) for K = 5 Gaussian components. The standard devi- ations are constrained to be positive using the ELU activation plus a stability constant (1 + ε). 4.7 Baselines To rigorously evaluate the proposed framework, we compare against four distinct base- lines representing the spectrum from naive deterministic regression to advanced prob- abilistic modeling. All distinct architectures share a common Multi-Layer Perceptron (MLP) backbone to ensure fair comparison of the loss functions. • MDN (NLL) [Primary Competitor]: A Mixture Density Network [4] that out- puts parameters (π,μ,σ) for a Gaussian Mixture Model (GMM). Optimized via Negative Log-Likelihood (NLL), this is the theoretical gold standard for multimodal regression but is often prone to optimization instability. • HMLP (GaussianNLL) [Probabilistic Baseline]: A Heteroscedastic MLP that outputs both a mean μ(x) and a variance σ 2 (x). It is trained using the Gaussian Neg- ative Log-Likelihood (GaussianNLL) loss. This represents the standard approach for modeling aleatoric uncertainty under a unimodal Gaussian assumption. • MLPQ (Quantile) [Robust Non-Parametric Baseline]: A Quantile Regres- sion MLP that outputs fixed quantiles (e.g., τ ∈ 0.1, 0.5, 0.9) instead of distribution parameters. Optimized via the Pinball Loss (QuantileLoss), it pro- vides a distribution-free estimate of prediction intervals and serves as a robust non-parametric competitor. • MLP (MSE) [Naive Baseline]: A standard MLP optimized with Mean Squared Error (MSE). This serves as the deterministic lower bound, predicting the condi- tional mean E[y|x]. It uses the same architecture as our proposed method but lacks distributional awareness. 4.8 Evaluation Metrics To comprehensively assess the performance, we employ a dual-metric strategy that evaluates both sample-level accuracy and distributional alignment. We utilize Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) to measure the precision of individual predictions. While our primary goal is distributional fidelity, these metrics serve as a control to ensure that the model 11 retains predictive power and does not simply output random samples from a correct distribution. Since standard regression metrics fail to penalize mode collapse, we introduce three distribution-aware metrics: 1. Jensen-Shannon (JS) Divergence: To quantify the overlap between the pre- dicted and ground-truth error distributions, we utilize the JS Divergence. Unlike KL-Divergence, JS is symmetric and bounded in [0, 1], providing a stable metric for comparing the probability distributions derived from the scalar predictions. 2. Wasserstein Distance: While our training objective incorporates a learned approximation of the Wasserstein distance, we explicitly calculate the exact dis- tance (Earth Mover’s Distance) during evaluation. This serves two purposes: first, it acts as an objective, critic-independent verification of convergence, ensuring that the generator is not merely fooling a sub-optimal critic. Second, unlike JS diver- gence, which relies on density overlap, W 1 captures the geometry of the error space, quantifying the physical ”effort” required to transport the predicted distribution to the ground truth. 3. Structural Fidelity (∆ BC ): To explicitly evaluate whether the model captures the structural complexity of the error surface (e.g., the bimodal nature of Easy vs. Hard samples), we first compute the Bimodality Coefficient (BC) [40] of the predictions: BC = γ 2 + 1 κ + 3(n−1) 2 (n−2)(n−3) (9) where γ is skewness, κ is excess kurtosis, and n is sample size. Values > 0.555 indicate bimodality [40–42]. We define Structural Fidelity as the absolute deviation of the predicted BC from the target: ∆ BC =|BC target − BC pred |(10) A ∆ BC close to 0 indicates that the model successfully recovers the distribution structure, whereas high values quantify the severity of the regression to the mean (unimodal collapse). 5 Results We structure our analysis hierarchically, moving from controlled synthetic environ- ments to the motivating problem of deep learning complexity assessment. We analyze performance across four distinct dataset types: Synthetic, Unimodal Real-World Tab- ular, Bimodal Real-World Tabular, and Image Error Prediction. There is also one Trimodal Real-World Tabular for comparison. 5.1 Stage I: Validation on Synthetic Distributions To verify if the proposed loss functions can successfully capture non-Gaussian, bimodal structures in a controlled environment, and to quantify the trade-off between pointwise 12 fidelity (RMSE) and distributional fidelity (JS Divergence). We evaluate the models on the Inverse Square and Two Path synthetic datasets. Stability and Pointwise Fidelity: The primary challenge in bimodal regression is maintaining the training stability of deterministic models while introducing prob- abilistic flexibility. As shown in Table 2, the Wasser-Simple configuration achieves an RMSE of 2.62, exactly matching the HMLP baseline (2.62) and outperforming the standard MSE baseline (2.70). This confirms that our proposed loss preserves the robust convergence properties of the MLP backbone. In contrast, the MDN base- line, while theoretically capable of modeling mixtures, suffers from severe optimization instability, resulting in a significantly degraded RMSE of 3.82. Distributional Trade-off Analysis: Beyond pointwise accuracy, we analyze the trade-off between stability and distributional fidelity using Jensen-Shannon (JS) Divergence. Improving Fidelity while maintaining Stability: The HMLP baseline, constrained by its unimodal Gaussian assumption, yields a high JS Divergence of 0.73. By incorporating our transport-based loss, the Wasser-Simple model reduces this diver- gence to 0.63 without any penalty to RMSE (remaining at 2.62). This represents an improvement in distributional alignment. Maximizing Fidelity: For applications prioritizing pure distribution matching, the Cramer-Range variation reduces JS Divergence further to 0.54, effectively matching the best probabilistic baseline (MDN, 0.54). However, unlike the MDN which sacrifices stability (RMSE 3.82), the Cramer-Range maintains a competitive RMSE of 2.94. Additionally, among the tuned configurations, Cramer-Range achieves the minimum Wasserstein distance of 1.38. Shape Recovery: In terms of capturing the specific bimodal shape, the Cramer- Simple variation excels, achieving a ∆ BC of 0.14, marginally outperforming the MDN (0.15) and significantly improving upon the unimodal baselines (HMLP 0.24, MSE 0.34). Comparison with Default Configuration: We also evaluate the Wasser-Simple (Def ) configuration (α = 1), representing an unweighted, out-of-the-box setting. This model highlights a critical limitation of pointwise metrics in bimodal regression: while base- lines like HMLP and Quantile achieve lower RMSE, they do so by collapsing to a unimodal average. In contrast, the default configuration preserves the bimodal geom- etry, achieving exceptional fidelity (∆ BC of 0.03, Wasserstein of 0.30, JS Divergence of 0.37). It drastically outperforms the MDN (∆ BC of 0.15), which, despite being probabilistic, suffers from optimization instability, manifesting as shifted modes and exaggerated peak amplitudes (as observed in the Two Path dataset in Figure 2), demonstrating that our method provides superior structural recovery even without hyperparameter tuning. Qualitative Shape Recovery: To visualize the capacity for mode recovery, Figure 2 presents the prediction densities on the Two Path dataset. Visual inspection reveals that while the HMLP and Quantile baselines suffer from unimodal collapse, the MDN recovers the bimodal structure but exhibits optimization instability, evi- denced by shifted modes and exaggerated peak amplitudes. In contrast, the proposed 13 method accurately splits the probability mass, aligning closely with the ground truth geometry. (a) Proposed (α = 1)(b) MDN (c) Quantile Regression(d) Heteroscedastic MLP (HMLP) Fig. 2: Qualitative comparison of predictive distributions (Orange) versus Ground Truth (Blue) on the Two Path dataset. In summary, the Stage I results confirm that the proposed framework successfully bridges the gap between the stability of deterministic baselines and the expressive- ness of mixture models. Specifically, Wasser-Simple offers the stability of HMLP with improved divergence metrics, while Cramer-Range matches the distributional fidelity of MDNs with significantly better stability. Overall, the proposed methods achieve a balance between stability (RMSE) and shape recovery (Wasserstein/∆ BC ). Table 2 summarizes the performance of the proposed loss functions (averaged across both datasets) against the baselines. 14 Table 2: Performance comparison on synthetic datasets. ConfigurationTest Loss RMSE ↓ Wasserstein ↓ JS Div ↓ ∆ BC ↓ Wasser-Simple0.5352.621.800.630.29 Wasser-Range0.5812.901.710.590.34 Cramer-Simple0.5422.771.730.620.14 Cramer-Range0.5762.941.380.540.39 Wasser-Simple (Def)0.7773.650.300.370.03 MDN (NLL)0.0713.821.610.540.15 HMLP (Gauss)0.2802.622.150.730.24 MLPQ (Quant)0.2142.871.560.550.31 MSE (Standard)1.0502.702.200.720.34 5.2 Stage I: Continuous Mode Separation Analysis To evaluate the stability of the loss landscape during the critical phase transition from unimodal to bimodal distributions. We utilize a controlled separation parameter S ∈ [0, 1] applied to the synthetic datasets. The bimodal targets are generated via a split-and-shift mechanism: data points below the K-Means midpoint are linearly interpolated towards the domain minimum, while those above are shifted towards the maximum, governed by S. At S = 0, the distribution is strictly unimodal; as S → 1, the modes diverge into a clear bimodal structure. Table 3 presents the performance metrics at critical transition points, addressing RQ1 regarding structural recovery and RQ2 regarding optimization stability during distribution shifts. Table 3: Performance snapshots at critical separation intervals (S) in synthetic datasets. Unimodal (S = 0) Transition (S = 0.5)Bimodal (S = 1) ConfigurationRMSEJS DivRMSEJS DivRMSE JS Div ∆ BC Proposed Methods Wasser-Simple0.9390.4270.4670.3870.5090.4410.105 Wasser-Range0.9770.4300.4910.3660.5320.4360.171 Cramer-Simple0.9510.4570.4810.3920.5010.4720.129 Cramer-Range0.9880.4400.4800.3780.5180.4410.156 Wasser-Simple (Def) 1.0930.2840.4930.3220.5810.339 0.052 Baselines MDN1.1540.3140.5910.3680.6070.4110.071 HMLP0.9680.4030.4720.4040.5350.5250.143 MLPQ0.9670.3690.4600.3930.5030.4690.061 MSE0.9420.4680.4690.4090.4960.4830.140 The results highlight distinct behavioral classes among the distribution-aware losses: 15 • The Singularity Phenomenon (S = 0 → 0.1): A critical transition occurs at the onset of separation. At S = 0, the distribution is continuous and unimodal. However, the generation mechanism applies a hard split at the distribution mid- point, linearly interpolating data toward the domain boundaries as S increases. The sharp drop in error metrics at S = 0.1 marks the singularity where this split creates two statistically resolvable attractors. While Cramer and Wasserstein models adapt immediately to this structural break, the MDN exhibits instability at S = 0, failing to collapse to the single mean when no gap exists. • Cramer vs. Wasserstein (Structure): While Cramer loss is theoretically a dis- tance metric like Wasserstein, the results show it behaves closer to MSE in practice. At S = 1, Cramer-Simple yields a ∆ BC of 0.129 and Cramer-Range 0.156, both com- parable to the collapsed HMLP (0.143). This indicates that Cramer loss struggles to drive the gradient towards explicit mode separation. In contrast, Wasser-Simple (Def ) achieves a ∆ BC of 0.052, confirming superior geometric recovery. • Information Fidelity (JS Div): The Wasser-Simple (Def ) is the only configura- tion that maintains high distributional fidelity across the entire spectrum, achieving the lowest JS Divergence at both the unimodal singularity (0.284) and the bimodal peak (0.339). To validate these metrics qualitatively, Figure 3 visualizes the predictive density evolution on the Bike Sharing (Hour) dataset. As the separation parameter increases from 0.0 to 1.0, the target distribution (blue) splits into two distinct modes. While all three methods successfully identify the bimodal structure in this specific dataset, subtle differences in alignment are observable. The Wasser-Simple (Def ) configuration yields a density estimate that coincides more precisely with the target distribution peaks compared to the slightly looser approximations of MSE and MDN, confirming the superior information fidelity scores observed in Table 3. In summary, the Stage I analysis highlights a critical trade-off in existing distri- bution learning methods: the choice between optimization stability (MSE, Cramer) and structural flexibility (MDN). The observed ”Singularity Phenomenon” at S = 0.1 underscores how sensitive these models are to the emergence of resolvable signal. The proposed Wasserstein loss resolves this dichotomy. By treating the output as a contin- uous probability measure rather than a fixed set of parameters, it matches the stability of deterministic losses in the unimodal regime (S = 0) while achieving the superior shape recovery and distributional fidelity of mixture models in the bimodal regime (S = 1), all without the need for discrete switching logic or component tuning. 5.3 Stage I: Validation on Real-World Bimodal Distributions To evaluate the robustness of the proposed framework on real-world tabular data specifically selected for their bimodal properties: Boston Housing, Protein Structure, and Energy Efficiency (Heating Load and Cooling Load). Table 4 summarizes the results across real-world bimodal datasets. These find- ings address RQ3, confirming that while MDNs often fail in high-noise tabular environments, the Wasserstein approach maintains robustness comparable to MSE. 16 (a) MSE (S = 0.0)(b) MSE (S = 0.5)(c) MSE (S = 1.0) (d) MDN (S = 0.0)(e) MDN (S = 0.5)(f) MDN (S = 1.0) (g) Proposed (S = 0.0)(h) Proposed (S = 0.5)(i) Proposed (S = 1.0) Fig. 3: Evolution of predictive densities on Bike Sharing (Hour). Blue Area: Ground Truth. Orange Area: Predicted Distribution. 5.3.1 Stability and Fidelity Analysis The results in Table 4 provide empirical evidence for the two core claims of this work, positioning the proposed method against both parametric (MDN) and non-parametric (MLPQ) baselines: Overcoming MDN Instability: The MDN baseline exhibits a critical failure mode in real-world regression. While it achieves a competitive Test Loss (NLL), it yields the worst RMSE (0.525) among all models. This indicates that while the MDN fits the probability density locally, it struggles to place the probability mass correctly in the global output space. In contrast, the Wasser-Simple configuration maintains an RMSE of 0.330, effectively matching the stability of standard MSE. Outperforming Robust Baselines (MLPQ): Notably, Quantile Regression (MLPQ) emerges as the most competitive baseline, offering a much better balance 17 Table 4: Aggregated performance on Real-World Bimodal Datasets. ConfigurationTest Loss RMSE ↓ Wasserstein ↓ JS Div ↓ ∆ BC ↓ Wasser-Simple0.2320.3300.0830.3060.045 Wasser-Range0.2510.3420.0790.3000.112 Cramer-Simple0.2220.3240.0850.3060.046 Cramer-Range0.2490.3440.1010.3290.071 Wasser-Simple (Def)0.3420.3670.0810.2910.018 MDN (NLL)-0.2460.5250.0840.3000.033 HMLP (Gauss)-0.5270.4230.1830.3690.099 MLPQ (Quant)0.0710.3380.0840.3010.019 MSE (Standard)0.1670.3460.1260.3270.065 of stability (RMSE 0.338) and shape recovery than the MDN. However, the pro- posed framework consistently outperforms this strong baseline. The Wasser-Simple (Def ) configuration achieves a lower Wasserstein distance (0.081 vs 0.084) and a lower Jensen-Shannon Divergence (0.291 vs 0.301) than MLPQ. This suggests that opti- mizing the Wasserstein distance directly provides a more accurate global distribution match than minimizing quantile loss, which approximates the distribution via discrete checkpoints. Superior Shape Recovery: Regarding structural fidelity, the proposed method achieves a δ Bimodality Coefficient of 0.018, which is nearly 45% lower than the MDN (0.033). This confirms that explicitly minimizing the Wasserstein distance is more effective at recovering complex geometries than minimizing Negative Log-Likelihood, which is prone to local optima in the mixture weight parameters. To visualize these dynamics, Figure 4 presents the predictive densities on the Pro- tein Structure dataset. The standard MSE and HMLP models exhibit mode collapse, approximating the bimodal ground truth with a single average. While the MDN cap- tures the split, it tends to produce sharp, disjoint peaks. Quantile Regression (MLPQ) offers a significant improvement, yet the proposed Wasser-Simple method (especially with higher coverage at α = 7) yields the smoothest approximation of the manifold, corroborating the superior Wasserstein metrics observed in Table 4. 5.4 Case Study: Generalization to Trimodal Dynamics While our framework targets bimodality, we examine the Bike Sharing (Day) dataset to evaluate generalization to higher-order complexity (K=3) without explicit model reconfiguration. As shown in Table 5, the MDN achieves the best performance across all metrics. This is expected, as the flexibility of mixture models is ideal when the number of components aligns with the physical modes of the target. However, the proposed methods serve as robust alternatives that significantly outperform the unimodal baselines. The Wasser-Simple configuration reduces the Wasserstein distance by 48% compared to HMLP (0.052 vs. 0.100) and achieves an RMSE of 0.088, which is closer to the MDN (0.069) than the unimodal baselines. 18 (a) MSE (Standard)(b) HMLP (GaussianNLL)(c) MDN (NLL) (d) MLPQ (Quantile)(e) Proposed (Def α = 1) (f) Proposed (HighCov α = 7) Fig. 4: Qualitative comparison on the Protein Structure dataset. Top Row: MSE and HMLP collapse to unimodal averages in sparse regions, while MDN captures bimodality but exhibits disjoint peaks. Bottom Row: MLPQ and Proposed method (Center) provide a strong approximation. Notably, increasing to α = 7 (Right) results in high coverage of the bimodal support. Uniquely, the Wasser-Range configuration achieves the lowest δ Bimodality Coefficient (0.032) of all models, including the MDN (0.046). This suggests that incorporating range constraints allows the model to capture the complex spread of the trimodal dis- tribution more accurately than likelihood maximization alone, even without explicit tuning for the number of modes. 5.5 Stage IV: High-Dimensional Image Complexity Assessment The final stage validates the framework on the motivating problem: assessing the ”dif- ficulty” of an image for a pre-trained classifier. Using embeddings from seven computer vision benchmarks, we group datasets into ”Easy” (Unimodal error distribution, e.g., CIFAR-10) and ”Hard” (Bimodal error distribution, e.g., CIFAR-100). Table 6 presents the comprehensive performance, directly answering RQ4 by val- idating the framework’s ability to distinguish complexity regimes in high-dimensional embeddings and providing a final visual confirmation of the trade-off discussed in RQ2. Key observations from the comprehensive comparison include: 19 Table 5: Performance on the Trimodal Bike Sharing dataset. ConfigurationTest Loss RMSE ↓ Wasserstein ↓ JS Div ↓ ∆ BC ↓ Wasser-Simple0.0920.0880.0520.4020.060 Wasser-Range0.1090.0880.0580.4340.032 Cramer-Simple0.0990.1080.0710.4370.061 Cramer-Range0.1310.1180.0860.4130.040 Wasser-Simple (Def)0.2370.1530.1120.4520.044 MDN (NLL)-2.1050.0690.0240.3380.046 HMLP (Gauss)-0.9880.1430.1000.4420.092 MLPQ (Quant)0.0330.1010.0730.4290.051 MSE (Standard)0.0170.1310.0870.4290.061 Table 6: Comprehensive performance on Image Complexity Assess- ment. Complexity ConfigurationRMSE ↓ Wasserstein ↓ JS Div. ↓ ∆ BC ↓ Unimodal Wasser-Simple1.0500.2540.4410.151 Wasser-Range1.0600.2530.4710.134 Cramer-Simple1.0360.2800.5360.146 Cramer-Range1.0680.2510.5060.149 Wasser-Simple (Def )1.0890.1720.1570.147 MDN1.0660.2610.2180.040 HMLP1.0680.2280.2480.134 MLPQ1.0570.2570.2380.119 MSE1.0360.2870.5410.216 Bimodal Wasser-Simple1.1620.2580.5290.137 Wasser-Range1.1630.2680.5420.137 Cramer-Simple1.1500.2790.5710.134 Cramer-Range1.1460.2780.5600.128 Wasser-Simple (Def )1.2080.1700.3010.073 MDN1.2540.1850.3200.035 HMLP1.1600.2460.4460.083 MLPQ1.1410.3010.4760.118 MSE1.1140.3310.5520.167 The ”Tuning Trap”: Interestingly, the tuned versions of the proposed losses often achieve lower JS fidelity than the default version. This occurs because hyperparam- eter optimization minimizes the composite loss (RMSE + Distance). In noise-heavy regimes, the optimizer often reduces α (the weight of the distance term) to favor RMSE, causing the model to revert toward the mean. The Default configuration (α = 1) enforces a stricter adherence to distributional matching, resulting in the best JS Divergence across the board. Baseline Trade-offs: MSE is the stability king (lowest RMSE 1.114 in Bimodal) but fails structurally (Highest JS 0.552). MDN is the structure king (lowest ∆ BC 0.035) but fails optimization (Highest RMSE 1.254). 20 Fig. 5: Distributional Fidelity (JS Divergence) across all Image Datasets. The Default Advantage: The Wasser-Simple (Def ) offers the only viable middle ground. In the difficult Bimodal regime, it achieves a JS Divergence of 0.301 (outper- forming MDN’s 0.320) and a Wasserstein Distance of 0.170 (the lowest of all methods). Figure 5 visualizes this dominance across the entire spectrum, showing that while the gap is negligible for ”Easy” datasets (left), the proposed method consistently achieves the lowest information divergence on ”Hard” datasets (right). Figure 6 explicitly visualizes this trade-off. The results reveal a clear efficiency frontier: The Stability Extreme: MSE occupies the far right of the fidelity spectrum (JS: 0.552), maximizing stability (RMSE: 1.114) at the cost of total structural collapse. The Pareto Domination: Crucially, the Wasser-Simple (Def ) configuration strictly dominates the Mixture Density Network. It achieves superior distributional fidelity (JS: 0.301 vs. MDN’s 0.320) while simultaneously maintaining significantly better optimization stability (RMSE: 1.208 vs. MDN’s 1.254). This confirms that the proposed Wasserstein loss identifies the optimal operat- ing point for bimodal regression, resolving the instability of mixture models without reverting to the unimodal collapse of MSE. To validate these metrics qualitatively, Figure 7 visualizes the predictive density on CIFAR-100. The failure of MSE is evident (Fig. 7a), as it predicts a single skewed mode. While the MDN (Fig. 7b) successfully splits the distribution, it suffers from misalignment, placing the second peak in the sparse uniform region between modes. In contrast, the proposed method (Fig. 7c-d) demonstrates better alignment. Notably, while the automated tuning often selects lower weights to minimize RMSE, explicitly setting α = 1.5 (Fig. 7d) forces the model to prioritize shape, resulting in a near-perfect recovery of the bimodal target that other methods fail to capture. 21 Fig. 6: Pareto Efficiency Frontier (stage IV Bimodal Datasets). 6 Discussion 6.1 Hypothesis Testing: Constraints vs. Distance A central theoretical question of this study was determining the necessary conditions for recovering bimodal distributions in regression. • Hypothesis A (Constraints are Necessary): Posited that without explicit boundary enforcement (the ℓ Range term with β = α/2), the regression model would fail to maintain the separation between modes, collapsing back to a unimodal mean. • Hypothesis B (Distance is Sufficient): Posited that minimizing a distributional metric (Wasserstein/Cram ́er) against the batch target distribution is mathe- matically sufficient to drive mode separation, making explicit range constraints redundant. Our empirical results across all four stages strongly support Hypothesis B and reject Hypothesis A. In the ”Controlled Separation” analysis (Stage I), the Wasser-Simple configuration (β = 0) successfully recovered the bimodal structure with a ∆ BC of 0.052, matching or outperforming constrained variants. Furthermore, in real-world tabular tasks (Stage I), the Range variants frequently exhibited higher RMSE variance, suggesting that rigid boundary constraints destabilize optimization when data contains heavy-tailed noise. We conclude that ”soft” distributional pressure is not only sufficient but superior to ”hard” geometric constraints for robust bimodal regression. 22 (a) MSE Baseline(b) MDN Baseline (c) Proposed (Default α = 1)(d) Proposed (Tuned α = 1.5) Fig. 7: Qualitative comparison on CIFAR-100. 6.2 The Stability-Fidelity Trade-off The comparison with Mixture Density Networks (MDN) reveals a fundamental ”No Free Lunch” dynamic in distribution learning. • The MDN Regime: MDNs achieve the highest structural separation (lowest ∆ BC in stage IV), but at the cost of extreme optimization fragility. As seen in Stage I (S = 0), they fail to converge to the mean in unimodal regimes, and in stage IV, they yield the highest RMSE (1.254). • The Wasserstein Regime: The proposed Wasserstein framework occupies a ”Pareto Optimal” middle ground. It matches the stability of MSE in unimodal tasks (RMSE ≈ 1.05 in stage IV Easy) while capturing bimodal structure significantly better than MSE (JS Div 0.301 vs 0.552 in stage IV Hard). This suggests that for practical applications where ”worst-case” stability is as impor- tant as ”best-case” fidelity, distance-based loss functions offer a safer alternative to likelihood-based mixture models. 23 6.3 Implications for Complexity Assessment The findings in stage IV directly validate the motivating application of this work. Standard regression models (MSE) equate ”uncertainty” with ”mean error,” effectively smoothing over the distinction between ”confident” and ”confused” predictions. By successfully recovering the bimodal error distribution on datasets like CIFAR-100 and Food-101, our framework provides a more semantically meaningful signal. This allows downstream systems to distinguish between aleatoric noise (irreducible confusion, Mode 2) and epistemic confidence (clean predictions, Mode 1), a distinction that is mathematically invisible to mean-squared error. 7 Conclusion This study addressed the limitation of standard regression losses, such as Mean Squared Error (MSE), in modeling bimodal error distributions. Motivated by the observation that deep learning ”complexity” often manifests as a bimodal signal (confident vs. confused ), we proposed a distribution-aware framework integrating Wasserstein and Cram ́er distances. Through a rigorous four-stage experimental protocol, we arrived at three primary conclusions. First, we demonstrated that the transition from unimodal to bimodal distributions induces a Singularity Phenomenon, where standard models abruptly fail. While Mix- ture Density Networks (MDNs) can recover the bimodal structure, they exhibit severe optimization instability, often failing to converge even in simple unimodal baselines (S = 0). Second, we refuted the hypothesis that explicit range constraints are necessary for mode separation. Our results show that minimizing the Wasserstein distance alone is sufficient to drive the gradient towards the correct distributional shape. In fact, remov- ing the constraints (as in the Wasser-Simple configuration) yields better robustness in real-world tabular tasks by avoiding the optimization rigidity associated with fixed boundaries. Finally, in the motivating task of High-Dimensional Image Complexity Assessment, our framework established a new Pareto efficiency frontier. The default Wasserstein loss (α = 1) successfully recovered the bimodal confidence vs. confusion structure of prediction errors, achieving a 45% reduction in Jensen-Shannon Divergence compared to MSE while avoiding the training volatility of MDNs. This confirms that treating regression targets as continuous probability measures, rather than single points, pro- vides a more reliable and semantically rich path for aleatoric uncertainty estimation in trustworthy AI systems. Declarations • Funding:ThisworkhasbeenpartiallysupportedbyMCIN/AEI /10.13039/501100011033 under the Mar ́ıa de Maeztu Units of Excellence Program (CEX2021-001195-M). • Conflict of interest: The authors declare no competing interests. • Ethics approval: Not applicable. 24 • Availability of data and materials: The datasets used in this study are publicly available. • Code availability: Code will be made available upon publication. • Author contribution: Abolfazl Mohammadi-Seif: Conceptualization, Methodology, Software, Formal Analysis, Investigation, Writing – Original Draft, Visualization. Carlos Soares: Conceptualization, Methodology, Writing – Review & Editing, Supervision. Rita P. Ribeiro: Conceptualization, Methodology, Writing – Review & Editing, Supervision. Ricardo Baeza-Yates: Conceptualization, Supervision. All authors read and approved the final manuscript. References [1] Tan, M., Le, Q.: Efficientnet: Rethinking model scaling for convolutional neu- ral networks. In: International Conference on Machine Learning, p. 6105–6114 (2019). PMLR [2] Nix, D.A., Weigend, A.S.: Estimating the mean and variance of the target prob- ability distribution. In: Proceedings of 1994 Ieee International Conference on Neural Networks (ICNN’94), vol. 1, p. 55–60 (1994). IEEE [3] Kendall, A., Gal, Y.: What uncertainties do we need in bayesian deep learning for computer vision? Advances in neural information processing systems 30 (2017) [4] Bishop, C.M.: Mixture density networks (1994) [5] LeCun, Y., Chopra, S., Hadsell, R., Ranzato, M., Huang, F., et al.: A tutorial on energy-based learning. Predicting structured data 1(0) (2006) [6] Du, Y., Mordatch, I.: Implicit Generation and Generalization in Energy-Based Models (2020). https://arxiv.org/abs/1903.08689 [7] Makansi, O., Ilg, E., Cicek, O., Brox, T.: Overcoming limitations of mixture den- sity networks: A sampling and fitting framework for multimodal future prediction. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 7144–7153 (2019) [8] Yang, Y., Zha, K., Chen, Y., Wang, H., Katabi, D.: Delving into deep imbalanced regression. In: International Conference on Machine Learning, p. 11842–11851 (2021). PMLR [9] Ren, J., Zhang, M., Yu, C., Liu, Z.: Balanced mse for imbalanced visual regression. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 7926–7935 (2022) 25 [10] Tian, H., Tian, C., Li, K., Jia, W.: Unbalanced regression sample generation algorithm based on confrontation. Information Sciences 642, 119157 (2023) [11] Steininger, M., Kobs, K., Davidson, P., Krause, A., Hotho, A.: Density-based weighting for imbalanced regression. Machine Learning 110(8), 2187–2211 (2021) [12] Vasconcelos, J.C.S., Cordeiro, G.M., Ortega, E.M.M., Rezende, ́ E.M.d.: A new regression model for bimodal data and applications in agriculture. Journal of Applied Statistics 48(2), 349–372 (2021) [13] Koenker, R., Hallock, K.F.: Quantile regression. Journal of economic perspectives 15(4), 143–156 (2001) [14] Meinshausen, N., Ridgeway, G.: Quantile regression forests. Journal of machine learning research 7(6) (2006) [15] Parzen, E.: On estimation of a probability density function and mode. The annals of mathematical statistics 33(3), 1065–1076 (1962) [16] Frogner, C., Zhang, C., Mobahi, H., Araya, M., Poggio, T.A.: Learning with a wasserstein loss. Advances in neural information processing systems 28 (2015) [17] Bengio, Y., Louradour, J., Collobert, R., Weston, J.: Curriculum learning. In: Proceedings of the 26th Annual International Conference on Machine Learning, p. 41–48 (2009) [18] Hacohen, G., Weinshall, D.: On the power of curriculum learning in training deep networks. In: International Conference on Machine Learning, p. 2535–2544 (2019). PMLR [19] Toneva, M., Sordoni, A., Combes, R.T., Trischler, A., Bengio, Y., Gordon, G.J.: An empirical study of example forgetting during deep neural network learning. In: International Conference on Learning Representations [20] Baldock, R., Maennel, H., Neyshabur, B.: Deep learning through the lens of example difficulty. Advances in Neural Information Processing Systems 34, 10876–10889 (2021) [21] Lakshminarayanan, B., Pritzel, A., Blundell, C.: Simple and scalable predictive uncertainty estimation using deep ensembles. Advances in neural information processing systems 30 (2017) [22] Gal, Y., Ghahramani, Z.: Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In: International Conference on Machine Learning, p. 1050–1059 (2016). PMLR [23] Arjovsky, M., Chintala, S., Bottou, L.: Wasserstein generative adversarial net- works. In: International Conference on Machine Learning, p. 214–223 (2017). 26 PMLR [24] Villani, C.: Optimal Transport: Old and New vol. 338. Springer, (2008) [25] Bellemare, M.G., Danihelka, I., Dabney, W., Mohamed, S., Lakshminarayanan, B., Hoyer, S., Munos, R.: The cramer distance as a solution to biased wasserstein gradients. arXiv preprint arXiv:1705.10743 (2017) [26] Brooks, T.F., Pope, D.S., Marcolini, M.A.: Airfoil self-noise and prediction. Nasa reference publication 1218, National Aeronautics and Space Administration, Hampton, Virginia (1989) [27] Fanaee-T, H., Gama, J.: Event labeling combining ensemble detectors and background knowledge. Progress in Artificial Intelligence 2(2), 113–127 (2014) [28] Yeh, I.-C.: Concrete Compressive Strength. UCI Machine Learning Repository. DOI: https://doi.org/10.24432/C5PK67 (1998) [29] Pace, R.K., Barry, R.: Sparse spatial autoregressions. Statistics & Probability Letters 33(3), 291–297 (1997) [30] Pace, R.K., Barry, R.: California Housing Dataset. StatLib. Accessed via OpenML (1997). https://w.openml.org [31] Tsanas, A., Xifara, A.: Energy Efficiency. UCI Machine Learning Repository. DOI: https://doi.org/10.24432/C51307 (2012) [32] Krizhevsky, A., Hinton, G., et al.: Learning multiple layers of features from tiny images (2009) [33] Xiao, H., Rasul, K., Vollgraf, R.: Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms. arXiv preprint arXiv:1708.07747 (2017) [34] Netzer, Y., Wang, T., Coates, A., Bissacco, A., Wu, B., Ng, A.Y., et al.: Reading digits in natural images with unsupervised feature learning. In: NIPS Workshop on Deep Learning and Unsupervised Feature Learning, vol. 2011, p. 7 (2011). Granada [35] Nilsback, M.-E., Zisserman, A.: Automated flower classification over a large num- ber of classes. In: 2008 Sixth Indian Conference on Computer Vision, Graphics & Image Processing, p. 722–729 (2008). IEEE [36] Bossard, L., Guillaumin, M., Van Gool, L.: Food-101–mining discriminative com- ponents with random forests. In: European Conference on Computer Vision, p. 446–461 (2014). Springer [37] Griffin, G., Holub, A., Perona, P., et al.: Caltech-256 object category dataset. 27 Technical report, Technical Report 7694, California Institute of Technology Pasadena (2007) [38] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., Chin- tala, S.: Pytorch: An imperative style, high-performance deep learning library. In: Wallach, H., Larochelle, H., Beygelzimer, A., Alch ́e-Buc, F., Fox, E., Garnett, R. (eds.) Advances in Neural Information Processing Systems, vol. 32. Curran Associates, Inc., (2019) [39] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014) [40] Pfister, R., Schwarz, K.A., Janczyk, M., Dale, R., Freeman, J.B.: Good things peak in pairs: a note on the bimodality coefficient. Frontiers in psychology 4, 700 (2013) [41] Freeman, J.B., Dale, R.: Assessing bimodality to detect the presence of a dual cognitive process. Behavior research methods 45(1), 83–97 (2013) [42] Kang, Y.-J., Noh, Y.: Development of hartigan’s dip statistic with bimodality coefficient to assess multimodality of distributions. Mathematical Problems in Engineering 2019(1), 4819475 (2019) 28