code: def load_diagonal(cov, load=None):
          
          n, m = cov.shape
          assert n == m, "matrix must be square, but found shape {}".format((n, m))
  
          if load is None:
              load = np.sqrt(np.finfo(np.float64).eps) # machine epsilon
          return cov + np.eye(n) * load
docs: Return the given square matrix with a small amount added to the diagonal
          to make it positive semi-definite.