c++ - 从模板类创建对象时出错

Question

我一直在尝试找到一种方法，从 C++ 中的多元正态分布中抽取随机 vector ，同时具有均值 vector 和协方差矩阵，就像 Matlab 的 mvnrnd 函数的工作原理一样。我在 this page 上找到了实现此功能的类的相关代码，但我在编译它时遇到了一些问题。我已经创建了一个包含在 main.cpp 中的头文件，并且我正在尝试创建 EigenMultivariateNormal 类的对象:

MatrixXd MN(10,1);
MatrixXd CVM(10,10);

EigenMultivariateNormal <double,int> (&MN,&CVM) mvn;

问题是我在编译时遇到与模板相关的错误:

error: type/value mismatch at argument 2 in template parameter list for ‘template<class _Scalar, int _size> class EigenMultivariateNormal’    
error:   expected a constant of type ‘int’, got ‘int’    
error: expected ‘;’ before ‘mvn’

我对如何使用模板只有一个肤浅的想法，而且我绝不是cpp专家，所以我想知道我到底做错了什么？显然我应该在代码中的某个地方有一个 const 。

Answer 1

该代码有点旧。这是一个更新的、可能改进的版本。也许还有一些不好的事情。例如，我认为应该更改为使用 MatrixBase 而不是实际的 Matrix。这可能会让它优化并更好地决定何时需要分配存储空间。这也使用了命名空间 internal ，这可能会让人皱眉，但似乎有必要利用 Eigen 的 NullaryExpr ，这似乎是正确的做法。这里使用了可怕的 mutable 关键字。这是必要的，因为 Eigen 认为在 NullaryExpr 中使用时应该是 const。 依赖boost也有点烦人。似乎在 C++11 中 necessary functions已成为标准。在类代码下方，有一个简短的使用示例。

类eigenmultivariatenormal.hpp

#ifndef __EIGENMULTIVARIATENORMAL_HPP
#define __EIGENMULTIVARIATENORMAL_HPP

#include <Eigen/Dense>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>    

/*
  We need a functor that can pretend it's const,
  but to be a good random number generator 
  it needs mutable state.  The standard Eigen function 
  Random() just calls rand(), which changes a global
  variable.
*/
namespace Eigen {
namespace internal {
template<typename Scalar> 
struct scalar_normal_dist_op 
{
  static boost::mt19937 rng;                        // The uniform pseudo-random algorithm
  mutable boost::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> 
boost::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
/**
    Find the eigen-decomposition of the covariance matrix
    and then store it for sampling from a multi-variate normal 
*/
template<typename Scalar, int Size>
class EigenMultivariateNormal
{
  Matrix<Scalar,Size,Size> _covar;
  Matrix<Scalar,Size,Size> _transform;
  Matrix< Scalar, Size, 1> _mean;
  internal::scalar_normal_dist_op<Scalar> randN; // Gaussian functor


public:
  EigenMultivariateNormal(const Matrix<Scalar,Size,1>& mean,const Matrix<Scalar,Size,Size>& covar)
  {
    setMean(mean);
    setCovar(covar);
  }

  void setMean(const Matrix<Scalar,Size,1>& mean) { _mean = mean; }
  void setCovar(const Matrix<Scalar,Size,Size>& covar) 
  {
    _covar = covar;

    // Assuming that we'll be using this repeatedly,
    // compute the transformation matrix that will
    // be applied to unit-variance independent normals

    /*
    Eigen::LDLT<Eigen::Matrix<Scalar,Size,Size> > cholSolver(_covar);
    // 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
    if (cholSolver.info()==Eigen::Success) {
      // Use cholesky solver
      _transform = cholSolver.matrixL();
    } else {*/
      SelfAdjointEigenSolver<Matrix<Scalar,Size,Size> > eigenSolver(_covar);
      _transform = eigenSolver.eigenvectors()*eigenSolver.eigenvalues().cwiseMax(0).cwiseSqrt().asDiagonal();
    /*}*/

  }

  /// Draw nn samples from the gaussian and return them
  /// as columns in a Size by nn matrix
  Matrix<Scalar,Size,-1> samples(int nn)
  {
    return (_transform * Matrix<Scalar,Size,-1>::NullaryExpr(Size,nn,randN)).colwise() + _mean;
  }
}; // end class EigenMultivariateNormal
} // end namespace Eigen
#endif

这是一个使用它的简单程序:

#include <fstream>
#include "eigenmultivariatenormal.hpp"
#ifndef M_PI
#define M_PI REAL(3.1415926535897932384626433832795029)
#endif

/**
  Take a pair of un-correlated variances.
  Create a covariance matrix by correlating 
  them, sandwiching them in a rotation matrix.
*/
Eigen::Matrix2d genCovar(double v0,double v1,double theta)
{
  Eigen::Matrix2d rot = Eigen::Rotation2Dd(theta).matrix();
  return rot*Eigen::DiagonalMatrix<double,2,2>(v0,v1)*rot.transpose();
}

void main()
{
  Eigen::Vector2d mean;
  Eigen::Matrix2d covar;
  mean << -1,0.5; // Set the mean
  // Create a covariance matrix
  // Much wider than it is tall
  // and rotated clockwise by a bit
  covar = genCovar(3,0.1,M_PI/5.0);

  // Create a bivariate gaussian distribution of doubles.
  // with our chosen mean and covariance
  Eigen::EigenMultivariateNormal<double,2> normX(mean,covar);
  std::ofstream file("samples.txt");

  // Generate some samples and write them out to file 
  // for plotting
  file << normX.samples(1000).transpose() << std::endl;
}

这是显示结果的图。

使用SelfAdjointEigenSolver可能比Cholesky分解慢很多，但它是稳定的，即使协方差矩阵是奇异的。如果您知道协方差矩阵始终是满的，那么您可以使用它。但是，如果您很少创建发行版并经常从中采样，那么这可能不是什么大问题。