Minimal Variance Model Aggregation: A principled, non-intrusive, and versatile integration of black box models

We have made several changes to the manuscript to improve the article.

We have modified the title, abstract, and introduction to better highlight the contributions of our work. In particular, we have revised the related work section and added a clear contributions section. We have proved a theorem to address concerns about the well-foundedness of the method. This theorem justifies the superiority of MEVA over MEEA by showing it has a faster convergence rate, making it more robustly estimated. We identify our setting as a low data regime, making robust estimation key for accuracy. We have reworked the Boston housing dataset example. Its role as a sanity check for the method has been emphasized. We also noticed the small data set size made results highly variable. We now average results over many seeds, and after better hyperparameter choice, our method shows a 2% overall improvement.

To answer your questions:

Performance on the Boston Housing dataset: We have updated the Boston Housing dataset example with better hyperparameters and averaging over multiple seeds. Note that ensemble learning techniques such as random forests and gradient boosting are state-of-the-art for regression tasks and are optimized for this setting, it is unsurprising that they remain difficult to surpass with the same data. That said, the data science setting is not the primary focus of our method or where it is most warranted. The inclusion of this experiment was intended to validate our approach by demonstrating that aggregation improves model performance. However, in practice, there is little reason to use non-trainable models in this context. In data science settings, our method serves as an effective aggregation tool, but trainable ensemble techniques like random forests and boosting remain better suited for this domain.

Robusteness of approximation of $A$: Most methods in the literature employing pointwise linear aggregation rely on convex combinations as a heuristic, effectively approximating A as diagonal. In our work, we simplify the problem by defining $A=P Diag(\lambda_i)P^T$ , which keeps the regression tractable. This approach can be compared to specifying a matrix-valued kernel in multi-dimensional Gaussian Process (GP) regression, as outlined in Appendix E. In general, defining a covariance structure across the entire domain x is a complex task, akin to recovering the full correlation structure in our setting. Commonly, GPs are assumed to be either independent or separable, where the covariance is represented as the product of a 1D kernel and a positive semi-definite (PSD) matrix Q [1]. To introduce additional complexity, kernels may be combined as described in Theorem 3 [2]:

$K(x,y) =\sum_{i=1}^p k_i(x,y) Q_i$

Our method is analogous, specifying the covariance structure as:

$A(x)=\sum_{i=1}^n e^{\lambda_i(x)}P_iP_i^T\text{ where }P_i\text{ are the columns of }P$

Since A must remain PSD, one potential approach when the approximation does not hold is to regress the Cholesky decomposition of A . However, this would come at the cost of having to regress a quadratic number of coefficients, significantly increasing complexity. [1] Sindhwani, V., Quang, M. H., & Lozano, A. C. (2012, October). Scalable Matrix-valued Kernel Learning for High-dimensional Nonlinear Multivariate Regression and Granger Causality. Retrieved from http://arxiv.org/abs/1210.4792 [2] Kadri, H., Duflos, E., Preux, P., Canu, S., Rakotomamonjy, A., & Audiffren, J. (2016, November). Operator-valued Kernels for Learning from Functional Response Data. Retrieved from http://arxiv.org/abs/1510.08231

Thank you for your assessment. To comment on the weaknesses you pointed out:

Theoretical guarantees: In the new version of the manuscript, we prove a theorem showing that MEVA's loss converges much faster to its limiting value than MEEA. Thus, despite making an error (which we argue is small in model aggregation) compared to MEA by adding extra constraint, this regularizes the aggregation and provides more robust estimates. In our context where $N$ datapoints are used to extrapolate to all $x$, which we identify as low data regime in general, this faster convergence ensures MEVA's superiority in practice.

Impact of bias: Our theoretical results take into account this misspecification and show MEVA's fast convergence even when there is bias. Many existing methods that employ pointwise linear aggregation rely on the heuristic that weights sum to 1. In our work, this idea naturally stems from the unbiasedness assumption, providing both a theoretical justification and a framework for improving the method when this assumption is violated. Appendix D specifically details a bias correction procedure for such cases. In practice, accounting for bias has not significantly influenced performance in our experiments. We attribute this to the minimal assumptions placed on the error covariance, allowing the aggregation process to interpret all errors as variance. Additionally, the distinction between bias and variance may often be more of a theoretical construct than a practical issue. For instance, in the PDE example, all methods are deterministic, meaning that all errors are inherently biased. Despite this, our method performs robustly, even in this misspecified setting. Similar to how linear regression assumes Gaussian errors but remains effective across a wide range of scenarios, our bias assumption serves to justify the aggregation formula and provide a theoretical foundation. However, the method demonstrates robustness in practice and does not necessarily fail when this assumption is not strictly met.

