code: def _compute_B_statistics(self, K, W, log_concave, *args, **kwargs):
          
          #w = GPy.util.diag.view(W)
          #W[:] = np.where(w<1e-6, 1e-6, w)
  
          #B = I + KW
          B = np.eye(K.shape[0]) + np.dot(K, W)
          #Bi, L, Li, logdetB = pdinv(B)
          Bi = np.linalg.inv(B)
  
          #K_Wi_i = np.eye(K.shape[0]) - mdot(W, Bi, K)
          K_Wi_i = np.dot(W, Bi)
  
          #self.K_Wi_i_brute = np.linalg.inv(K + np.linalg.inv(W))
          #self.B = B
          #self.Bi = Bi
          Ki_W_i = np.dot(Bi, K)
  
          sign, logdetB = np.linalg.slogdet(B)
          return K_Wi_i, sign*logdetB, Bi, Ki_W_i
docs: Rasmussen suggests the use of a numerically stable positive definite matrix B
          Which has a positive diagonal element and can be easyily inverted
  
          :param K: Prior Covariance matrix evaluated at locations X
          :type K: NxN matrix
          :param W: Negative hessian at a point (diagonal matrix)
          :type W: Vector of diagonal values of hessian (1xN)
          :returns: (K_Wi_i, L_B, not_provided)