Return the given square matrix with a small amount added to the diagonal
        to make it positive semi-definite.

def load_diagonal(cov, load=None):
        """Return the given square matrix with a small amount added to the diagonal
        to make it positive semi-definite.
        """
        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