Since this question has garnered a lot of views, I thought I'd post code for the final answer that I found, in part, by posting to the Eigen forums. The code uses Boost for the univariate normal and Eigen for matrix handling. It feels rather unorthodox, since it involves using the "internal" namespace, but it works. I'm open to improving it if someone suggests a way.
#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.
*/
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
} // end namespace Eigen
/*
Draw nn samples from a size-dimensional normal distribution
with a specified mean and covariance
*/
void main()
{
int size = 2; // Dimensionality (rows)
int nn=5; // How many samples (columns) to draw
Eigen::internal::scalar_normal_dist_op<double> randN; // Gaussian functor
Eigen::internal::scalar_normal_dist_op<double>::rng.seed(1); // Seed the rng
// Define mean and covariance of the distribution
Eigen::VectorXd mean(size);
Eigen::MatrixXd covar(size,size);
mean << 0, 0;
covar << 1, .5,
.5, 1;
Eigen::MatrixXd normTransform(size,size);
Eigen::LLT<Eigen::MatrixXd> 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
normTransform = cholSolver.matrixL();
} else {
// Use eigen solver
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigenSolver(covar);
normTransform = eigenSolver.eigenvectors()
* eigenSolver.eigenvalues().cwiseSqrt().asDiagonal();
}
Eigen::MatrixXd samples = (normTransform
* Eigen::MatrixXd::NullaryExpr(size,nn,randN)).colwise()
+ mean;
std::cout << "Mean\n" << mean << std::endl;
std::cout << "Covar\n" << covar << std::endl;
std::cout << "Samples\n" << samples << std::endl;
}
