DoctoralThis is a two hours class on conditional Gaussian processes, i.e., kriging.  We attempt to strike a compromise between a good theoretical foundation on Gaussian processes and practical issues (e.g., how to sample a Gaussian process). Note also that the case of Gaussian Processes with trends is discussed. Finally, we try to link kriging to Bayesian regression and Support Vector Machines. Illustrations are based on the R package DiceKriging

Le Riche, Rodolphe

HAL-EMSE

Introduction to krigingRodolphe Le RicheTo cite this version:Rodolphe Le Riche. Introduction to kriging. Doctoral. France. 2014. <cel-01081304>HAL Id: cel-01081304https://hal.archives-ouvertes.fr/cel-01081304Submitted on 7 Nov 2014HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestine´e au de´poˆt et a` la diffusion de documentsscientifiques de niveau recherche, publie´s ou non,e´manant des e´tablissements d’enseignement et derecherche franc¸ais ou e´trangers, des laboratoirespublics ou prive´s.  Introduction to krigingRodolphe Le Riche11 CNRS and Ecole des Mines de St-Etienne, FRClass given as part of the “Modeling and Numerical Methods for Uncertainty Quantification” French-German summer school, Porquerolles, Sept. 2014   Course outline1. Introduction to kriging  (R. Le Riche)  1.1. Gaussian Processes  1.2. Covariance functions  1.3. Conditional Gaussian Processes (kriging)1.3.1. No trend (simple kriging)1.3.2. With trend (universal kriging)  1.4. Issues, links with other methods  Kriging : introductionWhat can be said about possible measures at any x using probabilities ? (Krige, 1951; Matheron, 1963)Here, kriging for regression.Kriging = a family of metamodels (surrogate) with embedded uncertainty model.x5x2x1y5y2y1...Context : scalar measurements ( y1,…, yn) at n positions (x1 ,…, xn) in a d -dimensional space X⊂ℝd  Random processesRandom variable, Yrandom event ( e.g., throw dice )get an instance y Expl : if dice ≤ 3 , y =13 < dice ≤ 5 , y = 2dice = 6  , y = 3 Random process, Y(x)ω∈Ωrandom event ω∈Ω( e.g., wheather )get a function y (x ):x ∈ X ⊂ ℝd → ℝA set of RVs indexed by x  Random processesRepeat the random eventxyEx : three events instances, three y(x)'s.They are different, yet bear strong similarities.  Gaussian processesEach Y(x) follows a Gaussian lawY (x) ∼ N (μ(x) , C (x , x) )∀X = (x1…xn) ∈ X n ⊂ℝd×n , Y = (Y (x1)…Y (xn)) = (Y 1…Y n) ∼ N (μ , C )C ij =C (xi , x j)and, p( y ) = 1(2π)n/2 det1/2(C )exp(−12 ( y−μ)T C−1( y−μ))with probability density function (multi-Normal law), Note : C is called Gram matrix in the SVM class (J.-M. Bourinet)  Gaussian processes  (illustration)x1 x2 x3yY1 Y2 Y3C12C13C23y1 y2y3y(x)µ2µ3µ1n = 3Cij's covariances between Yi and Yj  (linear couplings)Other possible illustration : contour lines of p(y) as ellipses with principal axes as eigen-vectors/values of C-1   Special case : no spatial covarianceY = (Y (x1)…Y (xn)) = (Y 1…Y n) ∼ N (μ , [σ12 0 00 ⋱ 00 0 σn2])Y (x) ∼ N (μ(x ) , σ (x)2)Y (x ) = μ(x ) + Ε(x )  where Ε(x) ∼ N (0 , σ(x )2 )i.e., GP generalize a trend with white noise  “          “          regressionAt n observation points,Example : μ(x) = x2 , σ(x ) = 1  Numerical sampling of a GPC = U D2U TTo plot GP trajectories (in 1 and 2Ds), perform the eigen analysisIn R,  and then notice Y =μ+U DΕ , Ε ∼ N (0, I )( prove that E Y =μ  and Cov (Y ) = C  ,    + a Gaussian vector is fully determined by its 2 first moments )Ceig <- eigen(C)y <- mu + Ceig$vectors %*% diag(sqrt(Ceig$values))  %*% matrix(rnorm(n))( Cf. previous illustrations, functions plotted with large n's)More efficient implementation : cf. Carsten Proppe class, sl. “Discretization of random processes”  Definition of covariance functionsHow do we build C and µ and account for the data points ( xi , yi ) ?→ Start with C .Cov (Y (x ), Y (x ' )) = C (x , x ')The covariance function, a.k.a. kernel, is only a function of x and x' , The kernel defines the covariance matrix throughC ij =C ( xi , x j)C : X × X →ℝ  Valid covariance functionsAll functions C : ℝd×ℝd→ℝ  are not  valid kernels.Kernels must yield positive semidefinite covariance matrices, C :∀u∈ℝn , uT C u⩾ 0Functional view (cf. Mercer's theorem) : let φi be square integrable eigenfunctions of x and λi ≥ 0 the associated eigenvalues, C (x , x ' ) =∑i=1Nλiφi (x)φi(x ') , N=∞  but finite for degenerate kernelsInterpretation : kernels actually work in an N - (possibly infinite) dimensional feature space where a point is (φ1(x ) , … , φN ( x))T( but C(x,x') < 0 may happen )  Stationary covariance functionsThe covariance function depends only on τ = x-x' , but not on the position in spaceCov (Y (x ), Y (x ' )) =C (x , x ') = C (x−x ' ) = C (τ )C (0) = σ2 , C (τ) = σ2 R (τ)   ( R  the correlation)Example :  the squared exponential covariance function (Gaussian)Cov (Y (x ), Y (x ' )) =C (x−x ') = σ2 exp(−∑i=1d ∣x i−xi '∣22θi2 )= σ2∏i=1dexp(−∣τi∣22θi2 )Note : θi  is a length scale ≡ bandwidth (SVM class) ≡ 1/ γ  Gaussian kernel, illustrationyx∥x−x '∥C (x , x ' )Cov (Y (x ), Y (x ' )) =C (x−x ') = σ2 exp(−∑i=1d ∣x i−xi '∣22θi2 )The regularity and frequency content of the trajectories is controlled by the covariance functions. The θi's act as length scales.  Regularity of covariance functionsFor stationary processes, the trajectories y(x) are p times differentiable (in the mean square sense) if C(τ) is 2p times differentiable at τ=0. →The properties of C(τ) at τ=0 define the regularity of the process. Expl : trajectories with Gaussian kernels are infinitely differentiable (very  – unrealistically ? – smooth)   Recycling covariance functionsThe product of kernels is a kernel, (expl, d>1 kernels like the Gaussian kernel)The sum of kernels is a kernel, Let W(x) = a(x) Y(x) , where a(x) is a deterministic function. Then,C (x , x ' ) =∏i=1NC i(x , x ' )C (x , x ' ) =∑i=1NC i (x , x ')Cov (W (x),W (x ')) = a(x)C (x , x ' )a(x ')  Examples of stationary kernelsGeneral form  ,   C (x , x ' ) = σ2∏i=1dR (∣xi−x ' i∣)Gaussian  ,  R (τ) = exp(− τ22θ2 )    (infinitely differentiable trajectories)Matérn ν=5/2  ,  R (τ) = (1+ √5∣τ∣θ + 5 τ23θ2 )exp(−√5∣τ∣θ )    (twice diff. tr.)Matérn ν=3/2   ,  R (τ) = (1+ √3∣τ∣θ )exp(−√3∣τ∣θ )    (once diff. tr.)Power-exponential  ,  R (τ) = exp(−∣τ∣θp) , 0< p≤2(tr. not diff. except for p=2 )(the ones implemented in the DiceKriging R package)Matérn ν=5/2 is the default choice.They are functions of d+1 hyperparameters : σ and the θi's to learn from data.  Tuning of hyperparametersOur first use of the observations ( X ,Y )where  X = (x1…xn)  and  Y = (Y (x1)…Y (xn))Three paths to selecting hyperparameters ( σ and the θi's ) :                              maximum likelihood, cross-validation, Bayesian.Current statistical model  :  Y (x ) ∼ N (μ(x ) , C (x , x))                    equivalently ,  Y ( x) = μ( x)+Z (x)  where  μ(x ) known (deterministic) and  Z (x ) ∼ N (0 , C (x , x ))(discussed here)  Maximum of likelihood estimate (1/2)L(σ ,θ)= p ( y∣σ ,θ) = 1(2π)n/2 det1 /2(C)exp(−12 ( y−μ)T C−1( y−μ))  where  C ij = C (xi , x j ; σ ,θ)= σ2 R (x i , x j ; θ)  and  μ i =μ(xi)Likelihood : the probability of observing the observations as a function of the hyperparametersmaxσ ,θL(σ ,θ) ⇔ minσ ,θ−log L(σ ,θ)⏞mLLNote : compare the likelihood L(σ ,θ) to the regularized loss functionof the SVM class,   L(σ ,θ)= reg_terms(σ ,θ) × exp(−loss_function(σ ,θ))  Maximum of likelihood estimate (2/2)However, the calculations cannot be carried out analytically for θ.So, numerically, minimize the “concentrated” likelihood, mLL(σ ,θ)= n2log (2π)+n log (σ)+12log (det (R))+σ−22( y−μ )T R−1 ( y−μ )∂mLL∂σ = 0 ⇒ σ̂2 = 1n( y−μ )T R(θ)−1 ( y−μ )minθ∈[θmin ,θmax ]mLL(σ̂ (θ) ,θ)   where  θmin > 0A nonlinear, multimodal optimization problem. In R / DiceKriging , solved by a mix of evolutionary – global – and BFGS – local – algorithms.  xy (x )μ(x)μ(x)+1.95 σ̂μ(x)−1.95 σ̂Where do we stand ?Further use the observations ( X , Y=y ). Make such a model interpolate them → last step to kriging.Y (x ) normal with known average μ(x )                   and we have learned the covariance C (x−x ' )  Conditioning of a Gaussian random vectorLet U  and V  be jointly Gaussian random vectors,[UV ]∼ N([μUμV ] , [CU CUVCUVT CV ])then the conditional distribution of U  knowing V=v  is U |V=v ∼ N ( μU+CUV CV−1(v−μV )⏟cond. mean, CU−CUV CV−1 CUVT⏟cond. covar. )✔ the conditional distribution is still Gaussian✔ the conditional covariance does not depend on the observations vand this is all we need ...  Simple krigingApply the conditioning result to the vector[Y (x *)Y (x1)…Y ( xn) ]= [Y ( x *)Y ] ∼ N ([μ( x *)μ ] , [ σ2 C (x * , X )C (x * , X )T C ])which directly yields the simple kriging* formulaY (x *)|Y= y ∼ N (mSK (x *) , vSK ( x *))mSK (x *) = μ( x *)+C (x * , X )C−1( y−μ)vSK ( x *) = σ2−C (x * , X )C−1C (x * , X )T* simple kriging : the trend ( µ(x*) , µ ) is assumed to be knownFor prediction at 1 point,  C (x * , x *) = σ2 R (x *−x *) = σ2For joined prediction at many points, same formula but x * → x * , σ2 → C (x * , x *) , v SK (x *) → CSK (x * , x *) = C *SK( identical to LS-SVR formula ? )  mSK (x)+1.95√vSK (x )mSK (x)−1.95√vSK (x)mSK (x)x1 sample,ySK (x )Simple kriging (illustration)Matérn 5/2 kernel, constant trend –  µ(x) = -0.2 –  θ≈0.02 by MLE. Not a good statistical model for these data.  xSimple kriging (other illustration)Matérn 5/2 kernel, constant trend –  µ(x) = -0.2 –  θ=0.5 fixed a priori. Since Y KS ∼ N (mKS , CKS )  , kriging trajectories can be sampled  Interpolation properties● The kriging mean and trajectories are interpolating the data points● The kriging variance is null at data points Proof for vSK  :    C (xi , X)T = C i  , i-th column of CC−1 C = I ⇒ C−1C i = e i  , i -th basis vector  vSK (xi) = σ2−C iT C−1C i = σ2−C iT ei = σ2−C (xi , xi) = 0Same idea with mSK (xi) = yiY (xi)|Y= y ∼ N ( yi , 0), thus the trajectories are interpolating  Other paths towards kriginglinear  ,  Ŷ (x) =∑i=1nλi( x)Y i = λ(x )T Yunbiased  ,  EŶ (x ) = λ (x)T E Y = EY (x ) = μ(x )best  ,  λ (x ) = arg minλ∈ℝnE (∥Ŷ (x)−Y (x )∥2 )We just used a Bayesian approach (GP conditioning) to justify kriging.Often, kriging is introduced as a best linear estimator :This constrained optimization problem is solved in λ(x) and the kriging equations are recovered throughBut the link with GP interpretation is typically not discussedmK ( x) = E (Ŷ (x )∣Y= y )  and  v K (x ) = E ((Ŷ (x )−Y (x ))2)E (Y (x )∣Y= y ) =?E (Ŷ (x)∣Y= y ) , E((Y (x )−μ(x ))2∣Y= y ) =?E ((Ŷ (x)−Y (x))2)RKHS view : kriging as minimum norm interpolator in the Reproducing Kernel Hilbert Space generated by C(.,.) .   Cf. Nicolas Durrande's PhD thesis.  Universal kriging (1/4)Finally, let us learn the trend within the same framework.The UK statistical model is, Y (x ) = ∑i=1pai (x )βi + Z (x) = a (x)T β + Z (x)where  Z (x) ∼ N (0,C (x , x ))   and  β ∼ N (b , B )i.e., Gaussian prior on the trend weightsConsider the Gaussian vector, new points {Y (x *1)…Y (x *m)}observation points {Y (x1)…Y (xn)}= {Y *Y } and note   [a(x1)T…a(xn)T ]= A[a(x *1)T…a(x *n)T ]= A *  Universal kriging (2/4){Y *Y } ∼ N (A * bA b , [ C *+A * B A *T C (X * , X )+A * B ATC (X , X *)+A B A *T C+A B AT ])and apply the Gaussian vector conditioning formula (see earlier slide). The kriging mean is the conditional average,m(X *) = A * b+ (C (X * , X )+A * B AT ) (C+A B AT )−1 ( y−A b )and the kriging covariance is the conditional covariance,v (X *) = C *+A * B A *T−  (C (X * , X )+A * B AT ) (C+A B AT )−1 (C (X , X *)+A B A *T )The universal kriging formula are obtained by taking the limit of m(X*) and v(X*) when the prior on the weights becomes non informativeλ∈ℝ  , B= λ B̄  , limλ→∞b=0m(X *) = mUK (X *)  ,  limλ→∞b=0v (X *) = vUK (X *)  Universal kriging (3/4)Proof for mUK(.)Two matrix inversion lemma are used :(A−1+B−1)−1 = A−A (A+B)−1 A (1)(Z+UWV T )−1 = Z−1−Z−1 U (W−1+V T Z−1 U )−1V T Z−1 (2)(C+A B AT)−1 =(2) C−1−C−1 A(B−1+AT C−1 A)−1 AT C−1=(1) C−1−C−1 A [(AT C−1 A)−1−(AT C−1 A)−1((AT C−1 A)−1+B)−1(AT C−1 A)−1]AT C−1⇒ A * B AT (C+A B AT)−1 y = A * B((AT C−1 A)−1+B)−1(AT C−1 A)−1 AT C−1 y→λ→∞A * β̂ (3)C (X * , X )(C+A B AT)−1 y = C (X * , X )[C−1−C−1 A(B−1+AT C−1 A)−1 AT C−1] y→λ→∞C (X * , X )C−1( y−A β̂) (4)(3) + (4) yield mUK(.) □  Universal kriging (4/4)mUK (X *) = A(X *)β̂⏟linearpart+ C (X * , X )C−1 ( y−A β̂ )⏟localcorrectingpartvUK (X *)= C *−C (X * , X )C−1C (X , X *)⏟ = vSK+ UT ( AT C−1 A )−1 U⏟addit. partdue to βestimationwhere   β̂ = ( AT C−1 A )−1 AT C−1 ywhere   U = AT C−1C (X , X *)−A(X *)TNote : Ordinary kriging = universal kriging with constant trend  mUK (x)+1.95√vUK (x)mUK (x )−1.95√vUK (x )mUK (x )x1 sample,yUK (x )Universal kriging (illustration)Matérn 5/2 kernel, quadratic trend –  a(x)T = ( 1 , x , x2 ) –  θ≈0.02 by MLE.   x xIllustrative comparison of krigingsComparison of kriging mean ± 1.95 std. dev. for simple, ordinary and universal kriging (SK, OK, UK, resp.)Matérn 5/2 kernel, linear trend –  a(x)T = ( 1 , x ) –  for UK, MLE estimationNote how SK and OK go back to average, how OK uncertainty slightly > that of SK , uncertainty of UK small because of small estimated σθ̂SK = 0.5  , v̂SK = 1.2  ; θ̂OK = 0.5  , v̂OK = 1.2 ; θ̂UK = 0.01  , v̂UK = 0.02  Kriging in more than 1DSample code from DiceKriging R package : library(DiceKriging)library(lhs)n <- 80d <- 6X <- optimumLHS(n, d)X <- data.frame(X)y <- apply(X, 1, hartman6)mlog <- km(design = X,          response = -log(-y))plot(mlog)  Rough sensitivity analysis from krigingEstimated length scales (the θ's) can serve for rapid sensitivity analysis. The larger θ, the least function variation, the least sensitivity.   Expl. in 4D.True function : f ( x) = 2+x12−3 x3+ x13 x32+sin (3 x1 x4)n <- 80d <- 4X <- optimumLHS(n, d)X <- data.frame(X)y <- apply(X, 1, truefunc)model <- km(design = X, response = y, upper=c(5,5,5,5))⇒θ̂1=1.2θ̂2=5θ̂3=2.9θ̂4=1.7one checks that the function is the least sensitive to x2, the most to x1  .  Links with other methods (1/2)Bayesian linear regressionY (x ) =∑i=1NW iφi (x ) + Ε = φ(x )T W + Ε   where   φ (x) = [φ1(x )…φN (x )] , W ∼ N (0 ,Σw ) , Ε ∼ N (0,σn2 )The posterior distribution of the model follows the SK equations, (φ (x )T W ∣Y= y) ∼ N (mSK (x) , v SK (x ))provided      C (x , x ') = φ( x)T Σw φ(x )⇒ kriging is equivalent to Bayesian linear regression in the N (possibly infinite) dimensional feature space (φ1(x ),…,φN (x ))( GPML, p.12 )SVRLS-SVR ~ mK, cf. previous comments, but regularization control (C SVR parameter) embedded in the likelihood.  Links with other methods (2/2)Case of the Gaussian kernelThe feature space associated to the Gaussian kernel is  an infinite series of gaussians with varying centersC (x , x ' ) = σ2 exp(−∥x−x '∥22θ2 )= limN→∞∑c=1Nσ2Nφc(x )φc( x ')              where φc(x ) = exp(−∥x−c∥2θ2 )     (see GPML, p.84)Bayesian linear regression (cont)  Kriging issuesC not invertible If there are linear dependencies between the covariances of subset of data points. This happens in particular if data points are close to each other according to the cov function ← frequent with the very smooth Gaussian kernel.Solution : regularization, either by adding a small >0 quantities on the diagonal or by replacing C-1  by the pseudo-inverse C† .( cf. Mohammadi et al., 2013)Too many data pointsIt is not standard to invert a matrix beyond 1000 data points.Solution : regional kriging, still a research issue.Choosing the covariance function Max. likelihood is a multimodal optimization problem. More generally, how to choose the trend and covariance model ?  BibliographyO'Hagan, A., Bayesian analysis of computer code outputs: a tutorial, Reliability Engineering and System Safety, 91, pp.1290-1300, 2006.Mohammadi, H., Le Riche, R., Touboul, E and Bay, X., An analytic comparison of regularization methods for Gaussian processes, submitted to J. of multivariate statistics, 2014.Rasmussen, C.E., and Williams, C.K.I., Gaussian Processes for Machine Learning, (GPML),  the MIT Press, 2006.Roustant, O., Ginsbourger, D. and Deville, Y., DiceKriging, DiceOptim: Two R Packages for the analysis of Computer Experiments by Kriging-based Metamodeling and Optimization, J. of Statistical Software, 51(1), pp.1-55, 2012.  Compatibility noteFrom now on, Bruno Sudret will write mK ( x) → μŶ (x)v K (x) → σŶ2 (x)

Introduction to kriging

HAL Clermont Université

https://hal.archives-ouvertes.fr/cel-01081304

Introduction to kriging

Abstract

Similar works

Full text

Available Versions

HAL-EMSE

HAL Clermont Université