從模板類創建對象時出錯

Question

我一直在試圖找到一種方法來從C++中的多元正態分佈中抽取隨機向量，它具有均值向量和協方差矩陣，很像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

template<class _Scalar, int _size> class EigenMultivariateNormal是專門的模板類。第一個class _Scalar要求一個類型，但int _size要求一個整數。

你應該用一個常量int來代替int類型來調用它。其次，你的語法來實例化一個新類EigenMultivariateNormal是錯誤的。 試試這個：

EigenMultivariateNormal<double, 10> mvn (&MN, &CVM); // with 10 is the size 

Answer 2

該代碼有點舊。這是一個更新的，可能改進的版本。可能還有一些不好的事情。例如，我認爲應該更改爲使用MatrixBase而不是實際的Matrix。這可能會讓它優化並更好地決定何時需要分配存儲空間。這也使用命名空間internal這可能是皺眉，但似乎有必要使用Eigen的NullaryExpr這似乎是正確的事情。有可怕的mutable關鍵字的用法。這是必要的，因爲當在NullaryExpr中使用時，Eigen認爲應該是const。 依靠提升也有點煩人。看來在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分解慢了許多，但它是穩定的，即使協方差矩陣是奇異的。如果你知道你的協方差矩陣總是滿的，那麼你可以用它來代替。但是，如果您很少創建分佈並從中抽樣，那麼這可能不是什麼大問題。