Complexity of examples: We acknowledge the limited inclusion of real-world examples in this work. However, this paper represents the first attempt to extend aggregation techniques beyond simple regression settings, making it inherently exploratory. Our focus was on developing and validating the method in controlled environments. We are actively pursuing applications in scientific computing, particularly in fields such as weather and climate prediction, to tackle real-world problems and demonstrate the broader applicability of our approach.

Benchmarking: Upon further review of existing methods, we recognize that the current state-of-the-art is well established in regression tasks involving surrogate models, where the goal is to recover known analytic functions. We have updated the related works section to reflect this observation. Notably, these methods can be interpreted as specific cases of MEVA, using a predefined regressor for the error. While the Boston Housing dataset aligns broadly with the regression setting, its inclusion in our work primarily serves as a sanity check to validate whether aggregation improves the performance of the models being combined. The core of our experimental contribution lies in the novel application of model aggregation to PDE solvers and operator learning. In these contexts, benchmarking is inherently challenging due to the lack of existing comparable work. These examples better capture the complexities of real-world scientific computing problems, which we aim to address in future applications.

Thank you for your constructive feedback, which helped us improve the paper, and thank you for the support!

(figure reference was updated)

Bias: Our theoretical results take into account this misspecification and show MEVA's fast convergence even when there is bias. The extension to biased models is addressed in Appendix D, and we have added a clearer pointer to this appendix in the presentation of MEVA within the paper. It is worth noting that the distinction between bias and variance may often be more of a theoretical construct than a practical concern. For instance, in the PDE example, all methods are deterministic, meaning that all errors are, by definition, attributed to bias. Despite this, our method still performs well. Similar to how linear regression assumes Gaussian errors but remains effective in a wide range of scenarios, our bias assumption serves to justify the aggregation formula and provide a theoretical foundation. However, the method is robust in both theory and practice and does not necessarily fail when this assumption is not strictly met.

Error correlations: Correlated model errors are treated as a generalization in Appendix B. In the examples we considered, this complexification did not yield better results.

Comparison with other techniques: In all our experiments, we compare our approach with both MEEA (which can be seen as a variant of the popular Nested-Kriging approach) and the simple average of models. Since we do not prescribe a specific regression method for the aggregation process, weighted averages (which are common in the literature) are naturally included as a possible instantiation of our framework. The choice of an appropriate regressor depends on the specific application, allowing flexibility in adapting the method to different scenarios. Furthermore, our work represents the first application of model aggregation in the context of operator learning. As such, there are no existing methods to directly compare against for the PDE example, which highlights the novelty of our approach in this domain.

Practical guidelines: While we describe the class of regressor used (FNO or Matern kernel), practical considerations are left in the GitHub of the paper, where the full detail of the implementation is detailed. We note that we used the classical ML approach to tackle this supervised learning approach, and in particular, the PDE examples required very little tuning.

Thank you for your assessment. To comment on the weaknesses you pointed out:

Impact of bias: Many existing methods that employ pointwise linear aggregation rely on the heuristic that the weights sum to 1. In our work, this natural idea is derived directly from the unbiasedness assumption, providing a theoretical justification and a framework for improvement if this assumption is violated. Specifically, Appendix D outlines a bias correction procedure for such scenarios. Notably, accounting for bias did not significantly affect performance in our experiments. We attribute this to:

• The fast convergence of the loss of MEVA to its limiting value, at rate $\frac{1}{N}$ (compared to the Monte Carlo rate $\frac{1}{\sqrt{N}}$). This discussion is part of the larger theoretical contributions we added to our manuscript, which prove the superiority of MEVA over MEEA.
• The minimal assumptions placed on the error covariance, which enables the aggregation process to interpret and adapt to all errors as variance. This observation reinforces the robustness of the method in practice. These considerations have been clarified in the section describing MEVA.

Complexity of examples: We acknowledge that the operator learning literature now includes very complex examples, particularly in higher-dimensional settings. However, this work represents the first attempt to extend model aggregation beyond simple regression tasks, making it inherently exploratory. Our focus was on establishing the method and demonstrating its potential in controlled settings. We are actively seeking applications in scientific computing, particularly in fields like weather and climate prediction, to address real-world challenges and further validate our approach.

The Data Science Setting: The data science setting is not the primary focus of our method, nor is it where our approach is most warranted. While the question of aggregating existing models is valid and was addressed in the paper as a means to validate our approach, there is little practical reason to use non-trainable models in this context. In data science, our method effectively aggregates models and improves upon their individual performance. However, treating the models as trainable would likely yield even better results, as ensemble learning techniques like boosting have already been thoroughly explored and perform exceptionally well in this domain. It is not our goal to compete with such established methods in this setting. On the other hand, structured problems with specific design constraints, such as those encountered in scientific computing, particularly PDE solving, are better suited for our method. In these settings, non-trainable models or models with distinct drawbacks are common, making aggregation critical. The data science experiment included in the paper serves as a sanity check using a known benchmark, but our primary contribution and the applicability of our methods lie in addressing problems of the PDE kind.

