Minimal Variance Model Aggregation: A principled, non-intrusive, and versatile integration of black box models
Authors: Theo Bourdais
Houman Owhadi
Venue: ICLR 2025 Poster (2025)
Keywords: Model aggregation; Ensemble Learning; Scientific Machine Learning; Gaussian Processes; Neural operators
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We present a general framework to combine existing models in a robust manner. Our method is flexible, places few assumptions on the models aggregated, and performs well in our experiments.

Abstract: Whether deterministic or stochastic, models can be viewed as functions designed to approximate a specific quantity of interest. 
We introduce Minimal Empirical Variance Aggregation (MEVA), a data-driven framework that integrates predictions from various models, enhancing overall accuracy by leveraging the individual strengths of each. This non-intrusive, model-agnostic approach treats the contributing models as black boxes and accommodates outputs from diverse methodologies, including machine learning algorithms and traditional numerical solvers.
We advocate for a point-wise linear aggregation process and consider two methods for optimizing this aggregate: Minimal Error Aggregation (MEA), which minimizes the prediction error, and Minimal Variance Aggregation (MVA), which focuses on reducing variance. We prove a theorem showing that MVA can be more robustly estimated from data than MEA, making MEVA superior to Minimal Empirical Error Aggregation (MEEA). Unlike MEEA, which interpolates target values directly, MEVA formulates aggregation as an error estimation problem, which can be performed
using any backbone learning paradigm. We demonstrate the versatility and effectiveness of our framework across various applications, including data science and partial differential equations, illustrating its ability to significantly enhance both robustness and accuracy.

Published as a conference paper at ICLR 2025
MINIMAL VARIANCE MODEL AGGREGATION: A PRIN- CIPLED, NON-INTRUSIVE, AND VERSATILE INTEGRA- TION OF BLACK BOX MODELS
Théo Bourdais California Institute of Technology www..bouudaisscaaccco.
Houman Owhadi
California Institute of Technology
owhadi@caltech.ed.
ABSTRACT
Whether deterministic or stochastic, models can be viewed as functions designed approximate specific quantity of interest. We introduce Minimal Empiri- cal Variance Aggregation (MEVA), a data-driven framework that integrates pre- dictions from various models, enhancing overall accuracy by leveraging the in- dividual strengths of each. This non-intrusive, model-agnostic approach treats the contributing models as black boxes and accommodates outputs from diverse methodologies, including machine learning algorithms and traditional numerical solvers. We advocate for point-wise linear aggregation process and consider two methods for optimizing this aggregate: Minimal Error Aggregation (MEA), which minimizes the prediction error, and Minimal Variance Aggregation (MVA), which focuses reducing variance. We prove theorem showing that MVA can be more robustly estimated from data than MEA, making MEVA superior to Minimal Empirical Error Aggregation (MEEA). Unlike MEEA, which interpo- lates target values directly, MEVA formulates aggregation as an егrоr estimation problem, which can be performed using any backbone learning paradigm. We demonstrate the versatility and effectiveness of our framework across various ap- plications, including data science and partial differential equations, illustrating its ability to significantly enhance both robustness and accuracy.
INTRODUCTION
Many challenges in scientific computing and Machine Learning (ML) involve approximating func- tions using models. These models can vary significantly, ranging from numerical solvers that in- tegrate differential equations to various ML methods and other approximation techniques. Their performance can differ widely across different benchmarks, (Russakovsky et al. 2015). In some cases, no single model outperforms others across all scenarios; each has its own strengths and weaknesses. For instance, the Intergovernmental Panel Climate Change (IPCC) report on model evaluation (Flato et al.) 2014) highlights numerous evaluation met- rics, showing that different models excel in different areas. When multiple models are available to predict the same quantity of interest, a key challenge arises: how can their predictions be combined most effectively? Ensembling MI. models is well-established practice with many successful appli- cations, such bagging, boosting, and Mixture of Experts (MoE) (Zhang & Ma 2012). Gradient boosting, particular, is widely regarded for its effectiveness, combining simplicity and accuracy in data analysis (Bentéjac et al. 2021). Most ensembling methods focus on training strategies that build and combine multiple models (Zhang & Ma 2012). However, methods for aggregating pre- dictions from existing models are less developed and broadly used, and practitioners often rely simplistic techniques like averaging (Eberi) (001). This highlights the need for a simple and general framework aggregate predictions and better exploit diverse models' unique strengths.
in the Imagenet leaderboards
Corresponding author.
Department of Computing and Mathematical Sciences, Califoria Institute of Technology, Pasadena, CA 911 125, USA.

---

