Eigen LLT (Cholesky) fails, while SVD works

Question

I'm trying to reproduce some numpy code on Gaussian Processes (from here) using Eigen. Basically, I need to sample from a multivariate normal distribution:

samples = np.random.multivariate_normal(mu.ravel(), cov, 1)

The mean vector is currently zero, while the covariance matrix is a square matrix generated via the isotropic squared exponential kernel:

sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
return sigma_f**2 * np.exp(-0.5 / l**2 * sqdist)

I can generate the covariance matrix just fine for now (it can probably be cleaned but for now it's a POC):

Matrix2D kernel(const Matrix2D & x1, const Matrix2D & x2, double l = 1.0, double sigma = 1.0) {
    auto dists = ((- 2.0 * (x1 * x2.transpose())).colwise()
                 + x1.rowwise().squaredNorm()).rowwise() +
                 + x2.rowwise().squaredNorm().transpose();
    return std::pow(sigma, 2) * ((-0.5 / std::pow(l, 2)) * dists).array().exp();
}

However, my problems start when I need to sample the multivariate normal.

I've tried using the solution proposed in this accepted answer; however, the decomposition only works with covariance matrices of size up to 30x30; more than that and LLT fails to decompose the matrix. The alternative version provided in the answer also does not work, and creates NaNs. I tried LDLT as well but it also breaks (D contains negative values, so sqrt gives NaN).

Then, I got curious, and I looked into how numpy does this. Turns out the numpy implementation uses SVD decomposition (with LAPACK), rather than Cholesky. So I tried copying their implementation:

// SVD on the covariance matrix generated via kernel function
Eigen::BDCSVD<Matrix2D> solver(covs, Eigen::ComputeFullV);
normTransform = (-solver.matrixV().transpose()).array().colwise() * solver.singularValues().array().sqrt();
// Generate gaussian samples, "randN" is from the multivariate StackOverflow answer
Matrix2D gaussianSamples = Eigen::MatrixXd::NullaryExpr(1, means.size(), randN);

Eigen::MatrixXd samples = (gaussianSamples * normTransform).rowwise() + means.transpose();

The various minuses are me trying to exactly reproduce numpy's results.

In any case, this works perfectly fine, even with large dimensions. I was wondering why Eigen is not able to do LLT, but SVD works. The covariance matrix I use is the same. Is there something I can do to simply use LLT?

EDIT: Here is my full example:

#include <iostream>
#include <random>
#include <Eigen/Cholesky>
#include <Eigen/SVD>
#include <Eigen/Eigenvalues>

using Matrix2D = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor | Eigen::AutoAlign>;
using Vector = Eigen::Matrix<double, Eigen::Dynamic, 1>;

/*
  We need a functor that can pretend it's const,
  but to be a good random number generator
  it needs mutable state.
*/
namespace Eigen {
    namespace internal {
        template<typename Scalar>
        struct scalar_normal_dist_op
        {
            static std::mt19937 rng;    // The uniform pseudo-random algorithm
            mutable std::normal_distribution<Scalar> norm;  // The gaussian combinator

            EIGEN_EMPTY_STRUCT_CTOR(scalar_normal_dist_op)

            template<typename Index>
            inline const Scalar operator() (Index, Index = 0) const { return norm(rng); }
        };

        template<typename Scalar> std::mt19937 scalar_normal_dist_op<Scalar>::rng;

        template<typename Scalar>
        struct functor_traits<scalar_normal_dist_op<Scalar> >
        { enum { Cost = 50 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = false }; };
    } // end namespace internal
} // end namespace Eigen

Matrix2D kernel(const Matrix2D & x1, const Matrix2D & x2, double l = 1.0, double sigma = 1.0) {
    auto dists = ((- 2.0 * (x1 * x2.transpose())).colwise() + x1.rowwise().squaredNorm()).rowwise() + x2.rowwise().squaredNorm().transpose();
    return std::pow(sigma, 2) * ((-0.5 / std::pow(l, 2)) * dists).array().exp();
}

int main() {
    unsigned size = 50;
    unsigned seed = 1;

    Matrix2D X = Vector::LinSpaced(size, -5.0, 4.8);
    Eigen::internal::scalar_normal_dist_op<double> randN; // Gaussian functor
    Eigen::internal::scalar_normal_dist_op<double>::rng.seed(seed); // Seed the rng

    Vector means = Vector::Zero(X.rows());

    auto covs = kernel(X, X);

    Eigen::LLT<Matrix2D> cholSolver(covs);

    // We can only use the cholesky decomposition if
    // the covariance matrix is symmetric, pos-definite.
    // But a covariance matrix might be pos-semi-definite.
    // In that case, we'll go to an EigenSolver
    Eigen::MatrixXd normTransform;
    if (cholSolver.info()==Eigen::Success) {
        std::cout << "Used LLT\n";
        // Use cholesky solver
        normTransform = cholSolver.matrixL();
    } else {
        std::cout << "Broken\n";
        Eigen::BDCSVD<Matrix2D> solver(covs, Eigen::ComputeFullV);
        normTransform = (-solver.matrixV().transpose()).array().colwise() * solver.singularValues().array().sqrt();
    }
    Matrix2D gaussianSamples = Eigen::MatrixXd::NullaryExpr(1, means.size(), randN);

    Eigen::MatrixXd samples = (gaussianSamples * normTransform).rowwise() + means.transpose();

    return 0;
}