We have made several changes to the manuscript to improve the article.

We have modified the title, abstract, and introduction to better highlight the contributions of our work. In particular, we have revised the related work section and added a clear contributions section. We have proved a theorem to address concerns about the well-foundedness of the method. This theorem justifies the superiority of MEVA over MEEA by showing it has a faster convergence rate, making it more robustly estimated. We identify our setting as a low data regime, making robust estimation key for accuracy. We have reworked the Boston housing dataset example. Its role as a sanity check for the method has been emphasized. We also noticed the small data set size made results highly variable. We now average results over many seeds, and after better hyperparameter choice, our method shows a 2% overall improvement.

Thank you for your detailed responses and the thoughtful revisions to the manuscript. The clarifications on bias robustness and the focus on MEVA’s strengths in PDE applications are appreciated and add depth to the paper. I believe that the scores given earlier accurately reflect the strengths and contributions of the work as it stands.

We have made several changes to the manuscript to improve the article.

We have modified the title, abstract, and introduction to better highlight the contributions of our work. In particular, we have revised the related work section and added a clear contributions section. We have proved a theorem to address concerns about the well-foundedness of the method. This theorem justifies the superiority of MEVA over MEEA by showing it has a faster convergence rate, making it more robustly estimated. We identify our setting as a low data regime, making robust estimation key for accuracy. We have reworked the Boston housing dataset example. Its role as a sanity check for the method has been emphasized. We also noticed the small data set size made results highly variable. We now average results over many seeds, and after better hyperparameter choice, our method shows a 2% overall improvement.

• In equation 13, we use for the Boston housing dataset a Gaussian process with Matern kernel (equation 13), and for the PDE example, the Fourier Neural Operator (lines 427 to 431)

• Numerical solvers for PDEs, such as Laplace and Burger’s equations, often achieve machine precision within reasonable computational timeframes. In contrast, machine learning approaches, including operator learning methods like FNO, DeepONet, and kernel-based methods, typically struggle to surpass 1% accuracy in general settings. This creates an interesting challenge for machine learning methods to aggregate models that operate with orders of magnitude greater precision than direct approximation methods. Moreover, operator learning problems are inherently infinite-dimensional tasks, making them much closer to scientific computing challenges seen in fields like climate modeling or material science. This highlights the relevance of our work in addressing such problems, where the need to combine highly precise models with varying strengths becomes critical.

• Obtaining multiple predictions for a single quantity of interest is undoubtedly more computationally expensive than relying on the best individual model. However, in many practical settings, such as weather prediction, ensembles are widely used because accuracy is prioritized over computational cost. The trade-off between computational cost and accuracy is context-dependent and can only be evaluated in practice, where specific constraints and requirements are known. In our experiments, most examples were completed within a few seconds and demonstrated an order-of-magnitude improvement in error, which justifies the additional computational expense in these cases.

Thank you for your assessment. To comment on the weaknesses you pointed out:

Benchmarking: The method referred to as MEEA can be understood as the empirical variant of nested kriging, a commonly used approach in surrogate modeling and regression tasks. Current model aggregation techniques are largely focused on regression scenarios involving surrogate models, which is why the Boston Housing dataset was included as a sanity check to demonstrate that our method can improve performance in such settings. The primary contribution of our work, however, lies in the seamless, non-intrusive, and versatile applicability of the proposed approach, which enables novel applications of model aggregation to PDE solvers and operator learning. These domains present unique challenges, such as aggregating pre-trained, heterogeneous, or black-box models, and lack established benchmarks or comparable methods. By addressing these complexities, our method provides a robust solution for real-world scientific computing tasks, demonstrating its capability to unify and enhance predictions across diverse legacy models for complex systems. As observed by referee 3, `` The capacity to aggregate a variety of models, including those derived from different methodologies and domains, underscores the versatility of this approach. This has important implications for numerous fields where the integration of multiple predictive models is advantageous, such as climate modeling, where ensemble methods are employed but often rely on simplistic aggregation techniques.’’

Theoretical guarantees : MEVA can be seen as MEEA with the added assumption of unbiased models, which introduces a level of regularization to the aggregation process. This trade-off between regularization and generalization performance is a well-known principle in machine learning. In our revised version, we prove a theorem showing that MEVA converges faster than MEEA, making it a more robust aggregation. As we operate in a low data regime (because N datapoints need to be extrapolated to all x), this faster convergence rate is key in making MEVA superior to MEEA. Importantly, in all our experiments, MEVA consistently outperformed MEEA, providing strong empirical evidence for its superiority.

Remark: simple averaging is also an unbiased estimator (as coefficients sum to 1), thus it is automatically inferior to MVA.