Published as a conference paper at ICLR 2025
Related work Ensemble learning methods combine multiple learnable models improve accu- racy. While these works are influential training paradigms, our approach assumes the models are not trainable, and focuses solely on their aggregation, yielding different task. We emphasize that our work is not comparable to popular ensembles such as Gradient Boosting (Schapire 2003) as it combines black box models.
Extensive research has been made on model averaging, especially in cases where the models are as- sumed to be linear (Zhang & Liu 2023 Hansen 2007 Liang et al. (024). Combination strategies with constant factors include Viana et al. (2009); Pawlikowski & Chorowska (2020). Other studies have explored the pointwise linear aggregation of surrogate models, with various heuristics using weighted averages. Zerpa et al. (2005) uses the variance predicted by Gaussian processes, Lee & Choi (2014): Liang et al. (202); Ye et al. (2020) try to predict the cross-validation error at each point and use it to combine models while Liu et al. (2016) trains aggregation to minimize the Mean Squared Error of the aggregate. Crucially, these methods are heuristics, which do not justify why the aggregation should be convex combination, why it should be weighted using cross-validatee Secondly, they do not formulate the aggregation as a learning problem and instead use specific ap- proaches, such K-nearest neighbours interpolation through weighted averages. This limits the scope of the aggregation, in particular, these studies only aggregate one-dimensional regression, and prevent the methods from benefiting from advances in other areas of Machine Learning.
A growing trend in scientific computing involves applying statistical inference and machine learning (ML) methods (Raissi al.) 2019 Batlle al. 2024; Chen al. 2021] Bourdais et al. (0024). These problems are often framed as operator learning tasks, which can be seen as
re- gression problems. In this work, we adopt statistical/ML. framework to unify and aggregate diverse methods, even those not inherently inference or ML-based. This approach is particularly well-suited to structured problems with inherent design constraints, where heterogeneous models naturally arise. Such contexts often involve a combination of non-trainable numerical solvers, legacy simulation codes, and specialized ML models with distinct trade-offs, which are challenging to integrate or retrain using standard ensemble techniques. Our work focuses addressing these challenges by reconciling and leveraging the strengths of pre-trained or specialized models in constraint-heave heavy, structured problems.
Contributions We introduce Minimal Empirical Variance Aggregation (MEVA), a versatile method for aggregating predictions from diverse models, treating them as black-box onput-output functions (deterministic or stochastic). By focusing on variance minimization, MEVA provides a robust, non-intrusive framework, enhancing reliability in data-scarce or noisy environments where traditional methods often fall short. While the approach is broadly applicable, its primary strength lies in addressing challenges in scientific computing and other structured domains. Key contribu- tions include:
Contrary to previous heuristics, we present a probabilistic model for aggregation, grounded in clear assumptions that justify using convex combina- tions and approximations of cross-validation error (Sections 2 3 Future works may build upon this model and adapt it to constraints specific the application hand. We also for- mulate MEVA as a learning problem (sec. (.. making an extension to other settings such as graphs or computer vision simple using popular learning techniques of these fields.
Justification and formulation:
Ruling out MEEA We assess Minimal Empirical Error Aggregation (MEEA), which min- imizes aggregate empirical error directly, and identify its limitations via pathological example (Section 2.2). MEVA, introduced in Section 3 addresses these challenges and consistently outperforms MEEA in experiments. We rigorously prove that MEVA con- verges faster than MEEA, making it more robust in scarce data regime.
Extension to operator learning Our approach extends model aggregation beyond machine learning to include scientific computing, enabling the integration of diverse methods for es- timating shared target quantities. We demonstrate this by expanding from one-dimensional problems to operator learning using established techniques (Section 4.2). To our knowl- edge, this is the first work to aggregate outputs across machine learning and classical PDE- solving methods.
The remainder of the paper is structured as follows: MEEA'S limitations. Section 3 introduces the probabilistic foundation of MEVA, theorem proving
Section 2 discusses pointwise aggregation and
2

---

Published as a conference paper at ICLR 2025
its superiority, and its computation. Section 4 validates MEVA with two experiments: sanity check in data science and examples from Scientific Machine Learning (SciML) involving PDE solving.
THE MINIMAL ERROR AGGREGATION
2.1 BEST LINEAR AGGREGATE WHEN MODEL CORRELATIONS ARE KNOWN
An effective way to derive aggregation method is to define it as optimal with respect to specific loss function. If the target function (x) and the models M2 are viewed as random variables on the same probability space, a commonly used loss function is the Mean Square Error (MSE), In this context, an aggregated model MA (1) can defined as any measurable function of the models that minimizes the expected squared error:
(x = argmin E ((M (, (x)x)) TE (x)) (M), f measarable
(I)
where the second equality follows from the Le characterization of conditional expectation. While computing the conditional expectation can intractable (Ackerman et al. 2017). this com- putation reduces to solving a linear system in the Gaussian case. Specifically, if the vector (x)) Mz(z). ...).................. (x)) Gaussian with mean 0 (see remark below), ) then the conditional ex- pectation in the equation above is a linear combination of the models' outputs. Rullière et al. (2017) uses this assumption to aggregate Gaussian processes. Bajgiran Owhadil (202T) demonstrates that any consistent and rational aggregation of models must be weighted average that is point-wise lin- car. Finally, many ensembling methods (random forest (Breiman 2001). MoE (Yuksel et al. 2012)) also linear combination of models. Motivated by this, restrict aggregation approach to a pointwise linear combination the models:
Ma(x) = o Mxx)
With this simplification, aggregation can be achieved by determining the Minimal Error Aggregation (MEA) a" (x:) such that:
Ox (x) argmin E [(Y (())) = (M(+)
(M)
(3)
oER"
By definition, the MEA is the best possible point-wise linear aggregation. This equation reveals that the MEA requires knowledge of both the correlation matrix of the models, E [M (=)) as well as the correlation between the models and the target E (x))) We demonstrate in appendix B.2 that the aggregated model can achieve exceptional performance when these correlations are perfectly known.
Remark In the general case with non-zero mean, the aggregation of a Gaussian vector would be affine. Since we do not assume anything on the models, affine coefficients can be recovered by adding extra constant model, i.e., by a linear aggregation of M (z) (M)),, 1). In section 2.2 we will see that this can lead catastrophic failure for empirical error minimization. Appendix 1 shows that empirical variance minimization provides more robust approach for estimating bias.
2.2 DATA-DRIVEN AGGREGATION
We will now consider the situation where the quantities E [Y and E (x) re- quired for the best linear aggregate equation B may be ill-defined (e.g., because the underlying models are not stochastic) or difficult to estimate. In this situation, a natural alternative approach is
estimate the aggregation coefficients a" appearing in qquation[dirretty from data by minimizer of (possibly regularized) empirical version of the loss (3), as seen in Liu et al. 2016):
āg = argmin GEH
Σ y (V-aa(X)TM((>))) +|lal1 1=1
(4)
where the (X*, (X3,Yi:= ; Y(Xi) are N data points; H. is set of functions used for the approximation; 2 O and ||a / is a regularizing norm. We call this aggregation Minimal Empirical Error Aggre- gation (MEEA). While this transition from expected loss to empirical loss may seem rell-founded, it can introduce a significant loss information, which yields to the aggregation completely failing, showcased in the following sections.
3

---

Published as a conference paper at ICLR 2025
Pathological example: A dubious trend In this example, we introduce basic aggregation with linear coefficients, where the MEEA fails. This, along with appendix [[ G gives intuition as to why MEEA is not the appropriate strategy. Figure represents Figure[ pick Y, a good model Mg and bad models MB s.t.:
the results of this experiment. We
(()) = 2x + c05()) Mc(x) MB(²) =
((x) +e()) where e(x) 0.2)
(5)
We also choose a linear aggregation, which means that we will find aG. 4g by. by s.t. 6(x) = (agæ be, ag bø) RR? Finally, pick some trick data X (0.8 1/38-2/3-0 4, 0.0. - 2/3 3,.., 0.8, 0.8, -0.8 + 2/3 1, -0.8 + 2/3 4) so that the Y(X*) Y form a line. We may see that Me has lower error than Mn each X by evaluating the models at these data points. Thus, we expect an aggregation method to use the good model primarily in the aggregation. However, the opposite behaviour is observed. Minimizing the loss with X on the good model is O, so Mc is completely ignored, which is countcrintuitive. The aggregate uses the bad model, MB, to directly estimate Y by linear regression. This shows that the MEEA simply employs models as features to interpolate the data points instead of weighting each model according to its estimated accuracy.
we obtain that the coefficient
larqet
good madel 14200
Indiroudd Mels
X
Appropett
Figure 1: Pathological example (sec. 2.2
AGGREGATION THROUGH ERROR ESTIMATION: THE MINIMAL VARIANCE AGGREGATE (MVA)
3.1 MODELLING THE ERROR
By placing no assumption on our model, MEEA is vulnerable to pathological examples, and in experiments does not generalize well. In particular, there is nothing representing the fact that models are approximations of the target. To encode this, we add the assumption that the models, represented in the vector M, are unbiased estimators of Y. Thus, the errors of each model interpreted as follows:
Y(r)
{
[[Z(x)] (Cov[Z(z) A(x)
M(x)
+ Z(x) where
(6)
(())
Under this assumption, any linear combination with weights (2) E (x(x) = 1 is unbiased. Thus, we define the Minimal Variance Aggregate (MVA) as the Best Linear Unbiased Estimate (BLUE) of Y()):
Ma = (v = 17 ATA(x)-11
where 1
(1, 1)T the all ones vector.
Note that av
argmin,r...
E (Y())) a (x))) regardless of whether E [2(x) 0. In previous heuristics
LLe & Choi
2014
Liang al. 2023 Liu et al. 2016 Viana al. 2009 Σ (x(x)) is
4

---

Published as a conference paper at ICLR 2025
assumed, while we justify it as the set of unbiased linear aggregations. This leads to natural generalization for biased models in appendix [F] while in practice bias correction did not improve results. The following section proves that despite MVA being inferior to MEA in MSE, MVA can be more robustly estimated, making MEVA better performing than MEEA.
3.2 A THEOREM ON THE SUPERIORITY OF MEVA OVER MEEA
To investigate the theoretical properties of MEVA and MEEA, we fix and focus on estimating a and (vy Details and proofs are given in appendix A] Let Y be random real-valued variable and the random vector of models M E R Suppose we have i.i.d observations (M1, y (M YN), collected in the matrices M E R and Y E I. Define M = Y + <Z, where Z is the scaled error, with chosen so that |Zz ;= has Frobenius norm : and b == E [Y 2]. We find that satisfy:
the MEA a defined in equation and the MVA ay (equation
( A
a' = dov + (1 AaRR where
(8)
an =
with 1-A-11, Dt 1- and b" b assumed to be non quantities. Note that À E [0, 1], making a" convex combination of ay and aR- For the empirical variants of these aggregations, we use = 2M 1)T) y17)). Using the central limit theorem on Å 1--4 A VN 5A, we find that the scaled error 5A converges to Gaussian Matrix of finite covariance. Similarly, E [Y²] and b are estimated with estimators converging at vN rates, with 5V and ōb their scaled enror. Write åg == (MMM))-M'') for the empirical MEA, áv (A-11)/(1" A-' 1) the empirical MVA, and L(a) E [(Y M)]] the Mean Squared Error. We prove the following theorem:
I
Theorem 1. Assume 8, u, 1 + O, then there exist two sequences of random variables, K& and K each of which converges in distribution to a distinct finite random variable, such that
1 ((a²) + KV and L(āg) = L(aa)+++))))))))))))) + N
1
[(āv)
VN K
X
MEVA incurs a constant error relative to MEEA, given by 3. However, MEVA is less sensitive to perturbations, causing its error to converge more rapidly its limiting value than MEEA. Thus, in the model aggregation setting where À is close to O the smaller difference in limiting loss values is outweighed by MEVA's faster convergence making it outperform MEEA
While this theorem applies to a fixed x and N observations, we are interested in the case where these N data points are used to extrapolate across all x. such cases, model aggregation inherently operates in the scarce data regime where À the theorem, further emphasizing the practical advantages of MEVA.
large. This amplifies the behavior described in
3.3 COMPUTING THE AGGREGATION: THE MINIMAL EMPIRICAL VARIANCE AGGREGATE (MEVA)
To compute our empirical aggregation dy (z), must compute the covariance matrix of the error for all x E X' from the data. This poses several challenges when A(x) is unknown. First, we must approximate a symmetric matrix everywhere gregate n models. Moreover, we must ensure this matrix is always positive-definite To simplify the approximation and to easily enforee the positive-definte condition, we suppose that for all æ, A(x) can be diagonalized using fixed eigenbasis P. Then, its positive eigenvalues are approximated by (x), where the A are functions to be regressed from the data. Thus, we approximate A(z) with Á(x) pT Diag(e>1(x) (x). This formulation has the advantage of symmetrizing the co- variance and precision matrix estimation, as the big-eigenvslues of both matrices are opposite. This is desirable because while we are interested in the precision matrix, the covariance matrix is easier to estimate. As a further simplifying assumption, we will take P = In, which is equivalent to as- suming that models have independent error.. Appendix presents the general case P. y I,, which
the domain, needing n(n (1)22 coefficients ag-
See appendix A for more details
5

---

Published as a conference paper at ICLR 2025
in our experiments did not improve performance. Using this assumption on A(x) in equation 7 get:
A(z) (v(x)
Diag(e>1(z), ;..,
(d)( (()()))) Softmax(-M()))
where where
A are regressors
(10)
Softmar(-Xx))) = Σ -
Notice that using our model of the error (6., have re-discovered the popular softmax, used in particular in Mfxture-of-Experts (MoE). All that remains for our aggregation to be complete is to
train the regressors A₃. To do so, we will rely on covariance matrix estimation. Let Z!,., Z-V be i..dd samples of the random variable Z E R" such that E [2] = 0. An unbiased estimate of CoulZ) [Z]
N Σ
N
1
argmin Σ Cov = N Z'(ZT==
(11)
MEReX" i=1
where ||-I|= is the Frobenius norm. If we further assume that Cov Z must be diagonal (as we do with A(z)). then its coefficients are such that (Cov Z) = argmin.cs Ex= ||c (Z) We will now regularize this argmin formulation to approximate the model errors covariance matrix for all values of x. Writing e² for the vector with entries defined as the sample errors (e')x = Mx - Y(X), we identify the functions X;
n ΣΣ -c))) alli1
A argmin ICH
(12)
Here, denote the entries of the vector-valued function 1 by 4, and H represents the space of vector-valued functions selected to approximate the eigeeiins...es. This space could be a RKHS associated with particular kernel, a neural network with specific architecture, or any other normed function space. The notation |I-1| refers to the regularization norm intrinsic to that space, such as the RKHS norm or the L2 no-mo of the weights and biases in a neural network.
EXPERIMENTS
4.1 SANITY CHECK: AGGREGATION ON THE BOSTON HOUSING DATASET
The Boston housing dataset (Harrison & Rubinfeld 1978) is popular benchmarking dataset for regression. contains N = 506 samples of variables. We use this dataset as a sanity check that our approach can improve upon the predictions of the aggregated models, while our focus is on the PDE examples below. We split this dataset into train, validation, and test sets, and average results over many splits. To evaluate our method, we train several standard ML methods set (e.g. linear regression, gradient boosting, etc.). We then train our aggregation on the validation set before testing on the test set. To get fair comparison, the train and validation sets combined so that they see the same amount of data as the aggregate. We first implement our aggregation as described in section 3.3] using kernel method with kernel K. This is equivalent to minimizing equation [2] over A E HH where H2 is the RKHS defined by the kernel K. After experimenting with many kernels, the best we found uses the Matérn kernel
the train
also trained the ML methods on
Ma (u, U. p) = u-
exp /3luuool and is defined as P
<(x, 3) KM ( (M (1),
(13)
with adequate choice of p. This formulation uses the prediction (xx) to output the aggregation coefficients (z)).. The results of this experiment are shown in figure 2, In this situation, the mean performs well. The aggregate, denoted MEVA with kernel, performs better than all individual models being aggregated, their average and any ML model when fed the same data. This is a success for the MEVA, which difficult to obtain with MEEA. The aggregation improves upon the performance of the best ag- gregated model, Gradient Boosting, by Timain = overall by Tanl 2%.
A MSE S (Aggrggae)
7%, and the best model
(Gradiit Boosting
"See details in the paper's repository

---

Published as a conference paper at ICLR 2025
Model Performances on test set
Unbear Nigestaoon
Decibenn Sier
Flvest
Dioceet Beosting
Aleen montli
Asgropenc
MSSI with cNN
Aneigtioos
Model
Figure 2: Comparisons of the performance of the different models on the Boston housing dataset. Red bars are the performance of models trained using the training set and used in the aggregation. Blue bars show models trained on the training and validation set (Irainvval) to get fair comparison with the aggregations. The aggregation does not use the models trained on truin+val.
We conduct a second experiment to demonstrate the effectiveness of our proposed method compared to Minimal Error Aggregation, defined in equation 4 In this experiment, we use small, fully connected Deep Neural Network (DNN and train it using two different loss functions: the MEA described in section 2 and Minimal Variance Aggregate (MVA) using an end-to-end loss especially useful for neural networks, described in appendix [ The results, presented figure highlight the impact of the chosen loss function. Since the only difference between the two approaches is the loss function, the results clearly illustrate the advantages of variance minimization. Specifically, the minimal error aggregation (denoted MEEA with DNN in the figure) performs significantly worse compared to our minimal variance aggregation (denoted MEVA with DNN in the figure). It also performs worse than the two best models aggregated.
4.2 AGGREGATION OF PDE SOLVERS
This section introduces the approach for combining Partial Differential Equation (PDE) solvers through model aggregation. By incorporating classical PDE solvers into the aggregation frame- work, we lower the accuracy threshold necessary for the aggregated solution to outperform individ- ual models. Notably, this threshold is significantly below the typical 1 % relative error benchmark targeted by most ML-based PDE solvers. We first formulate the aggregation problem in the operator learning context, followed by two experimental case studies demonstrating the effectiveness of this approach.
4.2.1 OPERATOR LEARNING IN THE AGGREGATION CONTEXT
Consider the Laplace equation with zero Dirichlet boundary conditions on the domain =[0, as illustrative example of PDE:
Aut (w) f(w) for w E Ω W E Ue
(14)
(l(w)
for
where E F C. (2) and wl y C (22) n (T) is the solution. The solution opera- tor maps the source term the solution ut. In this context, the solution operator is defined as F + ut E y. This operator is often approximated using a PDE solver, which can be based on methods like finite element analysis (Bathe 2008) or spectral methods (Fornberg & Sloan 1994). Alternatively, the operator can be learned directly from pairs of input/output solutions, a supervised learning task known as operator learning (Kovachki al.) 2024). Since PDE solvers and machine learning methods approximate operators, we can naturally frame our aggregation task as an
f ::f
See details in the paper's repository

---

Published as a conference paper at ICLR 2025
aggregation of multiple operators (or PDE solvers), denoted Mi,..., Mn. The aggregated operator MA is then expressed using a, an operator that maps F to L (2, ()), so that:
Ma(f) = Σ {F
(15)
Remark: In the operator context, one may consider many aggregations that are pointwise lin- car with respect the models, for instance, the aggregation of Fourier coefficients, MA(f)
T (E (s(f)T(M(()))), where T represents the Fourier transform.
4.2.2 THE OPERATOR AGGREGATION LOSS
To get the loss for our problem, we will first state the general loss in the operator learning setting and then modify it for increased performance. Since the aggregation loss (12) is defined for real-valuce objectives with an abstract input x: E X, we must specify X' in the context of operator learning. We use X' = {ww, f) E Ω x F), which means
(w, f' is choice of point in space and source term. We define ((x) == (w:) = (w) E R. Given grid (w (w",...w M1) E Ω and a set of input functions (f (f',.,, f we get the N1 X N2 datapoints X == (w", for our aggregation. These pairs of datapoints (X+,,y"s
= ((F)(w)) can then be used in equation 12
When solving PDEs, the focus is often on the order of magnitude of the error, i.e., its logarithm. Therefore, we may adapt the loss function (12) to account for this. Specifically, observe that if the regularization parameter a is small, then functions obtained from regressors that can interpolate arbitrary data (e.g., GPs with universal kernels) may result in a near-zero loss and eAx(X* is
(Y(X)))- -)(x)))) = & as - O. Thus, we can linearize the exponential around log[ek (J]
see that:
(e - c (((Xii)--ooge
(16)
This lincarization shows that samples where the models perform well are down-weighted, and the aggregation would foeus on instances where they perform poorly. Instead, we the following more sensitive loss (17) where the down-weighting is removed, which, in practice, gives good results on these problems:
A argnin
(x(+))) - log ((Y(Xtt)) Mc((((1)3)) + aall =1j=1k=]
(17)
ICH
We also recall the process for obtaining a and the aggregate MA from X is given in equations TO and [5] In both examples, we use the Fourier Neural Operator (FNO) (Li et al.] (022) to learn
operator L. We also explored kemnel-based approach for operator learning,
detailed in Batlle
ct al. (2024). Although operator learning using a Matérn kernel yielded good results, FNO offered better performance with minimal hyperparameter tuning for the examples discussed in this section.
4.3 THE LAPLACE EQUATION
Our first example is the Laplace equation 14.. We use six simple solvers as models: a finite dif- ference method, two finite differences with an inhomogeneous grid (one denser for x < 0.4, the other for x > 0.4), a spectral method, and a GP solver. Our data is generated by sampling") random parameters fman= Pos P1, R for functions of the form
(- J 1/ p1
T
14 1)
u(xx y) = sin(TT) sin(Tr) sin fmax exp
14
R
x
140
(18)
Six hundred pairs of functions (ft := Au, u.) are sampled, of which 500 are used for training
and 100 for testing. Figure 6 shows an example of such a pair. These functions and the models are evaluated on grid of 100 x 100. The aggregation is performed by FNO, which takes input and the models' outputs and outputs the aggregation coefficients' logits X. When training FNO with parameters E C RN and loss (7). we find the minimizer of [7 with H. {f
FNO M1) Mf)) for 0 E Q) with [FNO ||0||2-
*See details in the paper's repository
8

---

Published as a conference paper at ICLR 2025
30"
fair_ssym__
falm coym right falm
10
mian zquired cror
spectral
30--
30
10 16
10
sumpll
(a) Laplace equation
(b) Burger's equation
Figure 3: log MSE of the different methods for: (3a) the Laplace equation; (3b) Burger's equation.
Samples are sorted by the error of the aggregate
Apart from the GP solver, all other methods perform similarly, with the finite difference method being the best overall. The GP solver, implemented with fixed length scale, has two regimes. When the function u, is not evolving too fast (i.e. when fmaz is not too large), the kernel is well specified, and the GP solver obtains the best performance. On other occasions where fmes larger, the method fails and has a substantial error. Thus, the aggregation challenge here is to use the GP when performing well but ignore when fails. The result of this experiment is shown figure 3a Because of the significant error in the failure case, we observe that the mean is poor aggrega- tion that performs worse than all models aggregated. Our aggregate consistently among the best performers on each sample and the best-performing model for many samples. On average, ag- gregate performs one order of magnitude better than the models aggregated, as shown in table An example of model outputs, errors, and final aggregation is shown in figure 7 As comparison, also train FNO using a direct MSE loss of the aggregation, akin to (4). We observe that this method fails the aggregation converges best-performing model. We believe that this convergence to a fixed value, without combining the different models get better
constant, i.e., a(S,w) (1,0, 0, 0, V J, w. picking the
observed in method, is due the absence of the log term that we added in the sharp loss (17). As we observed in section
without this log is difficult for a loss to differentiate between 10 and 10 error. Furthermore, while learning the log-error instead of the error is simple modification from [2 to (17), it is difficult to see where this log would be added in loss like
4.4 BURGER'S 'S EQUATION
We now turn to a non-linear PDE, Burger's equation with periodic boundary condition and finite time interval:
O,u(t, +z
v(, u(0, x)) = f(x) for (lt,0) u(t, for
for I, [0, 1 x E 0, E 1]
(19)
We choose 3.10 and pick initial conditions ft samples from the Gaussian process (0, K Kexpsin²) where Kexp sin= y) exp [=ln=( and = 1.5. With this set of param- eters, the initial conditions are infinitely differentiable periodic functions that often, but not always, form shock within the time frame studied (t E 1).. We use seven different solvers with the same fixed discretization grid. While the first five methods implement the viscid version of Burger's equation directly, for instance, using explicit scheme uses a flux limiter. The last one solves the inviscid Burger's equation. With this variety of solvers, have methods that are accurate in smooth cases when there is large shock, intermediate method that may be more robust to shocks, and a method that will never diverge even the presence of shocks, but will be less accurate in general because it solves slightly different equation. A good
finite volume method, the penultimate

---

Published as a conference paper at ICLR 2025
aggregation must be able recognize artefacts and divergences that may arise around shocks but use the more precise, although more brittle, solvers in cases where they are accurate. To perform the aggregation, similarly the previous section, use FNO, to which we feed the models' prediction as a stack of 2-dimensional arrays (one for time and one for space). Specifically, we minimize loss (J)) for some 0.
17 with X(f) F (N)
The results are displayed in figure 3b Our method successfully leverages the predictions of the different models to give a robust and accurate prediction. The behaviors of our aggregation can be separated into a few cases. On the left of the graph, we find easy initial conditions that all models can approximate correctly. We may notice the Riemann method, which solves the inviscid Burger's equation, has lower accuracy. Thus, on this easy half of the graph, the aggregation mainly relies on the TVD solver, the most precise, to achieve excellent accuracy. In the right half, we find the more complex cases where the aggregation has an accuracy closer to the more robust Riemann solver. In this half, there are cases where the TVD method fails and diverges and is avoided by our aggrega- tion. The process can identify each solver's strengths and weaknesses and obtain better average error, as shown in table 2 An example of model outputs, errors, a, and final aggregation is shown in figure 8
Geometric Mean of MSE (log1g)
Method Aggregate
Geometric Mean of MSE (log10)
Method
-3.106
Aggregate FDM
-6.282
Riemann
-2.734
-5.523
TVD
-2.568
Spectral
-4.988
FVM
-1.228
GP
-4.739
Spectral
-0.625
FDM (Asymmetric Right)
-4.685
Implicit
-0.488
FDM (Asymmetric Left)
-4.699
Explicit
-0.455
Lax-Wendroff
-0.455
Table 1: Result of Laplace equation experi-
ment (sec. 4.3
Table Result of Burger's equation experi- ment (sec. 4.4)
DISCUSSION
We introduced a general data-driven framework for aggregating models with minimal assumptions. As demonstrated in the pathological example Section 2.2 and the theorem 3.2, directly training the aggregate model is unstable and may result in aggregator that fails improve upon the performance of individual models. To address this, we introduced a simple assumption of unbiased models in (6) and reformulated the aggregation task as variance minimization problem (Section 3). In one data science problem (..1tion[4.]) and two PDE-solving examples (Sections 4.3 and loss function consistently outperformed the minimal error aggregate and led to aggregation that surpassed the performance of the individual models. These latter two examples also illustrate the flexibility of our method, which can aggregate not only machine learning models but also other types of models, such as PDE solvers, in a data-driven manner. Notably, we used FNO to get an aggregate performance much higher than it could get when trying to solve the PDE directly.
Our method relies on the availability of unseen data to estimate the error of different models. There- fore, as seen the Boston Housing example (section [.J)). splitting the data further to create an unseen validation set in scenarios where data is limited may significantly reduce the amount of training data available for the models, potentially impairing overall performance. Additionally, in challenging aggregation tasks, there is guarantee that the aggregate will consistently outperform the individual models or simpler methods, such as averaging, despite the higher computational cost.
MEVA's flexibility makes it well-suited for integrating diverse methods, particularly in structured contexts like scientific computing, where non-intrusive solutions are often necessary. By requiring minimal assumptions about the models or data, MEVA can effectively aggregate legacy solvers that are difficult to modify, physics-based models built on differing assumptions, or models with incompatible software frameworks.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGEMENTS
TB and HO acknowledge support from the Air Force Office of Scientific Research under MURI awards number FA9550-20-1-0333 1-0358 (Machine Learning and Physics-Based Modeling and Simula- tion), FOA-AFRL-AFOSR-DDU --00--2223-000 (Mathematics of Digital Twins) and by the Department of Energy under award number DE-SC0023163 (SEA-CROGS: Scalable, Efficient and Accelerated Causal Reasoning Operators, Graphs and Spikes for Earth and Embedded Systems). Additionally HO acknowledge support from the DoD Vannevar Bush Faculty Fellowship Program.
REFERENCES
Nathanael I. Ackerman, Cameron E Freer, and Daniel M Roy. On computability and disintegration. Mathematical Structures Computer Science, 2788):1277-131, 2017.
Mauricio A Alvarez, Lorenzo Rosasco, Neil D Lawrence, et al. Kernels for vector-valued functions: A review. Foundations and Trends@ in Machine Learning, 4(3):195-266) 2012.
Hamed Hamze Bajgiran and Houman Owhadi. Aggregation of Pareto optimal models, December 2021. DRLhttpa / /arxiv.vvv..////// hwt. http///as//111// 2212. 04161 12ivv2112.041111 [econ, math, stat].
Klaus-Jürgen Bathe. Finite Element Method, pp. 1-12. John Wiley & Sons, Ltd, 2008. 159. URL https: ecse151 159
ISBN 9780470050118 18. doi:
//onlinglibrary wi Tey. http://////bbooaa//////////////0
Pau Batlle, Matthicu Darey, Bamdad Hosseini, and Houman Owhadi. Kernel methods are competi- tive for operator learning. Journal of Computational Physics, 496:1 12549, 2024.
Candice Bentéjac, Anna Csörgö, and Gonzalo Martinez-Muñoz. A comparative analysis of gradient boosting algorithms. Artificial Intelligence Review, 143))1993711977, March 2021. ISSN 1563-7462. doi: 10. 610462-020-09896-
Patrick Billingsley. Convergence of probability measures. Wiley series in probability and statistics. Probability and statistics section. Wiley, New York, 2nd ed edition, 1999. ISBN 978-0-471- 19745-4.
Théo Bourdais, Pau Batlle, Xianjin Yang, Ricardo Baptista, Nicolas Rouquette, and Houman Owhadi. Codiscovering graphical structure and functional relationships within data: A gaussian process framework for connecting the dots Proceedings of the National Academy of Sciences, 121(32):e2003449121, 2024.
Leo Breiman. Random Forests. Machine Learning, 45(1):5-32, October 2001. ISSN 1573-0565. doi: 10.102//A110009333444444 URL htttps://do..org//// 1023/A:
1010933404324
Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M. Stuart. Solving and Learning Nonlinear PDEs with Gaussian Processes, August 2021. URL http://rxvv.rr////
2103. . 12959
arXiv:2103.1 12959 [cs, math, stat].
Elizabeth F. Ebert. Ability of a Poor Man's Ensemble to Predict the Probability and Dis- tribution of Precipitation. Monthly Weather Review, 129400)2411-2480, October 2001.
ISSN 1520-0493, 0027-0644. doi: 10.
CO,2. URL //oo://////ra.s..../ amet comeetocc.rr//iveeeeeee
journal httphaaa////eeee///////a
1520-0493 2001 129 2461 aoapms 2. 2.0.co. 2. xml
Publisher: American Meteo-
rological Society Section: Monthly Weather Review.
Gregory Flato, Jochem Marotzke, Babatunde Abiodun, Pascale Braconnot, Sin Chan Chou, William Collins, Peter Cox, Fatima Driouech, Seita Emori, Veronika Eyring, et al. Evaluation of climate models. In Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, pp. 741-866. Cambridge University Press, 2014.
11

---

Published as a conference paper at ICLR 2025
Bengt Fornberg and David M. Sloan. A review of pseudospectral methods for solving partial differ- ential equations. Acta Numerica, 3:203-267, 1994. doi: 10.1017/5062442290002444
Bruce E. Hansen. Leats Squares Model Averaging. Econometrica, 75(4): 2007. ISSN 0012-9682. URL hctps: /www://wwww...... jstor. oettor..rg//stabee/400000 Econometric Societyl.
Publisher: [Wiley,
David Harrison and Daniel L Rubinfeld. Hedonic housing prices and the demand for clean air. Journal of Environmental Economics and Management, 5(1):81-102, March 1978. ISSN 0095- 0696. doi: 10. 10.101660005-00696659999006666 SC /ccience/arriicle///////00533333
URL/httpst http://wwwwiieeeeeeeeeeeeeee com/
Nikola B. Kovachki, Samuel Lanthaler, and Andrew M. Stuart. Operator Learning: Algo- rithms and Analysis, February 2024. URL
arXiv:2402.15715 [es, math].
Yongbin Lee and Dong-Hoon Choi. Pointwise ensemble of meta-models using nearest points cross-validation, Structural and Multidiccipllnnary Optimization, 50(3):383-3399) September 2014.
ISSN 1615-1488, doi: 10. 10.1007/>00158-014 1067-1. DRLL:taps///ai..or//II...
5 1:0C588014-1066
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Fourier Neural Operator for Parametric Partial Differential Equations, May 2021. URL http:://aaat wwtp://woggggbbbbbeeeeooeee 08895 arXiv:2010.08895
[cs, math].
Pengwei Liang, Shuai Zhang, Yong Pang. Jianji Li, and Xueguan Song. A Pointwise-Optimal Ensemble Surrogate Models. Journal of Mechanical Design, 145(111705), August 2023, ISSN 1050-0472. doi: 406/9999
www.ii.orrraao. 111001.111111 10611///...669999 4062979
Zhongqi Liang. Cai, Suojin Wang, and Qihua Wang. K-fold cross-validation based fre- quentist model averaging for linear models with nonignorable missing responses, Octo- ber 2024. URL wttp::::/wwweeeerccsuuaaeeeeeee https//uuuaaaeeeee/rrrrree///sa/
ISSN: 2693-5015.
Haitao Liu, Shengli Xu, Xiaofang Wang, Jigang Meng, and Shuhua Yang. Optimal weighted pointwise ensemble of radial basis functions with different basis functions. AIAA Journal, 54 (10):3117-3133, 2016. doi: 22511//1J054666 URL https: org/ 00./10..251///1
J054664
Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondřej Čertík, Sergey B. Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, Thilina Rath- nayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, Štěpán Roučka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Scopatz. Sympy: symbolic computing in python. PeerJ Computer Science, 3 Jamuary 2017. ISSN 2376-5992. doi: 10.7717/peerj--oo ctapppeer--ssseeeeee cs.103.. 103. httpttttt///ooooo.
103
Houman Owhadi and Clint Scovel. Operator-Adapted Wavelets, Fast Solvers, and Numerical Ho- mogenization: From Game Theoretic Approach to Numerical Approximation and Algorithm Design, volume 35. Cambridge University Press, 2019.
Maciej Pawlikowski and Agata Chorowska. Weighted ensemble of statistical models. International Journal of Forecasting, 36(1):93-97, January 2020. ISSN 0169-2070. doi: 10. 10.10ea//.jfoorecass 2019.03.019. URL https: htt::://wwwww... www. tttp://wwwwiiiiiieeettttteeeeeee com http://rrt....oceeeeeeerrriaaaaa
20169207019301190 501692070192011
S016
M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equa- tions. Journal of Computational Physics, 378:686-700 707, February 2019. ISSN 00219991. doi:
20.1016/j.cc.2218...o 10.045. URL https: htti:////nnnnnnhbb/ eseevier com/ zetrieve/ pi 001/S00219000 0900002299999919999
12

---

Published as a conference paper at ICLR 2025
Didier Rullière, Nicolas Durrande, François Bachoc, and Clément Chevalier. Nested Kriging predic- tions for datasets with large number of observations, July 2017. URL
abs/ /207/160...m 05432 1XXiv:107.00444 [stat].
Olga Russakovsky, Jia Deng. Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhi- heng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li Fei-Fei. ImageNet Large Scale Visual Recognition Challenge. International Jour- nal of Computer Vision, 115(3):211-25)) December 2015. ISSN 0920-5691, 1573-1405. doi: 10. 10.1007711
URL http://ling spri wtt.iiinrrr....ooeoooommmmmoomme c0..m...mm..m.mmm....mm...mm.... 10mmm/10.10077/
s11263- 015-0816-1
Robert E. Schapire. The Boosting Approach Machine Learning: An Overview. In P. Bickel, P. Diggle, S. Fienberg, K. Krickeberg, Olkin, N. Wermuth, S. Zeger, David D. Denison, Mark H. Hansen, Christopher C. Holmes, Bani Mallick. and Bin Yu (eds.). Nonlinear Es- timation and Classification, volume 171, pp. 149-171. Springer New York, New York, NY, doi: 000000999770038--1157777 9. com wttp://......coo//....0.///63331
2003. ISBN 978-0-387-95444 -9471-4 978-0-387-21579-2 URL http:://img sarkssprrnng
Series Ti-
tle: Lecture Notes in Statistics.
Bernhard Schölkopf. Ralf Herbrich, and Alex J. Smola. A Generalized Representer Theorem. In David Helmbold and Bob Williamson (eds.), Computational Learning Theory, pp. 416-426, Berlin, Heidelberg, 2001. Springer. ISBN 978-3-540-44581-4, doi: 10.107733-540-44581-1177
Jack Sherman and Winifred J. Morrison. Adjustment of an inverse matrix corresponding to change in one element of given matrix. The Annals of Mathematical Statistics, 21(1):124-177, 1950.
ISSN 00034851. U : www. jstorsstooo.oooo wttp::/bbtt/eee5t53333333
Felipe A. C. Viana, Raphael T. Haftka, and Valder Steffen. Multiple surrogates: how cross-validation errors can help us to obtain the best predictor. Structural and Multidisciplinary Optimization, 39 (4):439-457, October 2009, ISSN 1615-1485 doi: 10.
URL https:
//doi.org//il http:////11110000001111000
008 0338-8
Yong Wang.
Gauss method.
WIREs Computational doi: 10.1002/wics.120..
Statistics, 4(4):
415-420, 2012.
ISSN 1939-0068.
URL https:
//nline ibrrar. wi ley com doi abs/ /10
_eprint:
Yifan Ye, Zhanxue Wang, and Xiaobo Zhang. An optimal pointwise weighted ensemble of surrogates based on minimization of local mean square error. Structural and Mul- tidisciplinary Optimization, 62(2))529-5444
August 2020. ISSN 1615-147X, 1615-1488. URL http://link springer wwt.eo...mmoo.......
doi: 10.1007/s00158-0200 020-00208---044 s0C158 020-020000000
Seniha Esen Yuksel, Joseph N. Wilson, and Paul D. Gader. Twenty Years of Mixture of Experts. JEEE Transactions on Neural Networks and Learning Systems, 23(8):1177-1199 1193,
August 2012. ISSN 2162-2388. doi: 10.1109/TNNIS.2011 2012..20999.. URL https:
lire.ieee..e.............oe..e.e ieee. org httg://ocumett./m /21eent/2111 621 5C56/ raramneuuerrrr..r...a httpuuubber6211555..55 Conference Name: IEEE Transactions on Neural Networks and Learning Systems.
7/ieeexp lore.
Luis E. Zerpa, Nestor V. Queipo, Salvador Pintos, and Jean-Louis Salager. An optimization methodology of alkaline-surfacttatt-oo-- ulation and multiple surrogates. Journal of Petroleum Science and Engineering, 47(3))197~ 208. June 2005. ISSN 0920-4105.
flooding processes using field scale numerical sim-
doi: 10.1016/) 03. 03.002 URL https: //www.
sci scienced rect. com/ sc www.e//rticeeeneeee htt.:///cc......a.aa////////////
Cha Zhang and Yunqian Ma (eds.). Ensemble Machine Learning: Methods and Applica- tions. Springer New York, New York, NY, 2012. ISBN 978 1778-14440 1-44199332-99--0 978- 178-1-4419-9326 7. doi: 10.10077978-1-444
URL https::/ wntk///nnnnn.///// www..n..sprri.rrr.....mm. com com.10.100......11
978-1-4419-9326-
Xinyu Zhang and Chu-An Liu. Model averaging prediction by K fold cross-validation. Journal of Econometrics, 235(1):2803011, July 2023. ISSN 03044076. doi: 10.1016/j.
13

---

Published as a conference paper at ICLR 2025
aecomom.2022.04.000 URL $0304407622000975
A THEOREM ON THE STABILITY OF THE MEVA
We consider the setting without a dependence on x to prove theorem on the robustness of MEVA compared to MEEA. Let Y be random real-valued variable and the random vector of models M E R", Suppose we have i.i.d observations (M )"...... (M1 Y collected in the matrices M E g& and y E RRN We wish to study the Mean Squared Error (MSE) of an aggregation method of the form MA = aTM. M. Define
((a) = E [(Y - - MM]]
(20)
Let C' = E [MM'] E [YM), have that the loss is minimized at a² C 1 the MEA. Let
M Y1 +Z, where Z is the scaled error, with chosen that IE [zZ"] = A has Frobenius norm 1 and b: == E [Y We have that
7 = E + cb
(21)
C - C= (E[Y²)1 + eb)1 + c1b
(22)
A. CHARACCERIZATION OF a
To facilitate our computations, we introduce = A (E|Y3) + eb) We define
s 11, { L..
(23)
(24)
u A²'A-"..
(25)
t R A' " + +bb)
E|
(26)
Using the Sherman & Morrison
1950) formula, we find that
D-
A-1
2 1 (F eb)1 )
(27)
2
D-
D'' D='16" D- where eb' D
(28)
1t
Let or, 3,5, S. such that
D o'1==A--1++2
(29)
+ +ÇA-
(30)
Thus, = D c(at + Bu) and 6T D- b=
öt + {u. The coefficients are:
2+d
8 = - रंं
(31)
+ Te++
The computation of C-1 is broken into two parts:
C c- -
(E 1 + eb)
+
(32)
The intermediate terms are:
C-1 1-11, 1,
(33)
+
C- b= = Dーb e(bT) D-16) D-11.
(34)
+
14

---

Published as a conference paper at ICLR 2025
Substitute D-1 1 and D
B A-b
C-1
A 1
(35)
1+ w
1+ w
8 c(ôt + {u)
A A=11+ + (s- (5-
<(ôt + [u)) 3) A bb
C-1b= =
(36)
1 w
1
Write C- xA-1 + 1yyAA11 b, then
2(öt + (u)
x = E
a
+ +
a
(37)
1 ++ w
1 w
B
(ôt + (u) + cc
y Y
(38)
++
1+ w
Using Sympy (Meurer et al.
2017) to simplify these expressions, we find that:
E|y²]
(39)
sE u) + (e + +t) y
+ c
(40)
8 ( u) + (e + t)2
Defining ov = A A-1 1-A-1 and aR =
c+1"A-16 we have that
A 1p
s(E [ u) s( u) + (e t)²
a' = aa= + (1 AlaR where ==
(41)
One may find that
E y - u = 16 12
E (M - Y1)) E [(M- -
Y1)) 20 0
using the Schur complement of the covariance matrix of the vector (Y.MT 7 YIT)T E Ra+1
Thus, À E [0, 11. and a² is convex combination of the MVA, av, and aR. Note that
- ey
(42)
and
-C = E [ (y * 1 + x 8) + (y u x 11)
(43)
1E 2 + 22
ct
E| su + 1| + Eu
(44)
s
sE Y u) + (e + t)²
E Y
= E 2 ?
(45)
YE u) + (e + 1)2
Thus, have that
((a") = de²
= xC(vv)
(46)
and
2
C(av)= = cart
T Aav =
(47)
A.2 PERTURATION OF ov
We want to study the finite sample properties of the empirical estimates, thus we define
A = N (MT - 1y XM y1T) == A
SA
(48)
i = eN (MT 13 = b
VN 5b
(49)
1 σv
E Y
|>>12 = F. [r
(50)
15

---

Published as a conference paper at ICLR 2025
The central limit theorem asserts that 8.A. 8b and 8V converge in distribution Gaussian random variable, with finite covariance matrices. If we study the perturbation of av when estimating it using
Á, we have that L(av) e² 21 1 AA '1 At the first order in VN
(1)-11)²
À- = = A 1 A
OAA
(51)
VN
-'A-1 = 1 1 11 A 8AA-11
(52)
VN
- -'A-AAA = 1 1 1
T A 81A
(53)
(54)
Using the first order expansion of a fraction,
M - n
M
MQ
0 Q 0
17
0
02
+ +017)
(55)
We find that first-order terms cancel out, meaning that
[(āv) = [(av) + N
(56)
=((AA,, where is continuous function of 8A such that [(6A)//1 0(1) as 6A 3--00 O. Since dA converges in distribution to a Gaussian with finite covariance, using the con- tinuous mapping theorem Billingsley 1999), we have that f(6A) converges to f(Z) in distribution, a finite random variable.
( (,AA),
A.3 PERTURBATION OF a'
A.3.1 PERTURBATIONS OF 1, u, AND 8
If (A + 5A)(# 55) = y VN Sy, then at the first order dx
VN
A- - ōAz):
5s - -7- -1A)
5t 10 16A-"---'
2b bf
A.3.2 PERTURBATION OF THE DENOMINATOR
The denominator of C-ty is:
D = s(E + (c+
The perturbation is:
VN( - = 8D - +50 sōu + 2(e + 1)81.
A.3.3 PERTURBATION OF THE NUMERATOR
The numerator of C- is:
N=
A.3.4 FIRST TERM: [ -
Using E| 2| + 5v. á Su, and the perturbation formula for à
(E) = - - - 1
- 11.
16

---

Published as a conference paper at ICLR 2025
A.3.5 SECOND TERM: (c + B
Using i = 1 őt, = b + 8b, and the perturbation formula for A-1ô:
(e + 1)--))1-1- ( (c+ VN 8IAA b (e +
()A 5b -
A.3.6 FINAL EXPRESSION FOR Ñ
Combining the two terms:
Ñ =
- - (6A-) (ōv òu)AA
VN
(57)
+5tA (e 1)A- (e ()A (64)))-(6AAAAA ((6A)--b]
(58)
A.3.7 FULL PERTURBATION OF C
The final perturbed C is:
C- N D
Where:
D D D+ JD, N 1 VN
and:
= - u) (E+ (e+))A
Now, we see that ' AÑ NT AN + - A&N at the first order. We have that
- VNASN = - -)
(ōv - 5u) + 8tb
(59)
- ( (c+
(60)
VNN (ETY9 u) (E A + (6v 8u)s - tôt
(61)
( -0"00"000000001--'00-001111
(62)
+
55AA- + (60 - 55)t
(63)
+uôt + b"
(64)
Then, we have that
-00= - - (8v 8u)
(65)
+ ust + (E+ () --A-AAAA-1
(66)
We have:
1 (c+ 1)b
1 1
(c + 1)b7 a 1
8th' 1
(e + O
ab
(67)
D
VN
VN
D
(e+1) /N
and 8²
Thus,
132
L(âx) =
E [ [[e + 1)2 26d(]
(68)
Da
2 2 t)bT
añ + atb' + (c t)b ON
(69)
+ N AN + ANN
(70)
We define
C. = E 2
(E + t)² - 2( t)6 aD +NT AN
(71)
C 2côt E 2 - 2 18tb"al aD + (E 1)6 &J + 2N A&N
(72)
17

---

Published as a conference paper at ICLR 2025
We have that
Lo
Li
2LodD)
E² +
[(āE) =
+
(73)
1)2
NN
D2
133
N
Sympy computations show that
22 VND2 (-1 A- 0 2 N
C(āx) L(a")+
(74)
Here again, 6A and 8b jointly converge Gaussian random variables, so using the continuous mapping theorem, EN = -1" 1--11---000 - 6d A converges to a Gaussian & with mean 0 and finite covariance (because is a linear combination of 5A and òb, which are jointly Gaussian in the limit).
A.4 COMPARISON BETWEEN MEVA AND MEEA
Building upon the previous sections, we derived the following results for the loss functions of the two methods:
[(åv) = x =(aa)) + 1 N C(ág) L(a') + 22
(75)
2 EN + O0
(76)
VND2
N
s(E Y? -u)
[(a*)= dc2
where X s(E u) + +c
(77)
8
We can reformulate these results as:
1 (c+ f)2
+ +0 (011)
[(âv) L(a")
IN)
(78)
Y uu)
Y EN + +0 +
((dE) L(a*)
1+ D(E Y u) VN
(79)
(80)
In the context of model aggregation, we expect errors to be small, leading being close 0.
While the MVA does not assume that b = E [Y M-Yi] = 0, small b makes À closer to 1. This corresponds to the errors being uncorrolated with the target. This setting of fixed I + 0 is one studied in the theorem. We may see that even for X = 0.5. and small N, say N = 100, the difference between L(av) and C(a*) ) is a factor while the difference between N vN is factor 10, making the latter effect much dominant.
and
However, t 0 makes the - term disappear, resulting in MEEA and MEVA having the same convergence rate. In this scenario: L(âv) L(a*) O(c* ). Thus, when t = 0, both methods perform well, with the performance gap being of the order e We now consider the more general case of a possibly small error correlation with the target, which leads to Ke. In this case.
(1+x)* + +01 x+ s(E Y 4)
L(âv) L(a*)
[1
(81)
2E Y KEN + +0
[(āg) = L(a")
1+
(82)
DEE Y2) - u)
(83)
This comparison reduces to evaluating the relative magnitudes of vi and e While it is not possible compare them a priori, our theorem is aimed at explaining the behavior the aggregation when extrapolating to all r, the regime of very scarce data. For models that are accurate, it expected that the model error c much less than the Monte-Carlo rate of the validation samples. / especially for N not very large. In this case, MEVA remains superior to MEEA, particularly for highly accurate models.
18

---

Published as a conference paper at ICLR 2025
B AGGREGATING PDE SOLVERS IN THE GAUSSIAN SETTING
This section describes the aggregation performance one may get with GP-based PDE solvers intro- duced in [Che al. 202I). Subsection B.I summarizes the method and subsection (.B.eectionBB.dsscriiiiii describes the performance of the aggregation of low-fidelity models obtained by constraining the GPs to satisfy the PDE at a small number of collocation points.
B.1 SOLVING PDES WITH GPS AND THEIR AGGREGATION
We now recall the GP-based PDE solver introduced in Chen et al. (2021) and describe the proposed method for aggregating such solvers. Consider the Laplace equation prototypical example:
(x) (x) = g(x)
for for
x E 1-1,11]
(84)
x x
Place Gaussian prior the solution, i.e., assume ut = ~ (0) where K is Radial Basis Function (RBF) kernel.
Given N interior collocation points X; E [-1, and N, boundary collocation points X3 0[[11, 1}3), we approximate the solution u with the Gaussian process estimator
i = E [S| - = S(X), {(X*) = o( s(X*)]
(85)
To express û, extend the action of distributions (linear functionals) described in Owhadi & Scovel (2019): for two distributions d1, d2, we write
(dd, d2) = / di(x))((, (dd2 dr dy
(86)
Note that (6,,, (,) = k(æ, y). while (6,, o A, 8,) = Ain(I, 8) where A is the laplacian, and A is the laplacian applied to the first variable.
Using this extended kernel, the Gaussian prior & ~ NN x) can be interpreted in terms of dis- tributions, allowing us to define {(x) (5) and Ag(z) 8₂ A(E). In particular, this Gaussian prior implies that E (d)(EE) = K ((1,, d2). Defining the stack of distributions (6x, Ax,,,xx A, xxp 5xx,) ), we can define the kernel matrix
K(o, = K( Øz)
(87)
and similarly, K(x Φ); = K(0,, Φ5). Using this, we have
(-E((X) E(X²)
E
E
A X), {(X*)
(88)
In scenarios with multiple GP approximations of w!, leveraging the GP prior allows computation of correlations. Let
(k) = E [ (k)) {(xb,)))
(89)
be the distinct approximations of the solution of the PDE obtained from different interior and bound- ary collocation points X and b.k, Write Φx (ōx² o ... ōxž o oA, ò xx... ..8g,.)) Since these
approximants (*),, interpreted as models, are also Gaussian processes, we can explicitly compute their correlations:
0 (k) (2)a = K K(x,, ( 2x)) - ((=A((xe) (~AEEr"))
(
77
K K(dr, - (o))--K(o)) K(hi,
K oh K (0k dzh) - K( Φ1) ()) -
K (0), 2). We then use these correlations in equation 3 to compute the aggregation analyzed in the following subsection.
Similarly, E [(( (x)) = K(x, Φx (( -
19

---

Published as a conference paper at ICLR 2025
B.2 PERFORMANCE THE AGGREGATION
Using the GP-based PDE-solver introduced in Chen et al. (2021) and described above, we generate 100 distinct low-fidelity models of the solution of the PDE equation 84] example of which is given figure [b As this method is GP-based, the correlations outlined in equation ( are known, and this approach reduces to nested kriging (Rullière al. 2017). In this context, equation is not only the best (in mean squared error) linear aggregate. It is also the best non-linear aggregate of all models. Figure f shows the results of this experiment. While the individual models are far from approximating the solution due to the limited number of collocation points, the aggregation shown Figure [4d closely matches the true solution and significantly surpasses any single model used in the aggregation. We also observe that mere averaging across all models does not yield the same accuracy level. This highlights that significant gains in aggregation performance can be achieved by optimally incorporating covariance information both within and across models relative the target.
4.35
80
65
-10
40
-23
10
-24
(a)
(b)
(c)
(d)
Figure 4: (a) Real solution of the PDE (b) One of the models aggregated (c) Uniform average of all
models (d) Proposed aggregate equation
C ANOTHER PATHOLOGICAL EXAMPLE: UNKNOWN REGION
In this second example illustrated in Figure we select
(()) 3 coss(2T2)
Mo(x) (()) + c(x) where e(x) N(0,033
(90)
Mn(r)
Mx(x) = -3
Write M = (Ma, MR, Mx) for the vector defined by the three models. It common for aggregation methods, such MoE, to constrain à be a convex combination, which we im- plement here by secking an aggregate of the form Ma(x) = (D(t))r(((())""> Mx)) where Soflmar(((x)))
Σ (xp (m)) To identify the functions ig we sample N 100 points
exp((x (x))
X; uniformly in [-1, 0.5) U 00.5, leaving the region [-0.5, 0.5) with data. We use kernel
method with an RBF kernel and loss 4 in which || 2 is the RKHS norm of k. Using the representer diiite-dimensionall non-linear optimization problem and solve it using Gauss-Newton iterations (Wang) 2012). In this example, the good model consistently better than the bad models over the entire dataset. Nevertheless, our aggregate fits the dataset perfectly by combining MB and MN instead of using Mc- In doing it fails to learn the quality of the models. Instead, it only interpolates the data, leading to poor performance outside the dataset despite having model with good performance everywhere. This pathological behavior arises because the aggregate fails to account for model errors in its approximation of Y even though it is convex combination.
theorem (Schölkopf et al.) (200)) reduce this loss
D NON-INDEPENDENT MODEL ERRORS
D.1 NEW FORMULAS FOR THE AGGREGATION
In section 3.3 we assume the models have independent errors. In general, this may not be true. Thus, present here extension to non-independent errors. Recall approximation of the
20

---

Published as a conference paper at ICLR 2025
2-
larget
-1.0
-0.05
-1.25
0.0
0.35
0.90
9.
1.00
Figure 5: Pathological example (sec. 2.2
covariance matrix A(x):
A1(2
( [Pp =
A(x) = PT
0
0
P where
(91)
(); are regressors
eAAa(u)
We may adapt the loss 12)
account for + Lns to get:
N
À argmin
p Diagle/b(X -
(92)
IEH
2m11
argmin
allil
(93)
IEH
With this formulation, we get more general aggregation. In particular, despite summing to 1, the coefficients a<z) are not necessarily positive, giving added flexibility in case models have nega- tive conditional correlations. The full Minimal Empirical Variance Aggregate (MEVA) with non- independent errors is given as:
2 (4(x9)((Pe))) -(Pe)) + all/112
argmin 1EH
(94)
alx)" 1TPT ( where D(x) = Diag 1TP"D(x)PP " MA(x) Σ (x(x)
I,..., n
(95)
(96)
D.2 CHOOSING THE MATRIX P
We have explored two options for the matrix P:
P P= I, the identity matrix, is the simplest choice.
as the orthonormal basis that diagonalizes the empirical covariance matrix C "
In the experiments presented below, we only show results for P = I. This choice is straightforward and delivers accuracy comparable to or better than using the eigenvectors of C.
We also observed that using the empirical covariance matrix render the method unstable if the models are highly correlated. Indeed, consider, for instance, two models whose errors are signifi- cantly correlated, leading to the following empirical covariance matrix:
1.1
0.9
C
(97)
0.9
21

---

Published as a conference paper at ICLR 2025
In this case, the second eigenvector of C, P2 = v2 (0.97. -1.02), results in P7 -0.04, which effectively subtracts the two models. This eigenvector associated with small eigenvalue (0.15), indicating that the difference between the two models has low variance average. Since equation 7 depends on the inverse of the estimated covariance matrix, a small estimated variance eigenvalue xz(x) leads (xx) having large coefficients with opposing signs. In this scenario, we observed that the aggregation method can become sensitive to small errors and thereby suffer from poor generalization.
E AN ALTERNATIVE, DIRECT MINIMAL EMPIRICAL VARIANCE AGGREGATION LOSS
Here, we propose an alternate end-to-end loss for the minimal empirical variance aggregation
(MEVA). In section B.I we observed that our approach can be interpreted as minimizing the aggregate's variance. using an unbiased combination s.t. 1"a = Under the model (ov(M(x)) = A(r), so our MVA is such that
axx)
argmin SH 1==1
vJ A(rvv
(98)
Instead of estimating A(x) and inverting this estimation to obtain the minimal variance aggregate as in 7 we propose using an empirical version of the above loss. Specifically, let the empirical covariance matrices A; be:
pd Diag((PP))) = I.., n.P Ding((Y' M& (x x k= I..,)
for general P if P In
A
(99)
Write Hi for normed subspace of the space of unbiased aggregators {u + u(x) E R" s.t. Ei= uz(x) = 1, Vz). Then, the proposed alternative approach identifies an aggregator by minimizing the sample variance of the aggregate:
ā = argmin N
u AX))TA,((X') + all allu/
(100)
1=1
We employ this strategy in the Boston dataset example (section 4. Here, H2 set of vector- valued functions parameterized by neural network, where the output is constrained to have positive entries summing to one via softmax layer. The term [|u] an 7-2moorm regularizer applied to the network's weights and biases. This method is both straightforward and practical.
F ADDING BIAS CORRECTION
Assuming the models M(x to be unbiased is necessary to obtain linear aggregate, as adding a bias correction will result in an affine aggregation of the form (x( = a (x) + (xx)- This assumption may, however, be violated in practice, and in this case, bias will be interpreted as variance. In a context with limited data and little knowledge of the models and target, would be difficult know if an are biased and interpreting as variance hurts accuracy, we may modify our method by accounting for the bias
be attributed to bias variance. In the case where is clear models
M(x) = Y (x) where B[ =(()) Cov [(x) A(x)
(101)
In such model, must, in addition to X, train μ. Using the definition of 2, we know the loss
for j is of the form E (xx(x(X)) A (x*) - e²). Taking P = 7, this bias-corrected
22

---

Published as a conference paper at ICLR 2025
aggregation would be:
n (culx)) - mx(x)))2) t (mx(X*) e²/
X, argmin LmEH, x
(102)
+ a + 611malliz
(103)
ā(x) = Softmaz( (-mmax(-((x)))
(104)
MA Σ āi
- (x()))
(105)
where H1, H2 are two spaces of functions used to approximate the log-eigenvalues of A(r) and the bias of the models.
In the Boston housing dataset example in section 4.1 we tried to implement this bias-corrected aggregation, and it did not yield better results than the simpler method MEVA. While we not prescribe its use in general, we did not consider this model in the rest of the paper.
G AGGREGATION USING VECTOR-VALURI GAUSSIAN PROCESSES AND MATRIX-VALUER KERNELS
We will now model aggregation as vector-valued Gaussian process regression to elucidate the pathologies observed in the abovementioned examples. We start with brief reminder this vector/matrix--vauue generalization of GPs/kernels (Alvarez et al. 2012), before computing the MEA in the GP context.
G.1 VECTOR- GAUSSIAN PROCESS AND INTRIX-VALURE KERNEL
This section is reminder vector/matrix-vallne valued GPs/kernels Alvarez et al. 2012) We define matrix-valued kernel as a function K X' x X + R such that,
Vz, x' E X, K(xxx)) = K(x',x) x)T
TN E X. Vyl,,., VN E R",
Σ ISUSN
of K(o, 20
(106)
This generalizes the positivity and symmetry conditions for defining a kernel. Each matrix-valued kernel uniquely defines zero-mean vector-valued Gaussian process & such that:
Vx: E ((x)) ~ (00 ((x,)))
Vn, E &, Cov (((x),(((')) K(r,x')
One interesting example of such a vector-valued Gaussian process is the case of independent Gaussian processes, corresponding to K being diagonal. One can then understand general vector- valued GP as adding correlation to the coordinates, each being its own GP.
G.2 MEA WITH MATRIX- MATRIX-VALUED KERNELS
We choose matrix-valued kernel K : X' x X' + Ruxn, and write HK for its associated Re- producing Kernel Hilbert Space (RKHS). We get the MEEA by solving equation (see Appendix
MA(x) = (x(()) M(x) = k(.. (Xx, x) +
(107)
where k(x, v) ((x) My). E [x(X) = + NO, AI)] with prior & ~ (N(0,)) E). This observation is crucial as it shows that the model values used feature for the regression, expected. However, it demonstrates that
This also the GP regressor MA
23

---

Published as a conference paper at ICLR 2025
200
0.75
0,50
100
0.25
0.00
0025
100
0.50
0.8
0.6
0.0
0.0
0.2
y1
0.2
y1
0.2
0.2
0.0
0.0
1.0
1.0
(a)
(b)
Figure 6: A pair (f,u) u) s.t. Au = f and u(8t)) = 0, sampled using
18)
the Mean Empirical Error Aggregate (MEEA) only tries to regress Y without considering the un- expect model aggregation to
derlying models' accuracy. This observation is counterintuitive leverage the individual models to approximate the target function better rather than merely regress- ing directly the target. This behavior is particularly evident in pathological examples, illustrating the limitations of the direct regression approach.
H MINIMIZING THE LOSS SECTION G
We aim to minimize the following optimization problem:
à argmin ogH
(x(x;) - Mxx)) + rllaa|
(108)
Applying the representer theorem, we deduce that the solution can be represented as (.). K. ((, where V RTN, Additionally, we define the matrix M E Nxn as
follows:
M
=
(109)
MN
where each M; (M. K(x,, 21), ..... (x;))² K(x;;, xN)). The optimal vector V then
given by:
V = argminn|YY - Meel + XuT Kuu
(110)
wERUN
where 10 is block matrix such that K.j. = K(o,ς1,) for 1 < i,j < N. Define the matrix D E X as:
M( O
0
M(r²)
D
(111)
[T
Observe that M = DK. Using the matrix identity:
(P-1 + Bd R
- R
PB [BPB RR)
(112)
24

---

Published as a conference paper at ICLR 2025
0.30
0.25
Model errors
0.002
4.000g
-3069
-0001
-000183
1bn, atrm
fmm
spectral
1000
33
02
033
3280
8:00
2004
Figure 7: Example of aggregation for the Laplace equation (section 4.3). Line 1: Outputs of different models for given input f. Line 2: Prediction errors for each model. Line 3: Values of a. Line 4: Aggregate and its error
where p-1 XK, B M, and R derive:
V = K-IMT (-MK-LMMP+ +1) y
(113)
= (DKDT +1I--Y
(114)
follows that (DKD"), tion at a new point æ is given by:
(2i, (I),))))) k(x,, Thus, the final model predic-
MA(E) M(x) X)V
(115)
where K(x, X) = (K( .... K (x, so that â(x) K(x, X)V. By identifying that K(z, = k(0, X), we can conclude:
Ma M(x,
-Y
(116)
25

---

Published as a conference paper at ICLR 2025
octpets
180
-20
Model WrTOrS
Depict
PMH
Impt ct
Lan-WWenrco
Spectra
TVD
Rarrant
130
308
300
800
300
300
100
15
Tackie
Pyee
Wireersf
0.49
Tipec
1V0
96.
58
100
300
440
6.1
196
56
E2
20
56
56
120
175
Figure 8: Example of aggregation for Burger's equation (section 4.4p. Line 1 Outputs of the dif- ferent models for given input Line 2: Prediction errors for each model. Line 3: Corresponding values of a. Line 4: Aggregate and its error. Note that the explicit, implicit, and Lax- -ennroff methods all diverged, while the spectral method exhibited spurious oscillations. contrast, the ag- gregated result avoids these issues.
26


Reviews:
[{"review_id": "NdNUodMylu", "reviewer": "Reviewer_E4o3", "summary": "In this paper, the authors study the problem of model aggregation where one intends to learn to combine the predictions given by different models. The authors study two variants, Minimal empirical error aggregation and Minimal empirical variance aggregation, both of which differ mainly with respect to the loss functions used. The former penalizes the error w.r.t. the aggretage function, while the latter takes into account the errors of the individual models. The authors study two applications: aggregating models in housing price prediction and learning aggregated solvers for PDE, the motive being illustrating that the MEVA learns to pick the model with low error for any given input.", "questions": "Would learning the functions in the pathological examples be valid if the input samples are chosen randomly rather than in the specific patterns ? Does this mean that the failure of learning the underlying model is partly due to the choice of the points which are presented for learning ?\n\nIt was mentioned that MEVA doesn't improve the performance in settings when P not equals Identity, (line 263), what is the implication of that in practical cases ? How practical is the assumption that expect multiple models' errors to be independent ?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "- The paper is written well, easy to follow. The key ideas are articulated nicely. \n- The contributed solutions seem significant more in the context of PDE solvers, where the metrics for the aggregated model outperform that of the baselines.", "weaknesses": "It is not very clear in the data-science settings when the aggregated models would be useful. Few more experiments migt add more clarity to that. For instance, the ensemble methods such as random forests achieve similar performance to that of the propsoed schemes. Would there be any advantage in using these aggregate schemes upon the existing benchmarks in data-science problems ?\n\nWhile the proposed formulations differ mainly with respect to the loss function, what may be the data assumptions which would help us understand the performance differences of one loss over other better ? Specifically, when do we expect the loss function (5) to be better than (13) and vice versa ? Some analysis w.r.t. the underlying data assumptions may help understanding them better.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ZnJFsaZ5hn", "reviewer": "Reviewer_k7LR", "summary": "This paper compares two different approaches to model aggregation: the first approach is called Minimal Empirical Error Aggregation (MEEA), which minimizes (over $a$): $\\sum_{i=1}^N\\left|Y^i-a\\left(X^i\\right)^T M\\left(X^i\\right)\\right|^2+\\lambda\\|a\\|_{\\mathcal{H}}^2$, where $(X^i, Y^i)$ are $N$ data points, $\\mathcal{H}$ is a set of functions used for the approximation, and $\\lambda \\geq 0$. In this approach, models $M\\left(X^i\\right)$ are treated as if they are features in regression. Section 3 gives pathological examples that show that this approach may produce an unreasonable aggregation. The second approach called Minimal Empirical Variance Aggregation (MEVA) has the same form of linear weights but $\\alpha(x)$ is obtained differently: i.e., $\\alpha(x)^T M(x)$, where the $i$-th element of $\\alpha(x)$ equals $\\exp(-\\lambda_i(x))/[\\sum \\exp(-\\lambda_k(x))]$, where summation is over $k$. To obtain $\\lambda_i$, first write $e^i$ for the vector with entries defined as the sample errors $\\left(e^i\\right)_k=M_k\\left(X^i\\right)-Y\\left(X^i\\right)$. Then, $\\lambda$ minimizes (over $l$): $\\sum \\sum \\left(\\exp(l_k\\left(X^i\\right))-\\left(e^i\\right)_k^2\\right)^2$ plus some penalty term, where double summation is over $i$ and $k$. The numerical experiments consider aggregation using the Boston housing dataset as well as aggregation of PDE solvers.", "questions": "- Perhaps I missed this, but it is unclear how to specify the $l(\\cdot)$ function in equation (13). It would be useful to provide a concrete example of this, after stating equation (13).\n- What are the most popular or best performing methods in PDE solvers? It might be good to indicate how well the propose method works relative to the most popular or best performing methods in practice.\n- Would it be possible to provide computational speed or difficulties as aggregation would take more time than using only one individual model?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- MEVA seems an interesting idea and can be potentially impactful when model aggregation is necessary. \n- It sounds neat to cast numerical PDE as combining Partial Differential Equation (PDE) solvers through model aggregation.\n- The paper is professionally edited and written.", "weaknesses": "- It is unclear what the benchmark is in the literature. The method called MEEA may not be commonly used in the literature and there is no clear sense that this paper improves the current state-of-the-art.\n- There is no theoretical result that guarantees the superior performance of MEVA relative to MEEA or simple averaging (that is, $\\alpha(x) = 1/n$ uniformly over $x$).\n- There is a hyper parameter to be tuned (e.g., the penalization parameter $a$ in equation (13)), but it is unclear how to choose it.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "JQZ0616TNs", "reviewer": "Reviewer_ppzU", "summary": "This paper presents a data-driven framework designed to consolidate predictions from various models, thereby improving overall accuracy by capitalizing on their unique strengths. This approach is non-intrusive and model-agnostic, treating the contributing models as black boxes and allowing for the integration of outputs from a wide range of methodologies, such as machine learning techniques and conventional numerical solvers. The authors introduce a point-wise linear aggregation method and investigate two strategies for optimizing this aggregation:\n\nThe first, Minimal Error Aggregation (MEA), aims to directly minimize prediction errors, while the second, Minimal Variance Aggregation (MVA), seeks to lower variance by estimating the errors of the models. Although MEA is theoretically optimal when the correlations between models and the target are perfectly understood, its empirical counterpart, Minimal Empirical Error Aggregation (MEEA), frequently encounters issues with overfitting due to inaccuracies in correlation estimates. The authors illustrate that the empirical variant of MVA, termed Minimal Empirical Variance Aggregation (MEVA), consistently surpasses MEEA by effectively estimating model errors and aggregating the results. This method mitigates the risk of overfitting and improves generalization capabilities. The efficacy of MEVA is validated through various applications, including the Boston Housing dataset and the resolution of partial differential equations (PDEs) such as the Laplace and Burgers' equations. In these tests, MEVA demonstrates superior performance compared to individual models and MEEA, leading to a more resilient and precise aggregate model.", "questions": "The inquiry into the handling of biased models raises significant questions regarding the assumption of individual model impartiality. It is essential to explore whether the MEVA method can be adapted to include biased models and what specific alterations would be required for such an extension. Additionally, it is important to consider how the existence of bias within the models might influence the overall performance of the aggregation process.\n\nThe presence of correlated model errors in practical applications necessitates an understanding of how MEVA addresses these correlations. It is pertinent to evaluate the sensitivity of aggregation performance to such error correlations and to consider whether incorporating off-diagonal elements in the covariance matrix could enhance the model's robustness.\n\nA comparative analysis of MEVA against other established aggregation methods, such as stacking, weighted averaging, or contemporary ensemble techniques, would provide valuable insights into its relative performance and practical benefits.\n\nLastly, the selection of function spaces, such as Reproducing Kernel Hilbert Spaces (RKHS), and the tuning of regularization parameters are critical components of the method. Practical guidelines or heuristics for these selections, as well as an assessment of how sensitive the method's performance is to these hyperparameters, would be beneficial for practitioners.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper introduces an innovative methodology for model aggregation that shifts the emphasis from the conventional goal of minimizing empirical error to the more nuanced objective of minimizing empirical variance. This approach is relatively uncommon in existing literature. Traditional ensemble techniques, such as bagging and boosting, primarily concentrate on reducing variance or bias through averaging or additive strategies. In contrast, this study highlights the importance of variance minimization via error estimation, a perspective that has not been extensively investigated.\n\nThe theoretical underpinnings of the proposed methods are robust, featuring comprehensive derivations and justifications. The authors conduct an in-depth examination of scenarios where the Model Error Estimation Aggregation (MEEA) falters while the Model Error Variance Aggregation (MEVA) excels, offering valuable insights into the benefits of focusing on variance minimization. This discussion is pertinent to the well-documented challenges of overfitting associated with empirical risk minimization, where models that prioritize empirical error reduction may struggle with generalization.\n\nThe organization and clarity of the paper are commendable, as it effectively presents complex ideas in an accessible manner. Necessary definitions are provided, and the step-by-step explanations aid in comprehension. The use of illustrative examples, figures, and comprehensive descriptions of experimental methodologies further enhances the work's readability and transparency.\n\nBy illustrating that a focus on minimizing empirical variance can yield superior aggregation performance compared to merely minimizing empirical error, this paper makes a good contribution to the field of ensemble learning. 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.", "weaknesses": "The methodology is predicated on the premise that individual models serve as unbiased estimators of the target variable. However, in real-world applications, models may demonstrate bias stemming from systematic errors or incorrect model specifications. This concern is well-documented in the context of model averaging, as highlighted by Claeskens and Hjort (2008), where the presence of biased models can detrimentally influence the overall aggregation. Although the paper acknowledges this limitation, it does not delve deeply into the implications of such assumption violations on the aggregation process or propose modifications to accommodate biased models.\n\n2. The experiments conducted, while illustrating the efficacy of MEVA in certain scenarios, exhibit a degree of limitation in their breadth. The Boston Housing dataset, being relatively small, may not adequately capture the complexities associated with larger and more intricate datasets. Furthermore, while the PDE examples are pertinent, they could be enhanced by incorporating additional real-world applications to better demonstrate the method's applicability across various contexts. Recent advancements in operator learning, as discussed by Kovachki et al. (2023), tackle complex PDEs in higher dimensions, and evaluating MEVA within these frameworks would bolster the validity of its claims.\n\n3. In the experiment utilizing the Boston Housing dataset, the enhancements offered by MEVA over alternative methods appear to be modest, particularly regarding mean squared error (MSE). This observation prompts a reevaluation of the practical significance of the method's purported benefits in conventional machine learning tasks. Established ensemble techniques, such as gradient boosting, have demonstrated substantial improvements in these areas, suggesting that a more comprehensive comparison of MEVA's performance against these well-regarded methods would be advantageous.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Fd4qHvRWYA", "reviewer": "Reviewer_aX4C", "summary": "This paper focuses on model aggregation by minimizing empirical variance rather than empirical bias, which can lead to better generalization. By making appropriate assumptions about the error model, the authors present a tractable aggregation scheme using nonparametric methods. They evaluate their approach on a housing dataset and two PDE solver problems, demonstrating its accuracy and effectiveness.", "questions": "1. Figure 2. It appears that MEVA does not outperform the random forest model in the housing data experiment. Can the authors explain the reasons for this outcome? And I have a minor comment on Figure 2. The labels in x-axis are too small to see them clearly. It is better to arrange this.\n2. Can the authors discuss the robustness of the approximation of $A(x)$ mentioned in Line 258? Is it reasonable to expect this approximation to hold in practice, and what are the potential impacts if the underlying assumptions are incorrect? A discussion on this issue, along with references to relevant literature, would be appreciated.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This paper highlights the drawbacks of current aggregation approaches that focus on minimizing prediction error, using extensive discussions and interpretable examples to provide valuable insights into model aggregation.\n2. The proposed Minimal Empirical Variance Aggregation method is intuitive and performs effectively in certain model aggregation tasks.\n3. The implementation on PDE solvers is both interesting and practically useful.", "weaknesses": "1. There is a lack of rigorous theoretical discussion on the superiority of MEVA over MEA. Including some theoretical investigations would be preferable. Besides, the implementation of MEVA requires certain approximations. Some theoretical guarantee for the final performance of MEVA towards MVA should be considered.\n2. The derivations of MVA and MEVA rely heavily on assumptions about the error model, which may not hold in practice. The authors should address the issue of misspecification in the paper.\n3. The experiments using real-world data are not comprehensive, as only a single housing dataset is evaluated. It is necessary for the authors to include more experiments to further demonstrate MEVA's performance if the paper is accepted.\n4. Although I am not familiar with the background of model aggregation, I noticed the absence of a literature review and comparison with relevant benchmarks. The cited literature in the related work section consists mainly of reviews and some unrelated papers in scientific computing, and there is no specific model aggregation approaches compared in the real-data experiments. More comparing benchmarks should be discussed and evaluated in the paper. For example, in econometrics, there is extensive literature on model averaging [1,2], which also uses linear combinations of models to improve predictions. The authors should discuss these approaches and include comparative experiments.\n\n[1] Hansen, B. E. (2007). Least squares model averaging. Econometrica, 75(4), 1175-1189.\n\n[2] Zhang, X., & Liu, C. A. (2023). Model averaging prediction by K-fold cross-validation. Journal of Econometrics, 235(1), 280-301.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)