# Best Linear Unbiased Prediction for Multifidelity Computer Experiments

Authors: Weiyan Mu; Qiuyue Wei; Dongli Cui; Shifeng Xiong
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Recently it becomes a growing trend to study complex systems which contain multiple computer codes with different levels of accuracy, and a number of hierarchical Gaussian process models are proposed to handle such multiple-fidelity codes. This paper derives the best linear unbiased prediction for three popular classes of multiple-level Gaussian process models. The predictors all have explicit expressions at each untried point. Empirical best linear unbiased predictors are also provided by plug-in methods with generalized maximum likelihood estimators of unknown parameters.


---

## Body

## 1. Introduction

With the rapid development of computer technology, computer experiments have been widely used in engineering and science [1]. Consider a computer code with inputs x∈[0,1]d. The common statistical modeling method is to model the response of the code as(1)Yx=∑j=1pfjxβj+Zx=f′xβ+Zx,where f1(·),…,fp(·) are known regression functions, β=β1,…,βp′ is a vector of unknown regression coefficients, and Z(·) is a stationary Gaussian process on [0,1]d having zero mean, variance σz2, and correlation function R(·). Let Y=y1,…,yn′ be the responses corresponding to the design points x1,…,xn∈[0,1]d. With a known correlation function R(·), the best linear unbiased predictor (BLUP) is given by(2)Y^x0=Y^0≡f0′β^+r0′R-1Yn-Fβ^for x0∈[0,1]d, where r0=Rx0-x1,…,Rx0-xn′, R=(R(xi-xj))i,j=1,…,n, F=f(x1),…,f(xn)′, and β^=(F′R-1F)-1F′R-1Y is the generalized least squares estimator of β [1].It is a growing trend to study the complex system which contains multifidelity computer codes with different levels of accuracy. For example, a Bayesian approach is described to predict and analyze complex computer codes which can be run at different levels of sophistication [2]. A novel approach is taken to integrate data from approximate and detailed simulations to build a surrogate model that describes the relationship between output and input parameters [3]. The Bayesian hierarchical Gaussian process (BHGP) models are introduced to integrate low-accuracy and high-accuracy [4]. A class of nonstationary Gaussian process models are proposed to link the computer outputs of different mesh densities [5]. However, there are few papers that provide BLUPs for multifidelity computer experiments.The purpose of this article is to find BLUPs for multifidelity computer experiments. The structure of the article is as follows. In Section2, BLUPs for two levels of accuracy [4] are discussed. In Section 3, BLUPs for general k-level cases in autoregressive model described by Kennedy and O’Hagan [2] are illustrated. In Section 4, BLUPs for continuous level in nonstationary Gaussian process model described by Tuo, Wu, and Yu [5] are demonstrated. We present a real application in Section 5. Concluding remarks are given in Section 6.


## 2. BLUPs for Two-Level Cases