Hyperparameter tuning: Our method is formulated as a supervised learning problem, and standard approaches for hyperparameter selection, such as cross-validation or evaluation on a separate test set, are applicable here. To aid reproducibility, we have also specified the hyperparameter choices used in our experiments in the code, which is publicly available on the paper’s GitHub repository.

We have made several changes to the manuscript to improve the article.

• We have modified the title, abstract, and introduction to better highlight the contributions of our work. In particular, we have revised the related work section and added a clear contributions section.
• We have proved a theorem to address concerns about the well-foundedness of the method. This theorem justifies the superiority of MEVA over MEEA by showing it has a faster convergence rate, making it more robustly estimated. We identify our setting as a low data regime, making robust estimation key for accuracy.
• We have reworked the Boston housing dataset example. Its role as a sanity check for the method has been emphasized. We also noticed the small data set size made results highly variable. We now average results over many seeds, and after better hyperparameter choice, our method shows a 2% overall improvement.

Pathological Examples and Input Distribution: The pathological example indeed relies on specific point placements to illustrate potential pitfalls of naive aggregation methods like MEEA. It is a deliberately extreme scenario designed to highlight how a method may overfit by interpolating data rather than learning a robust aggregation. While such specific patterns may not always occur in practice, the broader takeaway is that methods overly reliant on empirical error minimization (like MEEA) can fail in scenarios with structured or biased sampling. This extreme vulnerability is validated by our Gaussian Process (GP) computations, theorem, and experimental results, which consistently show MEEA performing poorly under less extreme conditions: when input samples are randomly distributed, the likelihood of extreme cases diminishes, but the underlying risks of overfitting or poor generalization remain for MEEA, particularly in data-scarce or heterogeneous settings.

Implications of Independence Assumptions : Allowing $P\neq I$ adds flexibility by enabling MEVA to account for correlated errors between models. However, our experiments show that this flexibility is often unnecessary, as independent errors are typically sufficient for robust performance. In practice, model independence is a reasonable assumption when models are developed using distinct assumptions, methods, or datasets. Even models trained on the same dataset (e.g., linear regression, SVR, boosting) often exhibit uncorrelated generalization errors due to their differing inductive biases and decision boundaries. When errors are significantly correlated, MEVA with $P\neq I$ can theoretically improve performance, but this comes at the cost of increased complexity and potential sensitivity to estimation inaccuracies in $P$. In our tests, the added flexibility did not lead to significant performance gains, making $P=I$ a practical and reliable default for most scenarios.

Thank you for your assessment. To comment on the weaknesses you pointed out:

The data science setting: The proposed method is most effective in scenarios involving heterogeneous models, where traditional ensemble methods like random forests or gradient boosting are not directly applicable. Unlike these approaches, which are designed for homogeneous, trainable models in structured data settings, MEVA operates seamlessly on pre-trained black-box models, requiring no access to their internal training processes or data. This flexibility makes MEVA particularly valuable in structured problems with inherent design constraints, such as scientific computing, multi-physics simulations, and domain adaptation. In these contexts, models often include non-trainable numerical solvers, legacy simulation codes, or specialized machine learning models with distinct trade-offs that cannot be easily retrained or integrated using standard ensemble techniques. By dynamically adapting to heterogeneous model errors and minimizing variance, MEVA ensures robustness to poor-performing models and excels in aggregating specialized or black-box models, offering a versatile solution for challenging aggregation tasks.

While data science tasks provide a useful validation context, they are not the primary domain where MEVA delivers the most practical value. Ensemble methods like boosting and random forests are highly effective in these settings, leveraging homogeneous, trainable models to achieve state-of-the-art performance. MEVA is not intended to compete with these established methods but instead addresses problems where standard approaches fail, such as reconciling the limitations of pre-trained or specialized models in constraint-heavy, structured problems. For instance, MEVA excels in aggregating heterogeneous models trained on different data subsets or domains, combining outputs from models with vastly different structures, or leveraging a mix of machine learning models and numerical solvers. While our experiment on the Boston Housing dataset serves as a sanity check to validate the method on a widely recognized benchmark, MEVA's core contributions and applicability lie in scientific computing and real-world aggregation scenarios where traditional ensemble techniques are insufficient.

Data Assumptions: Our work focuses on developing a robust aggregation method under the assumption that the data is given. While the quality and diversity of data undeniably influence performance, as in any data-driven computation, our primary contribution lies in addressing aggregation challenges rather than introducing novel insights into data requirements beyond standard machine learning principles. That said, MEVA’s regularization through the unbiasedness assumption makes it particularly effective in scarce data regimes, where overfitting or error compounding in traditional methods often leads to degraded performance. Our new theorem proves that MEVA converges faster than MEEA, making it a more robust aggregation.