CD-EXCALIBUR-FMS0063-M5.2_SelectionUncertaintyQuantification
============================================================

.. meta::
   :description: technical note

   :keywords: ExCALIBUR,Selection,of,Techniques,for,Uncertainty,Quan-,tification,M5.2,Version,1.00,Abstract,The,report,describes,work,for,ExCALIBUR,project,NEPTUNE,at,Mile-,stone,M5.2.,Since,contractual,issues,with,grantees,prevented,a,timely,,meaningful,exercise,in,code,integration,,acceptance,and,operation,,the,work,consists,of,a,discussion,of,alternative,approaches,to,the,production,of,surrogates.,Neural,networks,are,briefly,considered,,but,the,main,empha-,sis,is,on,a,comparison,between,adaptive,spline,fitting,and,Gaussian,pro-,cesses.,Expected,scaling,costs,with,number,of,datapoints,Nd,are,derived,mathematically,for,both,approaches,,and,verified,by,use,of,opensource,packages,PY-EARTH,(spline,fitting),and,SKLEARN,(Gaussian,processes).,UKAEA,REFERENCE,AND,APPROVAL,SHEET,Client,Reference:,UKAEA,Reference:,CD/EXCALIBUR-FMS/0063,Issue:,Date:,1.00,18,March,2022,Project,Name:,ExCALIBUR,Fusion,Modelling,System,Prepared,By:,Name,and,Department,Wayne,Arter,Joseph,Parker,Signature,N/A,N/A,Date,18,March,2022,18,March,2022,BD,Reviewed,By:,Wayne,Arter,18,March,2022,Project,Technical,Lead,1,1,Introduction,The,main,aim,of,this,activity,is,to,ensure,that,the,UQ,grantees,(successful,bidders,to,ITT,T/AW085/21),perform,their,work,such,as,to,bring,maximum,benefit,to,ExCALIBUR,project,NEPTUNE.,Unfor-,tunately,,despite,the,previous,contractual,agreement,between,UKAEA,and,UCL,around,the,earlier,UQ,and,MOR,ExCALIBUR,grants,,and,although,it,was,decided,that,UCL,were,the,successful,bidders,in,early,August,,contractual,negotiations,were,prolonged.,Thanks,to,the,delay,,at,the,time,of,writing,(Mid-March,2022),,grantee,tasks,have,not,been,completed,to,the,extent,that,many,code,integration,,acceptance,and,operation,tasks,have,been,possible.,It,has,been,previously,noted,[1],that,the,grantees,have,tended,to,focus,on,specific,aspects,of,UQ,,which,although,relevant,to,NEPTUNE,,may,not,answer,all,the,needs,of,the,project.,More-,over,,little,has,been,performed,in,the,way,of,comparison,of,the,proposed,Gaussian,Process,(GP),approach,with,older,methods,of,producing,surrogates,,such,as,spline,interpolation,and,neural,networks,(NNs).,Although,interest,in,the,latter,approach,has,revived,considerably,of,late,,other,UKAEA,and,cross-cutting,ExCALIBUR,projects,are,currently,examining,it.,If,NNs,can,be,shown,reliably,to,solve,the,basis-selection,problem,posed,in,[1,,§,5],,ie.,that,they,prove,capable,of,say,,without,supervision,,selecting,eigenfunction-type,bases,for,adjoint,operator,equations,and,pseu-,dospectral,bases,[2],for,non-normal,operators,,then,NNs,might,take,centre-stage,as,surrogates,for,NEPTUNE.,However,,to,be,fit,for,UQ,,currently,major,questions,need,to,be,answered,relating,to,the,stability,of,the,NN,approach,in,extrapolation,,and,to,the,determination,of,their,errors.,Herein,,discussion,of,NNs,is,relegated,to,a,brief,comparison,with,GPs,in,an,Annex,Section,A.,Thus,this,report,mainly,consists,of,a,comparison,between,spline-based,and,GP,approaches.,Note,that,as,foreshadowed,in,ref,[1,,§,1],,the,two,approaches,may,anyway,be,complementary,,in,that,splines,are,capable,of,fitting,large,data-sets,with,sizes,Nd,of,Petabytes,and,larger,,whereas,GPs,can,only,fit,relatively,small,data-sets.,As,preliminaries,,the,next,Section,2,establishes,mathe-,matically,why,GPs,are,so,costly,,scaling,as,N,3,d,,,and,Section,3,provides,an,introduction,to,the,mathematics,of,splines,,in,particular,of,Multivariate,adaptive,regression,splines,(MARS).,There-,after,,Section,4,contains,the,details,of,the,spline,and,GP,comparison,,and,Section,5,provides,a,summary.,2,2,Scaling,of,Production,of,Gaussian,Process,Surrogates,Fundamentally,the,Gaussian,Process,(GP),is,a,stochastic,interpolant,between,the,sample,data,points.,The,wikipedia,entry,[3],helpfully,shows,that,for,Gaussian,statistics,,the,GP,amounts,to,a,It,will,be,seen,why,there,is,an,optimal,Brownian,walk,which,takes,in,the,sample,points,in,turn.,mean,square,error,(MSE),,in,that,if,the,amplitude,of,the,walk,is,too,small,,then,the,interpolant,will,fail,to,change,enough,to,go,through,consecutive,points,that,have,very,different,sample,values,,too,large,an,MSE,and,it,is,becomes,increasingly,improbable,that,the,walk,will,‘hit’,the,points.,It,also,follows,that,it,is,vital,to,accurately,characterise,the,statistics,of,the,stochastic,process,as,they,are,integral,to,the,interpolation.,Ordinary,Kriging,assumes,that,the,data,has,a,constant,mean,at,least,locally.,There,follows,a,description,of,universal,Kriging,that,tries,to,fit,a,combination,of,stochastic,function,and,a,deterministic,functional,form,(typically,a,polynomial),through,the,sample,point.,The,most,comprehensive,work,appears,to,be,Sacks,et,al,[4],,who,suppose,that,,although,the,system,under,is,deterministic,(so,that,its,output,y,at,a,sample,point,x,should,always,be,the,same),,because,of,errors,y(x),can,be,regarded,as,resulting,from,a,random,function,aka,stochastic,process,as,follows:,Y,(x),=,ΣM,m=1cmfm(x),+,Z(x),=,c,·,f,(x),+,Z(x),(1),where,the,fm(x),are,a,basis,that,may,be,chosen,by,the,modeller,,cm,are,fitting,parameters,to,be,determined,,and,Z(x),is,a,(separate),random,process,at,every,x,with,zero,expected,value.,The,assumption,represented,by,Equation,(1),is,referred,to,as,universal,Kriging,,distinguishing,it,from,the,simple,Kriging,described,in,ref,[5,,§,2.3.1].,Further,suppose,that,the,experimental,design,,ie.,the,set,of,samples,,is,Experimental,design,S,=,(x1,,x2,,.,.,.,,,xN,),(2),leading,to,observed,values,y,=,(y(x1),,y(x2),,.,.,.,,,y(xN,)).,If,y,is,all,the,information,known,about,the,system,,then,it,is,reasonable,further,to,suppose,that,at,any,point,x,,y(x),may,be,predicted,as,a,linear,combination,of,the,observations,where,the,N,-vector,of,coefficients,β(x),is,chosen,to,minimise,the,mean,square-error,˜y(x),=,β(x),·,y,M,SE,(˜y(x)),=,E[(β(x),·,y,−,Y,(x))2],subject,to,the,constraint,E[(β(x),·,y,−,Y,(x))],=,0,Substituting,Equation,(1),in,the,argument,of,Equation,(5),gives,β(x),·,y,−,Y,(x),=,ΣN,n=1βn(x),(c,·,f,(xn),+,Z(xn)),−,(c,·,f,(x),+,Z(x)),giving,E[(β(x),·,y,−,Y,(x))],=,E[c,·,ΣN,n=1,(βn(x)f,(xn),−,f,(x))],since,E[Z(x)],=,0,at,each,value,of,x.,Thus,the,constraint,Equation,(5),implies,ΣN,n=1βn(x)f,(xn),=,f,(x),3,(3),(4),(5),(6),(7),(8),which,may,be,written,in,matrix-vector,form,as,F,T,β(x),=,f,(x),where,the,size,N,×,M,matrix,,,,,F,=,f1(x1),f1(x2),f2(x1),f2(x2),f1(xN,),f2(xN,),fM,(x1),fM,(x2),.,.,.,.,.,.,.,.,.,.,.,.,fM,(xN,),,,,,(9),(10),Substituting,Equation,(8),in,(β(x),·,y,−,Y,(x))2,,all,that,remains,are,the,terms,in,Z,,so,M,SE,(˜y(x)),=,E[(Σn,=,E[(cid:0)Σn,s=1βs(x)Z(xs),−,Z(x))2],s=1βs(x)Z(xs).Σn,q=1βq(x)Z(xq),−,2Σn,(11),s=1βs(x)Z(xs).Z(x),+,Z(x)2(cid:1)](12),Define,the,(normalised),correlation,matrix,by,E[Z(xi)Z(xj)],=,σ2R(xi,,xj),so,that,its,entries,Rij,are,symmetric,by,construction,,and,the,vector,r(x),=,(R(x1,,x),,R(x2,,x),,.,.,.,,,R(xn,,x)),then,M,SE,(˜y(x)),=,σ2,(1,+,β(x)Rβ(x),−,2β(x),·,r(x)),(13),(14),(15),which,must,be,a,minimum,as,a,function(al),of,β,subject,to,the,constraint,system,Equation,(9).,Introducing,λ(x),as,a,M,-vector,of,Lagrange,multipliers,corresponding,to,the,constraint,,it,follows,(since,varying,Equation,(9),gives,rise,simply,to,a,term,F,T,),that,thus,to,obtain,β(x),requires,solving,the,coupled,system,Rβ(x),−,r(x),+,F,λ(x),=,0,(cid:18),0,F,T,F,R,(cid:19),(cid:19),(cid:18)λ(x),β(x),(cid:19),(cid:18)f,(x),r(x),=,which,it,is,convenient,to,rearrange,as,(cid:18),R,F,F,T,0,(cid:19),(cid:19),(cid:18)β(x),λ(x),(cid:19),(cid:18)r(x),f,(x),=,(16),(17),(18),and,once,a,solution,(β(x),,λ(x)),has,been,obtained,then,substituting,for,Rβ,from,Equation,(16),in,Equation,(15),and,using,Equation,(9),gives,M,SE,(˜y(x)),=,σ2,(1,−,β(x),·,r(x),−,λ(x),·,f,(x)),(19),The,theory,of,block,matrices,in,Section,B,specifically,Equation,(47),,gives,the,solution,to,Equa-,tion,(18),in,terms,of,a,negative,inverse,Schur,complement,S,=,−S−1,C,=,(F,T,R−1F,)−1,as,(cid:19),(cid:18)β(x),λ(x),=,(cid:18)R−1,−,R−1F,SF,T,R−1,R−1F,S,SF,T,R−1,−S,(cid:19),(cid:19),(cid:18)r(x),f,(x),=,(cid:18)R−1,−,F,+T,S−1F,+,F,+T,−S,F,+,(cid:19),(cid:18)r(x),f,(x),(cid:19),(20),4,where,the,right,pseudo-inverse,of,F,,,F,+,=,SF,T,R−1,(note,SF,T,R−1.F,=,I),has,been,inroduced.,It,is,also,worth,noting,that,S,only,exists,if,N,≥,M,and,it,is,then,symmetric.,The,formulae,in,the,Sacks,et,al,paper,then,follow,from,the,solution,β(x),=,(I,−,R−1F,SF,T,)R−1r,+,R−1F,Sf,=,(R−1,−,F,+T,S−1F,+)r,+,F,+T,f,λ(x),=,SF,T,R−1r,−,Sf,=,F,+r,−,Sf,upon,introducing,˜c,=,SF,T,R−1y,=,F,+y,(21),(22),(23),and,using,the,fact,that,for,a,general,matrix,F,is,symmetric,,then,for,arbitrary,vectors,x,and,z,,F,T,x,·,z,=,x,·,F,z,=,xF,z,,so,that,Equation,(21),yields,˜y,=,β,·,y,=,˜c,·,f,+,R−1r,·,(y,−,F,˜c),Setting,˜c,=,0,recovers,simple,Kriging.,Another,formula,is,M,SE,(˜y(x)),=,σ2,(cid:0)1,−,r(R−1,−,F,+T,S−1F,+)r,−,2rF,+f,+,f,Sf,(cid:1),(24),(25),3,Spline,theory,This,section,consists,of,background,material,for,splines,based,on,wikipedia,entries,[6,,7,,8].,3.1,Spline,interpolation,Suppose,we,wish,to,fit,a,curve,piecewise,to,a,set,of,Nd,+,1,data,points,{(xi,,yi)}Nd,i=0,,called,knots,,using,Nd,functions,,called,splines,,defined,for,each,i,=,1,,.,.,.,,,Nd,as,y,=,qi(x),on,the,interval,[xi−1,,xi].,We,could,do,this,with,Nd,linear,functions,defined,on,each,of,the,intervals.,However,,this,will,not,generally,give,a,smooth,result,at,the,knots.,In,order,to,have,the,freedom,to,set,the,derivatives,at,the,knots,,we,need,to,consider,cubic,(or,higher-order),polynomials.,These,degrees,of,freedom,are,determined,by,specifying,that,the,splines,must,be,C2-differentiable,at,the,knots,,i(xi),=,q(cid:48),q(cid:48),i+1(xi),,i,(xi),=,q(cid:48)(cid:48),q(cid:48)(cid:48),i+1(xi),,(26),for,i,=,1,,.,.,.,,,Nd.,The,determination,of,spline,coefficients,in,general,is,covered,by,de,Boor,[9].,In,the,cubic,case,,we,have,Nd,intervals,each,containing,a,cubic,,we,have,4n,parameters,to,determine.,At,each,of,the,Nd,−,1,interior,knots,,we,have,4,equations,,qi(xi),=,yi,,qi+1(xi),=,yi,,i+1(xi),,i+1(xi),,q(cid:48),i(xi),=,q(cid:48),i,(xi),=,q(cid:48)(cid:48),q(cid:48)(cid:48),5,(27),(28),(29),(30),giving,4n,−,4,conditions,in,total.,At,the,edge,knots,,we,only,have,q0(x0),=,y0,,qNd(xNd),=,yNd,,(31),(32),and,so,require,only,a,further,conditions,to,be,able,to,determine,all,coefficients.,There,are,different,options,for,doing,this,,but,a,common,choice,is,the,“natural,conditions”,q(cid:48)(cid:48),0,(x0),=,0,,(xNd),=,0.,q(cid:48)(cid:48),Nd,(33),(34),Writing,the,system,symmetrically,,the,problem,becomes,one,of,inverting,a,tridiagonal,system,,an,O(Nd),complexity,procedure.,3.2,Spline,smoothing,This,determines,a,cubic,spline,fitting,function,to,observed,noisy,data,,Yi,=,f,(xi),+,(cid:15)i,,where,(cid:15)i,is,additive,noise,,assumed,to,have,zero,mean.,The,fitting,function,is,ˆf,is,the,function,that,minimizes,Yi,−,ˆf,(xi),(cid:90),(cid:17)2,+,λ,ˆf,(cid:48)(cid:48)(x)2,dx,,Nd(cid:88),(cid:16),i=1,(35),over,all,twice,differentiable,functions.,Detailed,inspection,of,ref,[7],suggests,that,once,again,a,tridiagonal,matrix,inversion,is,needed,to,find,the,vector,of,values,at,the,data,points/knots,,ˆf,(xi).,These,are,then,used,as,coefficient,in,the,expansion,in,spline,function,,ˆf,(x),=,Nd(cid:88),i=1,ˆf,(xi)fi(x),,(36),where,fi(x),are,the,set,of,spline,basis,functions.,3.3,Adaptive,spline,fitting,Multivariate,adaptive,regression,spline,is,also,known,as,MARS.,“It,is,a,non-parametric,regression,technique,and,can,be,seen,as,an,extension,of,linear,models,that,automatically,models,nonlinear-,ities,and,interactions,between,variables.”,The,basic,idea,is,to,replace,a,linear,regression,model,with,a,piecewise,linear,regression,model,,using,hinge,functions,,max(0,,x,−,c),or,max(0,,c,−,x),(37),where,constant,c,is,called,the,knot.,MARS,model,fits,are,sums,of,functions,of,three,types:,6,•,a,constant,,•,a,constant,times,a,hinge,function,,and,•,a,constant,times,a,product,of,hinge,functions.,The,overall,number,of,terms,,and,the,number,of,hinge,functions,that,may,be,multiplied,together,are,parameters,defined,by,the,user;,otherwise,,the,model,is,parameter-free.,(cid:88),ˆf,=,i,aiBi(x),=,a0,+,(cid:88),i,ai,max(0,,x,−,ci),(38),3.3.1,Implementation,The,python,package,SKLEARN,provides,an,implementation,of,MARS.,MARS,is,itself,a,trademark,,so,the,implementation,of,the,MARS,model,is,instead,called,EARTH,,in,Python,PY-EARTH,[10]..,3.3.2,Smoothness,MARS,differs,from,spline,fitting,in,that,MARS,splines,are,not,smooth.,This,is,potentially,attractive,,for,many,solutions,to,fluid,dynamical,problems,exhibit,shocks,or,at,least,internal/boundary,layer,behaviour,where,the,flow,solution,exhibits,local,rapid,variation,that,is,hard,to,capture,in,detail,,hence,resembles,discontinuity.,Nonetheless,,SKLEARN,offers,a,parameter,called,“smooth”,that,produces,a,differentiable,fit,,fol-,lowing,[11,,§,3.7].,That,paper,is,interesting,on,the,need,to,make,smooth,fits:,If,the,intent,is,accurate,estimation,of,the,function,(as,opposed,to,its,derivatives,of,various,orders),,then,there,is,little,to,be,gained,by,imposing,continuity,beyond,that,of,the,function,itself.,If,the,true,underlying,function,nowhere,has,a,very,large,local,second,derivative,,then,a,small,additional,increase,in,accuracy,can,be,achieved,by,imposing,continuous,first,derivatives.,Also,,continuous,(first),derivative,approximations,have,considerably,more,cosmetic,appeal.,There,is,,however,,little,to,be,gained,and,in,fact,much,to,lose,,by,imposing,continuity,beyond,that,of,the,first,derivative,,especially,in,high,dimensional,settings.,The,difficulty,with,higher,order,regression,splines,centers,on,so,called,end,effects.,The,largest,contribution,to,the,average,approximation,error,(2),,(3),emanates,from,locations,x,near,the,boundaries,of,the,domain.,This,phenomenon,is,well,known,even,in,univariate,smoothing,(Nd,=,1),and,is,especially,severe,in,higher,dimensions.,As,the,dimension,of,the,covariate,space,increases,,the,fraction,of,the,data,points,near,a,boundary,increases,rapidly.,Fitting,high,degree,polynomials,(associated,with,high,degree,regression,splines),in,these,regions,leads,to,very,high,variance,of,the,function,estimate,there.,This,is,mainly,due,to,the,lack,of,constraints,on,the,fit,at,the,boundaries,7,This,implies,that,piecewise-linear,models,should,be,adequate,in,higher,dimensions,,however,it,is,not,clear,“how,high,is,high”.,It,seems,likely,that,an,ordinary,2-D,or,3-D,flow,field,will,benefit,from,higher,order,,and,possible,that,a,kinetic,(5-D,or,6-D),representation,will,also,benefit.,A,smooth,fit,could,be,achieved,by,using,higher-order,hinge,functions,,bq(x,−,t),=,max(0,,(x,−,t)q),=,(x,−,t)q,+,,(39),with,the,lowest-order,smooth,fit,being,quadratic,,q,=,2.,However,,this,suffers,from,the,end,effects,mentioned,above.,Therefore,instead,,MARS,uses,a,modified,basis,set,with,formulae,that,appear,complicated,but,in,fact,which,are,essentially,hinge,functions,with,the,corner,smoothed,over,using,cubic,polynomials.,These,basis,functions,are,denoted,C(x|s,=,±1,,t−,,t,,t+),,where,s,corresponds,to,the,sidedness,of,the,function,(the,sign,in,usual,hinge,functions).,Instead,of,being,characterised,by,a,single,knot,,t,,these,are,characterised,by,three,knots,,t−,,t,,and,t+.,They,are,continuous,and,smooth,,but,have,discontinuities,in,the,second,derivative,at,t±.,Central,knots,t,are,placed,in,the,same,location,as,in,the,piecewise,linear,case.,The,side,knots,t±,are,positioned,midway,between,the,central,knots.,4,Comparison,of,MARS,and,GPs,Here,,the,quality,of,approximate,function,fits,to,data,from,MARS,and,Gaussian,Processes,are,compared.,Data,drawn,from,a,smooth,function,is,considered,in,Section,4.1,,while,data,drawn,from,the,same,function,with,additional,Gaussian,noise,is,considered,in,Section,4.2.,Non-smooth,continous,functions,and,discontinuous,functions,will,be,discussed,in,future,reports.,The,tests,described,are,implemented,in,Python,using,the,SKLEARN,packages,[12],for,Gaussian,Processes,and,py-earth,[10],(a,contributed,SKLEARN,package),for,MARS.,The,latter,implements,the,MARS,algorithm,as,described,by,Friedman,[11].,4.1,Fitting,a,noise-free,smooth,function,The,first,test,function,is,f,(x),=,x,sin(x),(40),The,function,f,(x),is,discretized,on,a,uniform,grid,of,1000,points,,and,MARS,and,GPs,are,used,to,create,fits,for,this,function,using,a,data,sample,of,Nd,(cid:28),1000,random,training,points.,The,seed,of,the,random,number,generator,is,fixed,so,that,both,methods,use,the,same,training,points,,and,moreover,that,the,set,of,Nd,+,1,training,points,is,the,same,as,the,set,of,Nd,training,points,with,one,additional,point.,4.1.1,MARS,MARS,fits,to,the,function,(40),are,shown,in,Figure,1,for,various,Nd.,These,are,generated,using,the,smoothed,hinge,basis,functions,C(x|s),,allowing,up,to,100,summed,terms,,and,for,each,term,to,8,(a),(b),(c),(d),Figure,1:,MARS,fits,to,the,function,f,(x),=,x,sin(x),for,training,sets,of,various,sizes,,n,=,Nd.,contain,the,product,of,up,to,100,basis,functions.,This,generates,smoother,fitting,functions,without,a,noticeable,increase,in,computational,cost.,The,fits,produced,by,MARS,require,at,least,Nd,=,10,data,points,to,be,reasonable,,but,many,more,points,(Nd,around,50),to,be,numerically,accurate.,Even,with,many,overall,data,points,,regions,of,the,domain,with,no,data,may,be,fitted,poorly:,see,around,x,=,2,in,the,Nd,=,50,plot,,Figure,1d.,The,error,in,the,fit,and,run,time,of,the,code,are,shown,for,MARS,as,the,blue,line,in,Figure,2.,To,pro-,duce,the,error,,the,MARS,model,is,generated,using,a,training,set,TNd,of,Nd,point,,fM,ARS(x|TNd),,and,evaluated,on,the,full,grid,of,1000,points.,The,absolute,error,|f,(x),−,fM,ARS(x|TNd)|,is,plotted,as,a,function,of,Nd.,The,relative,error,|f,(x),−,fM,ARS(x|TNd)|/|f,(x)|,is,not,used,,as,f,(x),is,near,zero,in,many,places,,making,this,metric,uninformative.,The,error,in,Figure,2a,decreases,with,increasing,Nd,until,reaching,a,plateau,at,around,Nd,=,50.,That,is,,using,a,larger,training,set,than,Nd,=,50,is,not,worthwhile,,as,the,model,fit,does,not,improve.,The,run,time,in,Figure,2b,is,generated,using,Python,internal,timers,to,measure,the,fitting,stage,9,0246810x5040302010010f(x)MARS,on,noise-free,dataset,,n=5f(x)=xsin(x)ObservationsMean,prediction0246810x10505f(x)MARS,on,noise-free,dataset,,n=12f(x)=xsin(x)ObservationsMean,prediction0246810x64202468f(x)MARS,on,noise-free,dataset,,n=20f(x)=xsin(x)ObservationsMean,prediction0246810x64202468f(x)MARS,on,noise-free,dataset,,n=50f(x)=xsin(x)ObservationsMean,prediction,(a),(b),Figure,2:,The,errors,(a),and,run,times,(b),for,MARS,(blue),,Gaussian,Processes,(red),and,Gauss-,ian,Processes,with,adaptive,sampling,(green).,(n,=,Nd,of,text.),of,the,algorithm,only,–,the,timings,do,not,include,sampling,or,function,evaluations.,Moreover,,each,Nd,is,treated,as,a,separate,case.,In,principle,,as,the,training,sets,are,nested,it,is,possible,to,leverage,information,from,previous,runs,to,reduce,the,cost,of,fitting,,but,that,is,not,done,here.,The,cost,of,MARS,increases,roughly,linearly,in,the,range,3,≤,Nd,≤,30,,before,increasing,roughly,like,N,1/3,for,Nd,≥,30.,This,compares,favourably,with,the,theoretical,O(Nd),scaling.,d,4.1.2,Gaussian,Processes,Gaussian,Processes,fits,to,the,function,(40),are,shown,in,Figure,3,for,various,n,=,Nd.,In,this,noise-free,case,these,are,generated,by,interpolating,data,points,using,the,RBF,kernel,,which,is,a,distribution,with,squared,exponential,form,(not,a,radial,basis,function).,Along,with,the,predicted,function,value,,the,GP,fit,also,generates,a,standard,deviation,for,the,fit,,from,which,confidence,intervals,can,be,produced:,95%,confidence,intervals,are,also,shown,as,the,shaded,regions,in,Figure,3.,GPs,produces,meaningful,results,with,as,few,as,three,data,points,(see,Figure,3b).,The,fits,become,extremely,accurate,for,Nd,≥,8,when,interpolating,within,the,training,set,(see,Figure,3c),,and,are,essentially,perfect,fits,by,eye,by,Nd,=,10,(Figure,3d).,The,absolute,error,of,the,GP,fit,is,plotted,as,the,red,line,in,Figure,2a.,The,error,drops,off,extremely,rapidly,in,the,range,8,≤,Nd,≤,15,before,settling,at,a,small,plateau.,The,decay,in,error,is,much,faster,than,that,achieved,by,MARS,(blue),,and,the,final,plateau,reached,is,five,orders,of,magnitude,smaller.,In,the,range,where,the,error,decays,,there,are,large,discrete,jumps,where,the,error,decreases,dramatically,from,Nd,to,Nd,+,1.,These,occur,when,the,new,training,point,added,is,outside,the,range,of,the,existing,training,set.,This,means,that,the,GPs,are,(accurately),interpolating,over,a,10,100101102103n106104102100102errormarsGPGP,adaptiveaveragemaximum100101102103n103102101100101GP,fitting,time,/,secmarsgpgp_adaptive,(a),(b),(c),(d),Figure,3:,GP,fits,to,the,function,f,(x),=,x,sin(x),for,training,sets,of,various,sizes,,n,=,Nd.,larger,range,of,the,domain,and,(less,accurately),extrapolating,over,a,smaller,range.,The,error,is,therefore,much,reduced.,The,run,time,for,GPs,is,plotted,as,the,red,line,in,Figure,2b.,Each,case,takes,around,ten,times,longer,than,a,MARS,fit,,but,for,Nd,(cid:46),200,requires,less,than,1,second.,Moreover,,the,observed,scaling,is,noticeably,more,favourable,than,the,theoretical,scaling,of,O(N,3,d,),,even,at,the,largest,values,of,Nd,considered.,4.1.3,Adaptive,sampling,with,GPs,So,far,,randomly,selected,training,sets,have,been,used.,However,,it,has,been,observed,that,some,regions,of,the,domain,are,fitted,better,than,others,–,in,particular,interpolation,within,the,training,set,is,more,accurate,than,extrapolation.,This,raises,the,prospect,that,faster,convergence,could,be,achieved,,if,the,training,points,could,be,placed,where,needed.,11,0246810x105051015f(x)Gaussian,process,regression,on,noise-free,dataset,,n=2f(x)=xsin(x)ObservationsMean,prediction95%,confidence,interval0246810x105051015f(x)Gaussian,process,regression,on,noise-free,dataset,,n=3f(x)=xsin(x)ObservationsMean,prediction95%,confidence,interval0246810x64202468f(x)Gaussian,process,regression,on,noise-free,dataset,,n=8f(x)=xsin(x)ObservationsMean,prediction95%,confidence,interval0246810x64202468f(x)Gaussian,process,regression,on,noise-free,dataset,,n=10f(x)=xsin(x)ObservationsMean,prediction95%,confidence,interval,(a),(b),(c),(d),Figure,4:,GP,fits,to,the,function,f,(x),=,x,sin(x),for,adaptively,sampled,training,sets,of,various,sizes,,Nd.,To,achieve,this,,the,confidence,intervals,produced,by,the,Gaussian,Process,routine,are,used.,To,go,from,the,training,set,TNd,to,the,set,TNd+1,,the,new,point,is,added,at,the,position,where,the,confidence,interval,is,widest.,T1,is,chosen,as,a,random,point.,An,example,of,this,is,shown,in,Figure,4,where,fits,for,the,cases,Nd,=,7,to,Nd,=,10,are,given.,In,each,plot,,the,red,circle,shows,the,function,value,at,the,position,where,the,confidence,interval,is,widest.,That,point,is,added,to,the,training,set,and,used,in,the,fit,for,Nd,+,1.,Note,that,the,new,data,point,dramatically,improves,the,function,fit,near,that,point,,but,also,generally,improves,the,fit,across,the,whole,domain.,For,example,,comparing,Figures,4a,and,4b,,the,new,point,at,x,=,4,dramatically,reduces,the,error,near,x,=,4,,but,also,reduces,the,error,around,x,=,2,and,x,∈,[5,,8].,Note,also,that,by,Nd,=,7,(Figure,4a),,the,extrema,of,the,grid,x,=,0,and,x,=,10,have,already,been,added,to,the,training,set.,One,property,of,this,sampling,algorithm,is,that,it,tends,to,immediately,add,the,grid,end-points,so,that,there,is,no,extrapolation.,12,0246810x7.55.02.50.02.55.07.510.0f(x)GP,on,noise-free,dataset,,n=7f(x)=xsin(x)ObservationsNew,observationMean,prediction95%,confidence,intervalNext,observation0246810x64202468f(x)GP,on,noise-free,dataset,,n=8f(x)=xsin(x)ObservationsNew,observationMean,prediction95%,confidence,intervalNext,observation0246810x64202468f(x)GP,on,noise-free,dataset,,n=9f(x)=xsin(x)ObservationsNew,observationMean,prediction95%,confidence,intervalNext,observation0246810x64202468f(x)GP,on,noise-free,dataset,,n=10f(x)=xsin(x)ObservationsNew,observationMean,prediction95%,confidence,intervalNext,observation,The,absolute,error,of,this,approach,is,plotted,as,the,green,line,in,Figure,2a.,After,some,initial,noise,,the,error,is,smaller,than,that,for,randomly-sampled,GPs,(red),and,reaches,the,minimum,error,plateau,sooner.,This,is,not,a,dramatic,improvement,however,,as,the,GP,fit,to,a,smooth,function,is,already,very,good.,Preliminary,work,suggests,that,adaptive,sampling,becomes,more,important,when,treating,non-smooth,or,discontinuous,functions.,The,run,time,is,plotted,as,the,green,line,in,Figure,2b.,As,the,time,only,includes,the,GP,fit,and,not,the,sampling,,the,behaviour,is,the,same,as,that,for,GPs,to,within,timing,noise.,Note,that,the,adaptive,sampling,itself,is,cheap,,amounting,to,finding,the,maximum,value,in,a,length,1000,vector,that,is,produced,automatically,by,the,GP,routine.,4.2,Fitting,a,smooth,function,with,noise,In,this,second,test,,the,same,function,f,(x),=,x,sin(x),is,considered,but,now,with,additional,super-,imposed,Gaussian,noise,with,unit,amplitude,and,standard,deviations,σ,∈,{0.01,,0.1,,1,,10}.,No,modification,is,needed,for,MARS,in,this,case,,as,it,does,not,interpolate,the,training,data.,The,GP,approach,however,must,be,modified:,in,the,preceding,case,,GPs,interpolated,the,data,which,leads,to,overfitting,and,inaccurate,models,when,the,data,is,noisy.,To,deal,with,this,,SKLEARN’s,GaussianProcessRegressor,function,is,provided,with,the,parameter,α,=,σ2,,which,it,interprets,as,the,variance,of,additional,Gaussian,measurement,noise,on,the,training,observations.,Note,that,in,a,problem,with,real,data,,this,step,would,entail,making,an,estimate,for,the,variance,of,the,noise,σ2.,An,example,of,a,GP,with,Nd,=,10,and,α,=,σ2,=,1,is,shown,in,Figure,5.,Note,that,unlike,in,the,noise-free,case,,the,mean,prediction,does,not,interpolate,the,observations,,and,the,confidence,intervals,do,not,vanish,at,the,data,points.,Otherwise,,the,behaviour,of,function,fits,is,very,similar,to,those,shown,in,Figures,1,,3,and,4,,and,so,these,are,not,reproduced.,Likewise,,the,timing,plot,is,very,similar,to,that,in,Figure,2b.,The,significant,difference,is,in,the,error,graph,,plotted,in,Figure,6.,The,absolute,error,is,plotted,against,training,set,size,n,=,Nd,for,MARS,(blue),,GPs,(red),and,adaptively-sampled,GPs,(green),,with,different,markers,showing,the,value,of,σ.,As,in,the,noise-free,case,(Figure,2a),,there,is,generally,a,reduction,in,error,with,Nd,,and,then,a,plateau,where,the,error,does,not,reduce,further.,As,before,,the,error,for,GPs,approaches,decreases,more,quickly,and,,in,the,cases,with,smaller,σ,,attain,a,lower,error,plateau,than,the,MARS,approach.,However,,there,is,now,the,significant,difference,that,by,the,very,nature,of,having,noisy,data,,there,is,a,lower,bound,of,the,size,of,σ,that,any,model,would,be,able,to,attain.,This,means,that,when,the,noise,level,is,σ,(cid:38),0.1,,both,GPs,and,MARS,produce,models,of,the,same,accuracy,,but,MARS,does,so,at,around,one-tenth,of,the,computational,cost.,In,cases,where,the,is,less,noise,σ,(cid:46),0.1,,GPs,still,produce,much,more,accurate,models.,13,Figure,5:,Fitting,GP,to,noisy,data,14,0246810x7.55.02.50.02.55.07.510.0f(x)n=10f(x)=xsin(x)ObservationsMean,prediction95%,confidence,interval,Figure,7:,GPs,and,MARS.,Fitting,error,for,noisy,data,,as,n,=,Nd,varies.,Figure,6,15,100101102103n107105103101101errormarsgpgp_adaptive0.00.010.11.010.0,5,Summary,This,report,has,compared,MARS,and,GP,fitting,for,smooth,function,both,with,and,without,added,Gaussian,noise.,In,the,smooth,,noise-free,case,,GPs,converge,very,quickly,to,the,test,function,,requiring,only,∼,10,training,points.,The,computation,is,not,expensive,with,each,case,running,in,less,than,a,sec-,ond,,and,the,observed,scaling,is,rather,weaker,than,the,theoretical,O(N,3,d,).,Using,an,adaptively-,sampled,training,set,does,increase,the,rate,of,convergence,,but,the,improvement,is,not,significant,for,smooth,functions,where,randomly-sampled,points,already,provides,very,quick,convergence.,In,contrast,to,GPs,,MARS,converges,more,slowly,,requiring,∼,50,training,points,,and,has,quite,a,high,error,plateau,so,that,it,never,converges,exactly,to,the,test,function.,The,MARS,computation,is,however,around,ten,times,faster,than,the,GP,computation,,and,the,observed,scaling,is,weaker,than,the,theoretical,O(Nd).,The,behaviour,of,the,fits,is,not,significantly,altered,by,the,addition,of,Gaussian,noise,to,the,test,function.,The,main,difference,is,that,the,standard,deviation,of,the,random,noise,σ,sets,a,lower,bound,on,the,absolute,error,that,any,fitting,algorithm,can,achieve.,When,the,noise,level,is,large,σ,(cid:38),0.1,,neither,algorithm,can,attain,a,good,fit,,and,so,MARS,is,preferable,as,the,cheaper,ap-,proach.,However,,when,the,noise,is,small,σ,(cid:46),0.1,,GPs,are,preferable,as,more,accurate,and,still,relatively,inexpensive.,5.0.1,Future,work,This,report,has,compared,MARS,and,GPs,approaches,for,fitting,smooth,functions,with,and,without,imposed,Gaussian,noise.,These,functions,might,represent,mean,heat,transport,by,each,of,a,set,of,a,turbulent,calculations.,Regarding,the,structure,of,the,turbulent,flows,themselves,,it,is,of,interest,as,to,how,these,approaches,handle,rapidly,varying,,effectively,discontinuous,functions.,This,will,be,the,subject,of,future,reports,,but,preliminary,work,indicates,the,following:,•,The,quality,of,fits,from,all,approaches,deteriorates,around,discontinuities.,•,To,achieve,good,fits,it,is,necessary,to,both,use,adaptively-sampled,training,sets,,and,to,divide,the,domain,into,disjoint,regions,on,which,the,test,function,is,smooth,and,continuous.,•,As,such,,the,ability,to,identify,discontinuities,in,data,will,be,very,important.,The,last,point,implies,eg.,that,fluid,calculations,should,record,information,about,where,techniques,such,as,shock-fitting,have,been,employed,to,eliminate,spurious,oscillation,,since,these,may,in,effect,constitute,discontinuities.,Acknowledgement,The,support,of,the,UK,Meteorological,Office,and,Strategic,Priorities,Fund,is,acknowledged.,16,References,[1],W.,Arter,and,E.,Threlfall.,Management,of,external,curement.,Technical,Report,CD/EXCALIBUR-FMS/0040-M5.1,,UKAEA,,https://github.com/ExCALIBUR-NEPTUNE/Documents/blob/main/reports/ukaea_,reports/CD-EXCALIBUR-FMS0040-M5.1.pdf.,research.,Supports,UQ,Pro-,2021.,[2],L.N.,Trefethen,and,M.,Embree.,Spectra,And,Pseudospectra:,The,Behavior,of,Nonnormal,Matrices,And,Operators.,Princeton,University,Press,,2005.,[3],Kriging,wiki.,https://en.wikipedia.org/wiki/Kriging,,2020.,Accessed:,October,2020.,[4],J.,Sacks,,W.J.,Welch,,T.J.,Mitchell,,and,H.P.,Wynn.,Design,and,analysis,of,computer,experi-,ments.,Statistical,science,,pages,409–423,,1989.,[5],W.,Arter,,E.,Threlfall,,and,J.,Parker.,layer,design,for,Uncer-,Technical,Report,CD/EXCALIBUR-FMS/0024-M3.1.3,,UKAEA,,tainty,Quantification.,2020.,https://github.com/ExCALIBUR-NEPTUNE/Documents/blob/main/reports/ukaea_,reports/CD-EXCALIBUR-FMS0024-M3.1.3.pdf.,Report,on,user,[6],Spline,interpolation.,https://en.wikipedia.org/wiki/Spline_interpolation,,2022.,Ac-,cessed:,March,2022.,[7],Smoothing,spline.,https://en.wikipedia.org/wiki/Smoothing_spline,,2022.,Accessed:,March,2022.,[8],Multivariate,adaptive,regression,spline.,https://en.wikipedia.org/wiki/Multivariate_,adaptive_regression_spline,,2022.,Accessed:,March,2022.,[9],C.,de,Boor.,A,practical,guide,to,splines.,Springer,,New,York,,1978.,[10],Multivariate,adaptive,regression,spline,(MARS),software.,http://contrib.scikit-learn.,org/py-earth/content.html#api,,2022.,Accessed:,March,2022.,[11],J.H.,Friedman.,Multivariate,adaptive,regression,splines.,The,Annals,of,Statistics,,pages,1–67,,1991.,[12],scikit-learn:,Machine,Learning,in,Python.,https://scikit-learn.org/stable/index.html,,2022.,Accessed:,March,2022.,[13],G.H.,Golub,and,C.F.,Van,Loan.,Matrix,computations.,2nd,Edition.,Johns,Hopkins,,Baltimore,,1989.,A,Annex,A:,Comparison,of,GP,and,NN,surrogates,Both,Gaussian,processes,(GP),and,neural,networks,(NN),are,basically,techniques,for,interpola-,tion,,which,may,be,shown,to,be,very,expensive,to,derive,relative,to,spline-fitting.,17,1.,GP,(a),At,each,position/time,x,there,is,a,Gaussian,random,variable.,However,(b),the,correlation,between,different,x,need,not,be,Gaussian,,and,(c),the,correlation,function,R(x,,x(cid:48)),essentially,defines,how,the,interpolant,is,produced,,at,N,3,cost,,N,is,number,of,samples.,(d),The,correlation,function,has,(hyper)parameters,which,may,be,calculated,by,optimisa-,tion,-,maximising,the,likelihood,function,over,relatively,small,number,of,parameters,,but,optimiser,may,not,converge,(e),May,do,amazingly,well,as,a,consequence,of,the,Central,Limit,Theorem,(f),Likely,to,be,stable,when,extrapolating,if,R,smooth,2.,NN,(a),Defined,by,localised,functions,at,nodes,,in,general,functions,of,functions,of,.,.,.,,,where,the,basic,function,usually,has,simple,form,(b),Fitting,by,optimisation,techniques,for,number,of,parameters,∝,number,of,nodes,-,tech-,niques,rely,on,having,many,samples,and,optimiser,may,be,slow,to,converge,or,overfit,(c),Theorems,that,NN,approximation,has,resistance,to,“curse,of,dimensionality”,(d),Usually,unstable,when,extrapolating,(e),Lack,of,routine,estimates,for,error,in,surrogate,B,Annex,B:,Advanced,Results,in,Matrix,Theory,The,inverse,of,a,block,matrix,,ie.,a,matrix,written,as,B,=,(cid:19),(cid:18)A11,A12,A21,A22,structure,(cid:18)N,×,N,N,×,M,M,×,N,M,×,M,(cid:19),(41),may,be,deduced,from,the,formulae,obtained,using,Gaussian,elimination,for,2,simultaneous,linear,equations,(expressed,using,the,2,×,2,matrix,A),(cid:19),(cid:18)x,y,A,=,(cid:19),(cid:18)b1,b2,where,A,=,(cid:19),(cid:18)a11,a12,a21,a22,(42),As,is,well-known,,Gaussian,elimination,proceeds,by,scaling,the,x-coefficient,of,the,first,equation,to,be,equal,to,that,a21,of,the,second,,ie.,by,the,factor,a21a−1,11,and,then,subtracting,the,two,equations.,This,may,be,expressed,as,a,matrix,multiply,,viz.,where,the,coefficient,(cid:18),1,−a21a−1,11,(cid:19),(cid:18)a11,a12,0,a21,a22,1,(cid:19),=,(cid:19),(cid:18)a11,a12,sc,0,sc,=,a22,−,a21a−1,11,a12,18,(43),(44),It,follows,that,the,inverse,of,A,,provided,sc,(cid:54)=,0,is,A−1,=,(cid:18)1,a−1,21,a12,1,0,(cid:19),(cid:18)a−1,11,0,0,s−1,c,(cid:19),(cid:18),1,−a21a−1,11,(cid:19),0,1,The,Schur,complement,is,the,analogue,of,the,sc,factor,for,block,matrices,,viz.,SC,=,A22,−,A21A−1,11,A12,(45),(46),when,writing,Equation,(45),with,aij,→,Aij,(sc,→,SC,,invertible),gives,the,block,factorisation.,Multiplying,out,the,product,gives,B−1,=,B.1,Additional,Results,(cid:18)A−1,11,+,A−1,−S−1,11,A12S−1,C,A21A−1,11,C,A21A−1,11,−A−1,11,A12S−1,C,S−1,C,(cid:19),(47),Using,the,generalised,Sherman-Morrison-Woodbury,formula,,the,first,entry,simplifies.,The,for-,mula,,from,Golub,&,van,Loan,[13,,Eq.,(2.1.4)],is,(T,+,U,V,)−1,=,T,−1,−,T,−1U,(I,+,V,T,−1U,)−1V,T,−1,(48),where,T,is,an,N,×,N,matrix,,U,is,N,×,M,and,V,is,M,×,N,(cf.,use,of,V,T,in,[13]).,To,apply,,write,U,=,−W,C−1,,then,(T,−,W,C−1V,)−1,=,T,−1,+,T,−1W,(cid:2)C−1(I,−,V,T,−1W,C−1)−1(cid:3),V,T,−1,(49),and,observe,that,the,term,is,square,brackets,simplifies,on,remembering,that,where,inverses,of,matrices,Sj,exist,,S−1,2,=,(S2S1)−1,,so,that,1,S−1,(T,−,W,C−1V,)−1,=,T,−1,+,T,−1W,(cid:2)C,−,V,T,−1W,)(cid:3)−1,V,T,−1,With,the,obvious,substitutions,,this,implies,(A11,−,A12A−1,22,A21)−1,=,A−1,11,+,A−1,11,A12[SC]−1A21A−1,11,Thus,B−1,=,There,are,results,of,the,form,(cid:18)(A11,−,A12A−1,−S−1,C,A21A−1,11,22,A21)−1,−A−1,(cid:19),11,A12S−1,C,S−1,C,(T,−,U,V,−1W,)−1U,V,−1,=,T,−1U,(V,−,W,T,−1U,)−1,perhaps,more,memorable,as,(T,−1,−,U,V,W,)−1U,V,−1,=,T,U,(V,−1,−,W,T,U,)−1,(50),(51),(52),(53),(54),These,are,easily,proven,by,multiplication,by,the,terms,in,parenthesis,,when,(in,the,second,case,Equation,(54),,two,terms,in,U,V,W,T,U,cancel.,The,importance,of,Equation,(53),is,that,the,off-,diagonal,terms,in,S−1,C,in,Equation,(52),may,be,replaced,by,terms,expressed,using,(A11−A12A−1,22,A21)−1.,19

:pdfembed:`src:_static/CD-EXCALIBUR-FMS0063-M5.2_SelectionUncertaintyQuantification.pdf, height:1600, width:1100, align:middle`