The two experiments considered in this section are named as the low-accuracy experiment (LE) and high-accuracy experiment (HE). Let GP(μ,σ2,θ) denote the Gaussian process with mean μ, variance σ2, and correlation parameters θ. Let Dl=(x1,…,xn) and Dh=(s1,…,sm) denote the design set for the LE and HE, respectively. Following Qian and Wu [4], for any xi∈Dl, the LE is described by(3)Ylxi=f′xiβ+Zxi,where Z(·)~GP(0,σZ2,θ1) and f(x)=f1x,…,fqx′ is a set of prespecified regression functions. For any si∈Dh, the HE can be described by(4)Yhsi=ρsiylsi+δsi+εsi,where the scale changes from LE to HE ρ(·)~GP(ρ0,σρ2,θ3), the location adjustment δ(·)~GP(δ0,σδ2,θ2), the measurement error ε~N(0,σε2) and Z(·), δ(·), and ε are jointly independent. Let Yl=Ylx1,…,Ylxn′ and Yh=Yls1,…,Ylsm′.Theorem 1.
Forx0∈[0,1]d, the BLUP of Yh(x0) is a′Y, where(5)a=B-1AA′B-1A-1v1+B-1-B-1AA′B-1A-1A′B-1v2,(6)A=F10ρ0F21m,v1=ρ0fx01,v2=ρ0σZ2rθ11σρ2vββ′fx0+σρ2σZ2v1+ρ02σZ2rθ12+σδ2rθ2,andB=σZ2Rθ11ρ0σZ2Rθ13ρ0σZ2Rθ13′σρ2β′Vβ+σρ2σZ2V1+ρ02σZ2Rθ12+σδ2Rθ2+σϵ2Im.Proof.
The linear predictorY^h(x0)=a1′Yl+a2′Yh+C based on training data YL and YH at an untried point x0 is unbiased for Yh(x0) provided(7)a1′F1+ρ0a2′F2=ρ0f′x0,a2′1m=1,C=0.
For any linear unbiased predictor (LUP) ofY0=Yh(x0), Y^h(x0)=a1′Yl+a2′Yh, the mean squared prediction error (MSPE) of a1′Yl+a2′Yh is E{M2}=E{(a1′Yl+a2′Yh-Y0)2}.
We have(8)EM2=σρ2a2′β′Vβa2+f′x0ββ′fx0-2a2′vββ′fx0+σZ2a1′Rθ11a1+σZ2σρ2a2′V1a2+ρ02σZ2a2′Rθ12a2+σρ2+ρ02σZ2+2ρ0σZ2a1′Rθ13a2-2ρ0σZ2a1′rθ11-2σZ2σρ2a2′v1-2ρ02σZ2a2′rθ12+σδ2a2′Rθ2a2-2a2′rθ2+1+σϵ2a2′a2+1,where V=(Rθ3si-sjfsif′(sj)), V1=(Rθ1(si-sj)Rθ3(si-sj)), Rθ12=(Rθ1(si-sj)), Rθ2=(Rθ2(si-sj)),1≤i,j≤m; v=(Rθ3x0-sif′(si)), v1=(Rθ1(x0-si)Rθ3(x0-si)), rθ12=Rθ1x0-si′, rθ2=Rθ2x0-si′, 1≤i≤m; Rθ11=(Rθ1(xi-xj)),1≤i,j≤n; Rθ13=(Rθ1(xi-sj)),1≤i≤m,1≤j≤n; rθ11=Rθ1x0-xi′, 1≤i≤n.
Thus the Lagrange multipliers can be used to solve BLUP corresponding toa1 and a2 that minimize (5) subject to(9)a1′F1+ρ0a2′F2=ρ0f′x0,a2′1m=1.
The Lagrange function is(10)L=σρ2a2′β′Vβa2+f′x0ββ′fx0-2a2′vββ′fx0+σZ2a1′Rθ11a1σZ2σρ2a2′V1a2+ρ02σZ2a2′Rθ12a2+σρ2+ρ02σZ2+2ρ0σZ2a1′Rθ13a2-2ρ0σZ2a1′rθ11-2σZ2σρ2a2′v1-2ρ02σZ2a2′rθ12+σδ2a2′Rθ2a2-2a2′rθ2+1+σϵ2a2′a2+1+2λ1′F1′a1+ρ0F2′a2-ρ0fx0+2λ21m′a2-1Let the gradient with respect to a1,a2,λ1,λ2 be zero, and we have(11)0A′ABλa=v1v2,where (12)A=F10ρ0F21mv1=ρ0fx01v2=ρ0σZ2rθ11σρ2vββ′fx0+σρ2σZ2v1+ρ02σZ2rθ12+σδ2rθ2B=σZ2Rθ11ρ0σZ2Rθ13ρ0σZ2Rθ13′σρ2β′Vβ+σρ2σZ2V1+ρ02σZ2Rθ12+σδ2Rθ2+σϵ2Im. After a small amount of algebra, the a solution gives(13)a=B-1AA′B-1A-1v1+B-1-B-1AA′B-1A-1A′B-1v2, Thus the BLUP of Yh(x0), the computer output at an untried point x0, is a′Y.Note that the parametersβ and ρ0 are often unknown in (10). We can use their estimators to replace them and get the empirical BLUP. For β, we use its maximum likelihood estimate (MLE) based on the data YL=F1β+Z1,(14)β^=F1′Rθ11-1F1-1F1′Rθ11-1YL.To estimateρ0, we first predict the computer output at DH=(s1,…,sm), say Y^L∗, based on the assumption Z1(·)~GP(h′(·)β1,σ12,b1) and training data (D1,Z1). We then estimate by MLE based on the model Yh(xi)=ρ0Y^l(xi)+δ(xi)+ϵ(xi).The log likelihood, up to an additive constant, can be written as(15)-12YH-ρ0Y^L∗-δ01m′σρ2V+σδ2Rθ2+σϵ2Im-1YH-ρ0Y^L∗-δ01m,where (16)V=Y^l12Rθ3s1-s2Y^l1Y^l2⋯Rθ3s1-smY^l1Y^lmRθ3s1-s2Y^l1Y^l2Y^l22⋯Rθ3s2-smY^l2Y^lm⋮⋮⋮Rθ3s1-smY^l1Y^lmRθ3s2-smY^l2Y^lm⋯Y^lm2,and then the MLE of ρ0 is (17)ρ^0=Q-1Y^L∗σρ2V+σδ2Rθ2+σϵ2Im-1YH-Q-1A1′B1-11m′σρ2V+σδ2Rθ2+σϵ2Im-1YH.


