code: def covariance(self):
          
          if self._covariance is None:
              # LiK, _ = dtrtrs(self.woodbury_chol, self._K, lower=1)
              self._covariance = (
              np.atleast_3d(self._K) - np.tensordot(np.dot(np.atleast_3d(self.woodbury_inv).T, self._K), self._K,
                                                    [1, 0]).T).squeeze()
              # self._covariance = self._K - self._K.dot(self.woodbury_inv).dot(self._K)
          return self._covariance
docs: Posterior covariance
          $$
          K_{xx} - K_{xx}W_{xx}^{-1}K_{xx}
          W_{xx} := \texttt{Woodbury inv}
          $$