kriging_loo.py

import openturns as ot
import numpy as np

def compute_kriging_virtual_loo(inputSample, outputSample, covariance_model, basis=None, transformation=None):
    """
     We write here the Virtual cross validation system

      # K   b + F c = y
      # F^t b       = 0

     which leads to :

      # b = K_ y
      # c = Z y

     with :
      # K_ = K^{-1} - K^{-1}.(F.(F^t.K^{-1}.F)^{-1}.F^t.K^{-1}
      # Z = (F^t K^{-1}F)^{-1} F^t K^{-1} y

     Using the Cholesky factor L such as L.L^t=K, and introducing:

     # M = L^{-1}
     # phi = L^{-1}.F = M.F

     It follows that:

     K_ = M^t.M - M^t.phi.(phi^t.phi)^{-1}.phi^t.M
        = M^t.( M - phi.(phi^t.phi)^{-1}.phi^t.M)

     One could also note that solving :

         phi.beta = M

     with 'phi' a rectangular matrix, leads to the normal equation :
     (phi^t.phi).beta = phi^t.M

     We get :

       beta = (phi^t.phi)^{-1}.phi^t.M

     Thus finally:

     K_ = M^t.(M - phi.beta)

     Note that in case basis_size = 0, K_ = K^{-1}.

    """
    input_sample = np.array(inputSample)
    output_sample = np.array(outputSample)
    size = input_sample.shape[0]
    # Basis & its size
    if basis is None:
        basis_size = 0
    else:
        basis_size = len(basis)

    # get input transformation
    if transformation is None:
        normalized_input_sample = input_sample
    else:
        normalized_input_sample = np.array(transformation(input_sample))
    # Covariance matrix
    #correlation_model = ot.CovarianceModel(covariance_model)
    #correlation_model.setAmplitude([1.0])
    K = covariance_model.discretize(normalized_input_sample)
    # Compute the Cholesky factor of K
    L = K.computeCholesky()
    # Define M as inverse of L
    M = L.solveLinearSystem(ot.IdentityMatrix(size))
    # Compute M^t * M
    K_ = M.computeGram()

    if basis_size > 0:
        F = np.ones((size, basis_size))
        for i in range(basis_size):
            out = np.array( basis.build(i)(normalized_input_sample) ).ravel()
            F[:, i] = out[:]
        # Compute phi such as L.phi = F
        phi = L.solveLinearSystem(ot.Matrix(F))
        # We compute beta as described above:
        beta = phi.solveLinearSystem(M)
        # Finally we get the last part of  K_:
        K_ -=  M.transpose() * phi * beta

    # Now we get all elements.
    # Using F.Bachoc, y_loo = (Id - K / sigma2).y where sigma2 = 1 / diag(K_)
    # and sigma_i^2 = 1. / K_{i,i}
    y = np.array(output_sample).ravel()
    sigma2 = 1.0 / np.diag(K_)
    output_loo = y - np.array(K_* y) * sigma2
    return output_loo.reshape(size, 1) ,  np.sqrt(sigma2).reshape(size, 1)

See https://tel.archives-ouvertes.fr/tel-00881002/document, p64 (eq 3.16 & 3.17)
See https://tel.archives-ouvertes.fr/tel-00697026v2/document (eq 1.131, 1.132, 1.133 p38)