ContinuAR: Continuous Autoregression For Infinite-Fidelity Fusion
Authors: WEI W. XING
Yuxin Wang
Zheng Xing
Venue: NeurIPS 2023 poster (2023)
Keywords: Gaussian process
autoregression
multi fidelity
nonparametric Bayesian
TL;DR: Turning autoregression into a tractable continuous model to deal with infinite fidelity fusion; outperforming the SOTA IFC with up to 4x in accuracy and  62,500x in training time.

Abstract: Multi-fidelity fusion has become an important surrogate technique, which provides insights into expensive computer simulations and effectively improves decision-making, e.g., optimization, with less computational cost. Multi-fidelity fusion is much more computationally efficient compared to traditional single-fidelity surrogates. Despite the fast advancement of multi-fidelity fusion techniques, they lack a systematic framework to make use of the fidelity indicator, deal with high-dimensional and arbitrary data structure, and scale well to infinite-fidelity problems. In this work, we first generalize the popular autoregression (AR) to derive a novel linear fidelity differential equation (FiDE), paving the way to tractable infinite-fidelity fusion. We generalize FiDE to a high-dimensional system, which also provides a unifying framework to seemly bridge the gap between many multi- and single-fidelity GP-based models. We then propose ContinuAR, a rank-1 approximation solution to FiDEs, which is tractable to train, compatible with arbitrary multi-fidelity data structure, linearly scalable to the output dimension, and most importantly, delivers consistent SOTA performance with a significant margin over the baseline methods. Compared to the SOTA infinite-fidelity fusion, IFC, ContinuAR achieves up to 4x improvement in accuracy and 62,500x speedup in training time.