## 3. BLUPs form-Level Cases (k>2)

For eacht=1,…,k, Dt denotes the design set consisting of the nt points x1(t),…,xnt(t) and Z′=(Z1′,…,Zk′) denotes the output data, where Zt′=(Zt(x1(t)),…,Zt(xnt(t))) is the vector of outputs for the level t code.The autoregressive model [2] can be described by(18)Ztx=ρt-1Zt-1x+δtxt=2,…,k,where Z1~GP(h(·)′β1,σ12,b1), δt~GP(h·′βt,σt2,bt), and δt(·) is independent of Zt-1(·), …,Z1(·).Theorem 2.
The BLUP ofZk(x0) is a′Z, where(19)a=Q-1HH′Q-1H-1h+Q-1-Q-1HH′Q-1H-1H′Q-1V2,where the notation is given in the proof.Proof.
The linear predictora′Z+C is unbiased for Zs(x0) provided(20)a′H=h′,C=0.where(21)a=a1a2⋮ak,h=ρ1⋯ρs-1hx0ρ2⋯ρs-1hx0⋮hx0′),Ht=hx1thx2t⋮hxntt,H=H10⋯0ρ1H2H2⋯0⋮⋮⋮ρ1⋯ρk-1Hkρ2⋯ρs-1Hs⋯Hk.For any LUP of Z^k(x0), say a′Z, the MSPE of a′Z is (22)EM2=Ea1′Z1+a2′Z2+⋯+ak′Zk-Zsx02=∑j=1kaj′Vjaj+V0-2∑j=1kaj′Vj,0+∑j=1k-1∑l=j+1kaj′Vj,lal+∑j=1k-1∑l=j+1kal′Vl,jaj.Let At(Dj,Dl) denote the matrix of correlations between points in Dj and Dl; At(Dj) is the abbreviation for At(Dj,Dj). The notation in (22) is as follows: (23)Vj=σk2AjDj+ρj-12σj-12Aj-1Dj+ρj-12ρj-22σj-22Aj-2Dj+⋯+ρ1⋯ρj-12σ12A1Dj;V0=σk2+ρk-12σk-12+ρk-12ρk-22σk-22+⋯+ρ1⋯ρk-12σ12;Vj,l=ρj⋯ρl-1σj2AjDj,Dl+ρj-1⋯ρl-1ρj-1σj-12Aj-1Dj,Dl+ρj-2⋯ρl-1ρj-2ρj-1σj-22Aj-2Dj,Dl+⋯+ρ1⋯ρj-1ρ1⋯ρl-1σ12A1Dj,Dl;Vj,0=ρ1⋯ρk-1ρ1⋯ρj-1σ12A1Dj,x0+ρ2⋯ρk-1ρ2⋯ρj-1σ22A2Dj,x0+⋯+ρj-1ρj-1⋯ρk-1σj-12Aj-1Dj,x0+ρj⋯ρk-1σj2AkDk,x0.
The Lagrange function is(24)L=∑j=1kaj′Vjaj-2∑j=1kaj′Vj,0+∑j=1k-1∑l=j+1kaj′Vj,lal+∑j=1k-1∑l=j+1kal′Vl,jaj+V0+2λ′H′a-h=a′V1a-2a′V2+a′V3a+a′V4a+2λ′H′a-h,where(25)V1=V10⋯00V2⋯0⋮⋮⋮00⋯Vk,V2=V1,0V2,0⋮Vk,0,V3=0V1,2⋯V1,k00⋯V2,k⋮⋮⋮00⋯0,V4=00⋯0V2,10⋯0⋮⋮⋮Vk,1Vk,2⋯0.Let the gradient with respect to a,λ be zero, and we have(26)0H′HV1+V3+V4λa=hV2 After a small amount of algebra, the a solution gives(27)a=Q-1HH′Q-1H-1h+Q-1-Q-1HH′Q-1H-1H′Q-1V2.Note that the parametersρi-1 are often unknown in (27). We estimate them as follows:(a) Predict the computer output atD2=(x1(2),…,xn2(2)), say Z^1(2), based on the assumption Z1(·)~GP(h′(·)β1,σ12,b1) and training data (D1,Z1).(b) Estimate the parameterρ1, say ρ^1, by MLE based on the model Z2=ρ1Z^1(2)+δ2.(c) Predict the computer output atD3=(x1(3),…,xn3(3)), say Z^2(3), based on the model Z2(·)=ρ^1Z1(·)+δ2(·) and training data (D2,Z2).(d) Estimate other parametersρ2,…,ρk-1 iteratively.By the multivariate normal assumption based on the modelZi=ρi-1Z^i-1(i)+δi, the log likelihood (up to an additive constant) is (28)-12Zi-ρi-1Z^i-1i-Hiβi′σi2AiDi-1Zi-ρi-1Z^i-1i-Hiβi,and then the MLE of ρi-1 is(29)ρ^i-1=Q1-1Z^i-1i′AiDi-1Zi-Q1-1A1′B1-1Hi′AiDi-1Zi.The empirical BLUP can be obtained by using (29).


