Eigen LLT (Cholesky) не работает, а SVD работает - PullRequest
Новая
       72

Eigen LLT (Cholesky) не работает, а SVD работает

0 голосов
/ 12 марта 2020

Я пытаюсь воспроизвести некоторый код numpy в гауссовых процессах (от здесь ), используя Eigen. По сути, мне нужно сделать выборку из многомерного нормального распределения: 

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

Средний вектор в настоящее время равен нулю, в то время как ковариационная матрица представляет собой квадратную матрицу, сгенерированную из изотропного c квадрата экспоненциального ядра: 

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)

Пока что я могу сгенерировать ковариационную матрицу просто отлично (возможно, ее можно почистить, но сейчас это PO C): 

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();
}

Однако мои проблемы начинаются, когда мне нужно для выборки многомерного нормального.

Я пытался использовать решение, предложенное в этом принятом ответе ; однако разложение работает только с ковариационными матрицами размером до 30х30; более того, и LLT не сможет разложить матрицу. Альтернативная версия, представленная в ответе, также не работает и создает NaN. Я также пробовал LDLT, но он также ломается (D содержит отрицательные значения, поэтому sqrt дает NaN).

Затем мне стало любопытно, и я посмотрел, как numpy делает это. Оказывается, реализация numpy использует декомпозицию SVD (с LAPACK), а не Cholesky. Поэтому я попытался скопировать их реализацию: 

// 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();

Различные минусы - это попытка точно воспроизвести результаты numpy.

В любом случае, это прекрасно работает, даже при больших Габаритные размеры. Мне было интересно, почему Эйген не может делать LLT, но SVD работает. Ковариационная матрица, которую я использую, такая же. Есть ли что-то, что я могу сделать, чтобы просто использовать LLT?

РЕДАКТИРОВАТЬ: Вот мой полный пример: 

#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;
}
...