Multivariate Gaussian process emulation for multifidelity computer models with high-dimensional spatial outputs

Multiv ariate Gaussian pro cess em ulation for m ultiﬁdelit y computer mo dels with high-dimensional spatial outputs Cyrus S. McCrimmon and Pulong Ma Departmen t of Statistics, Iow a State Universit y 2438 Osb orn Dr, Ames, IA 50011 Abstract Risk assessmen t of h urricane-driv en storm surge relies on deterministic computer mo dels that pro duce outputs o ver a large spatial domain. The surge mo dels can often b e run at a range of ﬁdelit y levels, with greater precision yielding more accurate sim ulations. Impro ved accuracy comes with a signiﬁcan t increase in computational exp ense, necessitating the developmen t of an emulator which lev erages information from the more plentiful low-ﬁdelit y outputs to provide fast and accurate predictions of high-ﬁdelit y sim ulations. T o prop erly assess the risk of storm surge o ver a geographic region at aggregated spatial resolution, an em ulator must account for spatial dep endence b et ween outputs y et remain computationally feasible for high-dimensional sim ulations. T o address this c hallenge, we exploit the autoregressiv e cokriging framew ork to dev elop tw o cross-cov ariance structures to account for spatial dep endence. One approac h uses a separable co v ariance structure with a sparse Cholesky prior for the inv erse of the cross-cov ariance matrix; the other inv olves a low-rank appro ximation via basis representations. W e demonstrate their predictive p erformance in the storm surge application and a testb ed example. Keywor ds: Autoregressiv e cokriging; Ba yesian inference; Cross-co v ariance; Storm surge; Uncertaint y quan- tiﬁcation 1 1 In tro duction In science and engineering, computer mo dels (also known as sim ulators) are us ed to simulate the behavior of ph ysical systems. Computer mo dels provide a digitized mathematical representation of a particular natural pro cess. Most often, these computer mo dels are very exp ensive to run to p erform uncertaint y quan tiﬁcation (UQ). There is extensive research in UQ that fo cuses on constructing surrogate mo dels or emulators - fast probabilistic approximations to computer mo dels ( Gramacy , 2020 ; Santner et al. , 2003 ). F or complex physical systems, sim ulators are often run at a range of ﬁdelity levels from low accuracy to high accuracy . Sim ulators at lo w-ﬁdelity level (or low accuracy) are faster to run but less accurate than sim ulators at high-ﬁdelity level (or high accuracy) b ecause of diﬀerences in the complexity of the mathematical models, sophistication of the n umerical solv ers used and the resolution of meshes ( Peherstorfer et al. , 2018 ). The motiv ating application in this pap er is a hurricane driven storm surge simulator featuring t wo ﬁdelity levels. The ADCIRC o cean circulation mo del ( Luettich and W esterink , 2004 ; W esterink et al. , 2008 ) is the primary computer mo del employ ed to study the abnormal rise in sea levels caused by storms, app earing in a num b er of applications, including the F ederal Emergency Management Agency (FEMA) ﬂo o d hazard map up dates (e.g., FEMA , 2008 ; Hesser et al. , 2013 ; Jensen et al. , 2012 ; Niedoro da et al. , 2010 ) and in supp ort of United States Army Corps of Engineers (USACE) pro jects (e.g., Cialone et al. , 2017 ; W amsley et al. , 2013 ). The ADCIRC sim ulator is more accurate when coupled with the SW AN mo del, whic h incorp orates wa ve eﬀects, but a single high-resolution sim ulation in Southw estern Florida takes around 2,000 core-hours on a high-p erformance sup ercomputer ( T o wns et al. , 2014 ). Risk assessment of storm surge requires a large num b er of mo del runs, making uncertaint y quantiﬁcation prohibitiv ely costly . As a result, coastal ﬂoo d hazard studies (e.g., FEMA , 2015 , 2017 ) often use lo w-ﬁdelit y models without wa ve eﬀects, although previous researc h (e.g., Dietrich et al. , 2010 ; Marsooli and Lin , 2018 ; Yin et al. , 2016 ) has indicated the imp ortance of using more comprehensiv e mo dels to estimate storm surge. T o address this scientiﬁc challenge, the goal of this pap er is to dev elop a multiﬁdelit y em ulator that syn thesizes information from b oth lo w-ﬁdelit y simulator (ADCIR C) and high-ﬁdelit y simulator (ADCIR C + SW AN). There is extensive research on emulating a single computer mo del using Gaussian pro cesses since the original work b y Sac ks et al. ( 1989 ); see Santner et al. ( 2003 ) and Gramacy ( 2020 ) for some background. The most widely adopted approach to em ulate m ultiﬁdelit y computer mo dels was proposed by Kennedy and O’Hagan ( 2000 ), who developed the so-called autoregressiv e cokriging mo del using a Marko v prop erty for the outputs at successive ﬁdelity levels. Several extensions to the autoregressive cokriging mo del hav e b een dev elop ed for v arious contexts ( Konomi and Karagiannis , 2021 ; Le Gratiet , 2013 ; Ma , 2020 ; Qian and W u , 2008 ). The Mark ov property is appropriate for computer models where the outputs exhibit a natural ordering b y ﬁdelit y lev el, as is the case for the aforemen tioned w ork. Important exceptions include non-rank ed outputs ( Ji et al. , 2024 ) and non-additive structures b etw een ﬁdelit y levels ( Heo and Sung , 2025 ), where Ji et al. ( 2024 ) emulate multiﬁdelit y computer mo dels that may not b e ranked from lo w to high ﬁdelit y levels, and Heo and Sung ( 2025 ) emulate the nonlinear relationship b etw een outputs at diﬀeren t ﬁdelity levels. F or the storm surge application, the ADCIRC + SW AN simulator has b een v alidated against historical storms with impro ved p erformance compared to the ADICRC simulator ( Cialone et al. , 2017 ; Dietrich et al. , 2010 , 2011 , 2012 ; FEMA , 2017 ). Our numerical stud ies also suggest that the Mark ov prop erty is appropriate for the storm surge computer mo dels. The ADCIR C and ADCIRC + SW AN mo dels generate high-dimensional outputs o ver a large spatial 2 region, producing the p eak surge elev ation at more than 9,000 spatial lo cations in the Cap e Coral region of Florida, which is part of the FEMA study region ( FEMA , 2017 ). The high-dimensional nature of the problem complicates the developmen t of an emulator. The ﬁrst c hallenge is to extend the aforementioned univ ariate approaches to accommo date high-dimensional outputs while main taining computational eﬃciency . While w e only consider a single spatial region with around 9 , 000 lo cations in the storm surge application, the full mesh of the ADCIRC computer mo del can include millions of meshes when detailed surge outputs are needed for coastal ﬂo o d hazard studies ( FEMA , 2015 , 2017 ). This high-dimensional c hallenge is also broadly applicable for a wide range of computer sim ulations in ph ysical sciences and engineering. The second c hallenge is to prop erly assess uncertaint y quantiﬁcation (UQ) at predicted outputs with a cross-cov ariance structure. Real-world storm surges often exhibit spatial dep endence, where the storm surge for one lo cation is similar to nearby lo cations. While incorp orating the cross-cov ariance structure in an emulator is not strictly necessary to calibrate UQ at individual lo cations (see Gu and Berger , 2016 ; Ma et al. , 2022a ), this do es not imply that predicted outputs at an aggregated spatial resolution will hav e prop er UQ. In fact, if the outputs exhibit p ositive correlation in a spatial region (which is the case for storm surge), ignoring such cross-co v ariance structure would lead to under co v erage for outputs at aggregated spatial resolutions; on the contrary , if the outputs exhibit negative correlation in a spatial region, ignoring the cross-cov ariance structure would lead to ov er-cov erage at aggregated spatial resolutions. In the storm surge application, decision making often needs to b e p erformed by join tly considering the surge predictions. F or instance, storm surge risk assessment o ver a spatial region could b e used to inform preparation for future storms, including the construction of defense structures ( FEMA , 2017 ). Scientiﬁc decisions sometimes need to b e made at coarser spatial locations in terms of FEMA ﬂo o d zones ev en though the computer mo del can predict storm surges at ﬁne-scale spatial lo cations. Thus, it is of great imp ortance to build an emulator that can accoun t for the output cross-cov ariance to prop erly quantify uncertain ty around predictions made at coarser spatial resolutions. In the con text of multiv ariate or high-dimensional spatial outputs, several modeling strategies hav e b een dev elop ed in the literature. The ﬁrst strategy is to ignore the diﬀerence b etw een the output domain and the physical input domain, and then adopt computationally eﬃcient Gaussian process (GP) mo dels (e.g., Gramacy and Apley , 2015 ; Katzfuss and Guinness , 2021 ) on an extended domain by augmenting input domain with the output domain. While these GP approximations hav e b een v ery successful for mo deling univ ariate outputs when the underlying cov ariance structure is parametrically sp eciﬁed, t ypically with a stationary correlation function, this strategy o versimpliﬁes the intrinsic nature of computer mo del outputs, and often there is no need for interpolation o ver the output domain. The second strategy is to assume indep endent mo dels at eac h spatial lo cation. The independence assumption allows mo del ﬁtting in parallel, but it requires estimating a large num b er of mo del parameters, which is undesirable for mo del assessment when there is a massiv e n umber of outputs. A v ariant of this strategy is to use the parallel partial em ulator ( Gu and Berger , 2016 ), which assumes that the outputs are conditionally indep endent with common correlation parameters. Ho wev er, this assumption may not b e desirable when the joint distribution of outputs o v er a spatial domain is of primary in terest. The third strategy is to consider a separable cov ariance structure ( Con ti and O’Hagan , 2010 ), where the output cross-cov ariance do es not change ov er the input space. The non-informative prior for the cross-co v ariance matrix in Conti and O’Hagan ( 2010 ) is infeasible when the num b er of mo del outputs is larger than the num b er of mo del runs, as the resulting p osterior is improp er. A remedy is to regularize the output cov ariance matrix, as w as done in Mak et al. ( 2018 ) via the graphical lasso ( F riedman et al. , 3 2008 ). The fourth strategy is to use a nonseparable mo del ( Bay arri et al. , 2007 ; Higdon et al. , 2008 ). This strategy treats outputs as a linear combination of basis functions and random weigh ts. This strategy is closely related to the idea of the linear mo del of coregionalization (LMC) ( Matheron , 1982 ; W ac kernagel , 2003 , pp. 175-176). Several v ariants for the construction of basis functions in the nonseparable mo del ha ve also b een prop osed (e.g., Guillas et al. , 2018 ; Ma et al. , 2022b ). Among the four strategies, the separable and nonseparable cross-cov ariance structures are the most app ealing for the storm surge application, where the output cross-cov ariance is of interest. Past work ( F ric ker et al. , 2013 ) has indicated that neither the separable nor th e nonseparable mo del is universally dominant, and their predictiv e p erformance can v ary by applications. In this paper, we prop ose a unifying Bay esian multiv ariate Gaussian pro cess mo deling framework for em ulating m ultiﬁdelity simulators with high-dimensional outputs using an autoregressiv e structure ( Kennedy and O’Hagan , 2000 ). The mo deling framework is implemented using b oth the separable and nonseparable approac hes describ ed ab ov e. The ﬁrst prop osed Bay esian em ulator features a separable cov ariance structure b et ween the input space and the output space, with a sparse prior on the inv erse cross-co v ariance matrix. W e refer to this separable autoregressive cokriging mo del as SEP AR Cokrig (or simply SEP) hereafter. The main idea in SEP ARCokrig is to express the joint distribution of outputs using a series of successive m ultiv ariate regression models and then enforce conditional indep endence to induce sparsit y in the inv erse cross-co v ariance (or precision) matrix (e.g., Bick el and Levina , 2008 ; Lee and Lee , 2021 ). A prior is then assigned to the elemen ts of the mo diﬁed Cholesky decomp osition of the precision matrix ( Lee and Lee , 2021 ), whic h app ear as parameters in the sequence of multiv ariate regression models, enabling computationally eﬃcient Ba yesian estimation of the inv erse cross-cov ariance matrix. Our exp erience has indicated that p osterior inference via the graphical lasso prior ( W ang , 2012 ) can b e very slow due to the p osterior sampling of the full precision matrix, and hence the graphical lasso prior is undesirable even when there is a moderate (e.g., a few thousand) n umber of outputs for a single mo del run. In the SEP ARCokrig mo del, we assign sparse Cholesky priors to the inv erse cross -co v ariance matrix at each ﬁdelity level. The resulting mo del thus exhibits a separable structure b et ween the input space and the output space and an autoregressiv e structure b etw een ﬁdelity lev el and input space. The second prop osed Bay esian emulator, which has a nonseparable c o v ariance structure, is based on the basis representation framework in Higdon et al. ( 2008 ). W e refer to this nonseparable autoregressive cokriging mo del as NONSEP ARCokrig (or simply NONSEP) hereafter. At each ﬁdelit y level, the outputs are mo deled as a linear com bination of basis vectors and random weigh ts. The random weigh ts, which are mo deled using Gaus sian pro cesses, are linked b etw een ﬁdelity levels via an autoregressive cokriging mo del so that nonseparability b etw een the input space and the output space is maintained not only within ﬁdelity lev els but also across ﬁdelity levels. This mo deling strategy uses a Marko v prop erty for the latent pro cesses across ﬁdelity levels to maintain the assumed ordering of outputs b y accuracy . Alternatively , if one jointly mo dels the random weigh ts across diﬀerent ﬁdelity levels via a multiv ariate GP mo del with a separable co v ariance structure, one actually obtains a separable cross-cov ariance in the output space, which may b e undesirable. A crucial adv antage in the prop osed autoregressive cokriging mo deling framew ork is that b oth SEP and NONSEP mo dels can allow separate p osterior inference for parameters at each ﬁdelity level. W e also show that at eac h ﬁdelity lev el the regression parameters and v ariance parameters can b e integrated out explicitly to allo w fast statistical inference. This further implies that b oth SEP and NONSEP mo dels can b e easily extended to multiple ﬁdelit y levels without incurring extra computational b ottlenecks since separate 4 mo del ﬁtting can b e done. Suc h extension could b e broadly applicable in v arious UQ applications. The rest of the paper is organized as follo ws. Section 2 presents the prop osed SEP and NONSEP mo dels and their p osterior inferences. Section 3 demonstrates the p erformance of these tw o mo dels in a testb ed example for mo deling the spread of c hemical concentration and compares their p erformance with sev eral existing mo dels. Section 4 performs the m ultiﬁdeit y em ulation for the storm surge outputs in the Cap e Coral region and compare their p erformance with existing mo dels. Section 5 concludes the pap er with a summary and additional mo deling p ersp ectives. The supplement contains technical deriv ations and additional numerical demonstrations. 2 Statistical Metho dology 2.1 Notation Setup Consider the collection of m deterministic multiﬁdelit y computer mo dels, sorted in ascending order b y the accuracy of their represen tation to the physical pro cess of interest. Eac h simulator generates N outputs for a single input x in an input space denoted by X ⊂ R d . In the storm surge application, there are m = 2 sim ulators with their outputs collected in a spatial domain in Southw est Florida, but the presented metho ds are applicable to more general settings. The computer mo del at ﬁdelity level t is run at a collection of n t inputs, denoted b y X t . W e assume that the input design is nested, meaning X 1 ⊃ X 2 ⊃ · · · ⊃ X m , as higher ﬁdelity lev el means greater computational costs for running the computer mo del. In the storm surge application, the computational demand for the low-ﬁdelit y mo del ADCIRC is several orders of magnitude smaller than the high-ﬁdelit y model ADCIRC + SW AN. The nested input design assumption can be satisﬁed b y selecting inputs using a conditional Latin hypercub e sample ( Minasn y and McBratney , 2006 ). T o simplify notation, let y t,j ( x ) denote the output ov er input x at spatial co ordinate j and ﬁdelit y level t . Let y t,j = ( y t,j ( x 1 ) , . . . , y t,j ( x n t )) ⊤ b e a vector of observed computer mo del outputs at co ordinate j and ﬁdelit y level t . Let y t ( x ) := ( y t, 1 ( x ) , y t, 2 ( x ) , . . . , y t,N ( x )) ⊤ denote a v ector of N outputs o v er input x at ﬁdelit y level t . Let Y t = [ y t, 1 , . . . , y t,N ] denote the n t × N matrix of outputs across all inputs and spatial lo cations at ﬁdelity lev el t . Let y D = { Y 1 , . . . , Y m } b e the collection of corresp onding outputs across all ﬁdelit y levels. Let z ( x ) := ( z 1 ( x ) , . . . , z N ( x )) ⊤ ∈ R N b e a random vector. W e denote z ( · ) := { z ( x ) ∈ R N : x ∈ X } ∼ G P ( µ ( · ) , r ( · , · )Σ) if z ( · ) follo ws a matrix-v ariate Gaussian pro cess in R N o ver input space X , where µ ( · ) denotes the mean function, r ( · , · ) denotes the input correlation function, and Σ denotes the N × N cross- co v ariance matrix. z ( · ) ∼ G P ( µ ( · ) , r ( · , · )Σ) if and only if for any n inputs { x 1 , . . . , x n } ⊂ X and any integer n > 0 , the n × N random matrix Z := [ z ( x 1 ) , . . . , z ( x n )] ⊤ ∼ MN n,N ( µ , R , Σ ) with µ := [ µ ( x 1 ) , . . . , µ ( x n )] ⊤ and R := [ r ( x i , x j )] i,j =1 ,...,n ; that is, Z follows a matrix-v ariate normal distribution ( Gupta and Nagar , 1999 ) with the n × N matrix of means µ and C ov ( vec ( Z )) = Σ ⊗ R . 5 2.2 Separable Autoregressiv e Cokriging 2.2.1 Mo del F orm ulation Based on the autoregressive cokriging framew ork, w e assume the following multiv ariate autoregressiv e model: y t ( x ) = γ t − 1 y t − 1 ( x ) + δ t ( x ) , t = 2 , . . . , m, δ t ( · ) ∼ G P ( h ⊤ t ( · ) β t , r ( · , · ; θ t ) Σ t ) y 1 ( · ) ∼ G P ( h ⊤ 1 ( · ) β 1 , r ( · , · ; θ 1 ) Σ 1 ) , (2.1) where δ t ( x ) accounts for unexplained v ariability in the level- t simulator outputs after linking successive mo dels via the scale discrepancy parameter γ t − 1 , and δ t ( x ) is indep endent of y 1 ( x ) , . . . , y t − 1 ( x ) . The mean functions in y 1 ( · ) , δ 2 ( · ) , . . . , δ m ( · ) are mo deled with a q t dimensional vector of presp eciﬁed basis functions, h t ( · ) , and a q t × N matrix of regression co eﬃcients, β t . The matrix-v ariate Gaussian pro cesses y 1 ( · ) and δ t ( · ) feature separable cov ariance structures ( Conti and O’Hagan , 2010 ) with an input correlation function r ( · , · ; θ t ) and N × N cross-cov ariance matrix Σ t that captures dep endence of outputs ov er the output space. The mo del ( 2.1 ) is an extension of the univ ariate autoregressive cokriging mo del ( Kennedy and O’Hagan , 2000 ) for mo de ling multiv ariate computer mo dels. W e call this mo del a separable autoregressive cokriging (SEP) mo del. A t each ﬁdelity level, w e assume that a pro duct form of correlation functions ( Sacks et al. , 1989 ) is used r ( x , x ′ ; θ ) = Q d i =1 r 0 ( x i , x ′ i | θ i ) , where θ denotes all the correlation parameters and r 0 ( · , ·| θ ) is chosen as the isotropic Matérn correlation function of the form r 0 ( u | θ ) = 2 1 − ν Γ( ν ) √ 2 ν u θ ! ν K ν √ 2 ν u θ ! , where u is the Euclidean distance, Γ( · ) is the gamma function, θ is the range parameter, K ν is the mo diﬁed Bessel function of the second kind, and ν is the smo othness parameter that con trols the mean-square diﬀer- en tiability of the random pro cess. A standard practice is to ﬁx ν at several v alues such as ν = 2 . 5 to allow closed-form formula. In our application, ν = 2 . 5 gives satisfactory results. A special case of the separable autoregressiv e mo del is to assume a diagonal structure for Σ t , resulting in the so-called parallel partial (PP) autoregressive cokriging mo del ( Ma et al. , 2022a ). This is only app ealing when the marginal predictive distribution at each spatial lo cation is needed for uncertaint y quan tiﬁcation and decision making. A graphical representation comparing the PP and SEP mo dels is given in Figure 1 . Capturing the dep endence structure b etw een outputs is critical when predictions made at multiple outputs are jointly used for decision-making suc h as in the storm surge application. When n t ≪ N , estimating the dep endence structure b etw een outputs requires regularization in either the cross-cov ariance matrix Σ t or the precision matrix Ω t := Σ − 1 t . Pourahmadi ( 1999 ) prop osed an autore- gressiv e reparameterization of the precision matrix for a multiv ariate normal random vector that provides a natural metho d to construct a sparse representation of Ω t . W e generalize this autoregressive reparametriza- tion to a matrix-v ariate normal distribution. Our form ulation is stated in the following result with its deriv ation given in Section S.1.1 of the Supplementary Material. Result 2.1. Let Z ∼ MN n,N ( M , R , Ω − 1 ) with mean M , p ositive-deﬁnite correlation matrix R and cross- 6 AAAB/nicbVDLSsNAFJ34rK2PqLhyE6xCBSlJEXVZdOOygn1AG8JkOmnHTiZhZlIMIeCvuHGhiDtx7x+480N07aTtQlsPDBzOuZd75rghJUKa5qc2N7+wuLScW8kXVtfWN/TNrYYIIo5wHQU04C0XCkwJw3VJJMWtkGPouxQ33cFF5jeHmAsSsGsZh9j2YY8RjyAoleToO0mcOknl6CYtdXwo+66X3KaHjl40y+YIxiyxJqRY3f96fR8WvmuO/tHpBijyMZOIQiHalhlKO4FcEkRxmu9EAocQDWAPtxVl0MfCTkbxU+NAKV3DC7h6TBoj9fdGAn0hYt9Vk1lEMe1l4n9eO5LemZ0QFkYSMzQ+5EXUkIGRdWF0CcdI0lgRiDhRWQ3UhxwiqRrLqxKs6S/PkkalbJ2Uj6+sYvUcjJEDu2APlIAFTkEVXIIaqAMEEnAPHsGTdqc9aM/ay3h0TpvsbIM/0N5+APMcmes= y 2 ,j ( x ) AAAB/nicbVDLSsNAFJ34rK2PqLhyE6xCBSlJEXVZdOOygn1AG8JkOmnHTiZhZlIMIeCvuHGhiDtx7x+480N07aTtQlsPDBzOuZd75rghJUKa5qc2N7+wuLScW8kXVtfWN/TNrYYIIo5wHQU04C0XCkwJw3VJJMWtkGPouxQ33cFF5jeHmAsSsGsZh9j2YY8RjyAoleToO0mcOknl6CYtdXwo+66X3KaHjl40y+YIxiyxJqRY3f96fR8WvmuO/tHpBijyMZOIQiHalhlKO4FcEkRxmu9EAocQDWAPtxVl0MfCTkbxU+NAKV3DC7h6TBoj9fdGAn0hYt9Vk1lEMe1l4n9eO5LemZ0QFkYSMzQ+5EXUkIGRdWF0CcdI0lgRiDhRWQ3UhxwiqRrLqxKs6S/PkkalbJ2Uj6+sYvUcjJEDu2APlIAFTkEVXIIaqAMEEnAPHsGTdqc9aM/ay3h0TpvsbIM/0N5+APMcmes= y 2 ,j ( x ) AAAB/HicbVDLSsNAFJ3UV219RLt0E6xCBSmJiLosunFZwT6gDWEynbRjZ5IwMymGEH/FjQtF3PYH/AN3foiunbRdaOuBgcM593LPHDekREjT/NRyS8srq2v59UJxY3NrW9/ZbYog4gg3UEAD3nahwJT4uCGJpLgdcgyZS3HLHV5lfmuEuSCBfyvjENsM9n3iEQSlkhy9FDsJO75LK10G5cD1kvv0yNHLZtWcwFgk1oyUawdf4/dR8bvu6B/dXoAihn2JKBSiY5mhtBPIJUEUp4VuJHAI0RD2cUdRHzIs7GQSPjUOldIzvICr50tjov7eSCATImaumswiinkvE//zOpH0LuyE+GEksY+mh7yIGjIwsiaMHuEYSRorAhEnKquBBpBDJFVfBVWCNf/lRdI8qVpn1dMbq1y7BFPkwR7YBxVggXNQA9egDhoAgRg8gmfwoj1oT9qr9jYdzWmznRL4A238A3rfmRo= y m,j ( x ) AAAB/HicbVDLSsNAFJ3UV219RLt0E6xCBSmJiLosunFZwT6gDWEynbRjZ5IwMymGEH/FjQtF3PYH/AN3foiunbRdaOuBgcM593LPHDekREjT/NRyS8srq2v59UJxY3NrW9/ZbYog4gg3UEAD3nahwJT4uCGJpLgdcgyZS3HLHV5lfmuEuSCBfyvjENsM9n3iEQSlkhy9FDsJO75LK10G5cD1kvv0yNHLZtWcwFgk1oyUawdf4/dR8bvu6B/dXoAihn2JKBSiY5mhtBPIJUEUp4VuJHAI0RD2cUdRHzIs7GQSPjUOldIzvICr50tjov7eSCATImaumswiinkvE//zOpH0LuyE+GEksY+mh7yIGjIwsiaMHuEYSRorAhEnKquBBpBDJFVfBVWCNf/lRdI8qVpn1dMbq1y7BFPkwR7YBxVggXNQA9egDhoAgRg8gmfwoj1oT9qr9jYdzWmznRL4A238A3rfmRo= y m,j ( x ) AAAB/nicbVDLSsNAFJ3UV219RMWVm2AVKkhJRNRl0Y3LCvYBbQiT6aQdO5mEmUkxhIC/4saFIu7EvX/gzg/RtZO2C209MHA4517umeOGlAhpmp9abm5+YXEpv1worqyuresbmw0RRBzhOgpowFsuFJgShuuSSIpbIcfQdyluuoOLzG8OMRckYNcyDrHtwx4jHkFQKsnRt5M4dRLr8CYtd3wo+66X3KYHjl4yK+YIxiyxJqRU3ft6fR8Wv2uO/tHpBijyMZOIQiHalhlKO4FcEkRxWuhEAocQDWAPtxVl0MfCTkbxU2NfKV3DC7h6TBoj9fdGAn0hYt9Vk1lEMe1l4n9eO5LemZ0QFkYSMzQ+5EXUkIGRdWF0CcdI0lgRiDhRWQ3UhxwiqRorqBKs6S/PksZRxTqpHF9Zpeo5GCMPdsAuKAMLnIIquAQ1UAcIJOAePIIn7U570J61l/FoTpvsbIE/0N5+APGJmeo= y 1 ,j ( x ) AAAB/nicbVDLSsNAFJ3UV219RMWVm2AVKkhJRNRl0Y3LCvYBbQiT6aQdO5mEmUkxhIC/4saFIu7EvX/gzg/RtZO2C209MHA4517umeOGlAhpmp9abm5+YXEpv1worqyuresbmw0RRBzhOgpowFsuFJgShuuSSIpbIcfQdyluuoOLzG8OMRckYNcyDrHtwx4jHkFQKsnRt5M4dRLr8CYtd3wo+66X3KYHjl4yK+YIxiyxJqRU3ft6fR8Wv2uO/tHpBijyMZOIQiHalhlKO4FcEkRxWuhEAocQDWAPtxVl0MfCTkbxU2NfKV3DC7h6TBoj9fdGAn0hYt9Vk1lEMe1l4n9eO5LemZ0QFkYSMzQ+5EXUkIGRdWF0CcdI0lgRiDhRWQ3UhxwiqRorqBKs6S/PksZRxTqpHF9Zpeo5GCMPdsAuKAMLnIIquAQ1UAcIJOAePIIn7U570J61l/FoTpvsbIE/0N5+APGJmeo= y 1 ,j ( x ) AAAB7XicbVA9SwNBEJ3zM8avqKXNYRCswp2I2hmwsYxgPiQ5wt7eXrJmb/fYnRNCyH+wsVDE1v9jZ+dPcXNJoYkPBh7vzTAzL0wFN+h5X87S8srq2npho7i5tb2zW9rbbxiVacrqVAmlWyExTHDJ6shRsFaqGUlCwZrh4HriNx+ZNlzJOxymLEhIT/KYU4JWanRopNB0S2Wv4uVwF4k/I+Wrb8hR65Y+O5GiWcIkUkGMafteisGIaORUsHGxkxmWEjogPda2VJKEmWCUXzt2j60SubHStiS6ufp7YkQSY4ZJaDsTgn0z703E/7x2hvFlMOIyzZBJOl0UZ8JF5U5edyOuGUUxtIRQze2tLu0TTSjagIo2BH/+5UXSOK3455WzW79cvZ+mAQU4hCM4AR8uoAo3UIM6UHiAJ3iBV0c5z86b8z5tXXJmMwfwB87HDzyikGs= ·· · AAAB7XicbVA9SwNBEJ3zM8avqKXNYRCswp2I2hmwsYxgPiQ5wt7eXrJmb/fYnRNCyH+wsVDE1v9jZ+dPcXNJoYkPBh7vzTAzL0wFN+h5X87S8srq2npho7i5tb2zW9rbbxiVacrqVAmlWyExTHDJ6shRsFaqGUlCwZrh4HriNx+ZNlzJOxymLEhIT/KYU4JWanRopNB0S2Wv4uVwF4k/I+Wrb8hR65Y+O5GiWcIkUkGMafteisGIaORUsHGxkxmWEjogPda2VJKEmWCUXzt2j60SubHStiS6ufp7YkQSY4ZJaDsTgn0z703E/7x2hvFlMOIyzZBJOl0UZ8JF5U5edyOuGUUxtIRQze2tLu0TTSjagIo2BH/+5UXSOK3455WzW79cvZ+mAQU4hCM4AR8uoAo3UIM6UHiAJ3iBV0c5z86b8z5tXXJmMwfwB87HDzyikGs= ·· · AAAB83icbVC7SgNBFL0bXzHxEbW0GYyChYRdEbUM2lhGMA9IljA7mU1GZ2aXmdlAWPIbNhaK2Nr6A/6BnR+itZNHoYkHLhzOuZd77wlizrRx3U8ns7C4tLySXc3l19Y3Ngtb2zUdJYrQKol4pBoB1pQzSauGGU4bsaJYBJzWg7vLkV/vU6VZJG/MIKa+wF3JQkawsVKr1cVC4HbqHd0O24WiW3LHQPPEm5Jief/r7b2f/660Cx+tTkQSQaUhHGvd9NzY+ClWhhFOh7lWommMyR3u0qalEguq/XR88xAdWKWDwkjZkgaN1d8TKRZaD0RgOwU2PT3rjcT/vGZiwnM/ZTJODJVksihMODIRGgWAOkxRYvjAEkwUs7ci0sMKE2NjytkQvNmX50ntuOSdlk6uvWL5AibIwi7swSF4cAZluIIKVIFADPfwCE9O4jw4z87LpDXjTGd24A+c1x9ytZW+  1 ,j AAAB83icbVC7SgNBFL0bXzHxEbW0GYyChYRdEbUM2lhGMA9IljA7mU1GZ2aXmdlAWPIbNhaK2Nr6A/6BnR+itZNHoYkHLhzOuZd77wlizrRx3U8ns7C4tLySXc3l19Y3Ngtb2zUdJYrQKol4pBoB1pQzSauGGU4bsaJYBJzWg7vLkV/vU6VZJG/MIKa+wF3JQkawsVKr1cVC4HbqHd0O24WiW3LHQPPEm5Jief/r7b2f/660Cx+tTkQSQaUhHGvd9NzY+ClWhhFOh7lWommMyR3u0qalEguq/XR88xAdWKWDwkjZkgaN1d8TKRZaD0RgOwU2PT3rjcT/vGZiwnM/ZTJODJVksihMODIRGgWAOkxRYvjAEkwUs7ci0sMKE2NjytkQvNmX50ntuOSdlk6uvWL5AibIwi7swSF4cAZluIIKVIFADPfwCE9O4jw4z87LpDXjTGd24A+c1x9ytZW+  1 ,j AAACA3icbVC7TsMwFHXKq7Q8AmywWBQkBqgShICxgoWxSPQhtVHkOE5raieR7VSqokos/AoLAwghNib+gI0PgRmn7QCUI1k+Oude3XuPFzMqlWV9GLmZ2bn5hfxiobi0vLJqrq3XZZQITGo4YpFoekgSRkNSU1Qx0owFQdxjpOH1zjO/0SdC0ii8UoOYOBx1QhpQjJSWXHOz7UXMlwOuP9juIM6Rm/IDe/966Jolq2yNAKeJPSGlys7ny1u/+FV1zfe2H+GEk1BhhqRs2VasnBQJRTEjw0I7kSRGuIc6pKVpiDiRTjq6YQh3teLDIBL6hQqO1J8dKeIyW1NXcqS68q+Xif95rUQFp05KwzhRJMTjQUHCoIpgFgj0qSBYsYEmCAuqd4W4iwTCSsdW0CHYf0+eJvXDsn1cPrq0S5UzMEYebIFtsAdscAIq4AJUQQ1gcAPuwAN4NG6Ne+PJeB6X5oxJzwb4BeP1G2GJm8g=  m  1 ,j AAACA3icbVC7TsMwFHXKq7Q8AmywWBQkBqgShICxgoWxSPQhtVHkOE5raieR7VSqokos/AoLAwghNib+gI0PgRmn7QCUI1k+Oude3XuPFzMqlWV9GLmZ2bn5hfxiobi0vLJqrq3XZZQITGo4YpFoekgSRkNSU1Qx0owFQdxjpOH1zjO/0SdC0ii8UoOYOBx1QhpQjJSWXHOz7UXMlwOuP9juIM6Rm/IDe/966Jolq2yNAKeJPSGlys7ny1u/+FV1zfe2H+GEk1BhhqRs2VasnBQJRTEjw0I7kSRGuIc6pKVpiDiRTjq6YQh3teLDIBL6hQqO1J8dKeIyW1NXcqS68q+Xif95rUQFp05KwzhRJMTjQUHCoIpgFgj0qSBYsYEmCAuqd4W4iwTCSsdW0CHYf0+eJvXDsn1cPrq0S5UzMEYebIFtsAdscAIq4AJUQQ1gcAPuwAN4NG6Ne+PJeB6X5oxJzwb4BeP1G2GJm8g=  m  1 ,j AAAB9HicbVDLSgMxFM3UV62vqks3oUWoKGVGRF0OunFZ0T6gM5RMmmnTJpkxyRSG0r8Q3LhQxK0f465/Y/pYaOuBC4dz7uXee4KYUaVte2xlVlbX1jeym7mt7Z3dvfz+QU1FicSkiiMWyUaAFGFUkKqmmpFGLAniASP1oH878esDIhWNxKNOY+Jz1BE0pBhpI/neA+1w1HJKvbPeSStftMv2FHCZOHNSdAve6fPYTSut/LfXjnDCidCYIaWajh1rf4ikppiRUc5LFIkR7qMOaRoqECfKH06PHsFjo7RhGElTQsOp+ntiiLhSKQ9MJ0e6qxa9ifif10x0eO0PqYgTTQSeLQoTBnUEJwnANpUEa5YagrCk5laIu0girE1OOROCs/jyMqmdl53L8sW9U3RvwAxZcAQKoAQccAVccAcqoAoweAIv4A28WwPr1fqwPmetGWs+cwj+wPr6AXsclD8= ⌃ 1 ( j, j ) AAAB9HicbVDLSgMxFM3UV62vqks3oUWoKGVGRF0OunFZ0T6gM5RMmmnTJpkxyRSG0r8Q3LhQxK0f465/Y/pYaOuBC4dz7uXee4KYUaVte2xlVlbX1jeym7mt7Z3dvfz+QU1FicSkiiMWyUaAFGFUkKqmmpFGLAniASP1oH878esDIhWNxKNOY+Jz1BE0pBhpI/neA+1w1HJKvbPeSStftMv2FHCZOHNSdAve6fPYTSut/LfXjnDCidCYIaWajh1rf4ikppiRUc5LFIkR7qMOaRoqECfKH06PHsFjo7RhGElTQsOp+ntiiLhSKQ9MJ0e6qxa9ifif10x0eO0PqYgTTQSeLQoTBnUEJwnANpUEa5YagrCk5laIu0girE1OOROCs/jyMqmdl53L8sW9U3RvwAxZcAQKoAQccAVccAcqoAoweAIv4A28WwPr1fqwPmetGWs+cwj+wPr6AXsclD8= ⌃ 1 ( j, j ) AAAB/3icbVC7TsMwFHV4lvIKIHVhsaiQmKoEIWCsysLYSvQhNVHkOG5r1XnIvkGqQgd+hYWBh1jZ+QI2Fr4Fp+0ALUeyfHTOvfLx8RPBFVjWl7G0vLK6tl7YKG5ube/smnv7LRWnkrImjUUsOz5RTPCINYGDYJ1EMhL6grX94VXut2+ZVDyObmCUMDck/Yj3OCWgJc8sOX4sAjUK9YUdGDAgXmaPPbNsVawJ8CKxZ6RcLTW++Uvto+6Zn04Q0zRkEVBBlOraVgJuRiRwKti46KSKJYQOSZ91NY1IyJSbTfKP8bFWAtyLpT4R4In6eyMjocoj6smQwEDNe7n4n9dNoXfpZjxKUmARnT7USwWGGOdl4IBLRkGMNCFUcp0V0wGRhIKurKhLsOe/vEhapxX7vHLWsMvVGpqigA7RETpBNrpAVXSN6qiJKLpDD+gJPRv3xqPxarxNR5eM2c4B+gPj/QcWcpnb ✓ 1 AAAB/3icbVC7TsMwFHV4lvIKIHVhsaiQmKoEIWCsysLYSvQhNVHkOG5r1XnIvkGqQgd+hYWBh1jZ+QI2Fr4Fp+0ALUeyfHTOvfLx8RPBFVjWl7G0vLK6tl7YKG5ube/smnv7LRWnkrImjUUsOz5RTPCINYGDYJ1EMhL6grX94VXut2+ZVDyObmCUMDck/Yj3OCWgJc8sOX4sAjUK9YUdGDAgXmaPPbNsVawJ8CKxZ6RcLTW++Uvto+6Zn04Q0zRkEVBBlOraVgJuRiRwKti46KSKJYQOSZ91NY1IyJSbTfKP8bFWAtyLpT4R4In6eyMjocoj6smQwEDNe7n4n9dNoXfpZjxKUmARnT7USwWGGOdl4IBLRkGMNCFUcp0V0wGRhIKurKhLsOe/vEhapxX7vHLWsMvVGpqigA7RETpBNrpAVXSN6qiJKLpDD+gJPRv3xqPxarxNR5eM2c4B+gPj/QcWcpnb ✓ 1 AAACAHicbVC7TsMwFHV4lpZHgIGBxaIgMaAqQQgYK1gYi0QfUhNFjuu0po4T2U6lKsrCr7AwgBALAyt/wMaHwIzTdoCWI1k+Oude3XuPHzMqlWV9GnPzC4tLy4WVYml1bX3D3NxqyCgRmNRxxCLR8pEkjHJSV1Qx0ooFQaHPSNPvX+Z+c0CEpBG/UcOYuCHqchpQjJSWPHPH8SPWkcNQf9DxiUJeah/dZp5ZtirWCHCW2BNSru5/vb4PSt81z/xwOhFOQsIVZkjKtm3Fyk2RUBQzkhWdRJIY4T7qkramHIVEuunogAweaKUDg0joxxUcqb87UhTKfEddGSLVk9NeLv7ntRMVnLsp5XGiCMfjQUHCoIpgngbsUEGwYkNNEBZU7wpxDwmElc6sqEOwp0+eJY3jin1aObm2y9ULMEYB7II9cAhscAaq4ArUQB1gkIF78AiejDvjwXg2Xsalc8akZxv8gfH2A1CMmqk=  1 ,j AAACAHicbVC7TsMwFHV4lpZHgIGBxaIgMaAqQQgYK1gYi0QfUhNFjuu0po4T2U6lKsrCr7AwgBALAyt/wMaHwIzTdoCWI1k+Oude3XuPHzMqlWV9GnPzC4tLy4WVYml1bX3D3NxqyCgRmNRxxCLR8pEkjHJSV1Qx0ooFQaHPSNPvX+Z+c0CEpBG/UcOYuCHqchpQjJSWPHPH8SPWkcNQf9DxiUJeah/dZp5ZtirWCHCW2BNSru5/vb4PSt81z/xwOhFOQsIVZkjKtm3Fyk2RUBQzkhWdRJIY4T7qkramHIVEuunogAweaKUDg0joxxUcqb87UhTKfEddGSLVk9NeLv7ntRMVnLsp5XGiCMfjQUHCoIpgngbsUEGwYkNNEBZU7wpxDwmElc6sqEOwp0+eJY3jin1aObm2y9ULMEYB7II9cAhscAaq4ArUQB1gkIF78AiejDvjwXg2Xsalc8akZxv8gfH2A1CMmqk=  1 ,j AAACAHicbVC7TsMwFHXKq7Q8AgwMLBYFiQFVSYWAsYKFsUj0ITVR5DhOa+o8ZDuVqigLv8LCAEIsDKz8ARsfAjNO2wFajmT56Jx7de89bsyokIbxqRUWFpeWV4qrpfLa+samvrXdElHCMWniiEW84yJBGA1JU1LJSCfmBAUuI213cJn77SHhgkbhjRzFxA5QL6Q+xUgqydF3LTdinhgF6oOWSyRy0trxbeboFaNqjAHniTkllfrB1+v7sPzdcPQPy4twEpBQYoaE6JpGLO0UcUkxI1nJSgSJER6gHukqGqKACDsdH5DBQ6V40I+4eqGEY/V3R4oCke+oKgMk+2LWy8X/vG4i/XM7pWGcSBLiySA/YVBGME8DepQTLNlIEYQ5VbtC3EccYakyK6kQzNmT50mrVjVPqyfXZqV+ASYogj2wD46ACc5AHVyBBmgCDDJwDx7Bk3anPWjP2suktKBNe3bAH2hvP1ITmqo=  2 ,j AAACAHicbVC7TsMwFHXKq7Q8AgwMLBYFiQFVSYWAsYKFsUj0ITVR5DhOa+o8ZDuVqigLv8LCAEIsDKz8ARsfAjNO2wFajmT56Jx7de89bsyokIbxqRUWFpeWV4qrpfLa+samvrXdElHCMWniiEW84yJBGA1JU1LJSCfmBAUuI213cJn77SHhgkbhjRzFxA5QL6Q+xUgqydF3LTdinhgF6oOWSyRy0trxbeboFaNqjAHniTkllfrB1+v7sPzdcPQPy4twEpBQYoaE6JpGLO0UcUkxI1nJSgSJER6gHukqGqKACDsdH5DBQ6V40I+4eqGEY/V3R4oCke+oKgMk+2LWy8X/vG4i/XM7pWGcSBLiySA/YVBGME8DepQTLNlIEYQ5VbtC3EccYakyK6kQzNmT50mrVjVPqyfXZqV+ASYogj2wD46ACc5AHVyBBmgCDDJwDx7Bk3anPWjP2suktKBNe3bAH2hvP1ITmqo=  2 ,j AAAB/3icbVDLSgMxFM3UV62vUaEbN8EiuCozRdRlqRuXLdgHtMOQSTNtaOZBckcoYxf+ihsXPnDr3i9w58ZvMdN2oa0HQg7n3EtOjhcLrsCyvozcyura+kZ+s7C1vbO7Z+4ftFSUSMqaNBKR7HhEMcFD1gQOgnViyUjgCdb2RleZ375lUvEovIFxzJyADELuc0pAS65Z7HmR6KtxoC/cgyED4qaViWuWrLI1BV4m9pyUqsXGN3+pfdRd87PXj2gSsBCoIEp1bSsGJyUSOBVsUuglisWEjsiAdTUNScCUk07zT/CJVvrYj6Q+IeCp+nsjJYHKIurJgMBQLXqZ+J/XTcC/dFIexgmwkM4e8hOBIcJZGbjPJaMgxpoQKrnOiumQSEJBV1bQJdiLX14mrUrZPi+fNexStYZmyKMjdIxOkY0uUBVdozpqIoru0AN6Qs/GvfFovBpvs9GcMd85RH9gvP8AF/eZ3A== ✓ 2 AAAB/3icbVDLSgMxFM3UV62vUaEbN8EiuCozRdRlqRuXLdgHtMOQSTNtaOZBckcoYxf+ihsXPnDr3i9w58ZvMdN2oa0HQg7n3EtOjhcLrsCyvozcyura+kZ+s7C1vbO7Z+4ftFSUSMqaNBKR7HhEMcFD1gQOgnViyUjgCdb2RleZ375lUvEovIFxzJyADELuc0pAS65Z7HmR6KtxoC/cgyED4qaViWuWrLI1BV4m9pyUqsXGN3+pfdRd87PXj2gSsBCoIEp1bSsGJyUSOBVsUuglisWEjsiAdTUNScCUk07zT/CJVvrYj6Q+IeCp+nsjJYHKIurJgMBQLXqZ+J/XTcC/dFIexgmwkM4e8hOBIcJZGbjPJaMgxpoQKrnOiumQSEJBV1bQJdiLX14mrUrZPi+fNexStYZmyKMjdIxOkY0uUBVdozpqIoru0AN6Qs/GvfFovBpvs9GcMd85RH9gvP8AF/eZ3A== ✓ 2 AAAB/XicbVDLSsNAFJ3UV62v+Ni5GSyCq5JIUZdFNy4r2Ac0IUwmk3bo5MHMjVBD8VfcuFDErf/hzr9x0mahrQeGOZxzL3Pm+KngCizr26isrK6tb1Q3a1vbO7t75v5BVyWZpKxDE5HIvk8UEzxmHeAgWD+VjES+YD1/fFP4vQcmFU/ie5ikzI3IMOYhpwS05JlHjp+IQE0ifWEHRgyIF3lm3WpYM+BlYpekjkq0PfPLCRKaRSwGKohSA9tKwc2JBE4Fm9acTLGU0DEZsoGmMYmYcvNZ+ik+1UqAw0TqEwOeqb83chKpIqCejAiM1KJXiP95gwzCKzfncZoBi+n8oTATGBJcVIEDLhkFMdGEUMl1VkxHRBIKurCaLsFe/PIy6Z437ItG865Zb12XdVTRMTpBZ8hGl6iFblEbdRBFj+gZvaI348l4Md6Nj/loxSh3DtEfGJ8/v+mVbQ== ✓ m AAAB/XicbVDLSsNAFJ3UV62v+Ni5GSyCq5JIUZdFNy4r2Ac0IUwmk3bo5MHMjVBD8VfcuFDErf/hzr9x0mahrQeGOZxzL3Pm+KngCizr26isrK6tb1Q3a1vbO7t75v5BVyWZpKxDE5HIvk8UEzxmHeAgWD+VjES+YD1/fFP4vQcmFU/ie5ikzI3IMOYhpwS05JlHjp+IQE0ifWEHRgyIF3lm3WpYM+BlYpekjkq0PfPLCRKaRSwGKohSA9tKwc2JBE4Fm9acTLGU0DEZsoGmMYmYcvNZ+ik+1UqAw0TqEwOeqb83chKpIqCejAiM1KJXiP95gwzCKzfncZoBi+n8oTATGBJcVIEDLhkFMdGEUMl1VkxHRBIKurCaLsFe/PIy6Z437ItG865Zb12XdVTRMTpBZ8hGl6iFblEbdRBFj+gZvaI348l4Md6Nj/loxSh3DtEfGJ8/v+mVbQ== ✓ m AAACAHicbVC7TsMwFHV4lpZHgIGBxaIgMaAqQQgYK1gYi0QfUhNFjuu0pnYS2U6lKsrCr7AwgBALAyt/wMaHwIzTdoCWI1k+Oude3XuPHzMqlWV9GnPzC4tLy4WVYml1bX3D3NxqyCgRmNRxxCLR8pEkjIakrqhipBULgrjPSNPvX+Z+c0CEpFF4o4YxcTnqhjSgGCkteeaO40esI4dcf9DxiUJeyo9uM88sWxVrBDhL7AkpV/e/Xt8Hpe+aZ344nQgnnIQKMyRl27Zi5aZIKIoZyYpOIkmMcB91SVvTEHEi3XR0QAYPtNKBQST0CxUcqb87UsRlvqOu5Ej15LSXi/957UQF525KwzhRJMTjQUHCoIpgngbsUEGwYkNNEBZU7wpxDwmElc6sqEOwp0+eJY3jin1aObm2y9ULMEYB7II9cAhscAaq4ArUQB1gkIF78AiejDvjwXg2Xsalc8akZxv8gfH2A6wwmuU=  m,j AAACAHicbVC7TsMwFHV4lpZHgIGBxaIgMaAqQQgYK1gYi0QfUhNFjuu0pnYS2U6lKsrCr7AwgBALAyt/wMaHwIzTdoCWI1k+Oude3XuPHzMqlWV9GnPzC4tLy4WVYml1bX3D3NxqyCgRmNRxxCLR8pEkjIakrqhipBULgrjPSNPvX+Z+c0CEpFF4o4YxcTnqhjSgGCkteeaO40esI4dcf9DxiUJeyo9uM88sWxVrBDhL7AkpV/e/Xt8Hpe+aZ344nQgnnIQKMyRl27Zi5aZIKIoZyYpOIkmMcB91SVvTEHEi3XR0QAYPtNKBQST0CxUcqb87UsRlvqOu5Ej15LSXi/957UQF525KwzhRJMTjQUHCoIpgngbsUEGwYkNNEBZU7wpxDwmElc6sqEOwp0+eJY3jin1aObm2y9ULMEYB7II9cAhscAaq4ArUQB1gkIF78AiejDvjwXg2Xsalc8akZxv8gfH2A6wwmuU=  m,j AAACBXicbVC7TsMwFHXKq5RXAIkFBosKialKEALGqiyMrUQfUhtFjuu0Vh0nsh1EFWVh4VdYOoAQKwtfwMbCt+CkHaDlSJaOz7lX997jRYxKZVlfRmFpeWV1rbhe2tjc2t4xd/daMowFJk0cslB0PCQJo5w0FVWMdCJBUOAx0vZG15nfviNC0pDfqnFEnAANOPUpRkpLrnnUC5Aaen5yn8Ie5TD/YsSSTurarlm2KlYOuEjsGSlXDxrfdFL7qLvmZ68f4jggXGGGpOzaVqScBAlFMSNpqRdLEiE8QgPS1ZSjgEgnya9I4YlW+tAPhX5cwVz93ZGgQMpx4OnKbEk572Xif143Vv6Vk1AexYpwPB3kxwyqEGaRwD4VBCs21gRhQfWuEA+RQFjp4Eo6BHv+5EXSOqvYF5Xzhl2u1sAURXAIjsEpsMElqIIbUAdNgMEDeALP4MV4NCbGq/E2LS0Ys5598AfG+w9c45wx x 2 X 1 AAACBXicbVC7TsMwFHV4lvIKILHAYFEhMVVJhYCxKgtjK9GH1EaR4zqtVceJbAdRRVlY+BWWDiDEysIXsLHwLThpB2g5kqXjc+7Vvfd4EaNSWdaXsbS8srq2Xtgobm5t7+yae/stGcYCkyYOWSg6HpKEUU6aiipGOpEgKPAYaXuj68xv3xEhachv1TgiToAGnPoUI6Ul1zzuBUgNPT+5T2GPcph/MWJJJ3UrrlmyylYOuEjsGSlVDxvfdFL7qLvmZ68f4jggXGGGpOzaVqScBAlFMSNpsRdLEiE8QgPS1ZSjgEgnya9I4alW+tAPhX5cwVz93ZGgQMpx4OnKbEk572Xif143Vv6Vk1AexYpwPB3kxwyqEGaRwD4VBCs21gRhQfWuEA+RQFjp4Io6BHv+5EXSqpTti/J5wy5Va2CKAjgCJ+AM2OASVMENqIMmwOABPIFn8GI8GhPj1Xibli4Zs54D8AfG+w9eZ5wy x 2 X 2 AAAB9HicbVDLSgMxFM3UV62vqks3oUWoKGWmiLocdOOyon1AZyiZNNOmTTJjkikMpX8huHGhiFs/xl3/xrR1oa0HLhzOuZd77wliRpW27YmVWVldW9/Ibua2tnd29/L7B3UVJRKTGo5YJJsBUoRRQWqaakaasSSIB4w0gsHN1G8MiVQ0Eg86jYnPUVfQkGKkjeR797TLUbtS6p/1T9r5ol22Z4DLxPkhRbfgnT5N3LTazn95nQgnnAiNGVKq5dix9kdIaooZGee8RJEY4QHqkpahAnGi/NHs6DE8NkoHhpE0JTScqb8nRogrlfLAdHKke2rRm4r/ea1Eh1f+iIo40UTg+aIwYVBHcJoA7FBJsGapIQhLam6FuIckwtrklDMhOIsvL5N6pexclM/vnKJ7DebIgiNQACXggEvggltQBTWAwSN4Bq/gzRpaL9a79TFvzVg/M4fgD6zPb3yllEA= ⌃ 2 ( j, j ) AAAB9HicbVDLSgMxFM3UV62vqks3oUWoKGWmiLocdOOyon1AZyiZNNOmTTJjkikMpX8huHGhiFs/xl3/xrR1oa0HLhzOuZd77wliRpW27YmVWVldW9/Ibua2tnd29/L7B3UVJRKTGo5YJJsBUoRRQWqaakaasSSIB4w0gsHN1G8MiVQ0Eg86jYnPUVfQkGKkjeR797TLUbtS6p/1T9r5ol22Z4DLxPkhRbfgnT5N3LTazn95nQgnnAiNGVKq5dix9kdIaooZGee8RJEY4QHqkpahAnGi/NHs6DE8NkoHhpE0JTScqb8nRogrlfLAdHKke2rRm4r/ea1Eh1f+iIo40UTg+aIwYVBHcJoA7FBJsGapIQhLam6FuIckwtrklDMhOIsvL5N6pexclM/vnKJ7DebIgiNQACXggEvggltQBTWAwSN4Bq/gzRpaL9a79TFvzVg/M4fgD6zPb3yllEA= ⌃ 2 ( j, j ) AAACAnicbVDLSsNAFJ3UV62vqCtxM1iEClISKeqy4MZlRfuAJoTJZNJOO5OEmYlQQtGFv+LGhSJu/Qp3/o2TtgttPXC5h3PuMPceP2FUKsv6NgpLyyura8X10sbm1vaOubvXknEqMGnimMWi4yNJGI1IU1HFSCcRBHGfkbY/vMr99j0RksbRnRolxOWoF9GQYqS05JkHjh+zQI64btC5pT2OPF4ZnA5OPLNsVa0J4CKxZ6QMZmh45pcTxDjlJFKYISm7tpUoN0NCUczIuOSkkiQID1GPdDWNECfSzSYnjOGxVgIYxkJXpOBE/f0iQ1zmW+pJjlRfznu5+J/XTVV46WY0SlJFIjz9KEwZVDHM84ABFQQrNtIEYUH1rhD3kUBY6dRKOgR7/uRF0jqr2ufV2k2tXK8/TuMogkNwBCrABhegDq5BAzQBBg/gGbyCN+PJeDHejY/paMGYRbgP/sD4/AGO25dJ ⌃ m ( j, j ) AAACAnicbVDLSsNAFJ3UV62vqCtxM1iEClISKeqy4MZlRfuAJoTJZNJOO5OEmYlQQtGFv+LGhSJu/Qp3/o2TtgttPXC5h3PuMPceP2FUKsv6NgpLyyura8X10sbm1vaOubvXknEqMGnimMWi4yNJGI1IU1HFSCcRBHGfkbY/vMr99j0RksbRnRolxOWoF9GQYqS05JkHjh+zQI64btC5pT2OPF4ZnA5OPLNsVa0J4CKxZ6QMZmh45pcTxDjlJFKYISm7tpUoN0NCUczIuOSkkiQID1GPdDWNECfSzSYnjOGxVgIYxkJXpOBE/f0iQ1zmW+pJjlRfznu5+J/XTVV46WY0SlJFIjz9KEwZVDHM84ABFQQrNtIEYUH1rhD3kUBY6dRKOgR7/uRF0jqr2ufV2k2tXK8/TuMogkNwBCrABhegDq5BAzQBBg/gGbyCN+PJeDHejY/paMGYRbgP/sD4/AGO25dJ ⌃ m ( j, j ) AAAB9HicbVDLSgMxFM3UV62vqktFgkVwUcqMiLoRim5cSQv2Ae1QMplMG5tJxiRTKEOXfoMbF4q4dd3vcOc3+BOmj4W2HrhwOOde7r3HixhV2ra/rNTC4tLySno1s7a+sbmV3d6pKhFLTCpYMCHrHlKEUU4qmmpG6pEkKPQYqXnd65Ff6xGpqOB3uh8RN0RtTgOKkTaSe3/p5JvMF1rl4W0rm7ML9hhwnjhTkivuD8vfjwfDUiv72fQFjkPCNWZIqYZjR9pNkNQUMzLINGNFIoS7qE0ahnIUEuUm46MH8MgoPgyENMU1HKu/JxIUKtUPPdMZIt1Rs95I/M9rxDq4cBPKo1gTjieLgphBLeAoAehTSbBmfUMQltTcCnEHSYS1ySljQnBmX54n1ZOCc1Y4LTu54hWYIA32wCE4Bg44B0VwA0qgAjB4AE/gBbxaPevZerPeJ60pazqzC/7A+vgBNceUzQ== j =1 ,...,N AAACBXicbVC7TsMwFHXKq5RXAIkFBosKialKEALGqiyMrUQfUhtFjuu0Vh0nsh1EFWVh4VdYOoAQKwtfwMbCt+CkHaDlSJaOz7lX997jRYxKZVlfRmFpeWV1rbhe2tjc2t4xd/daMowFJk0cslB0PCQJo5w0FVWMdCJBUOAx0vZG15nfviNC0pDfqnFEnAANOPUpRkpLrnnUC5Aaen5yn8Ie5TD/YsSSTuoGrlm2KlYOuEjsGSlXDxrfdFL7qLvmZ68f4jggXGGGpOzaVqScBAlFMSNpqRdLEiE8QgPS1ZSjgEgnya9I4YlW+tAPhX5cwVz93ZGgQMpx4OnKbEk572Xif143Vv6Vk1AexYpwPB3kxwyqEGaRwD4VBCs21gRhQfWuEA+RQFjp4Eo6BHv+5EXSOqvYF5Xzhl2u1sAURXAIjsEpsMElqIIbUAdNgMEDeALP4MV4NCbGq/E2LS0Ys5598AfG+w+305xt x 2 X m (a) PP-ARCokrig AAAB/nicbVDLSsNAFJ34rK2PqLhyE6xCBSlJEXVZdOOygn1AG8JkOmnHTiZhZlIMIeCvuHGhiDtx7x+480N07aTtQlsPDBzOuZd75rghJUKa5qc2N7+wuLScW8kXVtfWN/TNrYYIIo5wHQU04C0XCkwJw3VJJMWtkGPouxQ33cFF5jeHmAsSsGsZh9j2YY8RjyAoleToO0mcOknl6CYtdXwo+66X3KaHjl40y+YIxiyxJqRY3f96fR8WvmuO/tHpBijyMZOIQiHalhlKO4FcEkRxmu9EAocQDWAPtxVl0MfCTkbxU+NAKV3DC7h6TBoj9fdGAn0hYt9Vk1lEMe1l4n9eO5LemZ0QFkYSMzQ+5EXUkIGRdWF0CcdI0lgRiDhRWQ3UhxwiqRrLqxKs6S/PkkalbJ2Uj6+sYvUcjJEDu2APlIAFTkEVXIIaqAMEEnAPHsGTdqc9aM/ay3h0TpvsbIM/0N5+APMcmes= y 2 ,j ( x ) AAAB/nicbVDLSsNAFJ34rK2PqLhyE6xCBSlJEXVZdOOygn1AG8JkOmnHTiZhZlIMIeCvuHGhiDtx7x+480N07aTtQlsPDBzOuZd75rghJUKa5qc2N7+wuLScW8kXVtfWN/TNrYYIIo5wHQU04C0XCkwJw3VJJMWtkGPouxQ33cFF5jeHmAsSsGsZh9j2YY8RjyAoleToO0mcOknl6CYtdXwo+66X3KaHjl40y+YIxiyxJqRY3f96fR8WvmuO/tHpBijyMZOIQiHalhlKO4FcEkRxmu9EAocQDWAPtxVl0MfCTkbxU+NAKV3DC7h6TBoj9fdGAn0hYt9Vk1lEMe1l4n9eO5LemZ0QFkYSMzQ+5EXUkIGRdWF0CcdI0lgRiDhRWQ3UhxwiqRrLqxKs6S/PkkalbJ2Uj6+sYvUcjJEDu2APlIAFTkEVXIIaqAMEEnAPHsGTdqc9aM/ay3h0TpvsbIM/0N5+APMcmes= y 2 ,j ( x ) AAAB/HicbVDLSsNAFJ3UV219RLt0E6xCBSmJiLosunFZwT6gDWEynbRjZ5IwMymGEH/FjQtF3PYH/AN3foiunbRdaOuBgcM593LPHDekREjT/NRyS8srq2v59UJxY3NrW9/ZbYog4gg3UEAD3nahwJT4uCGJpLgdcgyZS3HLHV5lfmuEuSCBfyvjENsM9n3iEQSlkhy9FDsJO75LK10G5cD1kvv0yNHLZtWcwFgk1oyUawdf4/dR8bvu6B/dXoAihn2JKBSiY5mhtBPIJUEUp4VuJHAI0RD2cUdRHzIs7GQSPjUOldIzvICr50tjov7eSCATImaumswiinkvE//zOpH0LuyE+GEksY+mh7yIGjIwsiaMHuEYSRorAhEnKquBBpBDJFVfBVWCNf/lRdI8qVpn1dMbq1y7BFPkwR7YBxVggXNQA9egDhoAgRg8gmfwoj1oT9qr9jYdzWmznRL4A238A3rfmRo= y m,j ( x ) AAAB/HicbVDLSsNAFJ3UV219RLt0E6xCBSmJiLosunFZwT6gDWEynbRjZ5IwMymGEH/FjQtF3PYH/AN3foiunbRdaOuBgcM593LPHDekREjT/NRyS8srq2v59UJxY3NrW9/ZbYog4gg3UEAD3nahwJT4uCGJpLgdcgyZS3HLHV5lfmuEuSCBfyvjENsM9n3iEQSlkhy9FDsJO75LK10G5cD1kvv0yNHLZtWcwFgk1oyUawdf4/dR8bvu6B/dXoAihn2JKBSiY5mhtBPIJUEUp4VuJHAI0RD2cUdRHzIs7GQSPjUOldIzvICr50tjov7eSCATImaumswiinkvE//zOpH0LuyE+GEksY+mh7yIGjIwsiaMHuEYSRorAhEnKquBBpBDJFVfBVWCNf/lRdI8qVpn1dMbq1y7BFPkwR7YBxVggXNQA9egDhoAgRg8gmfwoj1oT9qr9jYdzWmznRL4A238A3rfmRo= y m,j ( x ) AAAB/nicbVDLSsNAFJ3UV219RMWVm2AVKkhJRNRl0Y3LCvYBbQiT6aQdO5mEmUkxhIC/4saFIu7EvX/gzg/RtZO2C209MHA4517umeOGlAhpmp9abm5+YXEpv1worqyuresbmw0RRBzhOgpowFsuFJgShuuSSIpbIcfQdyluuoOLzG8OMRckYNcyDrHtwx4jHkFQKsnRt5M4dRLr8CYtd3wo+66X3KYHjl4yK+YIxiyxJqRU3ft6fR8Wv2uO/tHpBijyMZOIQiHalhlKO4FcEkRxWuhEAocQDWAPtxVl0MfCTkbxU2NfKV3DC7h6TBoj9fdGAn0hYt9Vk1lEMe1l4n9eO5LemZ0QFkYSMzQ+5EXUkIGRdWF0CcdI0lgRiDhRWQ3UhxwiqRorqBKs6S/PksZRxTqpHF9Zpeo5GCMPdsAuKAMLnIIquAQ1UAcIJOAePIIn7U570J61l/FoTpvsbIE/0N5+APGJmeo= y 1 ,j ( x ) AAAB/nicbVDLSsNAFJ3UV219RMWVm2AVKkhJRNRl0Y3LCvYBbQiT6aQdO5mEmUkxhIC/4saFIu7EvX/gzg/RtZO2C209MHA4517umeOGlAhpmp9abm5+YXEpv1worqyuresbmw0RRBzhOgpowFsuFJgShuuSSIpbIcfQdyluuoOLzG8OMRckYNcyDrHtwx4jHkFQKsnRt5M4dRLr8CYtd3wo+66X3KYHjl4yK+YIxiyxJqRU3ft6fR8Wv2uO/tHpBijyMZOIQiHalhlKO4FcEkRxWuhEAocQDWAPtxVl0MfCTkbxU2NfKV3DC7h6TBoj9fdGAn0hYt9Vk1lEMe1l4n9eO5LemZ0QFkYSMzQ+5EXUkIGRdWF0CcdI0lgRiDhRWQ3UhxwiqRorqBKs6S/PksZRxTqpHF9Zpeo5GCMPdsAuKAMLnIIquAQ1UAcIJOAePIIn7U570J61l/FoTpvsbIE/0N5+APGJmeo= y 1 ,j ( x ) AAAB7XicbVA9SwNBEJ3zM8avqKXNYRCswp2I2hmwsYxgPiQ5wt7eXrJmb/fYnRNCyH+wsVDE1v9jZ+dPcXNJoYkPBh7vzTAzL0wFN+h5X87S8srq2npho7i5tb2zW9rbbxiVacrqVAmlWyExTHDJ6shRsFaqGUlCwZrh4HriNx+ZNlzJOxymLEhIT/KYU4JWanRopNB0S2Wv4uVwF4k/I+Wrb8hR65Y+O5GiWcIkUkGMafteisGIaORUsHGxkxmWEjogPda2VJKEmWCUXzt2j60SubHStiS6ufp7YkQSY4ZJaDsTgn0z703E/7x2hvFlMOIyzZBJOl0UZ8JF5U5edyOuGUUxtIRQze2tLu0TTSjagIo2BH/+5UXSOK3455WzW79cvZ+mAQU4hCM4AR8uoAo3UIM6UHiAJ3iBV0c5z86b8z5tXXJmMwfwB87HDzyikGs= ·· · AAAB7XicbVA9SwNBEJ3zM8avqKXNYRCswp2I2hmwsYxgPiQ5wt7eXrJmb/fYnRNCyH+wsVDE1v9jZ+dPcXNJoYkPBh7vzTAzL0wFN+h5X87S8srq2npho7i5tb2zW9rbbxiVacrqVAmlWyExTHDJ6shRsFaqGUlCwZrh4HriNx+ZNlzJOxymLEhIT/KYU4JWanRopNB0S2Wv4uVwF4k/I+Wrb8hR65Y+O5GiWcIkUkGMafteisGIaORUsHGxkxmWEjogPda2VJKEmWCUXzt2j60SubHStiS6ufp7YkQSY4ZJaDsTgn0z703E/7x2hvFlMOIyzZBJOl0UZ8JF5U5edyOuGUUxtIRQze2tLu0TTSjagIo2BH/+5UXSOK3455WzW79cvZ+mAQU4hCM4AR8uoAo3UIM6UHiAJ3iBV0c5z86b8z5tXXJmMwfwB87HDzyikGs= ·· · AAAB73icbZDLSgMxFIYz9VbrrepSkWARXJUZEXVZdOOyBXuBdihn0kwbmmTGJCOUoUtfwI0LRdz6Bn0Odz6DL2F6WWjrD4GP/z+HnHOCmDNtXPfLySwtr6yuZddzG5tb2zv53b2ajhJFaJVEPFKNADTlTNKqYYbTRqwoiIDTetC/Gef1B6o0i+SdGcTUF9CVLGQEjLUarS4IAW2vnS+4RXcivAjeDAqlw1Hl+/FoVG7nP1udiCSCSkM4aN303Nj4KSjDCKfDXCvRNAbShy5tWpQgqPbTybxDfGKdDg4jZZ80eOL+7khBaD0Qga0UYHp6Phub/2XNxIRXfspknBgqyfSjMOHYRHi8PO4wRYnhAwtAFLOzYtIDBcTYE+XsEbz5lRehdlb0LornFa9QukZTZdEBOkanyEOXqIRuURlVEUEcPaEX9OrcO8/Om/M+Lc04s5599EfOxw+gcpNq  1 AAAB73icbZDLSgMxFIYz9VbrrepSkWARXJUZEXVZdOOyBXuBdihn0kwbmmTGJCOUoUtfwI0LRdz6Bn0Odz6DL2F6WWjrD4GP/z+HnHOCmDNtXPfLySwtr6yuZddzG5tb2zv53b2ajhJFaJVEPFKNADTlTNKqYYbTRqwoiIDTetC/Gef1B6o0i+SdGcTUF9CVLGQEjLUarS4IAW2vnS+4RXcivAjeDAqlw1Hl+/FoVG7nP1udiCSCSkM4aN303Nj4KSjDCKfDXCvRNAbShy5tWpQgqPbTybxDfGKdDg4jZZ80eOL+7khBaD0Qga0UYHp6Phub/2XNxIRXfspknBgqyfSjMOHYRHi8PO4wRYnhAwtAFLOzYtIDBcTYE+XsEbz5lRehdlb0LornFa9QukZTZdEBOkanyEOXqIRuURlVEUEcPaEX9OrcO8/Om/M+Lc04s5599EfOxw+gcpNq  1 AAAB83icbVC7SgNBFL0bXzG+opY2o0GIhWFXRC2DNpYRzAOyS5idzCZDZmaXmVkhLPkGOxsLRWz9CT/Bzg+xd/IoNHrgwuGce7n3njDhTBvX/XRyC4tLyyv51cLa+sbmVnF7p6HjVBFaJzGPVSvEmnImad0ww2krURSLkNNmOLga+807qjSL5a0ZJjQQuCdZxAg2VvL9HhYCdzJx7I06xZJbcSdAf4k3I6Xqfvnr/d4/qnWKH343Jqmg0hCOtW57bmKCDCvDCKejgp9qmmAywD3atlRiQXWQTW4eoUOrdFEUK1vSoIn6cyLDQuuhCG2nwKav572x+J/XTk10EWRMJqmhkkwXRSlHJkbjAFCXKUoMH1qCiWL2VkT6WGFibEwFG4I3//Jf0jipeGeV0xuvVL2EKfKwBwdQBg/OoQrXUIM6EEjgAZ7g2UmdR+fFeZ225pzZzC78gvP2DTBDlM0=  m  1 AAAB83icbVC7SgNBFL0bXzG+opY2o0GIhWFXRC2DNpYRzAOyS5idzCZDZmaXmVkhLPkGOxsLRWz9CT/Bzg+xd/IoNHrgwuGce7n3njDhTBvX/XRyC4tLyyv51cLa+sbmVnF7p6HjVBFaJzGPVSvEmnImad0ww2krURSLkNNmOLga+807qjSL5a0ZJjQQuCdZxAg2VvL9HhYCdzJx7I06xZJbcSdAf4k3I6Xqfvnr/d4/qnWKH343Jqmg0hCOtW57bmKCDCvDCKejgp9qmmAywD3atlRiQXWQTW4eoUOrdFEUK1vSoIn6cyLDQuuhCG2nwKav572x+J/XTk10EWRMJqmhkkwXRSlHJkbjAFCXKUoMH1qCiWL2VkT6WGFibEwFG4I3//Jf0jipeGeV0xuvVL2EKfKwBwdQBg/OoQrXUIM6EEjgAZ7g2UmdR+fFeZ225pzZzC78gvP2DTBDlM0=  m  1 AAAB73icbZDLSgMxFIbP1Futt6rgxk2wCK7KjIi6LHXjskV7gXYomTTThiaZMckIZehLuHGhiFsXvoVP4M6Nz2J6WWjrD4GP/z+HnHOCmDNtXPfLySwtr6yuZddzG5tb2zv53b26jhJFaI1EPFLNAGvKmaQ1wwynzVhRLAJOG8Hgapw37qnSLJK3ZhhTX+CeZCEj2Fir2b5hPYE7XidfcIvuRGgRvBkUSgfVb/Ze/qh08p/tbkQSQaUhHGvd8tzY+ClWhhFOR7l2ommMyQD3aMuixIJqP53MO0LH1umiMFL2SYMm7u+OFAuthyKwlQKbvp7PxuZ/WSsx4aWfMhknhkoy/ShMODIRGi+PukxRYvjQAiaK2VkR6WOFibEnytkjePMrL0L9tOidF8+qXqFUhqmycAhHcAIeXEAJrqECNSDA4QGe4Nm5cx6dF+d1WppxZj378EfO2w92a5NO ⌃ 1 AAAB73icbZDLSgMxFIbP1Futt6rgxk2wCK7KjIi6LHXjskV7gXYomTTThiaZMckIZehLuHGhiFsXvoVP4M6Nz2J6WWjrD4GP/z+HnHOCmDNtXPfLySwtr6yuZddzG5tb2zv53b26jhJFaI1EPFLNAGvKmaQ1wwynzVhRLAJOG8Hgapw37qnSLJK3ZhhTX+CeZCEj2Fir2b5hPYE7XidfcIvuRGgRvBkUSgfVb/Ze/qh08p/tbkQSQaUhHGvd8tzY+ClWhhFOR7l2ommMyQD3aMuixIJqP53MO0LH1umiMFL2SYMm7u+OFAuthyKwlQKbvp7PxuZ/WSsx4aWfMhknhkoy/ShMODIRGi+PukxRYvjQAiaK2VkR6WOFibEnytkjePMrL0L9tOidF8+qXqFUhqmycAhHcAIeXEAJrqECNSDA4QGe4Nm5cx6dF+d1WppxZj378EfO2w92a5NO ⌃ 1 AAAB/3icbVC7TsMwFHV4lvIKIHVhsaiQmKoEIWCsysLYSvQhNVHkOG5r1XnIvkGqQgd+hYWBh1jZ+QI2Fr4Fp+0ALUeyfHTOvfLx8RPBFVjWl7G0vLK6tl7YKG5ube/smnv7LRWnkrImjUUsOz5RTPCINYGDYJ1EMhL6grX94VXut2+ZVDyObmCUMDck/Yj3OCWgJc8sOX4sAjUK9YUdGDAgXmaPPbNsVawJ8CKxZ6RcLTW++Uvto+6Zn04Q0zRkEVBBlOraVgJuRiRwKti46KSKJYQOSZ91NY1IyJSbTfKP8bFWAtyLpT4R4In6eyMjocoj6smQwEDNe7n4n9dNoXfpZjxKUmARnT7USwWGGOdl4IBLRkGMNCFUcp0V0wGRhIKurKhLsOe/vEhapxX7vHLWsMvVGpqigA7RETpBNrpAVXSN6qiJKLpDD+gJPRv3xqPxarxNR5eM2c4B+gPj/QcWcpnb ✓ 1 AAAB/3icbVC7TsMwFHV4lvIKIHVhsaiQmKoEIWCsysLYSvQhNVHkOG5r1XnIvkGqQgd+hYWBh1jZ+QI2Fr4Fp+0ALUeyfHTOvfLx8RPBFVjWl7G0vLK6tl7YKG5ube/smnv7LRWnkrImjUUsOz5RTPCINYGDYJ1EMhL6grX94VXut2+ZVDyObmCUMDck/Yj3OCWgJc8sOX4sAjUK9YUdGDAgXmaPPbNsVawJ8CKxZ6RcLTW++Uvto+6Zn04Q0zRkEVBBlOraVgJuRiRwKti46KSKJYQOSZ91NY1IyJSbTfKP8bFWAtyLpT4R4In6eyMjocoj6smQwEDNe7n4n9dNoXfpZjxKUmARnT7USwWGGOdl4IBLRkGMNCFUcp0V0wGRhIKurKhLsOe/vEhapxX7vHLWsMvVGpqigA7RETpBNrpAVXSN6qiJKLpDD+gJPRv3xqPxarxNR5eM2c4B+gPj/QcWcpnb ✓ 1 AAACAHicbVC7TsMwFHV4lpZHgIGBxaIgMaAqQQgYK1gYi0QfUhNFjuu0po4T2U6lKsrCr7AwgBALAyt/wMaHwIzTdoCWI1k+Oude3XuPHzMqlWV9GnPzC4tLy4WVYml1bX3D3NxqyCgRmNRxxCLR8pEkjHJSV1Qx0ooFQaHPSNPvX+Z+c0CEpBG/UcOYuCHqchpQjJSWPHPH8SPWkcNQf9DxiUJeah/dZp5ZtirWCHCW2BNSru5/vb4PSt81z/xwOhFOQsIVZkjKtm3Fyk2RUBQzkhWdRJIY4T7qkramHIVEuunogAweaKUDg0joxxUcqb87UhTKfEddGSLVk9NeLv7ntRMVnLsp5XGiCMfjQUHCoIpgngbsUEGwYkNNEBZU7wpxDwmElc6sqEOwp0+eJY3jin1aObm2y9ULMEYB7II9cAhscAaq4ArUQB1gkIF78AiejDvjwXg2Xsalc8akZxv8gfH2A1CMmqk=  1 ,j AAACAHicbVC7TsMwFHV4lpZHgIGBxaIgMaAqQQgYK1gYi0QfUhNFjuu0po4T2U6lKsrCr7AwgBALAyt/wMaHwIzTdoCWI1k+Oude3XuPHzMqlWV9GnPzC4tLy4WVYml1bX3D3NxqyCgRmNRxxCLR8pEkjHJSV1Qx0ooFQaHPSNPvX+Z+c0CEpBG/UcOYuCHqchpQjJSWPHPH8SPWkcNQf9DxiUJeah/dZp5ZtirWCHCW2BNSru5/vb4PSt81z/xwOhFOQsIVZkjKtm3Fyk2RUBQzkhWdRJIY4T7qkramHIVEuunogAweaKUDg0joxxUcqb87UhTKfEddGSLVk9NeLv7ntRMVnLsp5XGiCMfjQUHCoIpgngbsUEGwYkNNEBZU7wpxDwmElc6sqEOwp0+eJY3jin1aObm2y9ULMEYB7II9cAhscAaq4ArUQB1gkIF78AiejDvjwXg2Xsalc8akZxv8gfH2A1CMmqk=  1 ,j AAACAHicbVC7TsMwFHXKq7Q8AgwMLBYFiQFVSYWAsYKFsUj0ITVR5DhOa+o8ZDuVqigLv8LCAEIsDKz8ARsfAjNO2wFajmT56Jx7de89bsyokIbxqRUWFpeWV4qrpfLa+samvrXdElHCMWniiEW84yJBGA1JU1LJSCfmBAUuI213cJn77SHhgkbhjRzFxA5QL6Q+xUgqydF3LTdinhgF6oOWSyRy0trxbeboFaNqjAHniTkllfrB1+v7sPzdcPQPy4twEpBQYoaE6JpGLO0UcUkxI1nJSgSJER6gHukqGqKACDsdH5DBQ6V40I+4eqGEY/V3R4oCke+oKgMk+2LWy8X/vG4i/XM7pWGcSBLiySA/YVBGME8DepQTLNlIEYQ5VbtC3EccYakyK6kQzNmT50mrVjVPqyfXZqV+ASYogj2wD46ACc5AHVyBBmgCDDJwDx7Bk3anPWjP2suktKBNe3bAH2hvP1ITmqo=  2 ,j AAACAHicbVC7TsMwFHXKq7Q8AgwMLBYFiQFVSYWAsYKFsUj0ITVR5DhOa+o8ZDuVqigLv8LCAEIsDKz8ARsfAjNO2wFajmT56Jx7de89bsyokIbxqRUWFpeWV4qrpfLa+samvrXdElHCMWniiEW84yJBGA1JU1LJSCfmBAUuI213cJn77SHhgkbhjRzFxA5QL6Q+xUgqydF3LTdinhgF6oOWSyRy0trxbeboFaNqjAHniTkllfrB1+v7sPzdcPQPy4twEpBQYoaE6JpGLO0UcUkxI1nJSgSJER6gHukqGqKACDsdH5DBQ6V40I+4eqGEY/V3R4oCke+oKgMk+2LWy8X/vG4i/XM7pWGcSBLiySA/YVBGME8DepQTLNlIEYQ5VbtC3EccYakyK6kQzNmT50mrVjVPqyfXZqV+ASYogj2wD46ACc5AHVyBBmgCDDJwDx7Bk3anPWjP2suktKBNe3bAH2hvP1ITmqo=  2 ,j AAAB/3icbVDLSgMxFM3UV62vUaEbN8EiuCozRdRlqRuXLdgHtMOQSTNtaOZBckcoYxf+ihsXPnDr3i9w58ZvMdN2oa0HQg7n3EtOjhcLrsCyvozcyura+kZ+s7C1vbO7Z+4ftFSUSMqaNBKR7HhEMcFD1gQOgnViyUjgCdb2RleZ375lUvEovIFxzJyADELuc0pAS65Z7HmR6KtxoC/cgyED4qaViWuWrLI1BV4m9pyUqsXGN3+pfdRd87PXj2gSsBCoIEp1bSsGJyUSOBVsUuglisWEjsiAdTUNScCUk07zT/CJVvrYj6Q+IeCp+nsjJYHKIurJgMBQLXqZ+J/XTcC/dFIexgmwkM4e8hOBIcJZGbjPJaMgxpoQKrnOiumQSEJBV1bQJdiLX14mrUrZPi+fNexStYZmyKMjdIxOkY0uUBVdozpqIoru0AN6Qs/GvfFovBpvs9GcMd85RH9gvP8AF/eZ3A== ✓ 2 AAAB/3icbVDLSgMxFM3UV62vUaEbN8EiuCozRdRlqRuXLdgHtMOQSTNtaOZBckcoYxf+ihsXPnDr3i9w58ZvMdN2oa0HQg7n3EtOjhcLrsCyvozcyura+kZ+s7C1vbO7Z+4ftFSUSMqaNBKR7HhEMcFD1gQOgnViyUjgCdb2RleZ375lUvEovIFxzJyADELuc0pAS65Z7HmR6KtxoC/cgyED4qaViWuWrLI1BV4m9pyUqsXGN3+pfdRd87PXj2gSsBCoIEp1bSsGJyUSOBVsUuglisWEjsiAdTUNScCUk07zT/CJVvrYj6Q+IeCp+nsjJYHKIurJgMBQLXqZ+J/XTcC/dFIexgmwkM4e8hOBIcJZGbjPJaMgxpoQKrnOiumQSEJBV1bQJdiLX14mrUrZPi+fNexStYZmyKMjdIxOkY0uUBVdozpqIoru0AN6Qs/GvfFovBpvs9GcMd85RH9gvP8AF/eZ3A== ✓ 2 AAAB/XicbVDLSsNAFJ3UV62v+Ni5GSyCq5JIUZdFNy4r2Ac0IUwmk3bo5MHMjVBD8VfcuFDErf/hzr9x0mahrQeGOZxzL3Pm+KngCizr26isrK6tb1Q3a1vbO7t75v5BVyWZpKxDE5HIvk8UEzxmHeAgWD+VjES+YD1/fFP4vQcmFU/ie5ikzI3IMOYhpwS05JlHjp+IQE0ifWEHRgyIF3lm3WpYM+BlYpekjkq0PfPLCRKaRSwGKohSA9tKwc2JBE4Fm9acTLGU0DEZsoGmMYmYcvNZ+ik+1UqAw0TqEwOeqb83chKpIqCejAiM1KJXiP95gwzCKzfncZoBi+n8oTATGBJcVIEDLhkFMdGEUMl1VkxHRBIKurCaLsFe/PIy6Z437ItG865Zb12XdVTRMTpBZ8hGl6iFblEbdRBFj+gZvaI348l4Md6Nj/loxSh3DtEfGJ8/v+mVbQ== ✓ m AAAB/XicbVDLSsNAFJ3UV62v+Ni5GSyCq5JIUZdFNy4r2Ac0IUwmk3bo5MHMjVBD8VfcuFDErf/hzr9x0mahrQeGOZxzL3Pm+KngCizr26isrK6tb1Q3a1vbO7t75v5BVyWZpKxDE5HIvk8UEzxmHeAgWD+VjES+YD1/fFP4vQcmFU/ie5ikzI3IMOYhpwS05JlHjp+IQE0ifWEHRgyIF3lm3WpYM+BlYpekjkq0PfPLCRKaRSwGKohSA9tKwc2JBE4Fm9acTLGU0DEZsoGmMYmYcvNZ+ik+1UqAw0TqEwOeqb83chKpIqCejAiM1KJXiP95gwzCKzfncZoBi+n8oTATGBJcVIEDLhkFMdGEUMl1VkxHRBIKurCaLsFe/PIy6Z437ItG865Zb12XdVTRMTpBZ8hGl6iFblEbdRBFj+gZvaI348l4Md6Nj/loxSh3DtEfGJ8/v+mVbQ== ✓ m AAACAHicbVC7TsMwFHV4lpZHgIGBxaIgMaAqQQgYK1gYi0QfUhNFjuu0pnYS2U6lKsrCr7AwgBALAyt/wMaHwIzTdoCWI1k+Oude3XuPHzMqlWV9GnPzC4tLy4WVYml1bX3D3NxqyCgRmNRxxCLR8pEkjIakrqhipBULgrjPSNPvX+Z+c0CEpFF4o4YxcTnqhjSgGCkteeaO40esI4dcf9DxiUJeyo9uM88sWxVrBDhL7AkpV/e/Xt8Hpe+aZ344nQgnnIQKMyRl27Zi5aZIKIoZyYpOIkmMcB91SVvTEHEi3XR0QAYPtNKBQST0CxUcqb87UsRlvqOu5Ej15LSXi/957UQF525KwzhRJMTjQUHCoIpgngbsUEGwYkNNEBZU7wpxDwmElc6sqEOwp0+eJY3jin1aObm2y9ULMEYB7II9cAhscAaq4ArUQB1gkIF78AiejDvjwXg2Xsalc8akZxv8gfH2A6wwmuU=  m,j AAACAHicbVC7TsMwFHV4lpZHgIGBxaIgMaAqQQgYK1gYi0QfUhNFjuu0pnYS2U6lKsrCr7AwgBALAyt/wMaHwIzTdoCWI1k+Oude3XuPHzMqlWV9GnPzC4tLy4WVYml1bX3D3NxqyCgRmNRxxCLR8pEkjIakrqhipBULgrjPSNPvX+Z+c0CEpFF4o4YxcTnqhjSgGCkteeaO40esI4dcf9DxiUJeyo9uM88sWxVrBDhL7AkpV/e/Xt8Hpe+aZ344nQgnnIQKMyRl27Zi5aZIKIoZyYpOIkmMcB91SVvTEHEi3XR0QAYPtNKBQST0CxUcqb87UsRlvqOu5Ej15LSXi/957UQF525KwzhRJMTjQUHCoIpgngbsUEGwYkNNEBZU7wpxDwmElc6sqEOwp0+eJY3jin1aObm2y9ULMEYB7II9cAhscAaq4ArUQB1gkIF78AiejDvjwXg2Xsalc8akZxv8gfH2A6wwmuU=  m,j AAACBXicbVC7TsMwFHXKq5RXAIkFBosKialKEALGqiyMrUQfUhtFjuu0Vh0nsh1EFWVh4VdYOoAQKwtfwMbCt+CkHaDlSJaOz7lX997jRYxKZVlfRmFpeWV1rbhe2tjc2t4xd/daMowFJk0cslB0PCQJo5w0FVWMdCJBUOAx0vZG15nfviNC0pDfqnFEnAANOPUpRkpLrnnUC5Aaen5yn8Ie5TD/YsSSTurarlm2KlYOuEjsGSlXDxrfdFL7qLvmZ68f4jggXGGGpOzaVqScBAlFMSNpqRdLEiE8QgPS1ZSjgEgnya9I4YlW+tAPhX5cwVz93ZGgQMpx4OnKbEk572Xif143Vv6Vk1AexYpwPB3kxwyqEGaRwD4VBCs21gRhQfWuEA+RQFjp4Eo6BHv+5EXSOqvYF5Xzhl2u1sAURXAIjsEpsMElqIIbUAdNgMEDeALP4MV4NCbGq/E2LS0Ys5598AfG+w9c45wx x 2 X 1 AAACBXicbVC7TsMwFHV4lvIKILHAYFEhMVVJhYCxKgtjK9GH1EaR4zqtVceJbAdRRVlY+BWWDiDEysIXsLHwLThpB2g5kqXjc+7Vvfd4EaNSWdaXsbS8srq2Xtgobm5t7+yae/stGcYCkyYOWSg6HpKEUU6aiipGOpEgKPAYaXuj68xv3xEhachv1TgiToAGnPoUI6Ul1zzuBUgNPT+5T2GPcph/MWJJJ3UrrlmyylYOuEjsGSlVDxvfdFL7qLvmZ68f4jggXGGGpOzaVqScBAlFMSNpsRdLEiE8QgPS1ZSjgEgnya9I4alW+tAPhX5cwVz93ZGgQMpx4OnKbEk572Xif143Vv6Vk1AexYpwPB3kxwyqEGaRwD4VBCs21gRhQfWuEA+RQFjp4Io6BHv+5EXSqpTti/J5wy5Va2CKAjgCJ+AM2OASVMENqIMmwOABPIFn8GI8GhPj1Xibli4Zs54D8AfG+w9eZ5wy x 2 X 2 AAAB73icbZDLSgMxFIbPeK31VhXcuAkWwVWZKaIuS924bNFeoB1KJs20oUlmTDJCGfoSblwo4taFb+ETuHPjs5heFtr6Q+Dj/88h55wg5kwb1/1ylpZXVtfWMxvZza3tnd3c3n5dR4kitEYiHqlmgDXlTNKaYYbTZqwoFgGnjWBwNc4b91RpFslbM4ypL3BPspARbKzVbN+wnsCdYieXdwvuRGgRvBnkS4fVb/Ze/qh0cp/tbkQSQaUhHGvd8tzY+ClWhhFOR9l2ommMyQD3aMuixIJqP53MO0In1umiMFL2SYMm7u+OFAuthyKwlQKbvp7PxuZ/WSsx4aWfMhknhkoy/ShMODIRGi+PukxRYvjQAiaK2VkR6WOFibEnytojePMrL0K9WPDOC2dVL18qw1QZOIJjOAUPLqAE11CBGhDg8ABP8OzcOY/Oi/M6LV1yZj0H8EfO2w9375NP ⌃ 2 AAAB73icbZDLSgMxFIbPeK31VhXcuAkWwVWZKaIuS924bNFeoB1KJs20oUlmTDJCGfoSblwo4taFb+ETuHPjs5heFtr6Q+Dj/88h55wg5kwb1/1ylpZXVtfWMxvZza3tnd3c3n5dR4kitEYiHqlmgDXlTNKaYYbTZqwoFgGnjWBwNc4b91RpFslbM4ypL3BPspARbKzVbN+wnsCdYieXdwvuRGgRvBnkS4fVb/Ze/qh0cp/tbkQSQaUhHGvd8tzY+ClWhhFOR9l2ommMyQD3aMuixIJqP53MO0In1umiMFL2SYMm7u+OFAuthyKwlQKbvp7PxuZ/WSsx4aWfMhknhkoy/ShMODIRGi+PukxRYvjQAiaK2VkR6WOFibEnytojePMrL0K9WPDOC2dVL18qw1QZOIJjOAUPLqAE11CBGhDg8ABP8OzcOY/Oi/M6LV1yZj0H8EfO2w9375NP ⌃ 2 AAAB/XicbVA7T8MwGHR4lvIKj43FokJiqhJUAWMFC2MR9CE1UeQ4TmvVdiLbQSpRxV9hYQAhVv4HG/8Gp80ALSdZPt19n3y+MGVUacf5tpaWV1bX1isb1c2t7Z1de2+/o5JMYtLGCUtkL0SKMCpIW1PNSC+VBPGQkW44ui787gORiibiXo9T4nM0EDSmGGkjBfahFyYsUmNuLujd0QFHAQ/smlN3poCLxC1JDZRoBfaXFyU440RozJBSfddJtZ8jqSlmZFL1MkVShEdoQPqGCsSJ8vNp+gk8MUoE40SaIzScqr83csRVEdBMcqSHat4rxP+8fqbjSz+nIs00EXj2UJwxqBNYVAEjKgnWbGwIwpKarBAPkURYm8KqpgR3/suLpHNWd8/rjdtGrXlV1lEBR+AYnAIXXIAmuAEt0AYYPIJn8ArerCfrxXq3PmajS1a5cwD+wPr8AYcHlUg= ⌃ m AAAB/XicbVA7T8MwGHR4lvIKj43FokJiqhJUAWMFC2MR9CE1UeQ4TmvVdiLbQSpRxV9hYQAhVv4HG/8Gp80ALSdZPt19n3y+MGVUacf5tpaWV1bX1isb1c2t7Z1de2+/o5JMYtLGCUtkL0SKMCpIW1PNSC+VBPGQkW44ui787gORiibiXo9T4nM0EDSmGGkjBfahFyYsUmNuLujd0QFHAQ/smlN3poCLxC1JDZRoBfaXFyU440RozJBSfddJtZ8jqSlmZFL1MkVShEdoQPqGCsSJ8vNp+gk8MUoE40SaIzScqr83csRVEdBMcqSHat4rxP+8fqbjSz+nIs00EXj2UJwxqBNYVAEjKgnWbGwIwpKarBAPkURYm8KqpgR3/suLpHNWd8/rjdtGrXlV1lEBR+AYnAIXXIAmuAEt0AYYPIJn8ArerCfrxXq3PmajS1a5cwD+wPr8AYcHlUg= ⌃ m AAAB9HicbVDLSgMxFM3UV62vqktFgkVwUcqMiLoRim5cSQv2Ae1QMplMG5tJxiRTKEOXfoMbF4q4dd3vcOc3+BOmj4W2HrhwOOde7r3HixhV2ra/rNTC4tLySno1s7a+sbmV3d6pKhFLTCpYMCHrHlKEUU4qmmpG6pEkKPQYqXnd65Ff6xGpqOB3uh8RN0RtTgOKkTaSe3/p5JvMF1rl4W0rm7ML9hhwnjhTkivuD8vfjwfDUiv72fQFjkPCNWZIqYZjR9pNkNQUMzLINGNFIoS7qE0ahnIUEuUm46MH8MgoPgyENMU1HKu/JxIUKtUPPdMZIt1Rs95I/M9rxDq4cBPKo1gTjieLgphBLeAoAehTSbBmfUMQltTcCnEHSYS1ySljQnBmX54n1ZOCc1Y4LTu54hWYIA32wCE4Bg44B0VwA0qgAjB4AE/gBbxaPevZerPeJ60pazqzC/7A+vgBNceUzQ== j =1 ,...,N AAACBXicbVC7TsMwFHXKq5RXAIkFBosKialKEALGqiyMrUQfUhtFjuu0Vh0nsh1EFWVh4VdYOoAQKwtfwMbCt+CkHaDlSJaOz7lX997jRYxKZVlfRmFpeWV1rbhe2tjc2t4xd/daMowFJk0cslB0PCQJo5w0FVWMdCJBUOAx0vZG15nfviNC0pDfqnFEnAANOPUpRkpLrnnUC5Aaen5yn8Ie5TD/YsSSTuoGrlm2KlYOuEjsGSlXDxrfdFL7qLvmZ68f4jggXGGGpOzaVqScBAlFMSNpqRdLEiE8QgPS1ZSjgEgnya9I4YlW+tAPhX5cwVz93ZGgQMpx4OnKbEk572Xif143Vv6Vk1AexYpwPB3kxwyqEGaRwD4VBCs21gRhQfWuEA+RQFjp4Eo6BHv+5EXSOqvYF5Xzhl2u1sAURXAIjsEpsMElqIIbUAdNgMEDeALP4MV4NCbGq/E2LS0Ys5598AfG+w+305xt x 2 X m (b) SEP Fig. 1. Graphical representation of tw o autoregressive cokriging mo dels. Here Σ t ( j, j ) in the left panel is the j th diagonal element of the cov ariance matrix Σ t . co v ariance matrix Σ . Supp ose that the matrix Ω admits a mo diﬁed Cholesky decomp osition: Ω = ( I − A ) ⊤ D − 1 ( I − A ) , where A and D are unique low er triangular and diagonal matrices, resp ectively . Let a j b e the vector of non-zero elements in row j of A with a j j = 0 and d j > 0 b e the j -th diagonal elemen t of D . Then Z 1 ∼ MN n, 1 ( M 1 , R , d j ) Z j | Z 1: j − 1 ∼ MN n, 1 ( M j + ( Z 1: j − 1 − M 1: j − 1 ) a j , R , d j ) (2.2) where Z j and M j denote the j -th column of matrices Z and M and Z 1: j − 1 and M 1: j − 1 denote the submatrices of Z and M containing columns 1 to j − 1 . Result 2.1 sho ws that the problem of estimating the in verse cross-co v ariance matrix for a matrix-v ariate normal random v ariable Z is transformed into a series of N autoregressive linear regression problems, where a j and d j are the co eﬃcients and v ariance, respectively , from regressing the centered elements of Z j on the centered elements of Z 1: j − 1 . If the columns of Z ha ve a natural ordering, the inv erse cross-co v ariance matrix can b e approximated by regressing Z j on the p preceding columns to co ordinate j of Z ( Bick el and Levina , 2008 ; Lee and Lee , 2021 ). This sparse representation of inv erse cross-cov ariance matrix also implies conditional independence. Indeed, if a j i = 0 for i < j , then Z j is conditionally indep endent of Z i giv en Z 1 , . . . , Z i − 1 , Z i +1 , . . . , Z j − 1 . F or the separable autoregressive mo del ( 2.1 ), we assume that Ω t := Σ − 1 t = ( I − A t ) ⊤ D − 1 t ( I − A t ) . Let a t,j b e the vector of non-zero elements in row j of A t and d t,j b e the j -th diagnoal element of D t . Because the N computer mo del outputs at each ﬁdelity lev el were recorded at the same set of spatial lo cations , w e can determine a sparse structure for A t based on the pro ximit y of one lo cation to another. Let s 1 , . . . , s N b e the spatial lo cations of the N simulator outputs, sorted using the maximum-minim um-distance ordering ( Guinness , 2018 ). The v ector of non-zero elements in ro w j of A t is deﬁned as a ( C t,j ) t,j , and it includes all elemen ts in the set { a t,j i : i ∈ C t,j } , where C t,j is a set comprised of the indices for the p nearest spatial lo cations to s j c hosen from s 1 , . . . , s j − 1 . When p > j − 1 , C t,j = { 1 , . . . , j − 1 } . In the storm surge application, the simulator outputs are very smo oth ov er the spatial lo cations. W e found that using a single nearest neighbor for each set C t,j w orks well. F or applications where there is more v ariability b et ween neigh b oring outputs, some regularization techniques (see Kidd and Katzfuss , 2022 ) could be used to allo w for 7 larger neighbor sets in a C t,j t,j . After constructing the neighbor sets, y 1 ,j ( x ) will b e conditionally indep endent of y 1 ,i ( x ) for i < j given y 1 , 1 ( x ) , . . . , y 1 ,i − 1 ( x ) , y t,i +1 ( x ) , . . . , y t,j − 1 ( x ) and the parameters in mo del ( 2.1 ). The same conditional indep endence assumption holds for levels t > 1 given the outputs from level t − 1 . The collections of parameters in mo del ( 2.1 ) across the m ﬁdelit y levels are denoted by β := { β 1 , . . . , β m } , γ := ( γ 1 , . . . , γ m − 1 ) ⊤ , θ := { θ 1 , . . . , θ m } , A := { A 1 , . . . , A m } and D := { D 1 , . . . , D m } . Because the input design is nested, it is easy to chec k that the marginal likelihoo d of the SEP mo del ( 2.1 ) allows a closed-form expression L ( y D | β , γ , θ , A , D ) = m Y t =1 MN n t ,N ( Y t | F t B t , R t , Ω − 1 t ) , where for t = 1 , B 1 := β 1 , F 1 = H 1 , where H t = h t ( X t ) is an n t × q t matrix of the basis functions ev aluated at n t inputs ov er X t . F or t > 1 , B t := [ β ⊤ t , γ t − 1 I N × N ] ⊤ , F t = [ H t , U t − 1 ] , where U t − 1 := [ y t − 1 , 1 ( X t ) , . . . , y t − 1 ,N ( X t )] is an n t × N matrix of outputs from the mo del at level t − 1 at n t inputs o ver X t . The level- t correlation matrix R t := r ( X t , X t ; θ t ) has the ( i, j ) element r ( x i , x j ; θ t ) for x i , x j ∈ X t . F or j = 1 , . . . , N and t = 1 , . . . , m , a ( C t,j ) t,j and d t,j are assigned the following priors a ( C t,j ) t,j | d t,j ind ∼ N |C t,j | ( 0 , d t,j τ 2 t I ) , d t,j ind ∼ IG ( η t 2 , λ t 2 ) , (2.3) where |C t,j | is the cardinality of the set C t,j , N ( · , · ) denotes a multiv ariate-normal distribution and IG ( · , · ) denotes an In verse-Gamma distribution. W e assume that the h yp erparameters τ 2 t , η t and λ t are ﬁxed. Similar versions of the prior in ( 2.3 ) hav e b een used to estimate the precision matrix of multiv ariate normal random v ariables in high-dimensional settings ( Kidd and Katzfuss , 2022 ; Lee and Lee , 2021 ). The joint prior distribution of β , γ , A , D and θ is π ( β , γ , A , D , θ ) ∝ π ( β , γ ) m Y t =1 { π ( A t | D t ) π ( D t ) d Y i =1 π ( θ t,i ) } , (2.4) where π ( β , γ ) ∝ 1 , π ( A t | D t ) and π ( D t ) are the priors deﬁned in ( 2.3 ), while θ t,i iid ∼ C + (0 , q t,i ) , the half- Cauc hy distribution with scale parameter q t,i , for i = 1 , . . . , d and t = 1 , . . . , m . In the numerical studies, q t,i is ﬁxed at half the range of input domain along the i -th axis. The scale-discrepancy parameter γ t − 1 is often ﬁxed at one due to identiﬁably issues. W e also notice that allowing estimation of γ do es not improv e the predictive accuracy in our n umerical studies. In the follo wing sections, we assume that γ is ﬁxed. 2.2.2 P osterior Inference The form of the join t prior distribution in ( 2.4 ) enables analytic integration of β , A and D in p ( β , A , D , θ | y D ) using standard conjugacy results, simplifying computation of the p osterior distribution for computer mo dels with high dimensional outputs. Let ˆ B 1 := ˆ β t and ˆ B t := [ ˆ β ⊤ t , γ t − 1 I N × N ] ⊤ , where ˆ β 1 := ( H ⊤ 1 R − 1 1 H 1 ) − 1 H ⊤ 1 R − 1 1 Y 1 and ˆ β t := ( H ⊤ t R − 1 t H t ) − 1 H ⊤ t R − 1 t ( Y t − γ t − 1 U t − 1 ) for t > 1 are the generalized least-squares estimates of β 1 and β t . W e state this result in the following prop osition with its pro of provided in Section S.1.2 of the Supplemen tary Material. 8 Pr oposition 2.1 . The p osterior distribution of θ is p ( θ | y D ) ∝ m Y t =1 | R t | − N 2 | H ⊤ t R − 1 t H t | − N 2 ( ˆ d t, 1 ) − ν t 2 N Y j =2 | ˆ V t,j | 1 2 ( ˆ d t,j ) − ν t 2 π ( θ t ) , (2.5) where ν t = n t + η t − q t , ˆ V t,j = ( S ⊤ t, C t,j S t, C t,j + τ − 2 t I ) − 1 , ˆ d t, 1 = ( S ⊤ t, 1 S t, 1 + λ t ) , ˆ d t,t = ( S ⊤ t,j ( I − S t, C t,j ˆ V t,j S ⊤ t, C t,j ) S t,j + λ t ) , S t = R − 1 2 t ( Y t − F t ˆ B t ) , and S t, C t,j is the submatrix of S t con taining all columns in the set C t,j . Prop osition 2.1 sho ws that the p osterior distributions of mo del parameters θ t are conditionally inde- p enden t across diﬀerent ﬁdelity levels, facilitating fast computation. Posterior samples of θ t are drawn separately for t = 1 , . . . , m using standard Metrop olis-Hastings algorithms. Sampling from p ( θ | y D ) remains computationally feasible for large N as long as the n um b er of neighbors in eac h set C t,j is small. The primary goal of multiﬁdelit y computer mo del em ulation is to determine the predictive distribu- tion of the outputs for the simulator at highest ﬁdelit y level, m . F or any new input x 0 , let y t ( x 0 ) := ( y t, 1 ( x 0 ) , . . . , y t,N ( x 0 )) ⊤ b e a vector of N outputs at ﬁdelit y level t ∈ { 1 , . . . , m } and new input x 0 . Let y ( x 0 ) := ( y 1 ( x 0 ) ⊤ , . . . , y m ( x 0 ) ⊤ ) ⊤ denote the collection of outputs at x 0 across m ﬁdelity levels. W e can establish that the posterior predictiv e distribution p ( y ( x 0 ) | θ , y D ) has a closed-form expression with its pro of giv en in Section S.1.3 of the Supplementary Material. Deﬁne M t ( x 0 ) =    ˆ β ⊤ 1 h 1 ( x 0 ) + ( Y 1 − F 1 ˆ B 1 ) ⊤ R − 1 1 r 1 ( x 0 ) if t = 1 ˆ β ⊤ t h t ( x 0 ) + γ t − 1 y t − 1 ( x 0 ) + ( Y t − F t ˆ B t ) ⊤ R − 1 t r t ( x 0 ) if t > 1 (2.6) R t ( x 0 ) = 1 − r ⊤ t ( x 0 ) R − 1 t r t ( x 0 ) + ( h t ( x 0 ) − H ⊤ t R − 1 t r t ( x 0 )) ⊤ ( H ⊤ t R − 1 t H t ) − 1 ( h t ( x 0 ) − H ⊤ t R − 1 t r t ( x 0 )) , (2.7) where r t ( x 0 ) := r ( x 0 , X t ; θ t ) is an n t × 1 vector of correlations b etw een x 0 and the inputs in X t . Pr oposition 2.2 . The distribution of y ( x 0 ) | θ , y D is p ( y ( x 0 ) | θ , y D ) = m Y t =1  t ν t ( y t, 1 ( x 0 ) | µ t, 1 ( x 0 ) , σ 2 t, 1 ( x 0 )) N Y j =2 t ν t ( y t,j ( x 0 ) | µ t,j ( x 0 ) , σ 2 t,j ( x 0 ))  , (2.8) where µ t, 1 ( x 0 ) = M t, 1 ( x 0 ) , σ 2 t, 1 ( x 0 ) = ν − 1 t R t ( x 0 ) ˆ d t, 1 , and for j > 1 , µ t,j ( x 0 ) = M t,j ( x 0 ) + ( y t, C t,j ( x 0 ) − M t, C t,j ( x 0 )) ⊤ ˆ a ( C t,j ) t,j , σ 2 t,j ( x 0 ) = ν − 1 t R t ( x 0 ) ˆ d t,j (( y t, C t,j ( x 0 ) − M t, C t,j ) ⊤ ˆ V t,j ( y t, C t,j ( x 0 ) − M t, C t,j ) + 1) , ˆ a ( C t,j ) t,j = ˆ V t,j S ⊤ t, C t,j S t,j , and ν t , ˆ d t,j , ˆ V t,j and S t are deﬁned in Prop osition 2.1 . Prop osition 2.2 shows that p osterior samples from p ( y ( x 0 ) | θ , y D ) can b e generated from a sequence of m distributions. More imp ortantly , b ecause of the conditional indep endence assumption for the outputs, at eac h ﬁdelit y lev el, the p osterior predictiv e distribution can be computed eﬃciently without factorizing a large cross-co v ariance matrix. Samples of y m ( x 0 ) given y D are obtained using typical Marko v chain Monte Carlo (MCMC) algorithms. In addition, plug-in estimates for prediction via p ( y ( x 0 ) | θ , y D ) is used to facilitate fast computation compared to a fully Bay es pro cedure. F or instance, prediction can b e obtained by plugging in θ with its p osterior mo de estimate ˆ θ MAP from the MCMC algorithm. 9 2.3 Nonseparable Autoregressiv e Cokriging 2.3.1 Mo del F orm ulation A separable cov ariance structure assumes the same cross-co v ariance in the output space, whic h may b e restrictiv e in some real-w orld applications (e.g., F rick er et al. , 2013 ). A nonseparable dep endence struc- ture is more ﬂexible, but allowing nonseparability is particularly challenging for simulators that pro duce high-dimensional outputs. F ollowing the idea in Higdon et al. ( 2008 ), we also develop a nonseparable autore- gressiv e mo del for multiﬁdelit y computer mo dels. Similar to Higdon et al. ( 2008 ), we ﬁrst p erform principal comp onen t analysis (PCA) on the outputs. In practice, this step is often p erformed on standardized outputs. F or notational simplicit y , w e also refer to Y t as an n t × N matrix of standardized outputs at ﬁdelity level t . Then we assume a basis represen tation for the level- t outputs y t ( x ) = p t X ℓ =1 k t,ℓ w t,ℓ ( x ) + ϵ t , (2.9) where k t,ℓ is the ℓ -th N × 1 principal comp onent (or basis vector) constructed from principal comp onent analysis of Y t , w t,ℓ ( x ) is the random w eight asso ciated with the ℓ -th basis vector, p t is the dimension of the reduced space and ϵ t is a random reconstruction error. The reconstruction error, which accounts for discrepancies b et w een the p t dimensional basis represen tation of the outputs and the actual v alues of Y t , is mo deled indep endently from the random weigh ts. The reconstruction error and the selection of p t are discussed at the end of the section. F or computer mo dels with a single ﬁdelit y level, the random weigh ts are often modeled using indep endent Gaussian pro cesses ( Higdon et al. , 2008 ). T o preserve the rank of output accuracy , we assume an autore- gressiv e cokriging mo del for the PC weigh ts. Let w t − 1 ( · ) = [ w t − 1 , 1 ( · ) , . . . , w t − 1 ,p t ( · )] ⊤ denote a vector of p t random functions whose partial realizations are given by PC scores. In particular, we consider indep enden t autoregressiv e cokriging mo dels for the random weigh ts across diﬀeren t PC weigh ts: w t,ℓ ( · ) | w t − 1 ( · ) , β t − 1 ,ℓ , σ 2 t,ℓ , θ t,ℓ ind ∼ G P ( β ⊤ t − 1 w t − 1 ( · ) , σ 2 t,ℓ r ( · , · , θ t,ℓ )) , ℓ = 1 , . . . , p t w 1 ,ℓ ( x ) | σ 2 1 ,ℓ , θ 1 ,ℓ ind ∼ G P (0 , σ 2 1 ,ℓ r ( · , · ; θ 1 ,ℓ )) , (2.10) where β t − 1 := [ β t − 1 , 1 , . . . , β t − 1 ,p t ] ⊤ denotes a vector of scale discrepancy parameters, σ 2 t,ℓ is a v ariance parameter and r ( · , · ; θ t,ℓ ) is a correlation function with parameter θ t,ℓ . F or t = 1 , . . . , m and ℓ = 1 , . . . , p t , w e used the same pro duct form for the correlation function as the SEP mo del in ( 2.1 ). The mean of the lev el- 1 weigh ts is set to 0 b ecause the outputs were centered b efore computing the principal comp onents. The term β ⊤ t − 1 w t − 1 ( · ) acts lik e as a regression mean function with predictors given b y the PC weigh t v ectors w t − 1 ( x 1 ) , . . . , w t − 1 ( x n t − 1 ) . In the numerical exp eriment and application, w e found that the PC w eights w t,ℓ ( x ) is often highly correlated with the PC weigh ts w t − 1 ,ℓ ( x ) . This implies that not all weigh ts { w t − 1 ,ℓ ( · ) , ℓ = 1 , . . . , p t − 1 } contribute to the v ariation of w t − 1 ,ℓ ( · ) . This actually can be used to improv e mo del iden tiﬁabilit y by only including those highly correlated weigh ts at low-ﬁdelit y levels to capture the v ariability in the weigh ts at high-ﬁdelity level since the autoregressive cokriging mo del is not identiﬁable without informative prior on parameters ( Kennedy and O’Hagan , 2001 ). Nev ertheless, prediction can still b e carried out legitimately . 10 2.3.2 P osterior Inference Let w t,ℓ = [ w t,ℓ ( x 1 ) , . . . , w t,ℓ ( x n t )] ⊤ . Given the PC weigh ts, the likelihoo d of mo del ( 2.10 ) can b e shown to ha ve closed-form expression: L ( w D | β w , σ 2 w , θ w ) = ( p 1 Y ℓ =1 N n 1 ( w 1 ,ℓ | 0 , σ 2 1 ,ℓ R 1 ,ℓ ) ) m Y t =2 p t Y ℓ =1 N n t  w t,ℓ | W t − 1 β t − 1 ,ℓ , σ 2 t,ℓ R t,ℓ )  , (2.11) where w D := { w 1 , 1 , . . . , w m, 1 , . . . , w m,p m } , β w = { β 1 , 1 , . . . , β m − 1 ,p m − 1 } , σ 2 w := { σ 2 1 , 1 , . . . , σ 2 m,p m } , θ w := { θ 1 , 1 , . . . , θ m,p m } , W t − 1 is an n t × p t − 1 matrix of PC w eights at the level t − 1 ov er the inputs in X t and R t,ℓ is an n t × n t correlation matrix with the ( i, j ) element r ( x i , x j ; θ t,ℓ ) . T o complete the Bay esian sp eciﬁcation, w e assume the following prior distribution for β w , σ 2 w , θ w π ( β w , σ 2 w , θ w ) ∝ m Y t =1 p t Y ℓ =1 1 σ 2 t,ℓ π ( θ t,ℓ ) , (2.12) where π ( θ t,ℓ ) is the density for the prior on correlation parameter θ t,ℓ . W e used indep endent half-Cauch y priors for the range parameters at t = 1 . F or t > 1 , we found that more informativ e prior distributions for the range parameters help stabilize p osterior sampling with indep endent half-normal prior distributions. Using standard conjugacy results, the lik eliho o d ( 2.11 ) and prior ( 2.12 ) enable analytic integration of β w and σ 2 w from the posterior p ( β w , σ 2 w , θ w | w D ) , yielding a closed-form expression for p ( θ w | w D ) . Samples from the p osterior distribution p ( θ w | w D ) are generated using MCMC (Mark ov chain Monte Carlo) algorithms . F or new input x 0 , we also deriv ed a closed-form expression for the predictiv e distribution p ( w ( x 0 ) | θ w , w D ) , where w ( x 0 ) := ( w 1 ( x 0 ) ⊤ , . . . , w m ( x 0 ) ⊤ ) ⊤ . The expressions for b oth p ( θ w | w D ) and p ( w ( x 0 ) | θ w , w D ) are pro vided in Section S.1.4 of the Supplemen tary Material. Similarly , we noticed that separate modeling ﬁtting and prediction can b e done b ecause of the autoregressiv e structure. T o predict the m -th cen tered and scaled computer mo del outputs for x 0 , p osterior samples drawn from p ( w m ( x 0 ) | w D ) are pro jected back to the original space using the basis represen tation in ( 2.9 ) with ϵ m = n m X ℓ = p m +1 k m,ℓ ϵ m,ℓ , where k m,ℓ is the ℓ th discarded principal comp onent direction and the ϵ m,ℓ are indep endent normal random v ariables with v ariance equal to the v ariance of Y m explained b y direction k m,ℓ ( Wilkinson , 2010 ). Predic- tions are then transformed to the original scale of Y m . While the primary goal of multiﬁdelit y computer mo del emulation is to predict the outputs at the highest ﬁdelity level, the basis representation ( 2.9 ) enables predictions at any ﬁdelity level. T o select the num b er of basis vectors, Higdon et al. ( 2008 ) found that c ho osing the ﬁrst few principal comp onents to explain greater than 95 or 99 p ercent of the total v ariation in the data works well. How ever, in some applications, the remaining principal comp onents that explain a very small p ercentage of total v ariability could add meaningful con tributions to predictions ( Wilkinson , 2010 ). T o guide the selection of basis vectors, we ev aluate the ro ot mean square prediction error (RMSPE) in the mo del based on diﬀerent num b er of basis vectors in the numerical studies. 11 3 T estBed Example W e study the predictiv e p erformance of the SEP and NONSEP em ulators using a simulated computer exp erimen t from Blizn yuk et al. ( 2008 ) As a comparison to the prop osed metho ds, we also ﬁt tw o mo dels that feature indep endent Gaussian pro cesses at eac h lo cation with shared correlation parameters: the parallel partial autoregressive cokriging (PP-AR Cokrig) mo del ( Ma et al. , 2022a ) via the R pack age ARCokrig ( Ma , 2022 ) and the parallel partial (PP) mo del ( Gu and Berger , 2016 ) via the R pack age RobustGaSP ( Gu et al. , 2019 ). The PP emulator is ﬁt to the high-ﬁdelit y outputs only . W e also ﬁt the sparse V ecchia appro ximation ( Katzfuss and Guinness , 2021 ) for the data ov er the pro duct space of the input space and spatial domain, but the results are omitted due to w orse predictiv e performance. T o assess the predictive p erformance of eac h emulator at multiple outputs, we also consider the a verage high-ﬁdelity output for inpu t x ¯ y 2 ( x ) = N − 1 N X j =1 y 2 ,j ( x ) , (3.1) whic h is useful when decision-making is needed at aggregated spatial resolutions. In what follo ws, w e illustrate the predictive p erformances among these four mo dels: PP , PP-ARCokrig, SEP , and NONSEP , where PP is ﬁtted to the high-ﬁdelit y data only . T o measure predictive p erformance, we use RMSPE and the co verage probability of 95% equal-tailed credible in terv als (CV G), and a v erage length of 95% equal-tailed credible interv als (ALCI). 3.1 Exp erimen tal Design In this testbed example, w e consider the computer mo del from Bliznyuk et al. ( 2008 ) that captures the spread of a contaminan t after a chemical spill ov er a space-time domain using the function y 2 ( x , s ) = √ 4 π C ( x , s ) , where C ( x , s ) , the concentration of the p ollutan t ev aluated at the input x := ( M , D , L, T ) ⊤ and space-time co ordinate s := ( s 1 , s 2 ) , is C ( x , s ) = M √ 4 π D s 2 exp  − s 2 1 4 D s 2  + M p 4 π D ( s 2 − T ) exp  − ( s 1 − L ) 2 4 D ( s 2 − T )  1 ( T < s 2 ) , (3.2) where the parameter M ∈ [7 , 13] is the mass of the p ollutant spilled at each lo cation, D ∈ [0 . 02 , 0 . 12] is the diﬀusion rate in the channel, L ∈ [0 . 01 , 3] is the lo cation of the second s pill and T ∈ [30 , 30 . 295] is the time of the second spill, s 1 is the spatial co ordinate and s 2 is the time. T o mimic a multiﬁdelit y computer mo del, the function y 2 ( x , s ) is approximated using a one-term T aylor series expansion of the tw o exp onen tial terms in Equation ( 3.2 ). In what follows, we refer to y 1 ( x , s ) , the approximation of y 2 ( x , s ) , as the low-ﬁdelit y mo del. F or eac h input setting, the low and high ﬁdelit y sim ulators were ev aluated at a grid of 20 evenly spaced-spatial co ordinates ov er [0 . 5 , 5] and 50 evenly-spaced time points ov er [35 , 60] , resp ectively . With the parameter T ﬁxed at 30 , 60 v alues of ( M , D , L ) were selected for the low-ﬁdelit y mo del using a Latin hypercub e sample ( Gramacy , 2020 ). The function y 2 ( x , s ) w as ev aluated at 30 inputs drawn from the collection of 60 low-ﬁdelit y inputs using a conditional Latin h yp ercub e sample ( Minasny and McBratney , 2006 ). T o assess the predictive p erformance, 30 test inputs of ( M , D , L ) were randomly sampled using a uniform distribution ov er [7 , 13] × [0 . 02 , 0 . 12] × [0 . 01 , 3] . 12 3.2 Em ulation Results Before w e ﬁt the models, it is often helpful to perform diagnostic c hecks to identify k ey data analytic c hallenges and rule out implausible mo dels. The empirical mean and v ariance of the outputs at each spatial lo cation indicate that the output space is nonstationary; see Figure S.1 of the Supplemen tary Material. Thus, a stationary output co v ariance structure is not appropriate for the application. Instead, we may wish to ﬁt PP , PP-ARCokrig, SEP , and NONSEP . After taking the diﬀerence b etw een the high-ﬁdelity outputs and lo w-ﬁdelity outputs, the residuals show m uch less v ariability than the original outputs but still exhibit a high degree of spatial correlation, suggesting that the autoregressive structure is reasonable for the application. T able 1. Predictiv e performance of diﬀerent mo dels at 30 test inputs and 1 , 000 spatial locations. Predictive p erformance is ev aluated for individual high-ﬁdelity outputs and for the aggregated high-ﬁdelity outputs. Marginal predictive p erformance of y 2 ,j ( · ) Predictiv e p erformance of ¯ y 2 ( · ) Mo del RMSPE CV G ( 95% ) ALCI ( 95% ) RMSPE CVG ( 95% ) ALCI ( 95% ) SEP 0.195 97.0 0.739 0.051 93.3 0.158 NONSEP 0.250 97.3 0.870 0.049 100 0.402 PP-AR Cokrig 0.208 95.7 0.626 0.045 23.3 0.023 PP 0.267 95.9 0.726 0.058 20 0.028 W e show the predictive p erformance for the four aforementioned mo dels in T able 1 . The metrics for marginal predictiv e p erformance in T able 1 are computed at individual spatial lo cations and then av eraged o ver the entire spatial domain. In terms of RMSPE, the SEP mo del achiev e the b est results. In terms of CV G and ALCI, all mo dels achiev e go o d p erformance. While PP and PP-ARCokrig do not account for output cross-cov ariance structure, it is not surprising that they can achiev e calibrated CVG and ALCI since these metrics ev aluate the p erformance of marginal predictive distributions at individual spatial lo cations. When comparing the p erformance at aggregated spatial lo cations, such as the spatial av erage deﬁned in Equation ( 3.1 ), we exp ect the p erformance of PP and PP-AR Cokrig to deteriorate substantially b ecause the t wo mo dels do not allow prop er UQ, which is conﬁrmed by the results in T able 1 . Among the t w o prop osed mo dels, SEP gives b etter results than NONSEP in terms of UQ. The detailed implementation of SEP is giv en in Section S.2.1 of the Supplementary Material. In the NONSEP mo del, we performed PCA separately for the outputs at eac h ﬁdelity level. Figure S.3 of the Supplementary Material shows that PC scores at diﬀerent ﬁdelit y lev els are only highly correlated when they share the same component n umber, indicating that an y PC weigh t function at the high-ﬁdelity level can just b e mo deled with the corresp onding weigh t function at the low-ﬁdelit y level. The reason for this is likely due to the go o d approximation of low-ﬁdelit y simulator to the high-ﬁdelity simulator. The ﬁrst four PCs explain more than 99% of the total output v ariability at each ﬁdelit y level, but the NONSEP mo del do es not give satisfactory results with only four PCs; in contrast, choosing the n um b er of PCs based on RMSPE w ould lead to substantial improv ement in prediction; see Figure 2 . Because the RMSPE plateaus after an eigh th principal comp onent is added to the mo del, w e c hose eigh t PCs for the basis representation ( 2.9 ) at b oth levels to ﬁt the autoregressive cokriging mo del ( 2.10 ). F or each PC w eight and ﬁdelity lev el, the Metrop olis-Hastings algorithm w as run for 60 , 000 iterations (with a burn-in of 6 , 000 samples) to obtain 13 0.25 0.30 0.35 0.40 5 6 7 8 # PCs RMSPE Cokriging Kriging Fig. 2. RMSPE as a function of n umber of principal comp onents. The black line represen ts the results under the NONSEP cokriging model while the blue-dashed line represen ts the results based on NONSEP kriging mo del ﬁt to the high-ﬁdelity simulator only . The RMSPE was computed using the av erage p osterior predictiv e weigh t. p osterior samples of the correlation parameters. W e noticed that the MCMC samples do not show any lack of conv ergence b ehaviors and p osterior distributions are not sensitive to diﬀerent hyperparameter settings. Figure 2 shows how the RMSPE for the NONSEP emulator (cokriging) compares to a mo del ﬁt using only the PC weigh ts from level-2 (kriging). The autoregressive structure for the PC weigh ts improv es predictions of the cokriging mo del relative to the kriging mo del for the seven th and eighth principal comp onents, where there is a weak er relationship b etw een the random w eigh ts and inputs. 4 Multiﬁdelit y Em ulation of Storm Surges 4.1 Storm Surge Sim ulators and Exp erimental Design ADCIR C is an o cean circulation mo del used to simulate hurricane driven storm surge ( Luettich and W est- erink , 2004 ). ADCIR C is more accurate when coupled with the SW AN mo del ( Bo oij et al. , 1999 ; Zijlema , 2010 ), whic h incorp orates wa ve eﬀects. This computer mo del was studied by ( Ma et al. , 2022a ), where do wnstream applications and v alidation of computer mo dels hav e b een in tro duced. W e consider the p eak surge elev ation, measured in meters (m), generated at 9 , 284 lo cations in the Cap e Coral region that is selected based on the sea levels that are of primary interest in coastal ﬂo o d hazard studies. The ADICRC and ADCIRC + SW AN sim ulators take six input parameters: ∆ P ∈ [30 , 70] (central pressure deﬁcit (m b)), R p ∈ [16 , 39] (scale pressure radius (nautical miles)), V f ∈ [3 , 10] (forw ard sp eed (m/s)), θ ∈ [15 , 75] (heading in degrees clo c kwise from north), B ∈ [0 . 9 , 1 . 4] (Holland’s B parameter), ℓ (landfall lo cation in latitude and longitude). The ranges of these parameters are selected based on the characteristics of landfalling hurricanes. F ollowing Ma et al. ( 2022a ), exp erimental design is based on a maximin Latin h yp ercub e design to select 50 unique combinations of (∆ P , R p , V f , θ , B ) that c haracterize storm characteristics. F or each storm parameter of (∆ P , R p , V f , θ , B ) , the landfall lo cation ℓ w as chosen with one-half of R p uniformly spaced ov er the coastline but with random initial landfall lo cations on the southw est p osition; see Ma et al. ( 2022a ) for additional details. In total, the ADCIRC and ADCIRC + SW AN computer mo dels were run at 226 inputs 14 after excluding unrealistic storm c haracteristics. Since the coastline in Cap e Coral has very low v ariation in latitude, ℓ was con verted to a distance measure, d ℓ , using the lo cation farthest to the north west along the coastline as a reference. In terms of training, random samples of 120 and 60 inputs were selected for the ADCIR C and ADCIRC + SW AN simulators, resp ectiv ely , using a nested design. The remaining 106 inputs are used to assess predictive p erformance at individual spatial lo cations and aggregated spatial resolutions. 4.2 Diagnostics As in the testbed example, we ﬁrst p erform exploratory data analysis to c heck model assumptions. Figure S.6 of the Supplementary Material sho ws that there is region-sp eciﬁc spatial heterogeneity in mean and v ariance for b oth lo w and high ﬁdelity outputs, indicating that mo dels like PP and PP-ARCokrig might b e suitable. The residuals in Figure S.6 indicates that after taking out the v ariability due to low-ﬁdelit y simulation, the high-ﬁdelit y simulation shows some region sp eciﬁc nonstationary b eha viors, which fav ors a nonstationary co v ariance structure as in the four mo dels under comparison: PP , PP-ARCokrig, SEP , and NONSEP . Because nearb y outputs tend to b e similar (see Figure S.6 ), the output cross-cov ariance must b e accounted for. As suc h, mo dels like SEP and NONSEP , which include explicit cross-cov ariance structures are likely to improv e mo del ﬁtting. Figure S.6 of the Supplementary Material shows spatial maps of the empirical mean and empirical v ariance in the residuals (diﬀerence b etw een high-ﬁdelity sim ulator and low ﬁdelity sim ulator) ov er nested inputs, indicating that sensible mo dels need to account for spatially-v arying patterns in the mean function and cross-cov ariance matrix in the output space. The residual plots in Figure S.6 of the Supplementary Material also imply that the autoregressive structure is a reasonable assumption because the v ariation in the high-ﬁdelit y simulator is substantially reduced after regressing on the low-ﬁdelit y simulator for ov erlapping inputs. As shown in Figure S.9 of the Supplemen tary Material, there is a very strong linear trend b etw een lo w-ﬁdelity PC scores and high-ﬁdelity PC scores at the same PC comp onent num b er, whic h allows us to sp ecify the mean functions in the autoregressive structure in mo del ( 2.10 ). 4.3 Em ulation Results The predictiv e performance o ver four diﬀeren t models are given in T able 2 . In terms of RMSPE, PP- AR Cokrig gives the b est results for assessment ov er individual lo cation and aggregated spatial lo cations. This is not surprising since past work ( F rick er et al. , 2013 ) found that including an output dep endence structure in a multiv ariate emulator may lead to slightly worse marginal predictive p erformance in terms of RMSE. How ever, we can see from T able 2 that b oth SEP and NONSEP based on seven PCs give b etter UQ in terms of CV G ( 95% ) and ALCI ( 95% ). This agrees with the ﬁndings in Mak et al. ( 2018 ). When UQ is assessed for predictive distributions ov er aggregated quan tities, the lac k of cross-co v ariance in models suc h as PP and PP-ARCokrig b ecomes a serious issue as indicated in T able 2 . Compared with SEP , NONSEP gives b etter UQ of predictive distributions at individual lo cations and aggregated spatial resolution. Figure 3 also conﬁrms that the NONSEP mo del can b etter capture the residuals and prediction uncertain t y is smaller in smo oth spatial regions. The detailed implemen tation of the SEP mo del is given in Section S.3 of the Supplemen tary Material. In what follows, w e fo cus on presen ting the results on the NONSEP model b ecause of its sup erior p erformance. In the NONSEP mo del, we categorize the errors in predicting the outputs in to t w o diﬀeren t t yp es: 15 T able 2. Predictiv e p erformance of ADCIRC + SW AN emulators at 106 test inputs and 9 , 284 spatial lo cations. Marginal predictive p erformance of y 2 ,j ( · ) Predictiv e p erformance of ¯ y 2 ( · ) Mo del RMSPE CV G ( 95% ) ALCI ( 95% ) RMSPE CVG ( 95% ) ALCI ( 95% ) SEP 0.111 91.50 0.409 0.079 90.57 0.284 NONSEP 0.102 94.65 0.373 0.060 96.23 0.248 PP-AR Cokrig 0.099 91.32 0.345 0.057 4.72 0.004 PP 0.170 84.07 0.452 0.118 0.94 0.005 m −0.2 −0.1 0.0 0.1 (a) Diﬀerence SEP . m −0.2 −0.1 0.0 0.1 (b) Diﬀerence NONSEP . m 0.1 0.2 0.3 0.4 0.5 0.6 (c) CI Len. SEP . m 0.1 0.2 0.3 0.4 0.5 0.6 (d) CI Len. NONSEP . Fig. 3. Predictiv e comparison under SEP and NONSEP ov er one testing input: ∆ P = 37 . 58 , R p = 29 . 37 , V f = 7 . 96 , θ = 70 . 42 , B = 1 . 38 , ℓ = ( − 82 . 25 , 26 . 83) . Panel (a) sho ws the diﬀerence b etw een actual high-ﬁdelit y output and predicted output under the SEP mo del; P anel (b) shows the diﬀerence b etw een high-ﬁdelit y output and predicted output under the NONSEP mo del. Panels (c) and (d) show the lengths of 95% equal-tail credible interv al under the SEP and NONSEP mo dels, resp ectively . (a) reconstruction error (RE) due to the choice of n umber of PCs and (b) the mo del (or em ulator) error in predicting the PC w eights. The smallest reconstruction error is obtained using a basis representation constructed b y p erforming PCA simultaneously on both the training and testing data and will b e referred to as the minimum r e c onstruction err or (MRE) in T able 3 . T o assess ho w the ﬁxed basis v ectors in ( 2.9 ), whic h are constructed using only training data, inﬂuence the total prediction error, the MRE can b e compared to the reconstruction error for a reduced-rank represen tation of the test outputs based on k 2 , 1 , . . . , k 2 ,p 2 in 16 ( 2.9 ). W e refer to the second reconstruction error as the testing r e c onstruction err or (TRE). Both MRE and TRE are computed by measuring the mean square prediction error (MSPE) from pro jecting a reduced-rank represen tation of the test outputs back to the original output space. The last tw o columns of T able 3 show the prediction error for the NONSEP emulator (cokriging) and a NONSEP mo del estimated using only the lev el-2 weigh ts (kriging). Comparing the MRE and TRE to the cokriging prediction error in T able 3 , the c hoice of basis functions explains the smallest p ercentage of the total prediction error of the three error sources. The autoregressive structure for the random weigh ts in ( 2.10 ) reduces the MSPE for the cokriging em ulator with seven PCs by more than 50% compared to the kriging mo del. Improv emen ts in the predictive p erformance of the NONSEP em ulator for more than sev en PCs is negligible. Figure 4 sho ws how the RMSPE for the NONSEP emulator with one PC changes after a second PC is added to the mo del. The lo cations with the largest reduction in RMSPE after a second PC is added to the mo del are closely related to the lo cations with the largest absolute elements of the second principal comp onent direction, k 2 , 2 , in Equation ( 2.9 ); see Panel (c) of Figure 4 . T able 3. Principal comp onen t analysis diagnostics and prediction p erformance under NONSEP cokriging and NONSEP kriging mo dels as a function of num b er of PCs. NONSEP cokriging is p erformed based on outputs at low-ﬁdelit y and high-ﬁdelity . NONSEP kriging is p erformed based on outputs at high-ﬁdelit y only . MRE = minimum reconstruction error. TRE = testing reconstruction error. % V ariance Explained 10 2 × MSPE # PCs ADCIR C ADCIRC + SW AN MRE TRE NONSEP cokriging kriging 1 91.26 92.60 6.83 7.00 7.45 8.35 2 95.93 96.87 2.64 2.79 3.33 4.46 3 97.50 98.20 1.55 1.96 2.57 3.73 4 98.60 99.07 0.91 1.00 1.62 2.85 5 99.32 99.50 0.44 0.49 1.17 2.50 6 99.53 99.67 0.32 0.37 1.08 2.43 7 99.65 99.76 0.27 0.33 1.04 2.40 5 Discussion In this pap er, we hav e developed an autoregressiv e cokriging framew ork with tw o multiﬁdelit y em ulators for mo deling computer mo dels with high-dimensional spatial outputs. T o eﬀectively account for the cross- co v ariance in the output space, we ha ve constructed a separable mo del and nonseparable mo del. In the SEP AR Cokrig emulator, w e mo del the cross-cov ariance structure as a Kronec ker pro duct of the input correlation and the output cross-cov ariance. T o estimate the cross-co v ariance matrix, we reformulate the problem as a series of successiv e autoregressive regressions, which allows us to interpret the choice of prior sp eciﬁcation and enables fast shrink age estimation via a sparse Cholesky prior, even for high-dimensional outputs. In the NONSEP AR Cokrig emulator, w e represent the outputs at eac h ﬁdelit y lev el as a linear combination of basis functions and random weigh ts. The random weigh ts are mo deled using an autoregressive cokriging structure across ﬁdelity levels to enforce the rank ordering of output accuracy . 17 0.1 0.2 0.3 0.4 (a) RMSPE (1 PC) 0.1 0.2 0.3 0.4 (b) RMSPE (2 PCs) 0.000 0.005 0.010 0.015 0.020 (c) | k 2 , 2 | . Fig. 4. Spatial maps of RMSPE and principal comp onent direction under the NONSEP mo del. Panels (a) and (b) sho w the RMSPE based on the ﬁrst PC and ﬁrst tw o PCs, resp ectively . Panel (c) shows the absolute v alue of the second principal comp onent direction at ﬁdelity level 2. W e ha v e derived closed-form formulas for p osterior predictive distributions up to input correlation pa- rameters for the tw o multiv ariate GP emulators. Marko v chain Monte Carlo algorithms are required only to sample from the p osterior distributions of the input correlation parameters, and can b e implemented sepa- rately for each ﬁdelit y level. This is similar to ﬁtting independent m ultiv ariate GP mo dels without accoun ting for the cross-cov ariance b etw een outputs. The prop osed framework th us enables fast Bay esian estimation and prediction. W e demonstrate the p erformance of the prop osed emulators in a testb ed example and a real application of storm surges. W e consistently ﬁnd that the t wo emulators with more complex cross-co v ariance structures can better quantify uncertainties for aggregated quantities ov er the output space without incurring extra computational burden. The improv ed uncertaint y quan tiﬁcation at aggregated spatial resolutions of the t wo proposed emulators would help practitioners in coastal ﬂo o d hazard studies p erform risk assessmen t. The prop osed framework exploits the practical implications of Marko v prop erty via autoregressiv e struc- tures, separable structure, and nonseparable structure for mo deling computer model outputs. These mo deling c hoices hav e b een empirically veriﬁed in the numerical example and storm surge application through exten- siv e diagnostics. The Marko v prop erty based autoregressive cokriging is suitable to link computer mo dels whose outputs exhibit a natural ordering in terms of accuracy . The mo deling choice b etw een a separable structure and a nonseparable structure has received lots of attentions in the literature. Our ﬁndings are consisten t with prior w ork ( F rick er et al. , 2013 ) that for highly rough outputs such as in the testb ed exam- ple, a separable cross-cov ariance may give b etter predictive p erformance than a PCA-based nonseparable structure. How ever, the SEP emulator may exhibit slight undercov erage due to its assumption of a constant cross-co v ariance across the input space. In the storm surge application, the NONSEP emulator achiev es sup erior predictive p erformance to the SEP emulator. The tw o prop osed emulators can b e further explored. In the SEP ARCokrig mo del, the selection of neigh b ors is informed by the smo othness of outputs ov er the output space. W e currently use the maximin design to sort spatial lo cations and then choose the num b er of neighbors to deﬁne the sparse Cholesky prior. Since there is no parametric assumption imposed, unlike V ecchia’s approximation in the purely spatial setting, further inv estigation is needed to study how the choice of neighbor sets aﬀects shrink age estimation of the output cross-cov ariance. In the NONSEP autoregressive cokriging, we use exploratory data analysis 18 to ols to sp ecify which low-ﬁdelit y PC scores should b e included to mo del the PC scores at high-ﬁdelity via the autoregressive cokriging model. The pro cedure could b e fully automated by emplo ying v ariable selection metho ds, such as spike-and-slab priors on the regression co eﬃcients. There are several p otential future researc h directions. First, in many physical systems, the outputs could b e generated at diﬀeren t spatial resolutions. Building a multiv ariate em ulator w ould need to account for the resolution diﬀerence. This could b e accounted for by adopting the change of supp ort strategy ( Ma et al. , 2019 ). Second, when there are more than t wo multiﬁdelit y computer mo dels, the SEP emulator can b e generalized in a straightforw ard manner; in contrast, there may b e further challenges in linking PC scores at m ul tiple ﬁdelity levels via the autoregressive cokriging framework. The selection of PC scores in the autoregressive cokriging mo del needs to be treated carefully via techniques such as Bay esian v ariable selection. Supplemen t Material The Supplement Material consists of: (i) technical details and additional results (.p df ﬁle) and (ii) R co de to repro duce the numerical results in the example and application (.zip ﬁle). References Ba yarri, M. J., Berger, J. O., Cafeo, J., Garcia-Donato, G., Liu, F., Palomo, J., Parthasarath y , R. J., Paulo, R., Sac ks, J., and W alsh, D. (2007). Computer mo del v alidation with functional output. The Annals of Statistics , 35(5):1874–1906. Publisher: Institute of Mathe matical Statistics. Bic kel, P . J. and Levina, E. (2008). Regularized estimation of large cov ariance matrices. The Annals of Statistics , 36(1):199–227. Blizn yuk, N., Rupp ert, D., Sho emaker, C., Regis, R., and Mugunthan, S. W. (2008). Bay esian calibration and uncertain ty analysis for computationally exp ensiv e mo dels using optimization and radial basis function appro ximation. Journal of Computation and Gr aphic al Statistics , 17(2):270–294. Bo oij, N., Ris, R. C., and Holthuijsen, L. H. (1999). A th ird-generation wa ve mo del for coastal regions: 1. Mo del description and v alidation. Journal of Ge ophysic al R ese ar ch: Oc e ans , 104(C4):7649–7666. Cialone, M. A., Grzegorzewski, A. S., Mark, D. J., Bryan t, M. A., and Massey , T. C. (2017). Coastal-storm mo del developmen t and w ater-level v alidation for the North Atlan tic coast comprehensive study . Journal of W aterway, Port, Co astal, and Oc e an Engine ering , 143(5):04017031. Con ti, S. and O’Hagan, A. (2010). Bay esian emulation of complex m ulti-output and dynamic computer mo dels. Journal of Statistic al Planning and Infer enc e , 140(3):640–651. Dietric h, J., Buny a, S., W esterink, J., Eb ersole, B., Smith, J., Atkinson, J., Jensen, R., Resio, D., Luettic h, R., Dawson, C., et al. (2010). A high-resolution coupled riv erine ﬂow, tide, wind, wind wa ve, and storm surge mo del for southern Louisiana and Mississippi. P art I I: Synoptic description and analysis of h urricanes Katrina and Rita. Monthly W e ather R eview , 138(2):378–404. Dietric h, J., Zijlema, M., W esterink, J., Holthuijsen, L., Dawson, C., Luettic h, R., Jensen, R., Smith, J., Stelling, G., and Stone, G. (2011). Mo deling hurricane wa ves and storm surge using in tegrally-coupled, scalable computations. Co astal Engine ering , 58(1):45 – 65. Dietric h, J. C., T anak a, S., W esterink, J. J., Da wson, C. N., Luettic h, R. A., Zijlema, M., Holth uijsen, 19 L. H., Smith, J. M., W esterink, L. G., and W esterink, H. J. (2012). P erformance of the unstructured- mesh, SW AN+ADCIR C mo del in computing hurricane w av es and surge. Journal of Scientiﬁc Computing , 52(2):468–497. FEMA (2008). Mississippi coastal analysis pro ject. Pro ject rep orts prepared by URS Group Inc., (Graithers- burg MD and T allahassee FL) under HMT AP Contract HSFEHQ-06-D-0162, T ask Order 06-J-0018. FEMA (2015). South west Florida Storm Surge Study Intermediate Data Submittal 1. Report Prepared by the Risk Assessment, Mapping, and Planning Partners (RAMPP) under FEMA IDIQ Contract HSFEHQ- 09-D-0369 and T ask Order HSFE04-13-J-0097. FEMA Region IV. F airfax, Virginia. FEMA (2017). South west Florida Storm Surge Study Intermediate Data Submittal 2. Report Prepared by the Risk Assessment, Mapping, and Planning Partners (RAMPP) under FEMA IDIQ Contract HSFEHQ- 09-D-0369 and T ask Order HSFE04-13-J-0097. FEMA Region IV. F rick er, T. E., Oakley , J. E., and Urban, N. M. (2013). Multiv ariate Gaussian pro cess emulators with nonseparable cov ariance structures. T e chnometrics , 55(1):47–56. F riedman, J., Hastie, T., and Tibshirani, R. (2008). Sparse inv erse cov ariance estimation with the graphical lasso. Biostatistics , 9(3):432–441. Gramacy , R. B. (2020). Surr o gates: Gaussian Pr o c ess Mo deling, Design and Optimization for the Applie d Scienc es . Chapman Hall/CRC, Bo ca Raton, Florida. Gramacy , R. B. and Apley , D. W. (2015). Lo cal Gaussian pro cess approximation for large computer exp er- imen ts. Journal of Computational and Gr aphic al Statistics , 24(2):561–578. Gu, M. (2019). Jointly robust prior for Gaussian sto c hastic pro cess in emulation, calibration and v ariable selection. Bayesian Analysis , 14(3):877–905. Gu, M. and Berger, J. O. (2016). P arallel partial Gaussian pro cess emulation for computer mo dels with massiv e output. The Annals of Applie d Statistics , 10(3):1317–1347. Gu, M., P alomo, J., and Berger, J. O. (2019). RobustGaSP: Robust Gaussian sto chastic pro cess emulation in R. The R Journal , 11(1):112–136. Guillas, S., Sarri, A., Day , S. J., Liu, X., and Dias, F. (2018). F unctional emulation of high resolution tsunami mo delling ov er Cascadia. The Annals of Applie d Statistics , 12(4):2023–2053. Guinness, J. (2018). P ermutation and grouping metho ds for sharp ening Gaussian pro cess approximations. T e chnometrics , 60:415–429. Gupta, A. K. and Nagar, D. K. (1999). Matrix variate distributions . Chapman and Hall/CRC, Boca Raton. Heo, J. and Sung, C.-L. (2025). Activ e learning for a recursive non-additive emulator for multi-ﬁdelit y computer exp eriments. T e chnometrics , 67(1):58–72. Hesser, T. J., Cialone, M. A., and Anderson, M. E. (2013). Lak e St. Clair: Storm w av e and w ater level mo deling. Higdon, D., Gattik er, J., and Williams, B. (2008). Computer model calibration using high-dimensional output. Journal of the Americ an Statistic al Asso ciation , 103(482):570–583. Jensen, R. E., Cialone, M. A., Chapman, R. S., Eb e rsole, B. A., Anderson, M., and Thomas, L. (2012). Lake Mic higan storm: W av e and water level mo deling. Ji, Y., Mak, S., Soeder, D., Paquet, J.-F., and Bass, S. A. (2024). A graphical m ulti-ﬁdelity Gaussian pro cess mo del, with application to emulation of heavy-ion collisions. T e chnometrics , 66(2):267–281. Katzfuss, M. and Guinness, J. (2021). A general framew ork for Vecc hia approximations of Gaussian pro cesses. Statistic al Scienc e , 36(1):124 – 141. 20 Kennedy , M. and O’Hagan, A. (2000). Predicting the output from a complex computer co de when fast appro ximations are av ailable. Biometrika , 87(1):1–13. Kennedy , M. C. and O’Hagan, A. (2001). Ba yesian calibration of computer mo dels. Journal of the R oyal Statistic al So ciety: Series B (Statistic al Metho dolo gy) , 63(3):425–464. Kidd, B. and Katzfuss, M. (2022). Bay esian nonstationary and nonparametric cov ariance estimation for large spatial data (with discussion). Bayesian Analysis , 17(1):291–351. K onomi, B. A. and Karagiannis, G. (2021). Ba yesian analysis of multiﬁdelit y computer mo dels with lo cal features and non-nested exp erimental designs: Application to the WRF mo del. T e chnometrics , 63(4):510– 522. Le Gratiet, L. (2013). Bay esian analysis of hierarchical multiﬁdelit y co des. SIAM/ASA Journal on Unc er- tainty Quantiﬁc ation , 1(1):244–269. Lee, K. and Lee, J. (2021). Estimating large precision matrices via mo diﬁed cholesky decomp osition. Sta- tistic a Sinic a , 31(1):173–196. Luettic h, R. and W esterink, J. (2004). F ormulation and numerical implementation of the 2D/3D ADCIRC ﬁnite element mo del version 44.XX. Ma, P . (2020). Ob jective Bay esian analysis of a cokriging mo del for hierarc hical multiﬁdelit y co des. SIAM/ASA Journal on Unc ertainty Quantiﬁc ation , 8(4):1358–1382. Ma, P . (2022). AR Cokrig: Autoregressive cokriging mo dels for multiﬁdelit y co des. R pack age version 0.1.3. h ttps://CRAN.R-pro ject.org/pac k age=ARCokrig. Ma, P ., Kang, E. L., Bra v erman, A. J., and Nguyen, H. M. (2019). Spatial statistical downscaling for con- structing high-resolution nature runs in global observing system simulation exp eriments. T e chnometrics , 61(3):322–340. Ma, P ., Karagiannis, G., Konomi, B. A., Asher, T. G., T oro, G. R., and Cox, A. (2022a). Multiﬁdelity computer mo del emulation with high-dimensional output: An application to storm surge. Journal of the R oyal Statistic al So ciety: Series C , 74(4):861–883. Ma, P ., Mondal, A., K onomi, B. A., Hobbs, J., Song, J. J., and Kang, E. L. (2022b). Computer mo del em- ulation with high-dimensional functional output in large-scale observing system uncertaint y exp eriments. T e chnometrics , 64(1):65–79. Mak, S., Sung, C.-L., W ang, X., Y eh, S.-T., Chang, Y.-H., Joseph, V. R., Y ang, V., and W u, C. F. J. (2018). An eﬃcient s urrogate mo del for emulation and physics extraction of large eddy simulations. Journal of the Americ an Statistic al Asso ciation , 113(524):1443–1456. Marso oli, R. and Lin, N. (2018). Numerical modeling of historical storm tides and w a v es and their in teractions along the US east and Gulf coasts. Journal of Ge ophysic al R ese ar ch: Oc e ans , 123(5):3844–3874. Matheron (1982). Pour une analyse krigeante des données régionaliées. T echnical Rep ort. Ecole Nationale Sup érieure des Mines de Paris. Minasn y , B. and McBratney , A. B. (2006). A conditioned latin hypercub e method for sampling in the presence of ancillary information. Computers and Ge oscienc e , 32:1378–1388. Niedoro da, A., Re sio, D., T oro, G., Div oky , D., Das, H., and Reed, C. (2010). Analysis of the coastal Mississippi storm surge hazard. Oc e an Engine ering , 37(1):82–90. P eherstorfer, B., Willcox, K., and Gunzburger, M. (2018). Survey of Multiﬁdelity Metho ds in Uncertaint y Propagation, Inference, and Optimization. Siam R eview , 60(3):550–591. P ourahmadi, M. (1999). Joint mean-co v ariance models with applications to longitudinal data: unconstrained parameterisation. Biometrika , 86(3):677–690. 21 Qian, P . Z. G. and W u, C. F. J. (2008). Bay esian hierarchical mo deling for integrating low-accuracy and high-accuracy exp eriments. T e chnometrics , 50(2):192–204. Sac ks, J., W elch, W. J., Mitc hell, T. J., and W ynn, H. P . (1989). Design and analysis of computer exp eri- men ts. Statistic al Scienc e , 4(4):409–435. San tner, T. J., Williams, B. J., and Notz, W. I. (2003). The Design and Analysis of Computer Exp eriments . Springer Series in Statistics. Springer Science+Business Media, New Y ork, NY. T owns, J., Co ck erill, T., Dahan, M., F oster, I., Gaither, K., Grimshaw, A., Hazlewoo d, V., Lathrop, S., Lifk a, D., P eterson, G. D., Roskies, R., Scott, J. R., and Wilkins-Diehr, N. (2014). XSEDE: Accelerating scien tiﬁc discov ery . Computing in Scienc e & Engine ering , 16(5):62–74. W ack ernagel, H. (2003). Multivariate Ge ostatistics: An Intr o duction with Applic ations . Springer-V erlag, New Y ork. W amsley , T., Go dsey , E., Bunch, B. W., Chapman, R. S., Grav ens, M. B., Grzegorzewski, A. S., Johnson, B. D., King, D. B., P ermenter, R. L., and Tillman, D. L. (2013). Mississippi coastal impro vemen t program, barrier island restoration numerical mo deling. W ang, H. (2012). Bay esian graphical lasso mo dels and eﬃcient p osterior computation. Bayesian Analysis , 7(4):867 – 886. W esterink, J. J., Luettich, R. A., F eyen, J. C., Atkinson, J. H., Dawson, C., Rob erts, H. J., Po well, M. D., Dunion, J. P ., Kubatko, E. J., and Pourtaheri, H. (2008). A basin-to channel-scale unstructured grid h urricane storm surge mo del applied to southern Louisiana. Monthly we ather r eview , 136(3):833–864. Wilkinson, R. D. (2010). Bayesian Calibr ation of Exp ensive Multivariate Computer Exp eriments , chapter 10, pages 195–215. John Wiley & Sons, Ltd. Yin, J., Lin, N., and Y u, D. (2016). Coupled mo deling of storm surge and coastal inundation: A case study in New Y ork City during Hurricane Sandy. W ater R esour c es R ese ar ch , 52(11):8685–8699. Zijlema, M. (2010). Computation of wind-w a ve sp ectra in coastal waters with SW AN on unstructured grids. Co astal Engine ering , 57(3):267 – 277. 22 Supplemen tary Material: Multiv ariate Gaussian pro cess em ulation for m ultiﬁdelit y computer mo dels with high-dimensional spatial outputs b y Cyrus S. McCrimmon and Pulong Ma Dep artment of Statistics, Iowa State University The Supplementary Material is organized as follo ws. Section S.1 contains deriv ations on the technical results from Section 2 of the main pap er including detailed deriv ation for Result 2.1 (see Sections S.1.1 ), Prop osition 2.1 (see Section S.1.2 ), Prop osition 2.2 (see Section S.1.3 ), and the p osterior distribution and p osterior predictiv e distributions of the random weigh ts in Section 2.3 (see Section S.1.4 ). Section S.2 includes additional details for the testb ed example. Section S.3 contains addition details for the storm surge application. S.1 T echnical Results S.1.1 Deriv ation for Result 2.1 In this section, w e provide a deriv ation of Result 2.1 , extending the reparameterization in Pourahmadi ( 1999 ) to a matrix normal random v ariable . T o simplify the notation, we consider the case where the mean is 0 . A t the end of the section, we describ e how the result can b e extended to a matrix normal random v ariable with a non-zero mean. Supp ose Z ∼ MN n × N ( 0 , R , Σ ) , where R is an n × n correlation matrix and Σ is an N × N cross-co v ariance matrix. Both R and Σ are p ositive deﬁnite. Consider the problem of ﬁnding b j := ( b j, 1 , . . . , b j,j − 1 ) ⊺ to minimize E [( Z j − Z 1: j − 1 b j ) ⊺ ( Z j − Z 1: j − 1 b j )] , where Z j and Z 1: j − 1 are the submatrices of Z con taining the j th and 1 , . . . , j − 1 columns of Z , resp ectively . Rewriting the quadratic form, the previous exp ectation is equiv alent to h − b ⊺ j 1 i E ( Z ⊺ 1: j Z 1: j ) " − b j 1 # , whic h, from Theorem 2.3.5 in Gupta and Nagar ( 1999 ), is n h 1 − b ⊺ j i Σ j " 1 − b j # , (S.1) where Σ j is a j × j matrix comprised of the ﬁrst j rows and columns of Σ . Setting the partial deriv ative of expression ( S.1 ) with resp ect to b j equal to 0 and solving for b j yields the solution Σ − 1 j − 1 σ j , where Σ j − 1 and σ j are the submatrices of Σ containing the ﬁrst j − 1 rows and columns and the j − 1 elements in row j , resp ectiv ely . Let a j := Σ − 1 j − 1 σ j . Then V ar ( Z j − Z 1: j − 1 a j ) = d j R , where d j = σ j,j − σ ⊺ j Σ − 1 j σ j , and σ j,j is the element in row j and column j of Σ . F rom Theorem 2.3.12 in Gupta and Nagar ( 1999 ), the conditional distribution of Z j | Z 1: j − 1 ∼ MN n, 1 ( Z 1: j − 1 a j , R , d j ) . The random vectors Z j − Z 1: j − 1 a j and Z i − Z 1: i − 1 a i are indep enden t for any integers 1 ≤ i < j ≤ N ; see Pourahmadi ( 1999 ) or Theorem 2.3.14 in Gupta and Nagar ( 1999 ). Let A be an N × N lo wer-triangular matrix with zero es along the diagonal and ( a j, 1 , . . . , a j,j − 1 ) = a ⊤ j . A dditionally , let D b e a diagonal matrix with elements d 1 , . . . , d N . Using indep endence b etw een Z j − Z 1: j − 1 a j and Z i − Z 1: i − 1 a i for any integers S.1 1 ≤ i < j ≤ N and the fact that V ar ( Z j − Z 1: j − 1 a j ) = d j R , Z ( I − A ) ⊤ d = MN n,N ( 0 , R , D ) . F rom Theorem 2.3.10 in Gupta and Nagar ( 1999 ), Z d = MN n,N ( 0 , R , ( I − A ) − 1 D ( I − A ) −⊤ ) , meaning Ω := Σ − 1 = ( I − A ) ⊺ D − 1 ( I − A ) . Using the parameterization Σ − 1 = ( I − A ) ⊺ D − 1 ( I − A ) and the deﬁnitions of a j and d j , the density of Z can b e written as p ( Z ) = p ( Z 1 ) N Y j =2 p ( Z j | Z 1: j − 1 ) , (S.2) where Z j | Z 1: j − 1 ∼ MN n, 1 ( Z 1: j − 1 a j , R , d j ) . If Z ∼ MN n,N ( M , R , Σ ) , the same solution to the minimiza- tion problem at the start of the section can b e obtained after subtracting M from Z . Then ( Z − M )( I − A ) ⊺ d = MN n,N ( 0 , R , D ) , whic h simpliﬁes to Z ∼ MN n,N ( M , R , ( I − A ) − 1 D ( I − A ) −⊤ ) . Finally , from Theorem 2.3.12 in Gupta and Nagar ( 1999 ), Z j | Z 1: j − 1 ∼ MN n, 1 ( M j + ( Z 1: j − 1 − M 1: j − 1 ) a j , R , d j ) . S.1.2 Deriv ation of Prop osition 2.1 In this section, we derive the density of p ( θ | y D ) deﬁned in Prop osition 2.1 . The joint p osterior distribution of β , A , D , θ is prop ortional to p ( β , A , D , θ | y D ) ∝ m Y t =1  | ( I − A t ) ⊺ D − 1 t ( I − A t ) | n t 2 | R t | − N 2 × exp  − 1 2 tr  ( I − A t ) ⊺ D − 1 t ( I − A t )( Y t − F t B t ) ⊺ R − 1 t ( Y t − F t B t )  × d Y l =1 2 π (1 + ( θ t,ℓ q t,ℓ ) 2 ) 1 { θ t,l > 0 } × d − η t 2 − 1 t, 1 exp  − λ t 2 d t, 1  × N Y j =2 ( τ 2 t ) −|C t,j | 2 ( d t,j ) −|C t,j | 2 exp  − 1 2 τ 2 t d t,j a ( C t,j ) ⊺ t,j a ( C t,j ) t,j  d − η t 2 − 1 t,j exp  − λ t 2 d t,j  . (S.3) After completing the square, exp  − 1 2 tr  Ω t ( Y t − F t B t ) ⊺ R − 1 t ( Y t − F t B t )  = exp  − 1 2 tr  Ω t ( Y t − F t ˆ B t ) ⊺ R − 1 t ( Y t − F t ˆ B t )  S.2 × exp  − 1 2 tr  Ω t ( β t − ˆ β t ) ⊺ ( H ⊺ t R − 1 t H t )( β t − ˆ β t )  (S.4) where ˆ β t and ˆ B t are deﬁned in Section 2.2.2 . Since the second term in Expression ( S.4 ) is equiv alent to the k ernel of a matrix normal distribution, the p osterior distribution ( S.3 ) can b e rewritten as p ( A , D , θ | y D ) ∝ m Y t =1  | Ω t | n t 2 | R t | − N 2 | H ⊺ t R − 1 t H t | − N 2 × exp  − 1 2 tr  Ω t ( Y t − F t ˆ B t ) ⊺ R − 1 t ( Y t − F t ˆ B t )  × d Y l =1 2 π (1 + ( θ t,ℓ q t,ℓ ) 2 ) × d − η t 2 − 1 t, 1 exp  − λ t 2 d t, 1  × N Y j =2 ( τ 2 t ) −|C t,j | 2 ( d t,j ) −|C t,j | 2 exp  − 1 2 τ 2 t d t,j a ( C t,j ) ⊺ t,j a ( C t,j ) t,j  d − η t 2 − 1 t,j exp  − λ t 2 d t,j  (S.5) where p ( β | A , D , θ , y D ) ∝ m Y t =1 MN q t ,N ( ˆ β t , ( H ⊺ t R − 1 t H t ) − 1 , Ω − 1 t ) . (S.6) In tegrating β out of expression ( S.5 ) yields the distribution of p ( A , D , θ | y D ) deﬁned as p ( A , D , θ | y D ) ∝ m Y t =1  | Ω t | n t 2 | R t | − N 2 | H ⊺ t R − 1 t H t | − N 2 × exp  − 1 2 tr  Ω t ( Y t − F t ˆ B t ) ⊺ R − 1 t ( Y t − F t ˆ B t )  × d Y l =1 2 π (1 + ( θ t,ℓ q t,ℓ ) 2 ) × d − η t 2 − 1 t, 1 exp  − λ t 2 d t, 1  × N Y j =2 ( τ 2 t ) −|C t,j | 2 ( d t,j ) −|C t,j | 2 exp  − 1 2 τ 2 t d t,j a ( C t,j ) ⊺ t,j a ( C t,j ) t,j  d − η t 2 − 1 t,j exp  − λ t 2 d t,j  . (S.7) Since ( I − A t ) and D t are a low er triangular and diagonal matrix, resp ectively , | Ω t | = | ( I − A t ) ⊺ D − 1 t ( I − A t ) | = N Y j =1 d − 1 2 t,j , where d t,j is the j th elemen t along the diagonal of D t . Recall that a ( C t,j ) j is a column vector comprised of the non-zero elements in row j of A t . Using the results in Section S.1.1 , | Ω t | n t 2 exp  − 1 2 tr  Ω t ( Y t − F t ˆ B t ) ⊺ R − 1 t ( Y t − F t ˆ B t )  S.3 = N Y j =1 d − n t 2 t,j exp  − 1 2 tr  Ω t S ⊺ t S t  = N Y j =1 d − n t 2 t,j exp  − 1 2 d t,j ( S t,j − S t, C t,j a ( C t,j ) t,j ) ⊺ ( S t,j − S t, C t,j a ( C t,j ) t,j )  (S.8) where S t is deﬁned in Prop osition 2.1 . Recall that S t,j and S t, C t,j are the submatrices of S t that include the j th column and each column in the set C t,j , resp ectiv ely . Using the factorization in equation ( S.8 ) and expression ( S.7 ), p ( A , D | θ , y D ) ∝ m Y t =1 { d − ( n t − q t ) 2 t, 1 exp {− 1 2 d t, 1 S ⊺ t, 1 S t, 1 } d − η t 2 − 1 t, 1 exp( − λ t 2 d t, 1 ) × N Y j =2 d − ( n t − q t ) 2 t,j exp  − 1 2 d t,j ( S t,j − S t, C t,j a ( C t,j ) t,j ) ⊺ ( S t,j − S t, C t,j a ( C t,j ) t,j )  × d −|C t,j | 2 t,j exp {− 1 2 τ 2 t d t,j a ( C t,j ) ⊺ t,j a ( C t,j ) t,j } d − η t 2 − 1 t,j exp( − λ t 2 d t,j ) } whic h, after completing the square, is prop ortional to p ( A , D | θ , y D ) ∝ m Y t =1 p ( d t, 1 | θ t , y D ) N Y j =2 p ( a ( C t,j ) t,j | d t,j , θ t , y D ) p ( d t,j | θ t , y D ) ˆ d − ν t 2 t, 1 N Y j =2 | ˆ V t,j | 1 2 ˆ d − ν t 2 t,j , (S.9) with a ( C t,j ) t,j | d t,j , θ t , y D ∼ N |C t,j | ( ˆ a ( C t,j ) t,j , d t,j ˆ V t,j ) (S.10) d t,j | θ t , y D ∼ IG ( ν t 2 , ˆ d t,j 2 ) , (S.11) for ˆ a ( C t,j ) t,j , ˆ V t,j , ˆ d t,j and ν t deﬁned in Prop osition 2.1 . Using expression ( S.9 ) to integrate A and D out of the density in ( S.7 ) yields the expression for p ( θ | y D ) deﬁned in Prop osition 2.1 . S.1.3 Deriv ation in Prop osition 2.2 In this section, we derive the distribution of p ( y ( x 0 ) | θ , y D ) deﬁned in Prop osition 2.2 . Using the separable autoregressiv e mo del ( 2.1 ), the joint distribution of y ( x 0 ) and y D giv en all mo del parameters is a pro duct of matrix-normal distributions across m ﬁdelit y levels: y ( x 0 ) y D !     β , A , D , θ ∼ m Y t =1 MN n t +1 ,N h ⊤ t ( x 0 ) β t + γ t − 1 y ⊤ t − 1 , 1: N ( x 0 ) F t B t ! , r ( x 0 , x 0 ; θ t ) r ⊤ t ( x 0 ) r t ( x 0 ) R t ! , Ω − 1 t ! , (S.12) S.4 where r t ( x 0 ) is deﬁned following expression ( 2.7 ). Since the distribution in ( S.12 ) features a pro duct of matrix normal distributions, standard calculations yield the predictive distribution ( Conti and O’Hagan , 2010 ; Gupta and Nagar , 1999 , Theorem 2.3.12) y ( x 0 ) | β , A , D , θ , y D ∼ m Y t =1 MN 1 ,N ( M t, ∗ , R t, ∗ , Ω − 1 t ) , (S.13) where M t, ∗ =    h ⊺ 1 ( x 0 ) β 1 + r ⊺ 1 ( x 0 ) R − 1 1 ( y 1 − F 1 B 1 ) if t = 1 h ⊺ t ( x 0 ) β t + γ t − 1 y t − 1 ( x 0 ) + r ⊺ t ( x 0 ) R − 1 t ( y t − F t B t ) if t > 1 R t, ∗ = 1 − r ⊺ t ( x 0 ) R − 1 t r t ( x 0 ) . A similar strategy is exploited to in tegrate β out of the distribution in ( S.13 ). Since the distribution of β | A , D , θ , y D deﬁned in ( S.6 ) is prop ortional to a pro duct of matrix normal densities, the density of y t ( x 0 ) | y t − 1 ( x 0 ) , A , D , θ , y D will b e prop ortional to a matrix normal density . T o obtain the moments for the distribution, note that for t > 1 , E [ y t ( x 0 ) | y t − 1 ( x 0 )] = E β t [ E [ y t ( x 0 ) | y t − 1 ( x 0 ) , β t ]] = E β t [ M t, ∗ ] = h ⊺ t ( x 0 ) ˆ β t + γ t − 1 y t − 1 ( x 0 ) + r ⊺ t ( x 0 ) R − 1 t ( y t − F t ˆ B t ) , (S.14) whic h is equal to M t ( x 0 ) ⊺ in ( 2.6 ). A similar calculation provides the expression for M 1 ( x 0 ) . Using the law of total v ariance, V ar ( y t ( x 0 ) | y t − 1 ( x 0 )) = E β t [ V ar ( y t ( x 0 ) | y t − 1 ( x 0 ) , β t )] + V ar β t ( E [ y t ( x 0 ) | y t − 1 ( x 0 ) , β t ]) = E β t [ R t, ∗ ⊗ Ω − 1 t ] + V ar β t ( h ⊺ t ( x 0 ) β t + γ t − 1 y t − 1 ( x 0 ) + r ⊺ t ( x 0 ) R − 1 t ( y t − F t B t )) , (S.15) where ⊗ is the Kroneck er pro duct. After applying the deﬁnition of B t and dropping constan ts, the v ariance term on the right-hand-side of ( S.15 ) is equiv alent to V ar β t ( h ⊺ t ( x 0 ) β t − r ⊺ t ( x 0 ) R − 1 t H t β t ) . F rom Theorem 2.3.10 in Gupta and Nagar ( 1999 ) and expression ( S.6 ), V ar β t ( h ⊺ t ( x 0 ) β t − r ⊺ t ( x 0 ) R − 1 t H t β t ) = ( h ⊺ t ( x 0 ) − r ⊺ t ( x 0 ) R − 1 t H t )( H t R − 1 t H t ) − 1 ( h ⊺ t ( x 0 ) − r ⊺ t ( x 0 ) R − 1 t H t ) ⊺ ⊗ Ω − 1 t . (S.16) After combining expressions ( S.15 ) and ( S.16 ) and using prop erties of the Kroneck er pro duct, V ar ( y t ( x 0 ) | y t − 1 ( x 0 )) = R t ( x 0 ) ⊗ Ω − 1 t , S.5 where R t ( x 0 ) is deﬁned in ( 2.7 ). Th us, y ( x 0 ) | A , D , θ , y D ∼ m Y t =1 MN 1 ,N ( M t ( x 0 ) , R t ( x 0 ) , Ω − 1 t ) . (S.17) Using the factorization in Result 2.1 , p ( y t ( x 0 ) | A , D , θ , y D ) ∼ p ( y t, 1 | A , D , θ , y D ) N Y j =2 p ( y t,j | y t, C t,j , A , D , θ , y D ) , where y t,j | y t, C t,j , A , D , θ , y D ∼ N 1 ( M t,j ( x 0 ) + ( y t, C t,j − M t, C t,j ) ⊺ a C t,j t,j , d t,j R t ( x 0 )) . Then for each ﬁdelity level, the parameters a C t,j t,j and d t,j can b e integrated out using the distributions in ( S.10 ) and ( S.11 ) and standard conjugacy results, yielding the expression for p ( y ( x 0 ) | θ , y D ) in Prop osition 2.2 . S.1.4 P osterior Predictive Distribution in Section 2.3 In this section, we provide expressions for p ( θ w | w D ) and p ( w ( x 0 ) | θ w , w D ) from Section 2.3 . Let ˆ β t − 1 ,ℓ := ( W ⊺ t − 1 R − 1 t,ℓ W t − 1 ) − 1 W ⊺ t − 1 R − 1 t,ℓ w t,ℓ b e the generalized least-squares estimator of β t − 1 ,ℓ . Using standard conjugacy results, it can b e sho wn that p ( θ w | w D ) ∝ p 1 Y ℓ =1 | R 1 ,ℓ | − 1 2 ( w ⊺ 1 ,ℓ R − 1 1 ,ℓ w 1 ,ℓ ) − n 1 2 π ( θ 1 ,ℓ ) m Y t =2 p t Y ℓ =1  | R t,ℓ | − 1 2 | W ⊺ t − 1 R − 1 t,ℓ W t − 1 | − 1 2 (( w t,ℓ − W t − 1 ˆ β t − 1 ,ℓ ) ⊺ R − 1 t,ℓ ( w t,ℓ − W t − 1 ˆ β t − 1 ,ℓ )) − ( n t − p t − 1 ) 2 π ( θ t,ℓ )  Let r t,ℓ ( x 0 ) := r ( x 0 , X t ; θ t,ℓ ) b e an n t × 1 vector of correlations b etw een x 0 and the inputs in X t . The p osterior predictive distribution for the random weigh ts is p ( w ( x 0 ) | θ w , w D ) ∝ p 1 Y ℓ =1 t n 1 ( w 1 ,ℓ ( x 0 ); m 1 ,ℓ ( x 0 ) , ˆ σ 2 1 ,ℓ r 1 ,ℓ ( x 0 )) m Y t =2 p m Y ℓ =1 t n t − p t − 1 ( w t,ℓ ( x 0 ) | w t − 1 ,ℓ ( x 0 ); m t,ℓ ( x 0 ) , ˆ σ 2 t,ℓ r t,ℓ ( x 0 )) , where m 1 ,ℓ ( x 0 ) = r 1 ,ℓ ( x 0 ) ⊺ R − 1 1 ,ℓ w 1 ,ℓ , r 1 ,ℓ ( x 0 ) = 1 − r 1 ,ℓ ( x 0 ) ⊺ R − 1 1 ,ℓ r 1 ,ℓ ( x 0 ) , ˆ σ 2 1 ,ℓ = n − 1 1 w ⊺ 1 ,ℓ R − 1 1 ,ℓ w 1 ,ℓ , S.6 and for t > 1 , m t,ℓ ( x 0 ) = w t − 1 ( x 0 ) ⊺ ˆ β t − 1 ,ℓ + r t,ℓ ( x 0 ) ⊺ R − 1 t,ℓ ( w t,ℓ − W t − 1 ˆ β t − 1 ,ℓ ) , r t,ℓ ( x 0 ) = 1 − r t,ℓ ( x 0 ) ⊺ R − 1 t,ℓ r t,ℓ ( x 0 ) + ( w t − 1 ,ℓ ( x 0 ) − r t,ℓ ( x 0 ) ⊺ R − 1 t,ℓ W t − 1 ) ⊺ ( W ⊺ t − 1 R − 1 t,ℓ W t − 1 ) − 1 ( w t − 1 ,ℓ ( x 0 ) − r t,ℓ ( x 0 ) ⊺ R − 1 t,ℓ W t − 1 ) , ˆ σ 2 1 ,ℓ = ( n t − p t − 1 ) − 1 ( w t,ℓ − W t − 1 ˆ β t − 1 ,ℓ ) ⊺ R − 1 t,ℓ ( w t,ℓ − W t − 1 ˆ β t − 1 ,ℓ ) . S.7 S.2 A dditional Results for the T estb ed Example In Section S.2.1 , we justify several mo deling decisions for the SEP emulator, including the nonstationary cross-co v ariance matrix, autoregressive structure for outputs at multiple ﬁdelity lev els and the choice of h yp erparameters for the prior on A t and D t . In Section S.2.2 , we describ e the implementation details of the NONSEP emulator. In Section S.2.3 , we include visual comparisons of the predictive p erformance for the SEP and NONSEP emulators. S.2.1 Diagnostics and Implemen tation Details for the SEP Emulator 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate 5 10 15 (a) Mean t = 1 . 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate −2.5 −2.0 −1.5 −1.0 −0.5 0.0 (b) Mean t = 2 . 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate 20 40 60 80 (c) V ariance t = 1 . 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate 2.5 5.0 7.5 10.0 (d) V ariance t = 2 . 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate 0 4 8 (e) Residual t = 1 . 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate −1 0 1 (f ) Residual t = 2 . Fig. S.1. Exploratory plots of the sim ulation in the testb ed example. Panels (a) and (b) sho w the empirical mean for level-1 and level-2 outputs, resp ectively; panels (c) and (d) sho w the empirical v ariance of level- 1 and lev el-2 outputs. Panels (e) and (f ) show the residuals at ﬁdelity level 1 and ﬁdelity level 2 after subtracting the mean functions ov er the input setting M = 12 . 49 , D = 0 . 07 , L = 2 . 03 and T = 30 . The ﬁrst tw o rows of Figure S.1 show the empirical mean and v ariance of the level-1 and level-2 outputs. A t t = 2 , we subtracted the lev el-1 outputs from the level-2 outputs b efore calculating the mean and v ariance. At b oth ﬁdelit y levels, the mean and v ariance v ary b y space/time location, with the largest v ariability o ccurring near time 35 . Based on the diﬀerence in v ariance b etw een panels (c) and (d) of Figure S.1 , the level-1 outputs help explain a lot of the v ariability in the lev el-2 outputs, justifying the autoregressive structure in mo del ( 2.1 ). After accoun ting for the lev el-1 outputs in the mean function for the level-2 outputs and removing the spatial trend, the residuals for nearby locations are still correlated (see Panel (f ) of Figure S.1 for an illustration at a single input). After sorting the spatial lo cations, we found that a single nearest neighbor explains most of the output S.8 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate R 2 0.0 0.2 0.4 0.6 0.8 1.0 (a) Level 1 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate R 2 0.0 0.2 0.4 0.6 0.8 1.0 (b) Level 2 Fig. S.2. Output v ariability at each space/time lo cation explained b y the nearest neigh b or’s outputs. Panel (a) sho ws the co eﬃcient of determination for the centered outputs at a lo cation and its nearest neighbor for the level-1 simulator. Panel (b) shows the co eﬃcient of determination for the cen tered diﬀerence b etw een the level 2 and level 1 outputs. v ariability at many lo cations. Figure S.2 shows the co eﬃcient of determination from regressing the centered outputs at a lo cation on the cen tered outputs of its nearest neighbor. The outputs w ere centered at the empirical mean for each lo cation. The results for level 2 in panel (b) of Figure S.2 were calculated using the diﬀerence b etw een the lev el-2 and level-1 outputs. The smallest v alues of R 2 in Figure S.2 tend to o ccur at earlier space/time lo cations, where a p oint’s nearest neighbor is typically farther aw ay than it is at later lo cations. In general, we found that at lo cations where R 2 is close to 1 , when the neigh b or set includes more than one p oint, the centered columns of neighbor outputs are nearly linearly dep endent, causing numerical problems when ev aluating the p osterior distribution in ( 2.5 ). The results for the SEP emulator in T able 1 were obtained after ﬁxing the hyperparameters for the prior on A t and D t deﬁned in ( 2.3 ) at τ 2 t = 1 , η t = 4 and λ t = 2 for t = 1 , 2 . T o mak e predictions, the Metrop olis-Hastings algorithm was run for 30 , 000 iterations (with a burn-in of 3 , 000 samples). T he MCMC samples did not show any lack of conv ergence b ehavior s. The predictive p erformance of the SEP mo del is not particularly sensitive to the choice of hyperparame- ters. W e calculated the RMSPE, CVG ( 95% ) and ALCI ( 95%) using the MAP estimator of θ for three prior distributions on d t,j . Ordered from least to most informative, the hyperparameters were ﬁxed at ( η t , λ t ) = (0 . 2 , 0 . 02) , (4 , 2) , (4 , 0 . 002) . The results are sho wn in T able S.1 with τ 2 t ﬁxed at 1 . The predictive p erformance S.9 at diﬀeren t v alues of τ 2 t similar to what is presented in T able S.1 . The RMSPE of ¯ y 2 ( · ) is omitted from T able S.1 b ecause it is similar for all hyperparameter conﬁgurations. The RMSPE hardly changes with η t and λ t . The a v erage length of a 95% equal-tail credible interv al changes more than the RMSPE for diﬀerent v alues of η t and λ t , but that v ariability is small relative to the ov erall v ariability of y 2 ( x , s ) . T able S.1. Predictiv e p erformance of SEP model for diﬀeren t v alues of the h yp erparameters η t and λ t . Predictions were computed using the MAP estimator of θ . Hyp erparameters Predictiv e p erformance of y 2 ,j ( · ) Predictiv e p erformance of ¯ y 2 ( · ) η 1 λ 1 η 2 λ 2 RMSPE CV G (95%) ALCI (95%) CVG (95%) ALCI (95%) 0.2 0.02 0.2 0.02 0.196 96.50 0.700 93.33 0.147 4 2 0.195 96.32 0.681 93.33 0.145 4 0.002 0.193 96.76 0.760 93.33 0.171 4 2 0.2 0.02 0.196 96.52 0.701 93.33 0.150 4 2 0.195 96.37 0.682 93.33 0.148 4 0.002 0.193 96.83 0.761 93.33 0.174 4 0.002 0.2 0.02 0.195 96.59 0.704 93.33 0.152 4 2 0.195 96.44 0.685 9.333 0.150 4 0.002 0.193 96.90 0.763 93.33 0.176 S.2.2 Diagnostics and Implemen tation Details for the NONSEP emulator In the mean function for the ℓ th level-2 PC w eight, w e included only the ℓ th w eigh t from t = 1 . A scatterplot sho wing the relationship b etw een eac h PC weigh t at the inputs in X 2 for t = 2 (y-axis) and t = 1 (x-axis) is pro vided in Figure S.3 . The plots along the diagonal of Figure S.3 show a strong relationship b etw een the ℓ th PC weigh ts at successive ﬁdelity levels. W e found that including all of the w eights from t = 1 in the mean function for w 2 ,ℓ ( · ) led to a negligible decrease in th e RMSPE relative to a mo del that includes only the ℓ th weigh t from level 1. The scale parameters for the prior assigned to θ w,t,ℓ,i w as ﬁxed at half input x ′ i s range. T o stabilize the p osterior distribution of θ w, 2 , 5 , the scale parameter of θ w, 2 , 5 ,i w as set to a quarter of v ariable x ′ i s range. S.2.3 Visual Comparison for SEP and NONSEP Figure S.4 shows y 2 ( x , s ) ov er the spatial domain for a single input. P anels (a) and (b) of Figure S.5 show the diﬀerence betw een the actual and predicted v alues for the SEP and NONSEP emulators at the test input in Figure S.4 . Both emulators p erform well across the entire space/time domain, but the SEP mo del has sligh tly b etter predictions than the NONSEP mo del, esp ecially near the region with the greatest output v ariability . The credible interv al lengths are similar for b oth the SEP and NONSEP emulator; see row 2 of Figure S.5 . S.10 Fig. S.3. Scatterplot of PC weigh ts for ﬁdelity lev el t = 1 (x-axis) and ﬁdelity lev el t = 2 (y-axis) ev aluated at the inputs in X 2 . S.3 A dditional Results for the Storm Surge Application In Section S.3.1 , we describ e several mo deling decisions for the SEP emulator, including the num b er of lo cations in each neighbor set and the choice of h yp erparameters for the inv erse cross-cov ariance matrix at S.11 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate 5 10 15 20 Fig. S.4. Plot of y 2 ( x , s ) ev aluated at the test input M = 11 . 79 , D = 0 . 12 , L = 2 . 18 and T = 30 for s ∈ [0 . 5 , 5] × [35 , 60] . 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate 0.0 0.1 0.2 (a) Diﬀerence SEP 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate −0.5 0.0 0.5 1.0 (b) Diﬀerence NONSEP 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate 1 2 3 4 5 (c) CI Len. SEP 1 2 3 4 5 35 40 45 50 55 60 Time Spatial coordinate 1 2 3 4 5 (d) CI Len. NONSEP Fig. S.5. Predictive p erformance of the SEP and NONSEP em ulators at the test input M = 11 . 79 , D = 0 . 12 , L = 2 . 18 and T = 30 . Panels (a) and (b) show the diﬀerence b etw een the actual and predicted v alues of y 2 ( x , s ) ; panels (c) and (d) show 95% credible interv al lengths. eac h ﬁdelity level. In Section S.3.2 , we discuss the relationship b etw een the observed PC weigh ts at the tw o ﬁdelit y levels and include several visual diagnostics that are useful for determining the num b er of principal comp onen ts to include in the mo del. S.3.1 Diagnostics and Implemen tation Details for the SEP Emulator P anels (a), (b), (d) and (e) of Figure S.6 show the empirical mean and v ariance of the ADCIR C and ADCIR C + SW AN outputs, indicating nonstationarity at b oth ﬁdelit y levels. Based on the diﬀerence in v ariance b et w een panels (b) and (e), the lev el-1 outputs explain muc h of the v ariability in the level-2 outputs, suggesting that an autoregressive structure is suitable for a mo del of the ADCIRC + SW AN outputs. After remo ving the spatial trend from the outputs, Panels (c) and (f ) indicate that the residuals at each ﬁdelity lev el still exhibit a high degree of spatial correlation. Based on the plots in P anels (a) and (c), the basis function h t ( x ) was ﬁxed at 1 for t = 1 , 2 . After S.12 m 1.0 1.5 2.0 2.5 (a) Mean ADCIRC 0.5 1.0 1.5 (b) V ariance ADCIRC m 0.2 0.4 0.6 (c) Residual ADCIRC m 0.06 0.09 0.12 0.15 0.18 (d) Mean ADCIRC + SW AN 0.0025 0.0050 0.0075 0.0100 0.0125 (e) V ariance ADCIRC + SW AN m −0.05 0.00 0.05 0.10 0.15 (f ) Residual ADCIR C + SW AN Fig. S.6. Comparison of the spatial trend and v ariabilit y for the ADCIRC and ADCRIC + SW AN outputs. P anels (a) and (b) show the empirical mean and v ariance of ADCIR C outputs ev aluated at the inputs in X 1 . Panels (d) and (e) sho w the empirical mean and v ariance for the diﬀerence betw een the ADCIRC + SW AN and ADCIRC outputs o ver X 2 . P anels (c) and (f ) show the residuals for the tw o ﬁdelity levels after remo ving the spatial trend at the input ∆ P = 49 . 82 , R p = 19 . 47 , V f = 5 . 86 , θ = 50 . 23 , B = 1 . 02 and ℓ = ( − 82 . 20 , 26 . 75) . sorting the spatial lo cations using the maximum-minim um distance ordering ( Guinness , 2018 ), the neighbor set C t,j for t = 1 and t = 2 was constructed using the index of the nearest p oint to s j from s 1 , . . . , s j − 1 . Figure S.7 shows the p ercentage of explained v ariability from regressing the centered outputs at a lo cation on the centered outputs of the lo cation’s nearest neighbor. A single neighbor explains nearly all the residual v ariability at most lo cations. Because neighboring outputs are highly correlated, when C t,j includes more than one neighbor, the cen tered columns of y t, C t,j are nearly linearly dep endent for man y lo cations, causing n umerical problems when ev aluating p ( θ | y D ) . Conditional on d t,j , a ( C t,j ) t,j w as assigned a weakly informative prior with the hyperparameter τ 2 t = 10 7 for t = 1 , 2 . The conditional v ariance parameter d t,j w as also given a v ague prior, with η t = 0 . 2 for t = 1 , 2 , λ 1 = 0 . 02 and λ 2 = 0 . 002 . T o obtain samples from the p osterior distribution, the Metrop olis-Hastings algorithm was run for 30 , 000 iterations (with a burn-in of 3 , 000 samples), with every tenth sample sa v ed for predictions. Because neighboring outputs are highly correlated at b oth ﬁdelity levels, the distribution of p ( θ | y D ) is particularly sensitive to the choice of λ t . Small λ t p enalizes large v alues of θ t in the p osterior distribution p ( θ t | y D ) . If λ t is to o small, most of the mass for p ( θ t | y D ) concentrates in a region close to zero, leading S.13 R 2 0.4 0.6 0.8 1.0 (a) ADICRC R 2 0.4 0.6 0.8 1.0 (b) ADCIRC + SW AN Fig. S.7. Output v ariability by lo cation explained by the outputs of a nearest neighbor. Panel (a) shows the the co eﬃcien t of determination ( R 2 ) for the cen tered outputs at a lo cation and its nearest neigh b or for the ADCIR C sim ulator. Panel (b) shows the co eﬃcient of determination for the centered ADCIR C + SW AN outputs after subtracting the ADCIRC outputs. to p o or predictiv e p erformance for the emulator. A t t = 2 , for example, when τ 2 2 = 10 7 , η 2 = 0 . 2 and λ 2 = 0 (these hyperparameters corresp ond to the improp er prior p ( d 2 ,j ) ∝ d − 1 . 1 2 ,j ) , p ( θ 2 | y D ) is maximized at θ 2 = 0 . As λ t increases, the p enalt y for large θ t fades, leading to a ridge in the p osterior distribution of θ t , making it diﬃcult to implement a sampler or p erform MAP estimation. T o p enalize large v alues of θ t , the h yp erparameter λ t m ust b e ﬁxed close to zero. It could b e interesting to inv estigate if the assignment of the join tly robust prior ( Gu , 2019 ) to the range parameters improv es MAP estimation of θ for a wider range of λ t . T o assess the sensitivit y of the SEP mo del to λ 1 , we computed predictions using the MAP estimator of θ for sev eral v alues of λ 1 (see T able S.2 ). The remaining h yp erparameters in the priors for the cross-cov ariance matrices were ﬁxed at τ 2 t = 10 7 , η t = 0 . 2 for t = 1 , 2 and λ 2 = 0 . 002 . The results in T able S.2 are similar for other small v alues of λ 2 b ecause the ADCIRC simulations explain a high p ercentage of the ADCIRC + SW AN output v ariability . As λ 1 increases, the MAP estimate of θ 1 increases for each input v ariable. The predictiv e v ariance in T able S.2 decreases as a function of λ 1 b ecause larger v alues of θ 1 shrink the term R 1 ( x 0 ) , deﬁned in ( 2.7 ), which quantiﬁes predictive uncertaint y for a new input at all lo cations. T able S.2. Predictiv e p erformance of the SEP emulator as a function of λ 1 . The SEP mo del was estimated using the MAP estimator of θ and assessed on 106 test inputs. Marginal predictive p erformance of y 2 ,j ( · ) Predictiv e p erformance of ¯ y 2 ( · ) λ 1 RMSPE CV G ( 95% ) ALCI ( 95% ) CV G (95% ) ALCI ( 95% ) 0.002 0.137 97.89 0.613 96.23 0.509 0.02 0.111 90.61 0.394 89.62 0.275 0.1 0.114 84.59 0.347 82.08 0.214 S.14 S.3.2 Diagnostics and Implemen tation Details for the NONSEP Emulator Figure S.9 shows the relationship b etw een the PC weigh ts for the ADCIRC and ADCIRC + SW AN simu- lators. Because the ADCIRC simulation is a go o d approximation to the ADCIRC + SW AN mo del, the ℓ th w eights at the t wo ﬁdelity levels are highly correlated. Along the diagonal of Figure S.9 , the level-2 and lev el-1 w eigh ts are p ositively correlated, except for PC three. The direction of correlation for the weigh ts at the third principal comp onen t is negative b ecause the basis v ectors in ( 2.9 ) are unique up to a sign. Based on the relationship b etw een the PC weigh ts ov er the tw o ﬁdelity levels shown in Figure S.9 , the mean function for w 2 ,ℓ ( · ) w as simpliﬁed to only include the ℓ th w eight from level 1. After specifying the mo del for the random weigh ts, the scale parameters for the prior assigned to θ w,t,ℓ,i w ere ﬁxed at half v ariable x i ’s range. F or each PC w eigh t and ﬁdelity level, the Metrop olis-Hastings algorithm was run for 60 , 000 iterations (with a burn-in of 6 , 000 samples), with every tenth sample sav ed for predictions. Figure S.8 shows the RMSPE by lo cation for mo dels that include six or seven principal comp onents. P anels (a) and (b) sho w that the sev enth principal comp onent hardly changes the RMSPE at most locations. But in a small region along the southw estern coast, where the absolute v alues of k 2 , 7 are greatest, the addition of the seven th PC w eight leads to a meaningful reduction in the RMSPE. S.15 0.05 0.10 0.15 0.20 (a) RMSPE 6 PCs. 0.05 0.10 0.15 0.20 (b) RMSPE 7 PCs. 0.00 0.01 0.02 0.03 0.04 0.05 (c) | k 2 , 7 | Fig. S.8. Maps of RMSPE for mo dels with diﬀerent num b ers of principal comp onen ts. P anels (a) and (b) sho w the RMSPE ov er testing inputs when the NONSEP mo del includes six and seven PCs, resp ectiv ely . P anel (c) shows the absolute v alues of the basis vector k 2 , 7 in Equation ( 2.9 ). S.16 Fig. S.9. Scatterplot of PC w eights at ﬁdelity level t = 1 (x-axis) and ﬁdelit y level t = 2 (y-axis) ev aluated at the inputs in X 2 . S.17

Multivariate Gaussian process emulation for multifidelity computer models with high-dimensional spatial outputs

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment