1

我正在尝试使用 Eigen在高斯过程（从这里）上重现一些 numpy 代码。基本上，我需要从多元正态分布中抽样：

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

均值向量当前为零，而协方差矩阵是通过各向同性平方指数核生成的方阵：

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)

我现在可以很好地生成协方差矩阵（它可能可以被清理，但现在它是一个 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();
}

但是，当我需要对多元正态进行采样时，我的问题就开始了。

我试过使用这个接受的答案中提出的解决方案；但是，分解仅适用于最大为 30x30 的协方差矩阵；不仅如此，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 的结果。

无论如何，这工作得很好，即使是大尺寸。我想知道为什么 Eigen 不能做 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;
}
4

0 回答 0