EigenRand  0.3.0
MvNormal.h
Go to the documentation of this file.
1 
12 #ifndef EIGENRAND_MVDISTS_MVNORMAL_H
13 #define EIGENRAND_MVDISTS_MVNORMAL_H
14 
15 namespace Eigen
16 {
17  namespace Rand
18  {
19  namespace detail
20  {
21  template<typename _Scalar, Index Dim, typename _Mat>
22  Matrix<_Scalar, Dim, Dim> get_lt(const MatrixBase<_Mat>& mat)
23  {
24  LLT<Matrix<_Scalar, Dim, Dim>> llt(mat);
25  if (llt.info() == Eigen::Success)
26  {
27  return llt.matrixL();
28  }
29  else
30  {
31  SelfAdjointEigenSolver<Matrix<_Scalar, Dim, Dim>> solver(mat);
32  if (solver.info() != Eigen::Success)
33  {
34  throw std::runtime_error{ "The matrix cannot be solved!" };
35  }
36  return solver.eigenvectors() * solver.eigenvalues().cwiseMax(0).cwiseSqrt().asDiagonal();
37  }
38  }
39 
40  class FullMatrix {};
41  class LowerTriangular {};
42  class InvLowerTriangular {};
43  }
44 
45  constexpr detail::FullMatrix full_matrix;
46  constexpr detail::LowerTriangular lower_triangular;
47  constexpr detail::InvLowerTriangular inv_lower_triangular;
48 
49 
56  template<typename _Scalar, Index Dim = -1>
57  class MvNormalGen : public MvVecGenBase<MvNormalGen<_Scalar, Dim>, _Scalar, Dim>
58  {
59  Matrix<_Scalar, Dim, 1> mean;
60  Matrix<_Scalar, Dim, Dim> lt;
61  StdNormalGen<_Scalar> stdnorm;
62 
63  public:
72  template<typename MeanTy, typename LTTy>
73  MvNormalGen(const MatrixBase<MeanTy>& _mean, const MatrixBase<LTTy>& _lt, detail::LowerTriangular)
74  : mean{ _mean }, lt{ _lt }
75  {
76  eigen_assert(_mean.cols() == 1 && _mean.rows() == _lt.rows() && _lt.rows() == _lt.cols());
77  }
78 
87  template<typename MeanTy, typename CovTy>
88  MvNormalGen(const MatrixBase<MeanTy>& _mean, const MatrixBase<CovTy>& _cov, detail::FullMatrix = {})
89  : MvNormalGen{ _mean, detail::template get_lt<_Scalar, Dim>(_cov), lower_triangular }
90  {
91  }
92 
93  MvNormalGen(const MvNormalGen&) = default;
94  MvNormalGen(MvNormalGen&&) = default;
95 
96  Index dims() const { return mean.rows(); }
97 
98  template<typename Urng>
99  inline auto generate(Urng&& urng, Index samples)
100  -> decltype((lt * stdnorm.template generate<Matrix<_Scalar, Dim, -1>>(mean.rows(), samples, std::forward<Urng>(urng))).colwise() + mean)
101  {
102  return (lt * stdnorm.template generate<Matrix<_Scalar, Dim, -1>>(mean.rows(), samples, std::forward<Urng>(urng))).colwise() + mean;
103  }
104 
105  template<typename Urng>
106  inline auto generate(Urng&& urng)
107  -> decltype((lt * stdnorm.template generate<Matrix<_Scalar, Dim, 1>>(mean.rows(), 1, std::forward<Urng>(urng))).colwise() + mean)
108  {
109  return (lt * stdnorm.template generate<Matrix<_Scalar, Dim, 1>>(mean.rows(), 1, std::forward<Urng>(urng))).colwise() + mean;
110  }
111  };
112 
121  template<typename MeanTy, typename CovTy>
122  inline auto makeMvNormGen(const MatrixBase<MeanTy>& mean, const MatrixBase<CovTy>& cov)
123  -> MvNormalGen<typename MatrixBase<MeanTy>::Scalar, MatrixBase<MeanTy>::RowsAtCompileTime>
124  {
125  static_assert(
126  std::is_same<typename MatrixBase<MeanTy>::Scalar, typename MatrixBase<CovTy>::Scalar>::value,
127  "Derived::Scalar must be the same with `mean` and `cov`'s Scalar."
128  );
129  static_assert(
130  MatrixBase<MeanTy>::RowsAtCompileTime == MatrixBase<CovTy>::RowsAtCompileTime &&
131  MatrixBase<CovTy>::RowsAtCompileTime == MatrixBase<CovTy>::ColsAtCompileTime,
132  "assert: mean.RowsAtCompileTime == cov.RowsAtCompileTime && cov.RowsAtCompileTime == cov.ColsAtCompileTime"
133  );
134  return { mean, cov };
135  }
136 
145  template<typename MeanTy, typename LTTy>
146  inline auto makeMvNormGenFromLt(const MatrixBase<MeanTy>& mean, const MatrixBase<LTTy>& lt)
147  -> MvNormalGen<typename MatrixBase<MeanTy>::Scalar, MatrixBase<MeanTy>::RowsAtCompileTime>
148  {
149  static_assert(
150  std::is_same<typename MatrixBase<MeanTy>::Scalar, typename MatrixBase<LTTy>::Scalar>::value,
151  "Derived::Scalar must be the same with `mean` and `lt`'s Scalar."
152  );
153  static_assert(
154  MatrixBase<MeanTy>::RowsAtCompileTime == MatrixBase<LTTy>::RowsAtCompileTime &&
155  MatrixBase<LTTy>::RowsAtCompileTime == MatrixBase<LTTy>::ColsAtCompileTime,
156  "assert: mean.RowsAtCompileTime == lt.RowsAtCompileTime && lt.RowsAtCompileTime == lt.ColsAtCompileTime"
157  );
158  return { mean, lt, lower_triangular };
159  }
160 
167  template<typename _Scalar, Index Dim>
168  class WishartGen : public MvMatGenBase<WishartGen<_Scalar, Dim>, _Scalar, Dim>
169  {
170  Index df;
171  Matrix<_Scalar, Dim, Dim> chol;
172  StdNormalGen<_Scalar> stdnorm;
173  std::vector<ChiSquaredGen<_Scalar>> chisqs;
174  public:
182  template<typename ScaleTy>
183  WishartGen(Index _df, const MatrixBase<ScaleTy>& _lt, detail::LowerTriangular)
184  : df{ _df }, chol{ _lt }
185  {
186  eigen_assert(df > chol.rows() - 1);
187  eigen_assert(chol.rows() == chol.cols());
188 
189  for (Index i = 0; i < chol.rows(); ++i)
190  {
191  chisqs.emplace_back(df - i);
192  }
193  }
194 
202  template<typename ScaleTy>
203  WishartGen(Index _df, const MatrixBase<ScaleTy>& _scale, detail::FullMatrix = {})
204  : WishartGen{ _df, detail::template get_lt<_Scalar, Dim>(_scale), lower_triangular }
205  {
206  eigen_assert(_scale.rows() == _scale.cols());
207  }
208 
209  WishartGen(const WishartGen&) = default;
210  WishartGen(WishartGen&&) = default;
211 
212  Index dims() const { return chol.rows(); }
213 
214  template<typename Urng>
215  inline Matrix<_Scalar, Dim, -1> generate(Urng&& urng, Index samples)
216  {
217  const Index dim = chol.rows();
218  const Index normSamples = samples * dim * (dim - 1) / 2;
219  using ArrayXs = Array<_Scalar, -1, 1>;
220  Matrix<_Scalar, Dim, -1> rand_mat(dim, dim * samples), tmp(dim, dim * samples);
221 
222  _Scalar* ptr = tmp.data();
223  Map<ArrayXs>{ ptr, normSamples } = stdnorm.template generate<ArrayXs>(normSamples, 1, urng);
224  for (Index j = 0; j < samples; ++j)
225  {
226  for (Index i = 0; i < dim - 1; ++i)
227  {
228  rand_mat.col(i + j * dim).tail(dim - 1 - i) = Map<ArrayXs>{ ptr, dim - 1 - i };
229  ptr += dim - 1 - i;
230  }
231  }
232 
233  for (Index i = 0; i < dim; ++i)
234  {
235  _Scalar* ptr = tmp.data();
236  Map<ArrayXs>{ ptr, samples } = chisqs[i].template generate<ArrayXs>(samples, 1, urng).sqrt();
237  for (Index j = 0; j < samples; ++j)
238  {
239  rand_mat(i, i + j * dim) = *ptr++;
240  }
241  }
242 
243  for (Index j = 0; j < samples; ++j)
244  {
245  rand_mat.middleCols(j * dim, dim).template triangularView<StrictlyUpper>().setZero();
246  }
247  tmp.noalias() = chol * rand_mat;
248 
249  for (Index j = 0; j < samples; ++j)
250  {
251  auto t = tmp.middleCols(j * dim, dim);
252  rand_mat.middleCols(j * dim, dim).noalias() = t * t.transpose();
253  }
254  return rand_mat;
255  }
256 
257  template<typename Urng>
258  inline Matrix<_Scalar, Dim, -1> generate(Urng&& urng)
259  {
260  const Index dim = chol.rows();
261  const Index normSamples = dim * (dim - 1) / 2;
262  using ArrayXs = Array<_Scalar, -1, 1>;
263  Matrix<_Scalar, Dim, Dim> rand_mat(dim, dim);
264  Map<ArrayXs>{ rand_mat.data(), normSamples } = stdnorm.template generate<ArrayXs>(normSamples, 1, urng);
265 
266  for (Index i = 0; i < dim / 2; ++i)
267  {
268  rand_mat.col(dim - 2 - i).tail(i + 1) = rand_mat.col(i).head(i + 1);
269  }
270 
271  for (Index i = 0; i < dim; ++i)
272  {
273  rand_mat(i, i) = chisqs[i].template generate<Array<_Scalar, 1, 1>>(1, 1, urng).sqrt()(0);
274  }
275  rand_mat.template triangularView<StrictlyUpper>().setZero();
276 
277  auto t = (chol * rand_mat).eval();
278  return (t * t.transpose()).eval();
279  }
280  };
281 
289  template<typename ScaleTy>
290  inline auto makeWishartGen(Index df, const MatrixBase<ScaleTy>& scale)
291  -> WishartGen<typename MatrixBase<ScaleTy>::Scalar, MatrixBase<ScaleTy>::RowsAtCompileTime>
292  {
293  static_assert(
294  MatrixBase<ScaleTy>::RowsAtCompileTime == MatrixBase<ScaleTy>::ColsAtCompileTime,
295  "assert: scale.RowsAtCompileTime == scale.ColsAtCompileTime"
296  );
297  return { df, scale };
298  }
299 
307  template<typename LTTy>
308  inline auto makeWishartGenFromLt(Index df, const MatrixBase<LTTy>& lt)
309  -> WishartGen<typename MatrixBase<LTTy>::Scalar, MatrixBase<LTTy>::RowsAtCompileTime>
310  {
311  static_assert(
312  MatrixBase<LTTy>::RowsAtCompileTime == MatrixBase<LTTy>::ColsAtCompileTime,
313  "assert: lt.RowsAtCompileTime == lt.ColsAtCompileTime"
314  );
315  return { df, lt, lower_triangular };
316  }
317 
324  template<typename _Scalar, Index Dim>
325  class InvWishartGen : public MvMatGenBase<InvWishartGen<_Scalar, Dim>, _Scalar, Dim>
326  {
327  Index df;
328  Matrix<_Scalar, Dim, Dim> chol;
329  StdNormalGen<_Scalar> stdnorm;
330  std::vector<ChiSquaredGen<_Scalar>> chisqs;
331  public:
339  template<typename ScaleTy>
340  InvWishartGen(Index _df, const MatrixBase<ScaleTy>& _ilt, detail::InvLowerTriangular)
341  : df{ _df }, chol{ _ilt }
342  {
343  eigen_assert(df > chol.rows() - 1);
344  eigen_assert(chol.rows() == chol.cols());
345 
346  for (Index i = 0; i < chol.rows(); ++i)
347  {
348  chisqs.emplace_back(df - i);
349  }
350  }
351 
359  template<typename ScaleTy>
360  InvWishartGen(Index _df, const MatrixBase<ScaleTy>& _scale, detail::FullMatrix = {})
361  : InvWishartGen{ _df, detail::template get_lt<_Scalar, Dim>(_scale.inverse()), inv_lower_triangular }
362  {
363  eigen_assert(_scale.rows() == _scale.cols());
364  }
365 
366  InvWishartGen(const InvWishartGen&) = default;
367  InvWishartGen(InvWishartGen&&) = default;
368 
369  Index dims() const { return chol.rows(); }
370 
371  template<typename Urng>
372  inline Matrix<_Scalar, Dim, -1> generate(Urng&& urng, Index samples)
373  {
374  const Index dim = chol.rows();
375  const Index normSamples = samples * dim * (dim - 1) / 2;
376  using ArrayXs = Array<_Scalar, -1, 1>;
377  Matrix<_Scalar, Dim, -1> rand_mat(dim, dim * samples), tmp(dim, dim * samples);
378 
379  _Scalar* ptr = tmp.data();
380  Map<ArrayXs>{ ptr, normSamples } = stdnorm.template generate<ArrayXs>(normSamples, 1, urng);
381  for (Index j = 0; j < samples; ++j)
382  {
383  for (Index i = 0; i < dim - 1; ++i)
384  {
385  rand_mat.col(i + j * dim).tail(dim - 1 - i) = Map<ArrayXs>{ ptr, dim - 1 - i };
386  ptr += dim - 1 - i;
387  }
388  }
389 
390  for (Index i = 0; i < dim; ++i)
391  {
392  _Scalar* ptr = tmp.data();
393  Map<ArrayXs>{ ptr, samples } = chisqs[i].template generate<ArrayXs>(samples, 1, urng).sqrt();
394  for (Index j = 0; j < samples; ++j)
395  {
396  rand_mat(i, i + j * dim) = *ptr++;
397  }
398  }
399 
400  for (Index j = 0; j < samples; ++j)
401  {
402  rand_mat.middleCols(j * dim, dim).template triangularView<StrictlyUpper>().setZero();
403  }
404  tmp.noalias() = chol * rand_mat;
405 
406  auto id = Eigen::Matrix<_Scalar, Dim, Dim>::Identity(dim, dim);
407  for (Index j = 0; j < samples; ++j)
408  {
409  auto t = tmp.middleCols(j * dim, dim);
410  auto u = rand_mat.middleCols(j * dim, dim);
411  u.noalias() = t.template triangularView<Lower>().solve(id);
412  t.noalias() = u.transpose() * u;
413  }
414  return tmp;
415  }
416 
417  template<typename Urng>
418  inline Matrix<_Scalar, Dim, -1> generate(Urng&& urng)
419  {
420  const Index dim = chol.rows();
421  const Index normSamples = dim * (dim - 1) / 2;
422  using ArrayXs = Array<_Scalar, -1, 1>;
423  Matrix<_Scalar, Dim, Dim> rand_mat(dim, dim);
424  Map<ArrayXs>{ rand_mat.data(), normSamples } = stdnorm.template generate<ArrayXs>(normSamples, 1, urng);
425 
426  for (Index i = 0; i < dim / 2; ++i)
427  {
428  rand_mat.col(dim - 2 - i).tail(i + 1) = rand_mat.col(i).head(i + 1);
429  }
430 
431  for (Index i = 0; i < dim; ++i)
432  {
433  rand_mat(i, i) = chisqs[i].template generate<Array<_Scalar, 1, 1>>(1, 1, urng).sqrt()(0);
434  }
435  rand_mat.template triangularView<StrictlyUpper>().setZero();
436 
437  auto t = (chol * rand_mat).eval();
438  auto id = Eigen::Matrix<_Scalar, Dim, Dim>::Identity(dim, dim);
439  rand_mat.noalias() = t.template triangularView<Lower>().solve(id);
440 
441  return (rand_mat.transpose() * rand_mat).eval();
442  }
443  };
444 
452  template<typename ScaleTy>
453  inline auto makeInvWishartGen(Index df, const MatrixBase<ScaleTy>& scale)
454  -> InvWishartGen<typename MatrixBase<ScaleTy>::Scalar, MatrixBase<ScaleTy>::RowsAtCompileTime>
455  {
456  static_assert(
457  MatrixBase<ScaleTy>::RowsAtCompileTime == MatrixBase<ScaleTy>::ColsAtCompileTime,
458  "assert: scale.RowsAtCompileTime == scale.ColsAtCompileTime"
459  );
460  return { df, scale };
461  }
462 
470  template<typename ILTTy>
471  inline auto makeInvWishartGenFromIlt(Index df, const MatrixBase<ILTTy>& ilt)
472  -> InvWishartGen<typename MatrixBase<ILTTy>::Scalar, MatrixBase<ILTTy>::RowsAtCompileTime>
473  {
474  static_assert(
475  MatrixBase<ILTTy>::RowsAtCompileTime == MatrixBase<ILTTy>::ColsAtCompileTime,
476  "assert: ilt.RowsAtCompileTime == ilt.ColsAtCompileTime"
477  );
478  return { df, ilt, inv_lower_triangular };
479  }
480  }
481 }
482 
483 #endif
Eigen::Rand::InvWishartGen::InvWishartGen
InvWishartGen(Index _df, const MatrixBase< ScaleTy > &_ilt, detail::InvLowerTriangular)
Construct a new inverse Wishart generator.
Definition: MvNormal.h:340
Eigen::Rand::MvNormalGen
Generator of real vectors on a multivariate normal distribution.
Definition: MvNormal.h:58
Eigen::Rand::makeMvNormGen
auto makeMvNormGen(const MatrixBase< MeanTy > &mean, const MatrixBase< CovTy > &cov) -> MvNormalGen< typename MatrixBase< MeanTy >::Scalar, MatrixBase< MeanTy >::RowsAtCompileTime >
helper function constructing Eigen::Rand::MvNormal
Definition: MvNormal.h:122
Eigen::Rand::makeWishartGenFromLt
auto makeWishartGenFromLt(Index df, const MatrixBase< LTTy > &lt) -> WishartGen< typename MatrixBase< LTTy >::Scalar, MatrixBase< LTTy >::RowsAtCompileTime >
helper function constructing Eigen::Rand::WishartGen
Definition: MvNormal.h:308
Eigen::Rand::InvWishartGen
Generator of real matrices on a inverse Wishart distribution.
Definition: MvNormal.h:326
Eigen::Rand::MvVecGenBase
Base class of all multivariate random vector generators.
Definition: Basic.h:84
Eigen::Rand::MvMatGenBase
Base class of all multivariate random matrix generators.
Definition: Basic.h:129
Eigen::Rand::makeMvNormGenFromLt
auto makeMvNormGenFromLt(const MatrixBase< MeanTy > &mean, const MatrixBase< LTTy > &lt) -> MvNormalGen< typename MatrixBase< MeanTy >::Scalar, MatrixBase< MeanTy >::RowsAtCompileTime >
helper function constructing Eigen::Rand::MvNormal
Definition: MvNormal.h:146
Eigen::Rand::WishartGen
Generator of real matrices on a Wishart distribution.
Definition: MvNormal.h:169
Eigen::Rand::makeWishartGen
auto makeWishartGen(Index df, const MatrixBase< ScaleTy > &scale) -> WishartGen< typename MatrixBase< ScaleTy >::Scalar, MatrixBase< ScaleTy >::RowsAtCompileTime >
helper function constructing Eigen::Rand::WishartGen
Definition: MvNormal.h:290
Eigen::Rand::makeInvWishartGenFromIlt
auto makeInvWishartGenFromIlt(Index df, const MatrixBase< ILTTy > &ilt) -> InvWishartGen< typename MatrixBase< ILTTy >::Scalar, MatrixBase< ILTTy >::RowsAtCompileTime >
helper function constructing Eigen::Rand::InvWishartGen
Definition: MvNormal.h:471
Eigen::Rand::makeInvWishartGen
auto makeInvWishartGen(Index df, const MatrixBase< ScaleTy > &scale) -> InvWishartGen< typename MatrixBase< ScaleTy >::Scalar, MatrixBase< ScaleTy >::RowsAtCompileTime >
helper function constructing Eigen::Rand::InvWishartGen
Definition: MvNormal.h:453
Eigen::Rand::MvNormalGen::MvNormalGen
MvNormalGen(const MatrixBase< MeanTy > &_mean, const MatrixBase< LTTy > &_lt, detail::LowerTriangular)
Construct a new multivariate normal generator from lower triangular matrix of decomposed covariance.
Definition: MvNormal.h:73
Eigen::Rand::WishartGen::WishartGen
WishartGen(Index _df, const MatrixBase< ScaleTy > &_lt, detail::LowerTriangular)
Construct a new Wishart generator from lower triangular matrix of decomposed scale.
Definition: MvNormal.h:183
Eigen::Rand::StdNormalGen
Generator of reals on the standard normal distribution.
Definition: NormalExp.h:27
Eigen::Rand::MvNormalGen::MvNormalGen
MvNormalGen(const MatrixBase< MeanTy > &_mean, const MatrixBase< CovTy > &_cov, detail::FullMatrix={})
Construct a new multivariate normal generator from covariance matrix.
Definition: MvNormal.h:88
Eigen::Rand::InvWishartGen::InvWishartGen
InvWishartGen(Index _df, const MatrixBase< ScaleTy > &_scale, detail::FullMatrix={})
Construct a new inverse Wishart generator.
Definition: MvNormal.h:360
Eigen::Rand::WishartGen::WishartGen
WishartGen(Index _df, const MatrixBase< ScaleTy > &_scale, detail::FullMatrix={})
Construct a new Wishart generator from scale matrix.
Definition: MvNormal.h:203