ContinuAR: Continuous Autoregression For Infinite-Fidelity Fusion
Wei W. Xing
School of Mathematics and Statistics, University of Sheffield Hicks Building, Hounsfield Rd, Sheffield, UK, S3 7RH ac.uk.
Yuxin Wang
Zheng Xing
Department of Statistics and Data Science National Univerisity of Singapore 21 Lower Kent Ridge Road, Singapore, 119077. yuxinwang@u. nus.edu edu
Graphics& Computing Department Rockchip Electronics Co., Ltd Fuzhou, China, 350003. zheng wxiggrxicincciiipeeeee
Abstract
Multi-fidelity fusion has become an important surrogate technique, which pro- vides insights into expensive computer simulations and effectively improves decision-making, e.g., optimization, with less computational cost. Multi-fidelity fusion much more computationally efficient compared traditional single- fidelity surrogates. Despite the fast advancement of multi-fidelity fusion tech- niques, they lack systematic framework make use of the fidelity indicator, deal with high-dimensional and arbitrary data structure, and scale well to infinite- fidelity problems. In this work, we first generalize the popular autoregression (AR) to derive a novel linear fidelity differential equation (FiDE), paving the way tractable infinite- infinite-fidelity fusion. We generalize FiDE tem, which also provides unifying framework to seemly bridge the gap between many multi- and single-fidelity GP-based models. We then propose ContinuAR, rank-1 I approximation solution to FiDEs, which is tractable to train, compatible with arbitrary multi-fidelitt data structure, linearly scalable to the output dimen- sion, and most importantly, delivers consistent SOTA performance with signif- icant margin over the baseline methods. Compared to the SOTA infinite-fidelit fusion, IFC, ContinuAR achieves up to 4x improvement in accuracy and 62,500x speedup in training time.
nigh-dimensional sys-
1 Introduction
Contemporary scientific and engineering endeavors depend significantly on analyzing highly com- plex systems, where the repeated execution of large-scale differential equation numerical simula- tions, often with intricate interconnections, on-chip (SoC), more than 80% of the design time is spent on analysis based on different types of simulations, c.g., timing analysis and yield optimization. Such an intense computational demand prompts the use of a data-driven surrogate model, which essentially acts functional approxima- tion for the input-output mapping of a simulation. To enhance convergence efficiency, especially in scenarios involving recurrent simulations, such as those encountered Bayesian optimization (BO) [M]] and uncertainty quantification (UQ) [2]]
essential. For instance, in the design of a system-
Corresponding author.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Traditional surrogates are trained on many high-fidelity simulation results, which are still compu- tationally expensive to generate. To make high- high-fidelity predictions [B] while further reducing the computational burden, it is possible way to combine low-fidelity results. More specifically, we can run low-fidelity simulations based on simplified equations (e.g., reducing the levels of physical detail) or coarse solver setups (e.g., using coarse mesh, a large time step, a lower order of approx- imating basis, and a higher error tolerance) to generate cheap but inaccurate results, which offer a chance to train rough model with low cost. The multi-fidelity fusion techniques then improve such a rough model by using only a few high-fidelity simulation samples. In total, the computational cost is dramatically reduced. Multi-fidelit fusion, owing to its efficiency, has garnered growing interest in BO [.L. E], UQ [B], and surrogate modeling [μ]].
While many state-of-ther arr (SOTA) fusion approaches have been rapidly developing, most focus on improving the model accuracy or scalability for large-scale problems. They normally assume a small number of fidelities (say, less than five), and a particular fidelity data structure (i.e., the output space well aligned and the high-fidelity samples' corresponding inputs must form a subset of the low-fidelity inputs). However, in practice, a natural multi-fidelity problem can be much more complicated. For instance, the data does not admit any particular structure. Furthermore, the number of fidelities is in general countably infinite, very few works actually utilize the fact that the fidelity levels are implicitly quantified by a continuous variable, e.g., the number of nodes in a mesh and there are infinite number of fidelities. These limitations seriously hinder the applications of multi- fidelity- beseed methods, e.g., multi-fidelity BO and multi-fidelity Monte Carlo.
Recently, Li et al. [8] propose the first infinite-fieelity fusion, IFC, which utilizes NeuralODE [B] to resolve these challenges. Despite its success, IFC is difficult to train and scale poorly to high- dimensional problems. To make further step towards practicality while preserving tractability and accuracy, we propose the first tractable infinite-fidelit fusion framework, linear fidelity differential equations (FiDEs), and its rank-1 solution, ContinuAR, to deliver a powerful yet tractable fusion model. The novelty of this work is as follows,
1. We propose the first linear fidelity differential equation (FiDE) and its general solution, paving the way tractable fiddiit-fidelity fusion.
2. We extend FiDE to FiDEs, which handle the common siggh-dmennional simulation prob- lems. Furthermore, FiDEs bridge the gap between multi- fidelityyy and single-fidelity surro- gates and serve as a unifying framework for many existing GP-based surrogates, revealing some future directions for classic surrogate models (a.k.a emulators) with multi-fidelity fusion.
3.
We propose ContinuAR, a tractable and efficient rank-1 solution to FiDEs. It is scalable to infinite fidelity, capable of handling high-dimension problems, compatible with arbitrary multi- fidelity data structure, and delivering SOTA accuracy with low computational cost.
2 Backgronud
2.1 Statement of the problem
Provided with multi-fidelity dataset Dí, where i = O,...,T each set comprising entries (4,, (1) N' yn J==1, with 4, representing the fidelity indicator (e.g.. the number of nodes in a mesh generation), where a larger t indicates a more accurate solution; x(i) E Rº denotes the system inputs (e.g., a vector containing parameters found in the system of equations or initial-boundary conditions for a simulation), where Q is the dimensionality for x; y(9) E R² indicates the vec- torized outputs associated with x, where D is the dimensionality for y; T is the total number of fidelities. In general, higher fidelity simulations, being closer to the ground truth and more costly to acquire, result in limited number of available samples, i.e., No > N1 >> > NT. In most works, e.g., [II0-IZ], the system inputs of higher-fidelity are chosen to be the subset of the lower- fidelity, i.c., xT' c... C X2 X We call this the subset structure for a multi-fidelity dataset, as opposed to arbitrary data structures, which we will resolve in Section H. Our objective is to estimate the function y-(T) (x) based on observations at various fidelities {D')T
2.2 Autoregression (AR)
The classic AR model [B] only considers the sealar problem (i.e., y( E IR) and imposes a Markov property for its multi-iidelity formulation. Considering that to and br are indicators (e.g., time steps) for the low- and high-fidelity data, respectively, AR defines
y (T) (x) py (x) + u(x).
(1)
2

---

where is factor transferring knowledge from the low-fidelity in a linear fashion, while 2(x) aims encapsulate the remaining information. Assuming zero mean Gaussian process (GP) prior
[D3]] for y00 (x) and u(x), denoted as y0 (0, k. (xx x')) and u(x) (O k" x')), the
high-fidelity function similarly adheres to GP. This results in an elegant joint GP applicable to the combined observations y = [y(0))
y:u) ~ ~N 0,
)
K (0)
K
(OT)
y(())
pko (T))
pK (T)
+ K!))
(T)
(2)
where [K0 lis (x(x K 1y = (x(x x
K
(0T) les
h2 (x
[ 1
(x
x)))
0
x (T) We can see that the joint likelihood is still Gaussian, which admits tractable solution for model training as in standard GP. Further- more, the predictive posterior is also standard GP posterior, which can effectively utilize low- and high-fidelity data. Moreover, by decomposing the likelihood and predictive posterior into two separate components, Le Gratiet [MM]] managed to reduce the complexity from O((NO + N7) to O((N9)3 (NT))) using a subset data structure, where X C C
(T) (T)) K lis = (" (T)
2.3 Fidelity Differential Equation
In scientific numerical simulations, the fidelity be determined by quite different factors, e.g.. using simplified/compllte equations to describe the target system or implementing a dense/sparse mesh discretize the domain. For problems where the fidelity factor cannot be easily quantified, the classic multi-fidelity methods, e.g., AR, and other methods should work just fine. However, more frequent to find that the fidelity indicators have valuable information in them. For instance, the fidelity is often indicated by the number of elements in finite element simulation or time steps in forward/backward Euler timing scheme. In these cases, should treat the fidelity indicators as continuous variables and utilize the information they carry to boost the predictive accuracy. To this end, Li al. [8] propose general formulation of infinite- fidelity models,
dy(xx, l)/d/ = oxxx I, y(x, 1))-
(3)
Although this formulation is general enough to cover all multi-fidelity models, it lacks an insightful interpretation- -t is not clear how to design the function o(x, t, yxx t)). A clever workaround is to use a neural network to approximate o(x, (, l, (x, ()) [S], which introduces many challenges such as 1) requiring large amount of training data, 2) expensive computational cost for backpropagation even with adjoint method, and 3) instability as we can see how easy a simple nonlinear system can lead to chaotic behaviors (c.g.. the Lorenz system [T4])..
3 Proposed Method
3.1 Linear Fidelity Differential Equation
To remedy the over-general formulation of Eq. (B), we first revisit the classic AR and reveal its connection to linear ODE's forward Euler solution. Rewrite Eq. (II)
y(x, LT) y(x, lo) Boy(x, to) (((x
(4)
to
and take the limit (tr- (0r-to) to) + 0, we derive Proposition (see Appendix B for detailed derivations). Proposition The linear fidelity differential equation (FiDE):
j(x. t) ((t))(x,t) u(x,t).
(5)
If we take forward Euler solution to solve Eq. (6) discretized each fidelity, we recover the classic
AR in Eq. (M). Let stick to the continuous formulation and utilize calculus techniques [ we
can derive the general solution to Eq. (E)
y(x, t) = BB()) B(A) ( ([BWWuls u(x, A)dA+ + C(xx)
(6)
where C(x) a function x, and B(t) 1 B(T) is the antiderivative of Bt)). In order to design 3(t) such that the solution converges to the ground truth as t + co, a detailed stability analysis is necessitated. General stability considerations for systems of this kind are often complex and may warrant the application of Lyapunov's direct method. For the sake of maintaining tractable model, we make the simplifying assumption that 3(1) 3 a constant value. Provided that > our system is guaranteed to converge to what we define as the "ultra-high- high-fideliity ground truth" as
3

---

t + oo. Setting to as the lowest-fidelity indicator, we derive general solution (see Appendix
y(x,t) t) to)e to)e- + [
u(x,
To design subtle model with enough flexibility while trying to maintain tractability, we place zero mean GP priors for yx, (o) and u(x, 1) (with all data being normalized), i.e.,
v(x,ta) (0, k2 (xxxx))) x')), u(x,t) N (0, (x, t, x'))))
(8)
Similarly the latent force model [TG)), the resulting general solution y(x, t) is also GP with zero mean and covariance with the analytical form:
Covly(x, t), vxx, ()] x* (x,x')e A(t-to)
[" t
E-B1t'-r) Tx,,, d'dd'd...
(9)
Alvarez et al. [[G]] demonstrate that a simple ODE system, devoid x, yiclds nonstationary output covariance with closed-form solution when employing equared-exponential (SE) kemnel for u(t). This result, however, is not directly applicable to the FiDE due to the inherent x dependency.
To improve model flexibility, we enable the usage of arbitrary kemels, e.g.. deep ker- nel [لαس) by implementing Monte-Carlo integration and reparameterization trick [[S]] to conduct efficient backpropagation for model training- integral approximated by M Enj e e- (xx,,, (x,)))))) where t; and ł are M are random samples from [to. t] and [oo d'), respectively.
3.2 FiDEs: Multi-Variate Extension For High- dimennsina Problems
In practice, the outputs of simulation generally high- high-ddmeensioal [[S]], i... y(t) R. We hereby generalize Eq. (B) for a more general formulation that describes multi-variate system:
Proposition The linear fidelity differential equations (FiDEs) for multi-variate system: fxx,t) + B(t) Su(x, t). where B(t) E DDD and E RxR are affine transformations; u(x, t) E F is a source function.
(10)
We can derive a general solution for Eq. (III) with some calculus as model design guideline (see Appendix DD. For instance, for constant B, all eigenvalues of B must have positive real parts to keep the system stable. To derive an efficient model, we define a constant diagonal matrix and derive tractable general solution
R
va(x, Cx(x)e- BBtt + where Ga(u,(x,t))
Σ SanGa(u, ()),
(11)
1 c- ux T)dT,
(12)
with Bg being the d-diagonal element of B. Due to the linearity of the integral Ga(u, (x, f)), we derive the output correlation with an analytical form:
Covlya(x,, f), yu (x', ,t')] = - -B(t-to) x Ha. + Σ Sa. Sw Cov [Ga (u,(x,t)), Ir (u,(x', f'))]
(13)
where
Cov [ga (u,(x, t)), , Ge (xx(x(++))
e-aleB((--r) 1" eB (x, (x,T,"' 7' dr.
-Bul(t-r)
(14)
We recognize that this model is generalization of semiparametric latent factor model (SLFM) [0]] which is special case by setting t = to or Ba = 0 to consider single-fidelity problem.
3.3 ContinuAR: Rank-1 Solution to FiDEs
special yet practical case SLFM is the intrinsic model of coregionalization (IMC, also a rank-1 approximation to LMC [[I]),, where all u, (x, t) share the same kernel function, i.c.
u(x, t) ~ NO (0, k" (x, tx,'t')) @ I).
(15)
Here @ means the Kronecker product. Similar IMC assumptions are made in ResGP, a popular modification of AR for high-dimensional output fusion [12], which shows promising results in

---

many physics applications. Based on their conclusions and our intention to keep our model sim- ple and efficient, we define B = SI following our solution to FiDE and place an IMC model for the lowest-fidelity model, i.e., yxx ~ NO0 k0 (x, x') 60 H), where H is the output cor- relations. Substituting the new u(x, Z) and back into Eq. 1a (x,t,x,t) = Cov (ov,y((,,))) 1), yxx',"))
we get the output correlation
[[ e (e-(t--)) {" e (xx r') @ SS
e 8it 1 09 H
(16)
to
which uniquely defines a rank-1 solution to FiDEs. We call it ContinuAR. Lemma 1. Autokrigeability in ContinuAR: the particular values of the spatial correlation matrix H and SS do not matter in the predictive mean as they will be canceled out.
The proof is given in Appendıx E for clarity by basically deriving the predictive mean of ContinuAR. Since the predictive mean is the main concern in high-dimensional problems [8. [ we simply set H = SS in the mean predictions. The computational complexity w.r.t the output dimension is reduced from O(D)) to O(D) and the memory consumption is reduced from O(D?) to (1).
I to significantly improve model efficiency without introducing any additional error
4 Efficient Training and Inference Through Subset Decomposition
We have implicitly marginalized out the underlying function u(x, t) based on its GP prior and de- rived the output covariance of Eq. (IGG). Given a set of observations Y (0);; y())...; Y(T).
we have the joint distribution
y (0)
K (00)
(OT)
N
: KTTT)
(17)
KTO))
where y7(0) = vec (Y( is the vectorization; [K by = kd (x, Xy, ty) is the shorthand nota- tion of output correlation Eq. (IGG). For a small number of total training data of all fidelity, we can simply opt for maximum likelihood estimation (MLE) for the joint likelihood
ND
L L=
log Σ
log((T))
(18)
2
where E is the whole covariance matrix in Eq. (L and y "y.")..... yT)JT However, this approach will soon become invalid, particularly in a multi-fidelity scenario where we expect many low-fidelity simulations.
Since the integration of Eq. (IGG) can be done by parts, ((k) admits an additive structure exactly as in AR. We can follow Gratiet and Cannamela [22]] decompose the joint likelihood Eq. (IB) into independent components provided that corresponding inputs strictly follow a subset structure. i... X
X X!. For problems with only a small number of fidelity (a small T), the subset structure may not be too difficult to satisfy. However, for the infinite (countable) fidelity setting, such a requirement is not practical. Here, we derive a decomposition by introducing virtual observations Y for each fidelity such that Y(") satisfies the subset requirement for the completed set (T-) (T-) (T) is the part of (T) that forms the subset of (T-1) (with a selec- tion formulation ( E(T) ((T-), where (T-) corresponds to the previous-fidelity outputs ((-))). The joint likelihood can then be decomposed, f = =logp(Y = log (isspY((u (0:T--) (T))
= =logp(
(0:T-11) =logp(Y(0:T-1)) + log
p [P(YC)) ((T- (( v(1- (
(D(T-1)
= log p( (o (0:I-11) - (T)
D
(x!r)) T (K!TTT) -1 (T)
log(2TT)
log
K
(TT)
2
2
2
(19)
where
K(TT) - K(T))
(B())T
(20)
is the updated additive kernel with the uncertainty of the predictive variance (T) of the virtual points, whose computation details are given later in Section T.
K (
(
K (T-ITT-)) 1/
(21)
5

---

is the additive/residuall part of the kernel from (T - 1) To
Y(T-1) Y
Y!
Y( - - x
(22)
is the additive part for the completed outputs from (T -1) to T; AT is the time interval between tr-1 and ty; log log ( (0:T- is the log likelihood for the previous T fidelities, which is obtained by calling Eq. (19) recursively. The detailed derivation is preserved in Appendix I due to the space limitation. Through the decomposition of Eq. (I9), the computation complexity is reduced from O(DE )o) to EI
+ N (x(++) ((-)))) (1))) Here, I- indicates the number of samples. Furthermore, when data at new fidelity is obtained, we only need to modify the joint likelihood slightly by adding new blocks, which will be handy for active learning or Bayesian optimization. Model training is conducted easily by maximizing the joint likelihood Eq. (I9) with respect to the hyperparameters using gradient-based optimizations.
4.1 Predictive Posterior
Since the joint model (LL) is a Gaussian, the predictive posterior for the highest fidelity can be derived as in standard GP with a large covariance matrix that requires inversion for once. Similar to the derivation of an efficient joint likelihood in the previous section, we drive an efficient predictive
posterior yxx.. T) ~ (F T", [xx., T)) (see Appendix HD.
y(x,, T) = -BAT x (x., T - 1) +
K (T
Y!7) a
(23)
T),, = -e- -222Ar x Exx.,
Exx,, 1) + (T)) + r ((T) (T((T
where
(tr) I (TT))
(K(TTT)
k
Σ((t)= =I (K!!!!)
(K(TT))
k (TT)
(24)
((r)
(T)T'
111-117-1
K (K(T-1,-11) -1 E
Kas
Kno
with two selection matrixes that follow: (tr)
(E(TT) (T) X()), ((rr) (E (T(T-), (()
(25)
In these equations, (K KDd (x., T, X,, T) kD (x., T - 1, x, T 1) for x2 E (x(()) is the kernel additional part between 7 and 1 for x. and X (tr), [((r) is the predictive variance based on (T-) X of T fidelity model; [R = ko (x. T, X, T) - k-d (x., T - 1, - 1) for Xj E (X(())) is the kernel additional part between between T and T 1 for x. and (tr); Σ (tr)
is the predictive variance based on (x))).
5 Related Work
As stated, FiDEs serve as unifying frame- work that unifies many SOTA multi-fidelity and single-fidelity surrogates. We made connec- tions to some popular models in Table M.
Model
Assumptions under FIDEs
SLFM
T
0:
IMC
T 0:
u(x) ~ (O &I
HOGP
T
0:
R=1;
[88 Tis (x,22)
GPRN
T =
0,8
10.5-00 CP(X)
Multi- Variate GPs are commonly used surro- gates that can be recovered by setting L = Const or T 0 for the FiDEs, which gives us the fundamental semiparametric latent fac- tor models (SLFM) [20]. SLFM a simpli- fied Linear model of coregionalization (LMC)
AR S
(::(())
= =8tt))
ResGP
1.B(() 0,5 I.
NAR
B(0) -000=0((((x S)
MF-BNN
B(t) -0.0(u(x.t)) + f
SC
B(t) -0.0(u(x + P.j.
DC
P3
B(() 0.9( f1) + GP(x.y,))
ResPCA( PCCA(()
[S]
( U. U. -
[2 p. Based on SLFM, intrinsic model of coregionalization (IMC) reduces the computa- tional complexity by setting rank-1 approxi-
ContinuAR
(() =I
Table 1: FiDEs unifying GP surrogates
mation u(x) ~ (O, k" kx, x') @ I): HOGP [2a improves IMC by letting [S = (2, 2j) with tensor decomposition, where z are latent vectors. Higdon et al. [ further simplify SLFM us- ing singular value decomposition (SVD) to obtain S. For full review of the GP-based multi-varian surrogates, the readers are referred to [ for an excellent review. To overcome the fixed bases of SLFM, GP regression network (GPRN
3 P BR)) introduces another GP to model flexible S.

---

Tractable Fusion is the name we give to multi- models where the joint outputs form GP. AR [B] is the most fundamental and tractable solution to FiDEs for a univariate problem, i... D 1. which allows it to take a (discrete) t-dependent B(() = 3(4)) without leading to an intractable solution. To deal with high-dimensional problems, ResGP [|c avoids the difficulty involved with
(()) by setting it to zero and uses a conditional independent GP for u(x, (4) (equivalent to R =1).
Intractable Fusion refers to methods where the joint output is no longer a GP. Normally, they as- sume Bt) 0 to avoid the integral and replace the antiderivative G(u(x. t)) in Eq. (L3) with a regression model that takes y (at the previous fidelity) as model input. For instance, the popular non- linear AR (NAR) [MD]] uses GP((x. y, t4) to replace G(u(x, t)): Deep coregionalization (DC [(X]) extends NAR with a residual PCA to capture S; Wang et al. [M]] further introduce a fidelity vari- ating S(1,) to increase model flexibility at the cost of significant growth in the number of model parameters and few simplifications in the approximated inference; Li et al. [ take the advances of recent advancement of deep learning neural network (NN) and place a Bayesian neural network
[B>]] who add physics regularization layer to a deep NN; Li et al. [[3]] propose a Bayesian network approach to multi-fidelity fusion with active learning techniques for efficiency improvement. To account for the
(BNN) to replace G (u(x (x, f)); A similar idea is proposed by Meng and Karniadakis
missing uncertainty propagation in intractable fusion methods, Cutajar et al. [| use approximation
inference at the price of overfitting and scalability to high- nigh-dimensional problems; In the UQ com-
munity, multi-fidelity fusion has been implemented using stochastic collocation (SC) [27]] which
essentially uses a polynomial chaos expansion (PCE) to approximate ((u(x, t)) under FiDEs. All the above methods do not generalize to fideliyy problems. The seminal work IFC [S] uses a NeuralODE to solve general FiDEs, leading to a great challenge in model training.
Algorithm. The solutions of FiDEs are closely related to the latent force models (LFMs) [TG [S]
where they focus on dynamic processes without For tractable fusion, Le Gratiet [00] improves
the efficiency of AR by decomposing the likelihood functions based on a subset multi-fidelity data structure, which is relaxed through our non-subset decomposition in Section !!
6 Experiment
To assess ContinuAR, we compare with (1) AR [B]. (2) ResGP [Z] (3) MF- [SS and
(4) IFC [8], which are the most closely related SOTA methods for high-dimensions multi-fidelity fusion, particularly with infinite fidelities. Note that AR and ResGP are modification versions accord- ing to [BG] instead of their original versions that cannot deal with non-subset or high-dimensional problems. ContinuAR, AR, and ResGP are implemented using Pytorch. All GP-based methods use the RBF kernel for a fair comparison. MF-BNN and IFC are conducted using default settings from their source codes. As presented in IFC's original paper, there are two variations of IFC, one with deep learning (IFC-ODE) and the other one with Gaussian process ODE (IFC-GPODE), which shows better results. Thus, we show the results of IFC-GPODE in our experiments. All GP-based models trained with 200 iterations whereas MF-BNN and IFC with 1000 iterations reach con- vergence. All experiments are run on a workstation with an AMD 5950x CPU, Nvidia RTX3090 GPU, and 32 GB RAM. All experiments are repeated five times with different random seeds, and the mean performance and its standard deviation are reported.
6.1 Multi-Fidelity Fusion for Benchmark PDEs
We first assess ContinuAR in three canonical PDE simulation benchmark problems, namely, Heat, Burgers', and Poisson's equations as in [28, HT- which produces nigh-dimensional spatial/spatial- ttmpap//ppptalltmmpraallllll fields as model outputs. The multi-fidelity results are generated by solving the corresponding PDEs using finite difference with mesh nodes of {42,82, 82, 162 322. 64"). which is also used as the fidelity index t for ContinuAR and IFC. We use interpolation to upscale the lower fidelity fields and record them at the high- fideliiyiitty grid nodes to provide uniform output across different fidelity. Equation parameters in the PDEs and parameterized initial or boundary condi- tions represent the corresponding inputs. Please refer to Appendix J for comprehensive details of simulation setups.
IFC has an extremely high computational cost which scales poorly output dimension (see Table 1 for experiments with No = 32 and = 0.5). For Poisson's equation, considering the 5x more iterations required, IFC requires about 62,500x more training time than ContinuAR on CPU.
2 https://gituub.coo/hhbbiiiii
http:///ghhubbcoomshbbbiiiiiiint
7

---

0.20
NE-2NN
300
0.15
0.30
RMSS
Our IFC IFC/GPODE
RMSE
RMSE
0.10
0.20
200
0.05
0.10
100
32
96
128
32
64
96
128
32
96
128
Training Samples No
# Training Sampies N?
Training Samples N
reeGP
MI-ONN
0.13
200
0.10
RMSS
Ours GPODE
RMSE
0.08
RMSE
150
100
0.05
0.05
0.03
50
0.00
0.00
32
64
96
128
32
64
96
128
32
96
128
Training Samples No
# Training Samples A??
Training Samples N2
Figure 2: Subset Evaluation with 7 0.5 (top row) and π) 0.75 (bottom row): RMSE against
number of training samples No for Heat (left), Burger's (middle), and Poisson's 's (right) equation.
will take 32 days just to finish the Poisson experiment in our experimental sctup. Thus, we do not consider IFC a practical method for high-dimensional problems because its training cost is almost as expensive as running a high-fidelity simulation; we only apply it to relatively low-dimensional problems, namely, Heat equation and the later real-world application experiments.
FC-G
IFC (CPU/GPU)
Ours
AR
ResGP
MF-BNN
$
Heat
3.0/3.6
0,016
0.012
0.010
0.014
Burgers
66.9/63.9
0.010
0.009
0.024
0.026
Poisson
137/67.1
0.011
0.001
0.001
0.025
Table Training time (seconds) per iteration
Testing RMSE against training time for Heat equation.
Figure 1:
Classic Subset Assessment. We follow the classic experiment setup where the training samples consistently inereased; the fiderfiideltyy data forms superset of the higher-fidee fidelity data for this experiment. To deliver concrete results, the high-fidelity training samples are reduced at rate 7, i.e., ro (t=+))),) = (x In while the removed samples are randomly selected. For each experi- ment, gradually increase the number of lowest-fidelity training data N from 32 to 128 and caleu- late predictive accuracy (RMSE). The statistical result for five repeated experiments (with different random seeds) is demonstrated in Fig. DL The superiority of IFC and ContinuAR indeed highlight the benefits of harnessing the useful information hidden in the fidelity indicators t,. All method ben- efits from larger 77; MF-BNN is unstable due to its model complexity. ContinuAR outperforms the competitors with large margin consistently in call cases (with up to 6.8x improvement). Averaging over all experiments, ContinuAR achieves 3.5x, 2.5x, and 8.7x accuracy improvements over the best competitor for Heat, Burger's, s, and Poisson's 's equation, respectively. Testing error against training time for Heat N 32 and 17 0.5 is shown in Fig. M. where ContinuAR achieves 2.8x and 603x improvement in RMSE and training time. Detailed mean error fields (see Appendix for the compu- tational details) are also provided in Fig. E. which clearly reveals the superiority of ContinuAR by producing minimal red regions (high error) and maximal blue regions (low error).
Non-subset Assessment. We then assess the compatibility of ContinuAR for non-subset training data, which is inevitable in sany-fidelity problems. The setup is similar to subset assessment with the same decreasing rate 77, except that the training data for each fidelity is randomly selected without forcing a subset structure. The results are reported in Fig. BL The results are consistent with the subset assessment, and ContinuAR outperforms the competitors with large margin (with up to 3x improvements). whereas ME-BNN performs poorly as usual. Averaging over all experiments, ContinuAR achieves 1.9x. 1.6x, and 1.5x accuracy improvements over the best competitor for Heat, Burger's, s. and Poisson's equation, respectively. The mean error fields are demonstrated in Fig. !! which draws the same conclusion the subset assessment.
8

---

0.40
resGP
0.30
200
RMSS
0.20
Ours IFC-GPODE
RMSE
150
0.10
0.10
100
32
96
128
32
64
96
128
32
96
128
Training Samples No
# Training Sampies N?
Training Samples N
BNN
0.15
200
RMSE
AR
Ours
0.10
150
GPODE
0.05
50
516
128
32
96
128
Training Samples
Training Samples
Training Samples
Figure 3: Non-Subset Evaluation with 1j - 0.5 (top row) and 0.75 (bottom row): RMSE against number of training samples No for Heat (left), Burger's (middle), and Poisson's (right) equation.
Ours error
All eTTEr
reSOP enror
VP-0WN errrr
OLrS enor
A1 enrer
resOP emro
Nr-chb, aroo
030
400
800
0.0
Figure 4: Average RMSE fields of the subset (left four columns) and non-subset (right four columns) evaluation with = 0.5 and 128 n training samples for Heat equation row), Poisson's 's equation (2nd row), Burger's equation (3rd row), and TopOP (4th row).
6.2 Multi-Fidelity Fusion for Real-World Applications
Next, we assess ContinuAR with real-world applications. Particularly. we look at two practical simulation applications: topology optimization (TopOP) and Plasmonic nanoparticle arrays (PNA) simulation, both of which are known for their high computational cost and render the need for multi- fidelity fusion. The detailed problem descriptions and simulation setups are given in Appendrx The TopOP data contains five fidelity solutions of dimensionality of 1600 based on simulations on a mesh concatenating (64, 256, 576, 1024, 1600} nodes, whereas the PNA data has five-fidelity outputs of two dimensions.
The same subset and non-subset assessments are conducted, and the results are reported in Fig. We can see that the advantage of ContinuAR is clear but less pronounced due to the complexity of the data. Averaging over all experiments, ContinuAR achieves 1. and 1.4x improvements over the best competitor in TopOP and PNA, respectively. Although taking averagely 192s for training (140x more than ResGP and AR) non-subset settings, but not in the subset settings. In contrast, the training time of ours is only 5.7s (4.4x more than ResGP and AR) and the improvements are significant in all cases. Detailed mean error field in Fig. 0 shows significant reduction in red region volume for ContinuAR.
PNA, IFC only shows small improvements over ResGP and AR in

---

HF-INN
0,35
0.50
MP
0.50
AR Ours
0.30
0.40
0.40
0.30
RADE
0.30
0.30
0.20
0.20
0.20
0.15
0.10
96
96
Training eamples
Training Samples
Training Samples
Training Samples
0.35
0.50
0.30
RASE
0.25
0.30
0.20
0.20
0.25
0.20
0.10
64
96
128
32
64
128
32
64
128
64
95
125
Training Samples n
Training Samples 11
Training Samples A
Training Samples N
(a) TopOP with 77 = 0.5 (left) and 7) 0.75 (right)
(b) PNA with = 0.5 (left) and ry = 0.75 (right)
Figure 5: Subset (first row) and Non-Subee (second row) Evaluation.
HF
0.15
0.15
restP
0.10
Curs
0.05
0.00
200
0.30
0.22
300
930
1586
1914
2570
3220
3554
1914
2342
2570
2868
3226
3554
Cost
0.80
2.0.0
930
3258
1586
1934
2242
2570
2898
3226
3556
Figure 6: Cost Evaluation for Heat equation, Burger's equation, Poisson's equation, TopOP, and PNA.
Cost Evaluation. Finally, to assess ContinuAR in more realistic setting, relax the requirement that low- wow-fideltty samples are more than high-fidelity samples and randomly pick samples for each fidelity (ensuring that each fidelity has minimal of four samples) for each models. We do not test under different seeds because the randomness of setting is already sufficient to assess the robustness of the models. The computational cost for generating the training data against accuracy results for all assessment datasets are reported in Fig. B. The superiority of ContinuAR over the competitors consistent with previous experiments. The average accuracy improvements over all settings of ContinuAR over the best competitors is 4.3x, 9.3x, 2.7x, 1.2x, and 1.3x for Heat equation, Burger's equation, Poisson's equation, TopOP, and PNA, respectively.
7 Conclusions
In this work, we propose a unifying framework FiDEs to pave the way for tractable infinite-fided fusion and novel solution ContinuAR, which shows a significant improvement over SOTA competitors. We expect FiDEs and ContinuAR can lead to the development of surrogate-assisted systems deliver cheaper, more efficient, and eco-friendly scientific computation. The limitation of this work includes the rank-1 assumptions and further examination in highly-compla real systems.
10

---

References
[1] Bohak Shahriari, Kevin Swersky, Ziyu Wang. Ryan Adams, and Nando de Freitas. Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proceedings of IEEE, 104(1):148-17) January 2016. ISSN 0018-9219, 1558-2256. doi: 10.1 10.1109/P00..2215.2222 2.44418....
[2] Apostolos F. Psaros, Xuhui Meng, Zongren Zou, Ling Guo, and George Em Karniadakis. Uncer- tainty quantification in scientific machine learning: Methods, metrics, and comparisons. Journal of Computational Physics, 477:111902, 2023. ISSN 0021-9991. doi: https://doi..rr/100.0001 10t j... 111902.
[3] M. Kennedy, Predicting the output from a complex computer code when fast approximations are available. Biometrika, 87(1);1-13, March 2000. ISSN 0006-3444, 1464-3510. doi: 10. 0..109//iomee/87...1
[4] Matthias Poloczek, Jialei Wang. and Peter Frazier. MultiInfformation Source Optimization. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
Jialin Song, Yuxin Chen, and Yisong Yue. General Framework for Multi-fidelity Bayesian Opti- mization with Gaussian Processes. In Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, pages 3 3158-3167. PMLR, April 2019.
[6] L. Parussini, D. Venturi, P. Perdikaris, and G.E. Karniadakis. Multi-fidelity Gaussian process regression for prediction of random fields. Journal of Computational Physics, 336(C):36-50 50, May 2017. ISSN 0021-99911 doi: 10. jop.//.cp...01...00
[7] Zheng Wang. Wei Xing, Robert Kirby, and Shandian Zhe. Multi-Fidelity High-Orde Order Gaussian Processes for Physical Simulation. In International Conference on Artificial Intelligence and Statistics, pages 847- 855. PMLR, March 2021.
[8] Shibo Li, Zheng Wang, Robert Kirby, and Shandian Zhe. Infinite-ficllty coregionalization for physical simulation. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 25965-25978. Curran Associates, Inc., 2022.
[9] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differen- tial equations. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
[10] Loic Le Gratiet. Bayesian analysis of hierarchical multifidelity codes. SIAM/ASA Journal on Uncertainty Quantification, 1 1()::44 269, 2013.
[11] Paris Perdikaris, Raissi, Andreas Damianou, D. Lawrence, and George Karniadakis. Nonlinear Information Fusion Algorithms for Data-Efficient Multi-Fidelity Modelling, volume 473. Royal Society. doi: 10. wttp9///ppa2011.07555.
[12] W. W. Xing. A. A. Shah, P. Wang, S. Zhe, Q. Fu, and R. M. Kirby. Residual gaussian process: A tractable nonparametric bayesian emulator for multi-fidclity simulations. Applicd Mathematical Modelling, 97: 66-56, September 2021. ISSN 0307-904X. doi: 10.
20210030...
[13] Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. Adaptive Computation and Machine Learning, MIT Press, Cambridge, Mass, 2006. ISBN 978-0-262- 18253-9.
[14] Edward N Lorenz. Deterministic nonperiodic flow. Journal of atmospheric sciences, 2002):130-141, 1963.
[15] William E Boyce, Richard DiPrima, and Douglas B Meade. Elementary differential equations and boundary value problems. John Wiley & Sons, 2021.
[16] Mauricio A Alvarez, David Luengo, and Neil D Lawrence. Linear latent force models using gaussian processes. IEEE transactions on pattern analysis and machine intelligence, 35(11):293327755 2013.
[17] Andrew G Wilson, Zhiting Hu, Russ Salakhutdinov, and Eric P Xing. Stochastic variational deep kernel learning. In Lee. M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, Advances Neural Information Pracessing Systems, volume 29. Curran Associates, Inc., 2016.
[18] Diederik P. Kingma and Max Welling. Auto-Encoding Variational Bayes arXiv:1312.6114 [cs, stat], December 2013.
11

---

[19] Dave Higdon, James Gattiker, Brian Williams, and Maria Rightley. Computer Model Calibration Using Dimhnsioal Output. Journal of the American Statistical Association. 103(482):570-583, June 2008. ISSN 0162-1459, 1537-274X. doi: 10.11900161455000000000000000000 4507000000888.
[20] Yee Whye Teh, Matthias Seeger, and Michael I. Jordan. Semiparametric latent factor models. In International Workshop on Artificial Intelligence and Statistics, pages 333-340. PMIR, January 2005.
[21] Mauricio A. Álvarez, Lorenzo Rosasco, and Neil D. Lawrence. Kernels for Vector- Valuud Functions: Review. Foundations and Trends@ in Machine Learning, 4(3):195 266, 2012. ISSN 933-8237, 1935- 8245, doi: 10.1561/2200000036
[22] Loic Le Gratiet and Claire Cannamela. Cokriging- sequential design strategies using fast cross- validation techniques for multi-fidelity computer codes. Technometrics, (1(3)):48-421)) 2015. doi: 10. 1080/00401706.200 201.9.22314
[23] Michel Goulard and Mare Voltz. Linear coregionalization model: tools for estimation and choice of cross-vario variogrrr matrix. Mathematical Geology, 24(3):269 286, 1992.
[24] Shandian Zhe, Wei Xing, and Robert M. Kirby. Scalable High- Gaussian Process Regression, In The 22nd International Conference on Artificial Intelligence and Statistics, pages 261 1-2620. PMLR, April 2019.
[25] Andrew Gordon Wilson, David A Knowles, and Zoubin Ghahramani. Gaussian process regression net- works. In Proccedings of the 29th International Coference on International Conference on Machine Learning, pages 1139 1146, Omnipress.
[26] Shiho Li, Wei Xing, Robert Kirby, and Shandian Zhe. Multi-Fidelity Bayesian Optimization via Deep Neural Networks. Advances in Neural Information Processing Systems, 33:8521-8533 2020.
[27] Akil Narayan, Claude Gittelson, and Dongbin Xiu. A Stochastic Collocation Algorithm with Multifidelity Models. SIAM Journal on Scientific Computing, 36(2):A495 A521, January 2014. ISSN 1064-8275. doi: 10. 10.1197/130929441.
[28] Wei W. Xing, Robert M. Kirby, and Shandian Zhe, Deep coregionalization for the emulation of simulation- based seatial-temporal fields. Journal of Computational Physics, 428:109984, March 2021. ISSN 00219991. doi: 10. je.100/j.ppp..220... 1099949
[29] Pierre Goovaerts et al. Geostatisties for natural resources evaluation. Oxford University Press on Demand, 1997.
[30] Shibo Li, Wei Xing, Rohert M. Kirby, and Shandian Zhe. Scalable Gaussian Process Regression Networks. In Twenty-Ninth International Joint Conference on Artificial Intelligence, volume 3, pages 2456-2462, July 2020. doi: 10.24963/iccai 20200//334ca.222/344.... 2020/340.
[31] Trung Nguyen and Edwin Bonilla, Efficient variational inference for gaussian process regression net- works, In Artificial Intelligence and Statistics, pages 472-480. PMLR, 2013.
[32] Xuhui Meng and George Em Kamiadakis. A composite neural network that learns from multi-fidelity fidelity data: Application to function approximation and inverse PDE problems. Journal of Computational Physies, 401: 109020, January 2020. ISSN 0021-9991. doi: 10.10166j.jcp.2011 2010000000.
[33] Shibo Li. Zheng Wang, Robert Kirby, and Shandian Zhe. Deep multi-fidelit active learning of high- dimensional outputs. In Gustau Camps- Valls, Francisco J. R. Ruiz, and Isabel Valera, editors, Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, volume 151 Proceedings of Machine Learning Research, pages 1694-1711 1711. PMLR, 28 38-30kM Mar 2022.
[34] Kurt Cutajar, Mark Pullin, Andreas Damianou, Neil Lawrence, and Javier González. Deep Gaussian Processes for Multi-fidelity Modeling. arXiv:1903.0733 073200 [cs, stat], March 2019.
[35] Simo Särkkä, Mauricio A. Álvarez, and Neil D. Lawrence. Gaussian process latent force models for learning and stochastic control of physical systems. JEEE Transactions on Automatic Control, 64(7): 2953-2960, 2019. doi:
[36] Yuxin Wang, Zheng Xing, and WEI W. XING. GAR: Generalized autoregression for multi-fidelity fusion. In Advances in Neural Information Processing Systems, 2022.
[37] Rui Tuo, C. F. Jeff Wu, and Dan Yu. Surrogate Modeling of Computer Experiments With Different Mesh Densities. Technometrics, 56(3):372-36) 380, July 2014. ISSN 0040-1706, 1537-2723, doi: 10800 00401706.2011 00000.00....293....
12

---

[38] Mehmet Onder Ele and Hitay Ozbay. Proper orthogonal decomposition for reduced order modeling: 2d heat flow. In Proceedings of 2003 IEEE Conference on Control Applications. 2003. CCA 2003., volume pages 1273-1277. IEEE.
[39] Maziar Raissi and George Em Kariaadakis. Machine Learning of Linear Differential Equations using Gaussian Processes. Journal of Computational Physics, 348:683-693, November 2017. ISSN 00219991. doi: 10/gbzfgr.
13

---

Supplementary Materials
A
Gaussian process
Due to its capacity to handle complex black-box functions and quantify uncertainty, the Gaussian process (GP) is often a popular selection for the surrogate model. For now, let's contemplate simplified scenario where we possess noise-influenced observations, denoted as ys = f(x.) + e, I. In a GP model, prior distribution is placed over the true underlying function of fxx,, whose input is x:
{xx) ~ GP (mo(x), k(x, x(1)),
(A.1)
employing mean and covariance functions.
mx
[x], xxxx'00) (f(((()) mo mo(x))(f((')) mo(x))),
(A.2)
where the hyperparameters control kernel function E[: the expectation. Centering the data allows us to assume the mean function as a constant, mo =\t mon Other options, like a linear function of x, are possible but are seldom employed unless there is prior knowledge of the function's shape. Various forms can be adopted for the covariance function, with the ARD kernel being the most widely favored.
k(x; x x0 Po exp (-(x (x) diag(0," .....
diag((
".... 0,,00 ")(x (xxxxxx))
(A.3)
Starting here, we remove the specific mention of 0's dependence on k(x, x The hyperparameters 0,.... O: are termed as length-scales in this case. When is constant parameter, fxx represents its random variable. In contrast, set of values, fx:). where 1,..., N. constitutes partial realization of the GP. Realizations of GPs are deterministic functions of x. The key characteristic GPs is that the joint distribution of (y(2;), for i i= ,,..., N. a multivariate Gaussian.
We can derive the model likelihood by assuming the Gaussian distribution of the model deficiency N00 )²), along with the utilization of the prior Il and the existing data.
C a p(y|x, 0) (f(x) +e)d = K(ylmol,K+of
(y mo (K (K+o -o -m
(A.4)
+ - log(2m).
where K = [Ka]] is the covariance matrix, in which K, = k(x,, x)), = 1,..., The hyperparameters 0 are frequently obtained through point estimations, utilizing the maximum likelihood (MIE) of Eq. (B with respect 0. The joint Gaussian distribution of y and fx) also possesses mean value of mel and covariance matrix.
K 1+
k(x)
K'
(A.5)
k' (x)
k(x,x) + 0F
where k(x) = x)( k(xN, ()) xT The conditional predictive distribution at x is obtained by condi- tioning
f(x)ly ~ N (((x),u(xxx) xxxxx
x(x) mn1 k(x)" (K + 021) (y - mal),
(A.6)
v(x) /² + x(x,x) k" (x) (
o²I) (2)) k(x).
Given by (x)),, the expected value is [x)] and exx) represents the predictive variance. The transition from Eq. (A to Eq. () is critical, as the prediction posterior of this wake relies on comparable block covariance matrix.
14

---

B From AR to LiFiDE
We first revisit the classic AR and reveal its connection to an ODE's forward Euler solution. Rewrite Eq. (ID follows,
as
v(x, (r) = (y(x, (x, to) + oxx to), v(x,tr) (r) v(x, (v 1))(x, to) + u(x, to), v(x, v(x fo) axy(x, to) + to),
(A.7)
where a = 7 1. We then divide both sides by the constant (1r (o) indicating the fidelity difference and absorb the constant into (p 1) and u(x, 1o) and write
y(x, (T) fa) tr to
ty to v(x, to) + tr vxx, to).
(A.8)
Since a and v(x, to) are values and function to be estimated, we can absorb the constant tr to into a and n(xx to) and write
y(x, (T) v(x, fo) = Bxy(x, + u(x, ta). tu to
(A.9)
where Bo = ax/br to and uxx, 1o) = =(xx tu)/tr tu. We recognize that is explicit solution of a different equation. If we take the limit of (tr - ta) + 0. we have
dy(x, t) 3(Ov(x,)) () + u(x, 1). di
(A.10)
which is Proposition Notice that, this equation, we the constant 3 into function ((t), which allows us to control differential level information transfer depending on fidelity level as general formulation.
Taking a forward Euler solution to solve Eq. (A ITI. we have
(x, + At) y(x, t)
+ ((0)(x, = u(x,
(A.11)
,(x,1 / + A0) = (1 AI3(1))x(xx 1) + Alu(x, 1).
We recognize that this formulation exactly the the classic AR with p = At;((t) with residual GP Atu(x, t).
C A General Solution to LiFiDE
We derive the general solution to the derived linear fidelity differential equation (LiFiDE)
dy(x, t) + 3(0)((w) 1) = u(x, 1),
(A.12)
which is a standard non non-homogencous first order differential equation. We know that for the homogeneous equation, i.e., u(x, 1) = 0, the general solution
y(x, t) = C(x)e- R(O)
(A.13)
where C(x) a functional x, and B(t) = 1 ((r)dd is the antiderivative of B(t). Thus, we assume the
general solution the non-eeoogesssus equation
form y(x, 1) v(x, B(1), where v(x, 1)
functional of x and Substituting this into Eq. (B 12), we have
du(x, 1) e B(t) u(xx t)e B(0) 3(t) + ((t))(x, t)e B(I) u(xx f),
(A.14)
du(x,t) u(x, t)e B(O)
dt
Integrating both sides, have
v(x,t) t) = =f u(xx r)e B(T) C(x).
(A.15)
Putting this back into the assumed solution, we have
(xx =--BB) -B(1) f" u(xx B(() -) C(x)e B(t)
(A.16)
15

---

If we assume constant for B()), ie., [(1) = 3, then have
y(x,t) t) =8 I n(x; T) C(x)e
(A.17)
-(()) R + I u(x,7) r)e --ee- dr.
For practical from where the integral starts from the lowest- fowest-fidelity index, we have
Plo + I.² u(x, r)e Pr
y(x, to) C(x)e PPoo
(A.18)
= =((x)e
where we have assumed u(x; for because we not interested < to. Substituting this into Eq. (7), we have
v(x, t) ((x)e + f a1xx r)e T)e -
v(x, tto s0 for u(x, 7)e si -1e ( The ((x))))--0--))) d7
(A.19)
v(x, to p(e-to) ( u(x, T)e - -(.
This formulation is equivalent to Eq. (A). However, it allows for more flexibility in practice we can give start the model with the lowest-fidelity to and then the model predict the higher-fidelity to.
D General Solutions to LiFiDEs
Consider LiFiDEs taking this general form
dy f)
di + B(Iy(xx,)) = Su(x, L);
(A.20)
where Btt) and are matrices, yx, f) is a vector of unknown functions, and u(x, t) is non-zero vector. The general solution is the sum of the homogencous solution (x, t) and particular solution yp t).
The homogeneous solution comes from the homogencous equation:
dx,(xxt) ((t))n(x
(A.21)
dt
This has the solution:
(x, /BOdtt
(A.22)
where e constant vector. The particular solution comes from the oot-homogeneous equation can be
obtained using the variation parameters. The method involves finding function v(x, t) such that:
(x, = e B() v(x, t)
(A.23)
is a solution to the non-homogeneous system. Suhstituting into the non-homogeneous system gives:
dv(x, B())dt Su(x, f),
(A.24)
dt
which be integrated find v(x,t). Thus, the particular solution is:
y, t) = - B(nja / [e 1)] dt.
(A.25)
The general solution to the non- homoogeneneu system is the of the homogeneous and particular solutions:
16

---

v(x,t) = yn f) yy(x, 8)
=e B(t)dt + e B(tddt
Blodt Su(x, t) dt,
(A.26)
which is similar to solution to LiFiDE Eq. by
Note that the matrix exponential computation is done
(B()))
(B())²
(B(t))
Bit) =I I+ B(I)
(A.27)
2
3!
4!
which is not easy to compute. However, matrix (() a constant matrix 3I, then the matrix exponential can be computed by
0
B(t) 81 Ie
(A.28)
which gives efficient way to derive solution to LiFiDEs,
e y(x, = CC C+ St f t)dt.
Su(x,
(A.29)
E Proof of Lemma 1
Lemma 1. Autokrigeability in ContinuAR: the particular values of the spatial correlation matrix H and SS do not matter in the predictive mean they will be canceled out.
Proof. autokrigeahility ContinuAR is easy to derive once we derive the predictive posterior for subset case in Appendix G and the non-suhset case Appendix HJ. We can see the the matrix H and SS are canceled out in the predictive mean.
F Joint Likelihood of FiDEs-1
an illustrating example. For clarity, we slightly abuse the notations by replacing the previous spatial correlation in the main paper with ( and SS with S(3, The joint probability for Y is
We derive the joint likelihood of FiDEs for observations Y
(1) ((1) ~ N (o.
R (X (0) & ())) -
(X x S 1n (X x s )
pll to 2 (1 (0) & (())
X (
(A.30)
where 4n (1) (()) = (h B(t 1 Bld (1) ( '')d'''' dr.
A ) reveals that the for a two fidelity problem, ContinuAR is equivalent to AR. Thus, the likelihood Eq.
decomposition for subset structure data [[[ also holds for ContinuAR, which lays out the foundation for fast inference algorithm for non-subset problems.
G Predictive Posterior With Subset Structure
In order conduct the proof of decomposition for the likelihood function with non-subset data, need derive the form of predictive posterior distribution with multi-fidelity with subset structure i.c., ( C, ...
(D), first. Here, we take the simplest situation, i.e. holds. To stay consistent with AR, we denote p e B(tz-ta for clarity. Following the format of the covariance
an example, the assumption (1) (0)
17

---

matrix in Eq. AA 30 and Eq.
( R) the mean function and covariance matrix have the following expression.
3(2,,1) 1)
(pk_0 8 S (x
@ (
+ k @ STh) x-
)
(
((a) )
v10)
(ak! )0))
K (00) (I0))
-1 y:th
pkk )
(E") (p"" (000 (B"") @ ((("))-)) 5i0)
(o'k.)
E
z (k.0) @ ) (ok! (B @ 8 -1 500)
(ok! )00)
(I)
(okk!
ya + _2 x ) (K!") @ (()) yo
(pkc & (()))
K (00) ((0) y(u) k @ ((0) (iK (I()) @ (8( ) (00)
+ (k!") @ )
(kk!")) (E(")) T @ (800) 5:0 + k ) (K!") & sss) y.
[kk (K)
& 1 5(o + k ( -1 @1 y.
(A.31)
x!!
= (p2k(00)) @ ( + k(ll' @ (
(pk!") (0h p (x @ (0 + k!! o (t)) xー (
(pk! @ sto) T
o² (k!" (x11) @ StO))
+ (k." @ sty)
(p'k!0) 00)) (0) + ( 8 ) - (ok!") @ )
(K @ (TO)
(pk! SSS
(K!)) @ ) (kK!L)) (1)) @ ) (ok! @ ((0))
- (k!t) s (ok @ Sw)) (ok!0) @ sto))
(Kll @ ) (K.")) & (T))
k0 @ (())
=p² k - (k:0)) T (K:00) k
k@l) - (Kk!) T (K!!))
k() s
(A.32)
where k.)) 1. (X,,, (())
k0 (00) (X X.). K (00) = 1. (x ( denotes the covariance in the lowest fidelity, and k = 1 (X., X ) k. 1nd (X.,X.) ) and K." (1) = 1n11 (x)))
Notice that Eq. (A 3 decomposes the predictive posterior into two parts one related to the low-fidelity data and the other the high-fidelit data. This lays the foundation for our later derivation for the non-subset structure. In the following section, we will derive the predictive posterior for non non-subset structure, and show that the autokrigcahility also holds.
H Predictive Posterior With Non-Subset Structure
In this section, we derive the mean function and covariance matrix in the predictive posterior for non-subset structure following the previous two-fidelity setup.
( y(U) / p( p( 2), (1)) (+)) Y (d)
fidelity-2 posterior
Eidelty-1 posterior
(A.33)
=
- ((Y((x)))) 2) (()) (1))
(1) ()
18

---

Subser part Non subset part
Vinteal complement
X2-1)
X
X
80-1 Xm
Figure 7: Illustration for the selection matrix and multi-fidelity data structure. Let ( be the available system inputs for fidelity. The inputs for t fidelity will be two parts: the subset part X contained in (((-1), and the part that is not contained in ((t denoted as ( (where the hat and superscript indicate that it is complement set for the 1 fidelity). To extract these two parts, we define (() x and (1-1) b(0 (()
As we know, once the { decided, the high-fidelity posterior part can be written as the subset posterior distribution, in the following way,
p 2 (1) ((())) 2x Ngo x
|
1 (x!
-201) (x/a) -1
(x:50) -
exp
2/2
(A.34)
where
{ E [xx,, 1) + $12)
(44) (k( (K (11)
-1
( ) + (k((z)))
Y())
(K(22) (22)
(()))
(A.35)
This depends the conclusion in Eq. (A function.
{ and (2) posterior covariance matrix and mean
At the same time, we know that (1) are samples from fidelity I model, therefore the probability of getting (1) can be written
E @ s(u) exp (5 p(0) T (I!") @ =ww) -1 (s - -
(((-11) 2T
(A.36)
19

---

where the E.Sh @ S means the posterior covariance mairix of sampling YO from lower fidelity model and y()) are predicted mean for the non-subset data. Therefore, the posterior distribution of non-subset data structure is
pl Y (x., 2) (
ndelley-2 posterior
Rddlity-1 pasterior
p( (X((x,)2) (2)
()) 501 )
((
d'
/ { 2m D x
x
2 (xx-) T (I!e)) -1 (5:K) -y(o))
x exp
x1 D
-
1 (5:0) - (1.1)-9))
-(l) 7
(i & (.sse 1 (5 (10)-xa0)) - (iaa
x2
@ S0)
dxx
2
(M
1/2 x 1 8 (00
1 マン (xx -1 " 1/2 (stte) _ (ix") e (isse -1 |000s
X 11
-2
x exp
(2(3)) -1 rT((1 + (stti)
( 83 ()
(stti)) T r (xx3a) -1 ((
x
exp
y
2
2 (x
1
T
( s (()) -1 ) ddr(2)
-2m
x EX @ |
x rr (
r+ x(1
2 (eis)) (>!
T rT (2
T+ E!!
r (x:3) 88
X exp
2 (Vnr))
E & )T() -1 -((xx²)
r (x!3) r >.!
(x!)) (= (
X exp
(D 2 (Tn) - (5-1))]
exp
8T (x)))
r
r E
(A.37)
where ý and I' is defined by the following equations,
y²())
|100 -190-(1 (0- 0 )}
(((1)
y =y!") p (k((2)) (K!
(11)
@ I
(K! (2)
I
0
r (2:2) (2.2)
(((2)
(1.1)
(1.1)
((2)
(A.38)
where E!! denotes the selection matrix which selects the non-subset parts between two fidelities.
X1)
(1
E x x(1] and EE?)) also denotes the selection matrix but which selects all inputs for the
second fidelity, X (()) (B!?)
x()] Please see Fig. for the illustration of the selection matrix for [x
the multi-fidelity data structure.
After that, we can simplify the determinant and the exponential parts by decomposing them into different parts and using the
formula to obtain conclusions.
For the determinant part, we can derive
-1
x
2
S
x
T
r
5 8 /:0
+
(I!S) )1=11
EI @
) -1-
x 1 + +r (I!" -1 'T -
( ) 'r 1
(A.39)
20

---

For the exponential parts in Eq.
( - ( r r E!?
r -
r [! = Σ(') +r - & 8S r
EX ) -1 - (x! )
-
(r" r [
-1 r =."
(x!") 2 )5(1) -1 =r
Σ + r
-Σ 2 ) -1
r) r
( " r TT ( -1 r + x -
(D.)) 8 (0) = =-
(x -(xx!²) + r ( 8 ) " r r.
(A.40)
Therefore, the likelihood of the posterior distribution is
p (xx,,))) () y(1)
NoD =2m
E +T
DD @ (SS) Trr
1
(A.41)
1 (v: (v.)) -y: -
E! + r (x & (tn))) 'r) ' (y!²)
"
-vx-- -
x exp
From the upper formula, for the non-subset data structure, the posetrior mean and covairance matrix are,
y (k. (1) (K(1I) [T] (
y(1) )+ + (k (22))) (22
yz(1)
@ T (2)
(A.42)
2) (2)
T
(I.' @ 8 -1 r.
In Eq. (AZ22) we also prove that matrix S does not take part in posterior mean function which means autokrige- ability still holds in the non-subset data structure. By recursively applying this conclusion, we can easily extend it to multi-fidelity problem and we show the detailed application in our main paper.
I Decomposition of Joint Likelihood With Non-Subset Structure
As we have shown that the autokrigeability holds in subset and non-subset data structure, we thus assume an 1. for easy computation acceleration (with the cost of losing accuracy in the predictive variance). Due to the identical spatial correlations, we no longer need vectorization. First, we decompose the joint likelihood L into several independent parts,
identical spatial correlation.
= log (ogg((((0)) = log p (0) (0) y()), (())
= logp(( (1)) + log Y
(1) (1)) (0)) + logp( (og (+) (()),, y(u)
(A.43)
part
part
+ log -l))
(o) (")) + log
(2) y y p Y (1)
Once the (1) is tixed, the probability of part in Eq.
4 can be written as,
p( (2) (1), (1)) = 2n
x(2) p
x
K
(22)
T (K(22)) -1 Y!22) (
K! (2))
1
())
)
K (22)
-1 -[(
Y(2)
YO
=2m
x exp
2
(2)
(
(2)
(
(A.44)
Then is the part 2 in Fq.
(A it is based on the posterior distribution of lower fidelity,
p( (1)) = 2n
(1)
[-2
1
(800 ($ -1 (xa)-raa) -
x exp
(A.45)
where we use the (2) to denotes the missing point in second fidelity corresponding with fidelity 2, which means the (2) parts, yellow part, in Fig. m
21

---

Therefore, we can combine the Eq.
A with the Eq. 45, the integral part in Eq. A44 can be written as
log / ((X()))/X))
log
2T
k
exp
[2
1 ((
)
800)) )) (k! -1 (((
Y(2)
)
Y
Y
(x(1))
- ( T ($.0. - (x p200)] dre
x2
exp
2
-~(i3)
+
(
2 log |K!
1
2 log E110
1
2
8(a))
)) (k! (k.
- ({((
8(2) )
Y
)
log
exp
2
4/(2)
8(0)
Y
1 (m
T (i!!s) - (Yrt) ( | map
2
(A.46)
After that, we try to decompose vector as,
(
( )
Y( (1) = 2 - 8 8
(A.47)
22

---

The Y is divided into L parts substituting
Y Y and non-subset part Y. Then, we consider the data fitting part pa by W ( calculable,
S Eg. Eq.08 (B into intoEir
Eq.
( 4.
which
gives t the method to make the Eq. (A ( whichging sirdgiivee
/ / , [os]
( Y
+ (Nve
1og(2
1 2 1og
1 -
Toul 1og 2
2
7 (iyiee
_ )
(an 1.
(xpe
{
{
1
-
(K??)))
-
- Doaae Meaay -
- - ( )
(
- Baakeel ( i
(
_
1/2 2 (8(1)) ((55)) - ) ( com
1
(oosg (oaa))
(
(
|
(
- (.:s))
(
+
+No +
(08(277)
1
1 2 log 1
bog
2
- -
2
- - (iikkiie
2 (Boa)
- ((lee) )
/
( (all)
(2)) (2)
- Ceccae
((aay))
i. " "( (
/ 2 (Bookee))
3 ((nn))
/
(log)
( 1
(alaae (xxo / (oo
- 1 (ay) aa
i. (iaa)) K.
{
1- ((eey) Sakkee +
1- (
|
T (
1 1
i (oisg) (oaa))
1 (ine) + (
(pt0))
iiN + tavy (
(iti) 1og
1
log(2
-
3
--
- -
2
- 
-. (eeayy
1
I
V...eee -
/ " (Bleryta)
(Ke (
T ((
2 - (Bo))(e)
1
(010)
(
1
=
log
1
1 (loaa
(eae))
(
-
(igsy)
+
cxp
2
+
(
( -a..keei
-. | + (po)) /
( (5:0))
) 11(x0) (00) } 1xx
[
()
(2))
+
N
2
D
1og((2))
1og(2n)
- (ulkaai
-
- log
1
3
8og
1
Bog
(
T
( - + ( 1-
|
1 B
-
( (ooyyiy
(F
3
(pw)) (11)) - 1 (ma))
- (imaeee
(iezz)
-
(
(eaa
-
1
2
1/ T
(K!?)))
)
(
(
( 1 - ( + (.!5)) (1 --
((80) (k! ( (o
2
- )eee (
((KK)) (2) - ( + +( (E) )) (()) -1 )())
+
T
(
1
(y
T
( ( (K( 2 2
+ ( -1) -1/ ((5) (()) -1 "r()) )
) ))
(4 (4.45)
23

---

Y (2), 2(2))
T 00) (y 0)²
where o= Sherman-Morrisee formula.
After that we try to simplity (A by different parts with
First, we consider the determenant part,
2 log
((1))
1
log
log (B()T [(K!2))
( + (5:U) -1
2B(2)
2
log
K
log
(2))
2 log K log
|
1 log (K! ( -1 +B++ E (D()) -1 (B(E))
T
2
2
2
log
(12))
1 + ( (DS))
-1
(E())
(A.49)
Then, we gather terms with ø and simplify them as,
(( - B (K!
T
(B
(Ba) T (K!
) T (k!22)) -1 (1le)1111
2 T
2
(BB:0)) (K ( ( + (S!!)) ( (B(t) (K!-e)) )))
(k(
-1 (K(
2 ( ( + (311)- ) (B(a))
2
(K! -1 + B(2) (S)L) -1
((()) Φ
2
(A.50)
After that, we consider the interaction part between and Y
K (22)
-1
E(,)
(E
- + (5:00 )
-1 (()) -1 )
=6- T
kk?31)
-1
(K!? -1 ( (1) ( - 2 (S51)) -1) -17
((1)
(A.51)
=6
T
(K! + 23 ($:U)= -1 (B1e)) ((0)
Finally, we simplify the terms related with (1)
1
((())
(2()) (v:))
2
+ (rr))
1
T
(E.)) (K!22))-1 2 ($11) -1) (8) -1 )(1)
Σ
2
1
(v()
T
($is) (8) -1 (( (B122) T ((k))((K))) ( - -p + (S!!!)
( -1 (v(:)
2
1
(B1o)) )
T
[K( (lae) -1 + ( (>)')² ( B (x
2
(A.52)
24

---

Putting all the parts logether, the joint likelihood for the R = 1 is.
L log p( (0:2)) - log (0) (1), Y (())
(oxp(Y())) + log Y ( log / (1/0)/11)) y (x)0)))))000)000))))))))))))))))
= =logp(Y ) + logg(((()))) (( + log / p M( y )( (
=logp(Y ( + log )
1 log(2T)
z log
E(2) ($") -2
[E(2)
2
Y!!) (x:)) (K!
( -1 (E(2) Y!
(A.53)
Therefore, we can see that the joint likelihood of three different fidelities with non-subset dataset can be decom- posed as three independent model to train. For problems with arbitrary number of fidelity, we can easily apply this conclusion recursively to decompose them into the summation structure, which is scalable to the number of training data and fidelity level and is easy code up.
J Experiment in Detail
J.1 Canonical PDEs
Three canonical PDEs are under consideration: Poisson's equation, Heat equation, and Burger's equation. These equations play pivotal roles in scientific and technological applications [АT)) They present typical simulation scenarios that include high-dimensional spatial-tempont field outputs, nonlincarities, and disconti- nuities. These scenarios frequently used as benehmark problems for surrogate models [28, 67)-B9], and y denote the spatial coordinates, with representing the time coordinate, a departure from the notation used the primary paper. The notation in the appendix is employed solely to improve comprehension and does not influence or correlate with the main text.
Let's start with Burger's equation, a renowned nonlinear hyperbolic PDE commonly used to depict various physical phenomena such as fluid dynamics [B nonlinear acousties [H8], and traffic flows [ It is often employed as a standard benchmark for multiple numerical solvers and surrogate models [6T [S]] due its ability to generate discontinuities (shock waves) based on a standard conservation equation. Here the equation in its viscous form,
du Øt
8u
8uu 0x²
8x
Within the context, represents volume, æ pertains to spatial location, t signifies time, and v stands for viscosity. The parameters are defined as follows: x ranges between 0 and (in meters), spans from 0 to (in seconds), and the initial condition is set as u(x, 0) fan /2) with homogeneous Dirichlet boundary conditions. We conducted uniform sampling of viscosities u within the interval of [0.001,0.] (in milliPascal seconds) as the input parameter for generating the solution field.
Solving the problem in the space and time domains involves the utilization of finite elements employing hat functions and backward Euler, respectively. To produce solutions various fidelities, the solver addresses the PDEs by employing discretized spatial-temporal domains consisting of regular rectangular meshes with grid points ranging from 8², 16 32 64 128 facilitating simulations from low to high fidelity.
Next, the focus shifts to Poisson's equation, common elliptic PDE employed in mechanical engineering and physics to simulate potential ficlds, including gravitational and electrostatic fields [HI]]. Expressed in the form of
0 u
8u
O.
Oy²
Encountered frequently in physics, Poisson's equation, despite its simplicity, serves as a fundamental test case for surrogate models, as highlighted various stud- ies [B7, bell. Within our experiment, we apply Dirichlet boundary conditions to a 2D spatial domain, where x ranges over [0. 1] x [0, 1]. The input parameters involve fixed values assigned to the four boundaries and the central point of the rectangular domain, each ranging from 0.1 to 0.9. To generate the matching potential fields as outputs, we sampled the input parameters uniformly. Employing first-order center differencing scheme and regular rectangular meshes, we solved the PDE using the finite difference approach. Simulations of five different fidelities were generated using meshes containing grid nodes 8², 162 32 642. 128².
represents an extension of Laplace's equation
|53
Heat equation, established in 1822 to depict time- dependent changes in heat fluxes, remains fundamental PDE. Its scope extends beyond heat fluxes, finding applications in diverse scientific domains such as probability
25

---

Figure Geometry, boundary conditions, and simulation parameters for cantilever beam
theory |55, [ and financial mathematics [58]]
As a result, it is widely adopted as a representative model. Here
is the heat equation:
JT dx
O7 k:
(k 0z
0T k 0=z
ØT
dr
dy
ay
+ qv
pcp ōt
where k is the materials conductivity qv the rate which energy is generated per unit volume of the medium i is the density and ep the specific heat capacity. The flux rate at the left boundary, ranging from to 1 at æ = 0, the flux rate at the right boundary, ranging from -1 to 0 at x = 1, and the thermal conductivity. ranging from 0.01 to 0.1, serve as the input parameters. Within 2D spatial-temporal domain x E [0. 1] and and 1. Additionally,
t 50, we set up the Neumann boundary condition at H(z 0.25) H(= 0.75), where H(-) represents the Hcaviside step function. Using finite difference in space and backward Euler time domains, the equation is solved. The spatial-temporal domain is discretized into a 16 x 16 regular rectangular mesh for the first (lowest) fidelity solver, while a refined solver employs a 32 32 mesh for the second fidelity. The computed result fields are 100 x 100 spatial-temporal grid. Solving the equation involves applying finite difference approach in the spatial domain and employing reverse Euler the temporal domain. We discretize the spatial-temp metiaal-emppral domain into a regular rectangular mesh using 82, 16?, 322, 64", 128² nodes to produce simulation results at five varying fidelity levels.
define u(x, 0) as
J.2 Multi-Fidelity Fusion for Topology Optimization
In topology structure optimization problem, we employ ContinuAR to determine the optimal topology struc- ture, maximizing mechanical metries such stiffness, for material layout comprising alloys and concrete. This is accomplished by considering various design parameters including external force and angle. Particularly with the advancements in 3D printing methods, where materials are deposited in small increments, topology structure optimization has become a crucial technique in mechanical design, applied in areas like airfoils and slab bridges. Notably, the computational intensity of topology optimization stems from the need for gradicnt- based optimization and mechanical simulations. The demanding nature of a high-fidelity solution, requiring extensive discretization mesh and imposing substantial computational burdens both space and time, exacer- bates the situation.
Gaining popularity is the utilization of data-driven methods that offer suitable structures to facilitate the pro- cess [БА]] In this study, examine the topology optimization of cantilever beam. Utilizing the efficient implementation [S8]] conduct density-based topology optimization, aiming to minimize the compliance C while adhering volume constraints V
Utilizing the SIMP scheme [[S]]. we convert continuous density measurements into discrete, optimal topologies. System inputs include the position of point load Pl, the angle of point load P2, and the filter radius P3 [[0]] The issue is addressed using standard mesh containing nodes in the set 32 40 482 562, 642..
J.3
Multi-Fidelity Fusion for Plasmonic nanoparticle arrays
Using the Coupled Dipole Approximation (CDA) approach, we compute the extinction and scattering effi- ciencies Qent and Qae for plasmonic systems of different scatterer quantitics the final illustration. CDA represents a technique for emulating the optical response of an ensemble of identical, non-magnetic metallic nanoparticles, each with dimensions significantly smaller than the wavelength of light (in this case, 25 nm). In this work, define Qext and Q.e as the Qols. Our proposed method allowed to construct surrogate models efficiently, incorporating up three fidelities. examined particle arrays generated by Vogel spirals, where the interaction of incident waves with particles significantly influences the magnetic field. Consequently, the local extinction field produced by plasmonic arrays is greatly impacted by the number of nanoparticles. The configurations of Vogel spirals, characterized by particle numbers in the set 2, 50, 200, 500, 1000, defined
26

---

the parameters for five-fidelity simulations. The parameter space was determined to be À E [200, 800] nm. due [0, 2T] rad, and des E (1, 1500). These parameters represent the incidence wavelength, the diver- gence angle, and the scaling factor, respectively. Inputs were chosen using Sobol sequence. As the number of nanoparticles increases, the computing time required for CDA exccution grows exponentially. Hence, the suggested sampling approach leads to substantial reductions in computational costs.
By solving the corresponding linear equation, one can compute the local field (r,)(r,) for cach nano-sphere, given N metallic particles with the same volumetric polarizability axw) situated at vector coordinates r;- This calculation allows for the determination of the response of plasmonic array electromagnetic radiation, as per the solution [|T]] of the local electric fields, (r;)
ak²
Etoc (r+) = En(r:)
Σ 5=1.f4i
Gi,,Eoo(TT)
(A.54)
where the incident field is denoted by (r:),, where & represents the wave number in the background medium, and 6o signifies the dielectric permiuivitt of vacuum (eo = 1 in the CGS unit system). Gij is derived from 3 x 3 blocks of the overall 3N x 3NV Green's matrices for the ith and jth particles. When = i, Gij becomes a zero matrix. Otherwise, it is computed as
G.y exp(ikr,))
ikres + (kroy))
(A.55)
Tíj
Expressed Fij. the unit position vector represents the distance from particle j i, denoted by rij
|r
Solving Eq A and A enables the computation of overall local fields (r.,(r;). which in turn determines the scattering and extinction cross-sections. Further information on the numerical approach is available in [[2]].
Derived from the normalization of the scattering and extinction cross- sess--cctions concerning the array's 's entire projected area Qext and Qse are considered. Our focus was on the Vogel spiral class of particle arrays, as
described in [3
Pn Vrid... and On TICK
(A.56)
Vogel spiral array. the n-th particle's radial distance and polar angle are denoted by Pa and O.. respec- tively. Hence the unique definition of the Vogel spiral configuration includes the incidence wavelength A, the divergence angle avs, the scaling factor avz. and the number of particles n.
J.4 Detailed Prediction Analysis
Defining the average RMSE field y (RF) allows us to thoroughly examine the prediction error.
(EF)
ly v.)²
(A.57)
N
The square root taken element-wisg wise, where y. represents prediction and y, is the true ground value.
The average RMSE field corresponding the Heat equation, Burger's 's equation, Poisson's 's equation, TopOP problem with a decreasing rate of = 0.5, lowest-fidelity training samples, No = 128 and 128 testing samples are shown in Fig. E (left) for the subset assessment and Fig. E (right) for the non-subset assessment, Plasmonic nanoparticle arrays (PNA) have only two output variables and thus we do not show the average RMSE field for it.
For the classic subset assessment in Fig. H (left), we see that clear that ContinuAR outperforms the competi- tors all cases with a large margin by showing more blue areas and only tiny red areas. For Heat equation, the error is significantly reduced in most areas except for the bottom area of the domain, where a tiny thin bar of red area is shown. For the Burger's 's equation, ContinuAR shows some checker board pattern in the error field, which is probability caused by the conditional independence assumption. Nonetheless, the error significantly reduced in most areas the largest error also reduced. For the Poisson's cquation, different method has own crror patterns. ContinuAR show a more bluc area in the left part of the domain, while AR has high-error areas everywhere except the center. For the TopOP problem. ContinuAR show a significantly reduced error by a significant reductions of deep red areas. The AR the deepest blue ares, indicating its good performance. However, it also has lar areas of red, indicating its poor performance overall.
For the classic non-subset assessment in fig: errorplot both (right), the overall conclusion is similar the subset assessment. ContinuAR outperforms the competitors in cases with large margin by showing more blue areas and tiny red areas. The overall error pattern for most methods are similar to the subset assessment with subtle difference. For instance, for the Burger's equation, checker board pattern disappear, indicating a successful improvement by using training data across the input domain. The error for Poisson's 's equation
27

---

-1.00
resop
-1.00
6.00
Our
5.00
-3.00
-3.00
-4.00
128
126
Training Samples
Training N
Training Samples
Figure 9: Subset Evaluation with 7j = 0.5: Test negative log- log-likehood against number of training samples No for Heat (left), Burger's (middle), and Poisson's ''ss (right) equation.
MF-BNN
50.0
0.03
resGP
AR
0.15
40.0
RMSE
0.02
Ours
RMSE
0.10
30.0
0.01
10.0
32
64
96
126
32
64
128
32
64
96
126
47 Training Samples N
Training Samples N
Training Samples No
Figure 10: Two-fidelity subset evaluation with 0.5: RMSE against the number of training samples No for Heat (left), Burger's (middle), and Poisson's 's (right) equation.
also significantly reduced particularly on the right part domain. TopOP turns out to the most stable problem as the error patterns for all methods are almost the same as the subset assessment.
K
Additional Experimental Results
K.1
Likelihood Evaluation
Here are some additional results in Figure EL which shows the negative log-likelihood (without the constant term and thus the result
negative) of ContinuAR, ResGP, and AR. Surprisingly, the NLL of our method better than ResGP and more significantly than comparisons RMSE. We believe that this because the log-likelihood is more sensitive to the uncertainty the prediction. Since our method a joint lcarning (as discussed in the previous question), it is able to capture the uncertainty a joint model whercas the other methods treat each fidelity separately.
K.2 Two-fidelity Evaluation
In order to understand the limitation of our method, we follow the classic subset experiment setting with v/ = 0.5 and reduce the number of fidelity from five to two. The results are shown in the Figure IM In this case, we often that ContinuAR is almost identical the classic AR, losing advantages method. Also, when the number of training data is searee, our method does not perform better compared other basclines. There are certainly other factors that may affect the performance, such as the choice fidelity and Bxx). We will investigate this in future work.
infinite-fidelity fusion
28

---

Supplementary References
[10] Loic Le Gratict. Bayesian analysis of hicrarchical multifidelity codes. SIAM/ASA Journal on Uncertainty Quantification, 1(1):244-269) 2013.
[41] Steven C Chapra, Raymond P Canale, et al. Numerical methods for engineers. Higher Education,
Boston: McGraw-Hill
[42] TJ Chung. Computational fluid dynamics. Cambridge university press.
[43] Krzysztof Burdzy, Zhen-Qing Chen, John Sylvester, et al. The heat equation and reflected brownian motion in time-dependent domains. The Annals of Probability, 32(1B):775-804
[28] Wei Xing. Robert M. Kirby, and Shandian Zhe. Deep coregionalization for the emulation of simulation- based seatial-temporal fields. Journal of Computational Physics, 428:109984, March 2021. ISSN 00219991. doi: 10. wttp:.pp...0999000
[37] Rui Tuo, C. Jeff Wu, and Dan Yu. Surrogate Modeling Computer Experiments With Different Mesh Densities. Technometrics, 56(3):372-380, July 2014. ISSN 0040-1706, 1537-2723. doi: 10.1080/ 700.2013..04444
[38] Mehmet Onder Efe and Hitay Ozbay. Proper orthogonal decomposition for reduced order modeling: 2d heat flow. In Proceedings of 2003 IEEE Conference on Control Applications, 2003. CCA 2003., volume 2, pages 1273 1277. IEEE,
[39] Maziar Raissi and George Em Kariaaakii. Machine Leamning of Linear Differential Equations using Gaussian Processes. Journal of Computational Physics, 348:683-693, November 2017. ISSN 00219991. doi: 10/gbzfgr.
[48] N Sugimoto. Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves. Journal of fluid mechanics, 225:631 1-653.
[49] Kai Nagel. Particle hopping models and traffic flow theory. Physical review E, 53(5):6655.
[50] Kutluay, AR Bahadir, and A Özdecs. Numerical solution of mie-dimensional burgers equation: explicit and exact-explicit finite difference methods. Journal of Computational and Applicd Mathematics, 103(2): 251-261.
[51] A. A. Shah, W. Xing, V. Triantafyllidis. Reduced-order modelling perameter-dependent, lin- ear and nonlinear dynamic partial differential equation models. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2200):20160809) April 2017. ISSN 1364-5021, 147 1-716. doi: 10. (099///pa.201..08855 0809..
[52] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations. arXiv preprint arXiv: 17771111166...
[53] S Persides. The laplace and poisson equations in schwarzschild's space-time. Journal of Mathematical Analysis and Applications, 43(3):571- 558.
[54] I.F. Lagaris, A. Likas, and D.I. Fotiadis. Artificial neural networks for solving ordinary and partial differ- ential equations. JEEE Transactions on Neural Networks, 95):98771000, Sept./1998. ISSN 10459227. doi: 10. 10.7109/7..712777.
[55] Frank Spitzer. Electrostatic capacity, heat flow, and brownian motion. fields, 3(2)1100-12.
Probability theory and related
[56] Fischer Black and Myron Scholes. The pricing of options and corporate liabilities. Journal of political economy, 81(3):637-654
[57] Wei Xing, Shireen Elhabian, Vahid Keshavarzzadeh, and Robert Kirby. Shared-Gaussian Pro- cess: Learning Interpretable Shared Hidden Structure Across Data Spaces for Design Space Analysis and Exploration. Journal of Mechanical Design, 142(8), August 2020. ISSN 1050-0472, 00.1115114044077
doi:
[58] Erik Andreassen, Anders Clausen, Mattias Schevenels, Boyan S. Lazarov, and Ole Sigmund. Efficient topology optimization in matlab using 88 lines of code. Structural and Multidisciplinary Optimization. 43(1):1-16, Jan 2.11. ISSN 1615-1488.
[59] Martin Philip Bendsoe and Ole Sigmund. Topology optimization: Theory, methods and applications. Springer, 2004.
29

---

[60] Tyler E. Bruns and Daniel A. Tortorelli. Topology optimization of non-linear elastic structures and com-- pliant mechanisms. Computer Methods in Applied Mechanics and Engineering, 190(26):3443 3459, 2001. ISSN 0045-7825.
[61] Charles-Antoine Guérin, Piere Mallet, and Anne Sentenac. Effectivv- mediuu theory for finite-size ag- gregates. JOSA A, 23(2):349-353 2006.
[62] Mani Razi, Ren Wang. Yanyan He, Robert M. Kirby, and Luca Dal Negro. Optimization of Large-Scale Vogel Spiral Arrays of Plasmonic Nanoparticles. Plasmonics, 14(1):253 261, February 2019. ISSN 1557-1955, 1557-1963. doi: 10.1007/s 1468 018-0799.y.
[63] Aristi Christofi, Felipe A Pinheiro, and Luca Dal Negro. Probing scattering resonances vogels spirals with the greens matrix spectral method. Optics letters, 1199):1933-1936, 2016.
30


Reviews:
[{"review_id": "OG66ph8EgP", "reviewer": "Reviewer_KCHb", "summary": "The paper proposes an implementable and concrete version of Li et al's infinite dimensional fidelity DE, an method of fusing simulations at different levels of fidelity/resolution, to trade off between computational tractability and statistical accuracy.", "questions": "equation (1) is central in the paper, and yet I don't understand it. Perhaps I am misreading, or perhaps a typo? in what sense is it _auto_-regressive? This seems to be a linear relation between two fidelities, $t=0,t=T$. Do you want it to relate two variable fidelities, not simply the two fidelities $i=0$ and $i=T$?\nThings look more autoregressive in eq (4) but I'm not sure I understand that either. What is $\\Delta$ doing here? Maybe this is a  notational quibble, but It is multiplying the difference between $t_T$ and $t_0$. If we let $\\Delta(t_T-t_0)\\to 0$, does that mean that $t_T\\to t_0$ *and* $\\Delta\\to 0$ simultaneously? I think it makes more sense if we delete $\\Delta$ and let $t_T\\to t_0$, esp in the light of (5)\n\nAlso, what is the relationship between times and fidelities? Does the system we are looking at have a time dimension, and each simulation run has a different time and spatial discretization? In which case, how do I interpret something like $t_T$ as in eq(4)? Shouldn't each fixed fidelity have a _sequence_ of different timesteps at which the entire simulation is evaluated?\n\nin (19) our virtual sites look a lot like the inducing points of sparse GPs. Is that how I should be interpreting them?\n\nl70:\"the system inputs of higher-fidelity are chosen to be the subset of the lower-fidelity, i.e., $\\mathbf{X}^T \\subset \\cdots \\subset \\mathbf{X}^2 \\subset \\mathbf{X}^1$.\" OK, I think I must be confused; what exactly is being subset here? If the higher-fidelity model is sampled over a denser mesh than the lower one, for example, then its data points should be, if anything, a superset of the mesh points of the lower-fidelity model. So I guess it is not mesh points; what is in a subset relation with what then?\n\nl177: I got lost trying to understand the subset selection here. Can you diagram it? I suspect this is very simple, but I just can't parse the sentence \"setting, such a requirement is not practical. Here, we derive a decomposition by introducing virtual observations $\\hat{\\mathbf{Y}}$ for each fidelity such that $\\mathbf{Y}^{(T)}$ satisfies the subset requirement for the completed set $\\left\\{\\mathbf{Y}^{(T-1)}, \\hat{\\mathbf{Y}}^{(T-1)}\\right\\} . \\check{\\mathbf{Y}}^{(T)}$ is the part of $\\mathbf{Y}^{(T)}$ that forms the subset of $\\mathbf{Y}^{(T-1)}$ (with a selection formulation $\\check{\\mathbf{X}}^{(T)}=\\mathbf{E}^{(T)} \\mathbf{X}^{(T-1)}$, where $\\mathbf{X}^{(T-1)}$ corresponds to the previous-fidelity outputs $\\left.\\mathbf{Y}^{(T-1)}\\right)$.\"\n\nGenerally are we missing something from the framing? What even is the fusion problem, if we already have a fixed library of simulations?  If I have run a high-fidelity simulation already then I have a maximally dense mesh of points for some version (with maximal $T$; Then what do I  gain by fusing it with lower-fidelity runs also? If the simulation is deterministic, which it seems to be, then we might as well just take that and go home. I thought that the fusion setting made sense in an _adaptive_ design-of-experiments setting where we might start from low fidelity model and up-sample as necessary to reduce overall uncertainty until we are satisfied. the model in this paper seems to assume a fixed set of simulation runs then pool them; but why would we do this, rather than simply throw out all but the highest-fidelity model? Is there some quantification of uncertainty we get from the lo-fi models which is not apparent at the high fidelity?\n\nBonus question about related research: The setting of the LiFiDEs, as one of the \"tractable fusion\" methods (l222) looks a lot like the \"probabilistic numerics\" setting -see e.g. https://www.probabilistic-numerics.org/research/pde/ where there are multi-resolution and meshless methods. Can you position this work in relation to that literature?", "limitations": "The authors are transparent about the limitations of the method. Could probably test in settings where the model likelihoods are poorly approximated by Gaussians. Possibly some of the example problem do this; I have not have time to check the appendices to confirm this, however.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "if i understand correctly the \"infinite-fidelity\" model of Li et al has many attractive data fusion properties for varaible-resolution simulation, but is not implementable.\nThe claim of this paper is that an this specific parameterisation to the infinite-fidelity approach can be implemented with the GP induced by a linear ODE with a GP prior over functional inputs (to the ODE), which induces a GP posterior. Some fancy work with inducing points is done to make this tractable in practice.", "weaknesses": "many small typos, and some odd phrasing that undermine my confidence in the results. See questions.\n\nThe paper seems not to be about design of experiments but rather sharing uncertainty between low and high-fidelity simulations that have already been performed, and yet lacks a justification for the informativeness of the low-fidelity simulations.\n\nI suspect this paper could be great with some typo- and bugfixes, but in the current form I hesitate recommend with confidence. There are too many confusing things to be sure I have understood the paper correctly.\n\nI think a simple diagram or two could have made this much clearer.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**“What even is the fusion problem”**\n\nWe appreciate the reviewer for providing such valuable feedback to us.\nPlease allow us to clarify the fusion problem here.\nThe goal for multi-fidelity fusion is to accurately predict the output of $f(\\mathbf{x},T)$.\nHere $f$ is the simulator, $\\mathbf{x}$ is the system input, and $T$ highest fidelity indicator. Note that $\\mathbf{x}$ is the system input (such as the attack angle for an airfoil) NOT the space locations; T is the fidelity indicator having nothing to do with time.\nTraditional methods approximate $f(\\mathbf{x},T)$ using many simulations for different $\\mathbf{x}$ at fidelity T, which is computationally expensive.\nFidelity Fusion methods, for instance AR, decomposes $f(\\mathbf{x},T)=\\beta f(\\mathbf{x},T')+u(\\mathbf{x},T)$ such that $f(\\mathbf{x},T)$ can be approximated using many simulations at low-fidelity $T'$ (to approximate $\\beta f(\\mathbf{x},T')$) and a few simulations at high-fidelity $T$ (to approximate $ u(\\mathbf{x},T)$).\n\n**“equation (1) is central…”**\n\nThank you for your question. We double-check to confirm that Eq (1) is correct. It is equivalent to Eq (4). It is called autoregression because the highest-fidelity (indicated as T) prediction relies on the prediction of the low-fidelity (indicated as 0) solution plus some residual. It is indeed a linear relationship. \n\n**“What is Δ doing here?”**\n\nThank you for your insight! Our initial thought is to use $\\Delta(t_T-t_0)$ to denote the difference between the solutions at two fidelity $t_T$ and $t_0$ as a function of the fidelity difference $(t_T-t_0)$. We now realize that ∆ is redundant here. We will remove ∆ in the revision. Thank you again!\n\n**“relationship between times and fidelities?”**\n\nWe apologize for the confusion. We denote the fidelity using factor $t$, which has nothing to do with time. We do not consider time as a variable in this work. Instead, the values at particular time stamps and space locations are recorded to form the output (QoI) vector $y$, which contains the key information of the entire simulation evolution. This is a common workaround for learning a spatial-temporal field from complex simulations [1].\n\n[1] S. Conti and A. O’Hagan, “Bayesian emulation of complex multi-output and dynamic computer models,” Journal of Statistical Planning and Inference, vol. 140, no. 3, pp. 640–651, Mar. 2010\n\n**“(19) look a lot like the inducing points of sparse GPs.”**\n\nYes, thank you for your insight! It is very similar to the sparse GPs but there are also some differences.\nThe difference: the inducing points of sparse GPs are introduced to reduce the size of the kernel matrix whereas the inducing points in our work are introduced to fulfill the subset requirement.\nThe similarities: if the inducing points of sparse GPs are assumed Gaussians, they can be integrated out as is shown in [2]; the inducing points naturally admit a Gaussian distribution because they are predictions of the low-fidelity GP. The challenge becomes how to integrate them out and also decompose the large kernel matrix using the subset structure. We will add this discussion to our revision. Thanks again for your comments.\n\n[2] M. Titsias, “Variational Learning of Inducing Variables in Sparse Gaussian Processes,” AISTATS, PMLR, Apr. 2009, pp. 567–574.\n\n**“the system inputs of higher-fidelity are chosen to be the subset of the lower-fidelity”**\n\nAs we try to clarify at the beginning, $x$ is the system inputs rather than the spatial locations. Thus, the subset setting means that the system inputs for conducting high-fidelity simulations are a subset for those for low-fidelity simulations. This is intuitive as one can normally run many low-fidelity simulations and fewer high-fidelity simulations that are crucial for some goals, such as optimization.\nAlso, note that we do not consider the mesh/resolution difference between different fidelities. This is achieved by recording values at some pre-defined spatial-temporal locations based on interpolations on the simulation results.\n\n**“I got lost trying to understand the subset selection”**\n\nWe apologize for the confusion. We will add a diagram of the subset and non-subset structure in the supplementary materials. Also see Fig.3 in the rebuttal PDF. The idea is intuitive. \nlet $\\mathbf{X}^{(t)}$, be the available system inputs for $t$ fidelity. The inputs for $t$ fidelity will be two parts: the subset part $\\check{\\mathbf{X}}^{(t)}$ contained in $\\mathbf{X}^{(t-1)}$, and the part that is not contained in $\\mathbf{X}^{(t-1)}$, denoted as $\\hat{\\mathbf{X}}^{(t-1)}$ (where the hat and superscript indicate that it is a complement set for the $t-1$ fidelity).\nTo extract these two parts, we define $\\check{\\mathbf{X}}^{(t)}= \\mathbf{E}^{(t)} \\mathbf{X}^{(t-1)}$ and $\\hat{\\mathbf{X}}^{(t-1)}= \\hat{\\mathbf{E}}^{(t)} \\mathbf{X}^{(t)}$.\n\n**\"Bonus question about related research...\"**\n\nWe are impressed by the reviewer’s knowledge. \nProbabilistic numerics and multi-fidelity fusion (as special types of surrogate models) both involve probabilistic methods to approximate or replace deterministic computations. \nHowever, the key difference is that probabilistic numerics views numerical problems as statistical inference problems and aims to provide uncertainty estimates along with the solution, whereas surrogate models provide a computationally efficient approximation to the expensive simulation.\n\nThe setting of the LiFiDEs looks like the probabilistic numerics but differs in that the differential operator is an assumed model whereas the differential operator is specified by the target PDE in probabilistic numerics.\n\nDespite that the probabilistic numerics is meshless, it still relies on solving a modified PDE at collocation points, which equivalently defines the fidelity factor t in our work. There is a potential to use our method to accelerate solving probabilistic numerics. We will add such discussions to the revision with related references of probabilistic numerics."}, {"review_id": "R5ozw6stim", "reviewer": "Reviewer_JqSS", "summary": "This paper presents a Gaussian process (GP) based multi-fidelity model that makes use of fidelity indicators. This paper extends the autoregression two-fidelity formulation to a linear fidelity differential equation. By assuming the lowest fidelity function and all the residual functions follow GP, a joint GP model of all the fidelity can be derived. To allow flexible kernel choices, the integral in the kernel function is approximated with Monte Carle samples. In the case of multi-dimensional observations, the proposed model assumes a coregionalization formulation. To further speed up inference, the proposed model relies on the assumption that the inputs of high fidelity are a subset of the inputs of low fidelity. The proposed method compared to state-of-the-art multi-fidelity methods on both synthetic and real data sets and shows significant improvement on mean prediction accuracy. \n", "questions": "- In the experiments, the authors attribute the performance improvement to the usage of fidelity indicator. However, compared to AR, two factors potentially contribute to the performance difference: the joint modeling of all the fidelity data (instead of pairwise modeling of two consecutive fidelities) and the explicit usage of fidelity indicator in the linear fidelity differential equation. I wonder which factor is more important to the performance difference.\n- A big benefit of using GP is uncertainty quantification. In the experiments, only the accuracy of mean prediction is compared. I wonder what is the performance of the proposed method in terms of test log-likehood compared to the methods like AR.\n", "limitations": "The limitation of the proposed method has not been sufficiently discussed.\n", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- This paper extends the common autoregression multi-fidelity formulation to a linear  linear fidelity differential equation, which results in a joint GP model over the observations of all the fidelity.\n- The proposed method makes an explicit assumption about the role of fidelity indicator in the model, which allows it to use this information for modeling.\n- With a sophisticated GP model, the proposed method requires less training time and works better with low data compared to the neural network based multi-fidelity method.\n- The proposed method significantly outperforms state-of-the-art multi-fidelity methods on both synthetic and real data.\n", "weaknesses": "- The modeling assumption in the linear fidelity differential equation formulation is quite restrictive, which may not be applicable for many real world problems. For example, this model is not very effective if low fidelity data is only good at certain area, i.e., the knowledge transferring factor needs to depend on input x.\n- For simplicity, the proposed method assumes that $\\beta(t)$ is a constant. This means that the knowledge transferring factor is full determined by the fidelity indicator, which may be too restrictive for the use case where the fidelity indicator only shows the order of fidelity not the relative quality.\n- The proposed method jointly models the observed data of all the fidelity under a single GP model, which does not scale well when a lot of data are available.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**\"which factor is more important to the performance difference.\"**\n\nThank you for your insightful comment! We did some investigation into the issue and we believe that the main contributing factor is joint modeling of all the fidelity data. \nWe find this by giving a 1st- and 2nd-order polynomial form for B(t) function. This can be understood as a simple transformation of the fidelity indicator t so that it can take more different values. We did not observe a significant improvement in the performance in most experiments. \nOne possible reason is that the fidelity indicator t is already a good representation of the fidelity information. We will investigate this in future work. We believe that such investigation will lead to a new understanding of fidelity information and the infinite-fidelity fusion, leading to a more effective model.\n\n**\"test log-likelihood compared to the methods like AR\"**\n\nThank you for your insightful comment! We agree that as a probabilistic model, the log-likelihood is an important metric to evaluate performance. We have included some additional results in rebuttal PDF Fig.2, which shows the negative log-likelihood (without the constant term and thus the result can be negative) of our method, ResGP, and AR.\nSurprisingly, the NLL of our method is better than ResGP and AR more significantly that the comparisons in RMSE. We believe that this is because the log-likelihood is more sensitive to the uncertainty of the prediction. Since our method is a joint learning (as discussed in the previous question), it is able to capture the uncertainty as a joint model whereas the other methods treat each fidelity separately. We will add this new finding to the supplementary materials along with some discussions.\n\n**\"this model is not very effective if low fidelity data is only good at certain areas...\"**\n\nThank you for your valuable insight! We agree that the proposed method is to some extent limited in terms of model capacity for the tradeoff for tractability. Most SOTA methods choose to sacrifice tractability by using nonlinear mapping from low-fidelity to high-fidelity. We are investigating the possibility of using more flexible non-linear mapping while maintaining tractability. One particular direction we are working on is to derive an explicit form of nonlinear mapping using equation discovery techniques (such as SINDy). \n\n**\"$\\beta(t)$ is a constant... is restrictive\"**\n\nThank you for your insightful comment! As mentioned earlier, we have tested a polynomial $\\beta(t)$ trying to improve the performance. However, the improvements are not significant. We believe this is because the fidelity indicator t is already a good representation of the fidelity information. We will investigate this in future work. \n\n**\"does not scale well when a lot of data are available.\"**\n\nIndeed, the proposed method will suffer from the scalability issue when the number of data is large. However, there are already good solutions to resolve this issue. For example, we can use inducing point-based sparse GP to reduce the computational complexity or use tensor algebra to reduce the computation provided that the data has a certain structure. We will add this discussion to the revision. Thank you for your insightful comment!\n"}, {"review_id": "QZRGByv74e", "reviewer": "Reviewer_Pa1c", "summary": "Multi-fidelity models are widely used for combining training data obtained from information sources with different degrees of precision or accuracy. More specifically, this allows for the combination of greater quantities of noisier but more cheaply-obtained examples with more faithful (but limited) data. In this work, the authors describe an extension to infinite-fidelity fusion that incorporates information contained within the fidelity indicator itself, while also mitigating issues relating to training time and complexity, as well as scalability to high-dimensional outputs. The authors also formulate a surrogate model than unifies a large selection of pre-existing multi- and single-fidelity models. Experiments on synthetic and real-world data indicate that the model obtains significant performance improvements over competing techniques, without incurring an unreasonably large speed penalty (compared to IFC).", "questions": "I have listed my concerns with the paper in the *Weaknesses* section. I encourage the authors to focus on these comments when continuing the discussion on the paper.\n\nThe contributions in this paper are insightful, and could inspire further research within the community. However, I currently have major concerns on the writing and novelty of the paper (which are very similar to a pre-existing paper that is only trivially referenced in the submission - I am personally not comfortable with this degree of overlap), as well as clarity and presentation. In view of the above, I am currently inclined towards rejecting this submission, but I look forward to reading the feedback from other reviewers as well as the author rebuttal.", "limitations": "There shouldn't be any immediate negative societal impact resulting from this work.", "rating": 3, "confidence": 3, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "- The problems investigated in this work, along with the associated solutions, are non-trivial, and the authors diligently include detailed derivations for all their contributions. \n- The improvements over IFC are well-motivated in this work, and I appreciated how there was a strong emphasis on computational complexity and training stability. Both of these are highly prized by practitioners, and I would expect the performance improvements reported here to be transferrable to other problem domains as long as the training process is stable.\n", "weaknesses": "- The paper is currently quite dense and difficult to follow at times. While I appreciate that the authors present several varied contributions here, I believe the presentation of the main paper could be improved further to highlight the key takeaways while deferring detailed derivations to the supplementary.\n- The paper bears very strong writing similarities to *GAR: Generalized Autoregression for Multi-Fidelity Fusion* by Wang et al., where some sentences are nearly copied in their entirety with only a single word replaced here and there. This is especially noticeable in the *Introduction* and *Background* sections of the paper, as well as some of the *Related Work*. The contributions themselves are different, although I am surprised that this 2022 paper is only given a cursory reference given the degree of similarity in the problem statement and experimental set-up.\n- Maybe I missed this while reading the paper, but while is the IFC method listed as IFC-GPT in the figures and tables?\n- A handful of limitations for this method are listed at the very end of the paper, but these currently come across as an afterthought. I would prefer to see additional ablation studies or synthetic examples showing specific situations where the proposed models may not work as well as expected.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**“main paper could be improved further to highlight the key takeaways”**\n\nWe agree with the reviewer. However, feel this to be a challenging task since our work is quite theoretical as it tried to revise the classic multi-fidelity autoregression and extend it to a tractable form for infinite fidelity fusion by introducing the concept of differential equations.\nIn the meanwhile, our method is highly application-oriented. It is difficult to balance technicality and accessibility. We endeavor to improve the presentation by:\n1).\tHighlighting the theoretical novelty of the proposed linear differential equation approach versus the non-linear approach at the introduction and the end of Section 3.1;\n2).\tPointing out the connection between the proposed method and GAR by the end of the introduction and Section 3.3;\n3).\tRevising our abstract to highlight the novelty of the infinite fidelity problem instead of mentioning many other less-significant aspects; and\n4).  Shortening the derivations in the main paper and moving the detailed derivations to the supplementary materials.\n\n**“The paper bears very strong writing similarities to GAR: Generalized Autoregression for Multi-Fidelity Fusion by Wang et al…”**\n\nWe appreciate the review for pointing this out and will revise our manuscript according for a better presentation. We did learn a lot from the GAR paper in formulating the motivation, problem definition, and some related work due to the very close connection between these two works---both works are essentially extensions of the foundation AR model but for different types of problems (one for infinite-fidelity problem and one for non-aligned high-dimensional problem).\nThere are also significant differences in the method and the novelty.\nTo resolve the reviewer’s concerns, we will make this very clear in the introduction and methodology section. \nWe will include discussions on the detailed connections including 1) the way to handle high-dimensional outputs and non-subset data structure and 2) the particular setting to turn GAR a special case of the proposed method.\n\n**“Maybe I missed this while reading the paper, but while is the IFC method listed as IFC-GPT in the figures and tables?”**\n\nAs presented in their original work, there are two variations of IFC, one with deep learning (IFC-ODE) and the other one with Gaussian process ODE (IFC-GPODE), which shows better results in the original manuscript. Thus, we show the results of IFC-GPODE. We will make this statement in the experimental section clearly.\nDue to our carelessness, we use the wrong name IFC-GPT (the name used in the original IFC codes) rather than IFC-GPODE. We will correct this in the revision.\n\n**\"additional ablation studies or synthetic examples\"**\n\nThank you for your suggestions.\nWe did not find particularly poor performance of our method in the many-fidelity setting when compared with other methods. \nTo understand the limitation, we follow the classic subset experiment setting with η=0.5 and reduce the number of fidelity from five to two. The results are shown in the extra PDF file Fig.1.\nIn this case, we often see that our approach is almost identical to the classic AR, losing its advantages as an infinite-fidelity fusion method.\nAlso, when the number of training data is scarce, our method does not perform better compared to other baselines.\nThere are certainly other factors that may affect the performance, such as the choice of fidelity and B(x). We will investigate this in future work. We believe that the current experimental results are sufficient to demonstrate the main advantages of the first tractable infinite-fidelity fusion method that is orders of magnitude faster than the only existing infinite-fidelity fusion method.\n\n\n**“The contributions in this paper are insightful… I have major concerns on the writing and novelty of the paper”**\n\nThank you for the valuable comments.\nWe appreciate the reviewer for seeing the actual contribution of this work, particularly on the novelty and the insight into the classic multi-fidelity fusion problem.\n\nThe multi-fidelity fusion has been an important topic in the surrogate modeling community with many real-world applications. This work contains significant novelty by proposing the first tractable infinite-fidelity fusion (with one benefit being orders of magnitude faster than the only existing intractable infinite-fidelity fusion). We believe that this is a significant contribution to the community. We will revise the manuscript to highlight the novelty of the proposed method.\nDue to the very close connection between this work and GAR, the presentation did show a certain level of overlap with GAR. We will revise the manuscript to make this very clear in the introduction and methodology sections. We will include discussions on the detailed connections.\n\nWe would like to humbly bring the reviewer’s attention to the work itself and its contribution. We believe this work is practical and useful for many research (such as uncertainty quantification and Bayesian optimization) and we will open-source the codes to benefit the community. We will do our best to improve the writing and presentation.\n"}, {"review_id": "jeBWbYq3c9", "reviewer": "Reviewer_5GR5", "summary": "The author proposes a general auto-regression model for multi-fidelity fusion. By simplifying the ODEs over the fidelity indicator in a linear form, close-form solutions can be derived. And the computational efficiency can be further improved using a rank-1 approximation. The experiment results also show superior performance of the method. ", "questions": "1. Could you also provide some fidelity interpolation and extrapolation results?\n2. Could you also give some results with much less high-fidelity samples.\n", "limitations": "The authors did not address the limitations in the paper.\n\nI do not believe there is any potential negative societal impact of their work.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. A general linear fidelity differential equation is proposed. It serves as an simplified version of IFC and generalized version of IMC.\n2. Close-form solutions provide a more efficient way to conduct multi-fidelity fusion with infinite fidelities, especially for high-dimensional problems.\n3. The effectiveness of this approach is demonstrated through the testing of both simulated and real-world data.\n4. The paper exhibits a clear structure and is straightforward to comprehend.", "weaknesses": "1. In experiment section, it seems the baseline methods are tested with default settings without fine tuning. The comparisons are not completely fair.\n2. One big benefit using the ODE formulation is to extrapolate since the ODE formulation can capture the underlying dynamics between different fidelities. In the paper, there is no discussion regarding this matter.\n3. \\eta is set 0.5 and 0.75 in all experiment settings. It might be a bit high.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Reviewer 4\n\n**Could you also provide some fidelity interpolation and extrapolation results?** \n\nWe have tested this functionality and found the interpolation results to be accurate as expected. We did not include these results due to the limitation of the space and the motivation of this work---proposing a new paradigm for infinite multi-fidelity fusion based on a tractable model.\nOur experiments aim to convey the advantages over the traditional approaches in terms of the main concern: accuracy and computational cost. We agree with the reviewer very much and will investigate more of this aspect, particularly for applications like Bayesian optimization.\n\n**\"Could you also give some results with much less high-fidelity samples.\"**\n\nThank you for your valuable suggestions.\nSuch an investigation is implicitly included in our Cost Evaluation experiment (Fig.5 in the main paper).\nThe number of training profiles (from low- to high-fidelity) is [10, 62, 33, 4, 1], [41, 19, 50, 6, 9], [37, 3, 26, 32, 6] for the first three points in Fig.5.\nWe can see that even with 1 highest-fidelity training sample, our method can still achieve a relatively good performance compared with other baselines. \nIt also highlights the importance of choosing a proper training data profile. For instance, the 3rd point always outperforms the 2nd point in RMSE despite that the 2nd point has more highest -fidelity training samples. We will add more discussions on this in the revision if the space allows.\n\n**\"\\eta is set 0.5 and 0.75 in all experiment settings. It might be a bit high.\"**\n\nThank you for the comments. We will investigate a smaller $\\eta$ in future work. \nOne challenge is that when $\\eta$ is small, the number of training data for the highest fidelity will be very small. For instance, for a 5-fidelity setting, if $\\eta$ is 0.25, the number of training data for the highest fidelity will be 1/1024 of that for the lowest fidelity. We are currently solving this issue by letting the model decide $\\eta$ automatically for each fidelity. However, we believe that is beyond the scope of this work and will leave it for future work."}]

Author rebuttal:
We sincerely appreciate the time and effort of the reviewers. The valuable comments will be absorbed into our revision. Here we supply some additional graphical information to address the reviewers' concerns. 

The first part is about additional experiment results where Fig 1 shows the classic subset experiment setting with η=0.5 with reducing the number of fidelity from five to two and Fig 2 shows the test negative log-likelihood (without the constant term and thus the result can be negative) of our method, ResGP, and AR with 5 fidelities subset experiment setting η=0.5.
The second part (Fig 3) illustrates the notation system for the subset data structure.

Meta-review: Most reviewers agreed this is a good contribution to the multi-fidelity regression literature and suggested acceptance. The authors promised a significant rewrite to address text overlap with an existing publication and to improve the exposition of the proposed method. It would be great if the authors could do this in time for the camera ready.
Decision: Accept (poster)