## 4. BLUPs for Continuous Levels

For continuous levels, we assume that as the tuning parametert gets closer to zero, the output of the computer experiment gets closer to the truth-value. Model [5] is(30)Wx,t=Vx+Zx,t,where V(x) is a stationary Gaussian process with mean function (31)EVx=f1′xβ1and covariance function(32)covVx1,Vx2=σ12Kϕ1x1,x2,Z(x,t) is a nonstationary Gaussian process with mean structure (33)EZx,t=f2′xβ2and covariance structures are (34)covZx1,t1,Zx2,t2=σ22Kϕ2x1,x2min⁡t1,t2l,and V and Z are mutually independent. It follows that the covariance function of W has the form (35)covWx1,t1,Wx2,t2=cov(Vx1,Vx2+covZx1,t1,Zx2,t2=σ12Kϕ1x1,x2+σ22Kϕ2x1,x2min⁡t1,t2l.Theorem 3.
The BLUP ofW0=W(x0,0) based on training data W=Wx1,t1,…,Wxn,tn′ is a′W, where(36)a=B-1AA′B-1A-1v+σ12B-1-B-1AA′B-1A-1A′B-1r0,where the notation is given in the proof.Proof.
Since the linear predictora′W+C is unbiased for W0, (37)∑i=1nai′f1′xi=f1′x0,∑i=1nai′f2′xi=0,C=0The Lagrange function is (38)L=σ12a′R1a-2σ12a′r0+σ12+σ22a′R2a+2λ1′F1′a-f0+2λ2′F2′a where Fi=fix1,fix2,…,fixn′,i=1,2; R1=(Kϕ1(xi,xj)), r0=(Kϕ1(xi,x0)), and R2=(Kϕ2(xi,xj)min(ti,tj)l), 1≤i,j≤n. Let the gradient with respect to a,λ1,λ2 be zero, and we have(39)00F1′00F2′F1F2σ12R1+σ22R2λ1λ2a=f00σ12r0 After a small amount of algebra, the a solution gives(40)a=B-1AA′B-1A-1v+σ12B-1-B-1AA′B-1A-1A′B-1r0,where A=(F1,F2), B=σ12R1+σ22R2, and v=f0′,0′. This completes the proof.


## 5. A Real Application

Casting is one of the most important manufacturing processes for making metal products. Nowadays computer simulations for casting processes are commonly used. In the project we consider a low-pressure die-casting process that produces a certain component of satellites and study how the process parameters influence the temperature difference at a key location on the casting product. In the simulations the input factors are four-process parameters shown in Table1. We conduct the experiment with two levels of accuracy and specify the run sizes of the high-accuracy and low-accuracy simulations as 16 and 32, respectively. The design and its corresponding responses are showed in Table 2. We then use the BLUP (5) to predict the high-accuracy responses at the points where we have only low-accuracy responses. The results are also reported in Table 2; see asterisked values. It can be seen that the prediction results are satisfactory.Table 1
Input factors in Section5.
Factor filling velocity initial pressure increase rate of pressure initial temperature ( x 1 ) ( x 2 ) ( x 3 ) ( x 4 ) Range [ 30,60 ] [ 20,40 ] [ 0.5,9.5 ] [ 725,745 ] Unit mm/s Kpa Kpa/s °CTable 2
Design, responses, and prediction (asterisked) in Section5.
Run ( x 1 ) ( x 2 ) ( x 3 ) ( x 4 ) y l y h 1 55.782 22.812 8.7971 740.938 52.55 49.40 2 48.282 20.312 1.2029 733.438 48.81 49.21 3 37.968 21.562 7.6721 734.688 39.52 40.06 4 46.407 38.438 4.0154 727.812 45.20 44.33 5 51.093 29.062 4.8596 735.938 37.54 39.10 6 57.657 37.188 8.2346 738.438 48.32 44.11 7 49.218 28.438 0.6404 743.438 50.09 49.34 8 59.532 30.312 7.1096 726.562 39.30 35.39 9 44.532 32.188 9.3596 731.562 59.11 57.23 10 37.032 35.938 2.3279 737.188 32.60 31.81 11 53.907 39.688 2.8904 739.688 37.27 33.39 12 32.343 25.312 1.7654 730.938 34.70 40.56 13 30.468 34.688 5.9846 729.062 42.56 41.52 14 35.157 24.062 3.4529 742.188 34.86 38.36 15 42.657 26.562 5.4221 725.312 36.43 40.02 16 40.782 33.438 6.5471 744.688 52.34 47.85 17 40.782 33.438 6.5471 744.688 46.45 47.21∗ 18 41.718 24.688 6.8279 742.812 43.43 41.98∗ 19 45.468 36.562 6.2654 735.312 46.24 45.88∗ 20 52.032 31.562 4.2971 744.062 36.21 37.14∗ 21 31.407 30.938 7.9529 740.312 46.43 46.41∗ 22 58.593 39.062 5.7029 732.188 39.33 40.04∗ 23 43.593 32.812 1.4846 730.312 34.98 34.34∗ 24 39.843 27.812 0.9221 737.812 45.92 45.12∗ 25 56.718 34.062 2.6096 732.812 30.39 31.35∗ 26 56.718 34.062 2.6096 732.812 40.87 41.09∗ 27 34.218 27.188 8.5154 727.188 48.60 47.36∗ 28 47.343 20.938 5.1404 736.562 41.54 42.08∗ 29 50.157 22.188 7.3904 729.688 48.09 50.98∗ 30 38.907 37.812 9.0779 739.062 44.60 45.67∗ 31 54.843 23.438 3.1721 741.562 45.55 43.55∗ 32 52.968 25.938 3.7346 728.438 43.93 44.03∗


## 6. Concluding Remarks

We have provided computable BLUPs for three classes of Gaussian process models, which are commonly used in modeling multifidelity computer experiments. A real engineering application is presented to show the effectiveness of our method. An advantage of our method is its explicit forms. The predictors by our method are consistent (or close to) with Bayesian predictors in the literature. Furthermore, our method does not require nested structure of the design [6, 7]. A limitation of our method lies in that its expression involves some unknown parameters which we need to estimate in advance. The MLE method or modified MLE methods can be used to solve this problem. Our method to derive BLUPs can also be used for other related cases, including sequential settings [8], computer experiments with both qualitative and quantitative factors [9], and sensitivity analysis [10].



---

*Source: 8525736-2018-06-07.xml*