MultiVariateNormalDistribution.h
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Implements a multi-variate normal distribution with zero mean.
5  *
6  *
7  *
8  * \author T.Voss, O.Krause
9  * \date 2016
10  *
11  *
12  * \par Copyright 1995-2015 Shark Development Team
13  *
14  * <BR><HR>
15  * This file is part of Shark.
16  * <http://image.diku.dk/shark/>
17  *
18  * Shark is free software: you can redistribute it and/or modify
19  * it under the terms of the GNU Lesser General Public License as published
20  * by the Free Software Foundation, either version 3 of the License, or
21  * (at your option) any later version.
22  *
23  * Shark is distributed in the hope that it will be useful,
24  * but WITHOUT ANY WARRANTY; without even the implied warranty of
25  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26  * GNU Lesser General Public License for more details.
27  *
28  * You should have received a copy of the GNU Lesser General Public License
29  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
30  *
31  */
32 #ifndef SHARK_STATISTICS_MULTIVARIATENORMALDISTRIBUTION_H
33 #define SHARK_STATISTICS_MULTIVARIATENORMALDISTRIBUTION_H
34 
36 #include <shark/LinAlg/Cholesky.h>
37 #include <shark/Rng/GlobalRng.h>
38 #include <shark/Core/OpenMP.h>
39 namespace shark {
40 
41 /// \brief Implements a multi-variate normal distribution with zero mean.
43 public:
44  ///\brief Result type of a sampling operation.
45  ///
46  /// The first element is the result of sampling this distribution, the
47  /// second element is the original standard-normally distributed vector drawn
48  /// for sampling purposes.
49  typedef std::pair<RealVector,RealVector> result_type;
50 
51  /// \brief Constructor
52  /// \param [in] Sigma covariance matrix
53  MultiVariateNormalDistribution(RealMatrix const& Sigma ) {
54  m_covarianceMatrix = Sigma;
55  update();
56  }
57 
58  /// \brief Constructor
60 
61  /// \brief Stores/Restores the distribution from the supplied archive.
62  /// \param [in,out] ar The archive to read from/write to.
63  /// \param [in] version Currently unused.
64  template<typename Archive>
65  void serialize( Archive & ar, const std::size_t version ) {
66  ar & BOOST_SERIALIZATION_NVP( m_covarianceMatrix );
67  ar & BOOST_SERIALIZATION_NVP( m_eigenVectors );
68  ar & BOOST_SERIALIZATION_NVP( m_eigenValues );
69  }
70 
71  /// \brief Resizes the distribution. Updates both eigenvectors and eigenvalues.
72  /// \param [in] size The new size of the distribution
73  void resize( std::size_t size ) {
74  m_covarianceMatrix = blas::identity_matrix<double>( size );
75  m_eigenValues = blas::repeat(1.0,size);
76  m_eigenVectors = blas::identity_matrix<double>( size );
77  }
78 
79  /// \brief Accesses the covariance matrix defining the distribution.
80  RealMatrix const& covarianceMatrix() const {
81  return m_covarianceMatrix;
82  }
83 
84  /// \brief Accesses a mutable reference to the covariance matrix
85  /// defining the distribution. Allows for l-value semantics.
86  ///
87  /// ATTENTION: If the reference is altered, update needs to be called manually.
88  RealMatrix& covarianceMatrix() {
89  return m_covarianceMatrix;
90  }
91 
92  /// \brief Sets the covariance matrix and updates the internal variables. This is expensive
93  void setCovarianceMatrix(RealMatrix const& matrix){
94  covarianceMatrix() = matrix;
95  update();
96  }
97 
98  /// \brief Accesses an immutable reference to the eigenvectors of the covariance matrix.
99  const RealMatrix & eigenVectors() const {
100  return m_eigenVectors;
101  }
102 
103  /// \brief Accesses an immutable reference to the eigenvalues of the covariance matrix.
104  RealVector const& eigenValues() const {
105  return m_eigenValues;
106  }
107 
108  /// \brief Accesses a reference to the eigenvalues of the covariance matrix.
109  RealVector& eigenValues(){
110  return m_eigenValues;
111  }
112 
113  /// \brief Samples the distribution.
114  template<class RngType>
115  result_type operator()(RngType& rng) const {
116  RealVector result( m_eigenValues.size(), 0. );
117  RealVector z( m_eigenValues.size() );
118 
119  for( std::size_t i = 0; i < result.size(); i++ ) {
120  z( i ) = gauss(rng, 0., 1. );
121  }
122  for( std::size_t i = 0; i < result.size(); i++ )
123  for( std::size_t j = 0; j < result.size(); j++ )
124  result( i ) += m_eigenVectors( i, j ) * std::sqrt( std::abs( m_eigenValues(j) ) ) * z( j );
125 
126  return( std::make_pair( result, z ) );
127  }
128 
129  /// \brief Calculates the evd of the current covariance matrix.
130  void update() {
131  eigensymm( m_covarianceMatrix, m_eigenVectors, m_eigenValues );
132  }
133 
134 private:
135  RealMatrix m_covarianceMatrix; ///< Covariance matrix of the mutation distribution.
136  RealMatrix m_eigenVectors; ///< Eigenvectors of the covariance matrix.
137  RealVector m_eigenValues; ///< Eigenvalues of the covariance matrix.
138 };
139 
140 /// \brief Multivariate normal distribution with zero mean using a cholesky decomposition
142 public:
143  /// \brief Result type of a sampling operation.
144  ///
145  /// The first element is the result of sampling this distribution, the
146  /// second element is the original standard-normally distributed vector drawn
147  /// for sampling purposes.
148  typedef std::pair<RealVector,RealVector> result_type;
149 
150  /// \brief Constructor
151  /// \param [in] rng the random number generator
152  /// \param [in] covariance covariance matrix
154  setCovarianceMatrix(covariance);
155  }
156 
158 
159  /// \brief Stores/Restores the distribution from the supplied archive.
160  ///\param [in,out] ar The archive to read from/write to.
161  ///\param [in] version Currently unused.
162  template<typename Archive>
163  void serialize( Archive & ar, const std::size_t version ) {
164  ar & BOOST_SERIALIZATION_NVP( m_lowerCholesky);
165  }
166 
167  /// \brief Resizes the distribution. Updates both eigenvectors and eigenvalues.
168  /// \param [in] size The new size of the distribution
169  void resize( std::size_t size ) {
170  m_lowerCholesky = blas::identity_matrix<double>( size );
171  }
172 
173  /// \brief Returns the size of the created vectors
174  std::size_t size()const{
175  return m_lowerCholesky.size1();
176  }
177 
178  /// \brief Sets the new covariance matrix by computing the new cholesky dcomposition
179  void setCovarianceMatrix(RealMatrix const& matrix){
180  choleskyDecomposition(matrix,m_lowerCholesky);
181  }
182 
183  /// \brief Returns the lower cholsky factor.
185  return m_lowerCholesky;
186  }
187 
188  /// \brief Returns the lower cholesky factor.
190  return m_lowerCholesky;
191  }
192 
193  void rankOneUpdate(double alpha, double beta, RealVector const& v){
194  choleskyUpdate(m_lowerCholesky,v,alpha,beta);
195  }
196 
197 
198  template<class RngType, class Vector1, class Vector2>
199  void generate(RngType& rng, Vector1& y, Vector2& z)const{
200  z.resize(size());
201  y.resize(size());
202  for( std::size_t i = 0; i != size(); i++ ) {
203  z( i ) = gauss(rng, 0, 1 );
204  }
205  if(size() > 400){
206  y = z;
207  blas::triangular_prod<blas::lower>(m_lowerCholesky,y);
208  }else{
209  noalias(y) = prod(m_lowerCholesky,z);
210  }
211  }
212 
213  /// \brief Samples the distribution.
214  ///
215  /// Returns a vector pair (y,z) where y=Lz and, L is the lower cholesky factor and z is a vector
216  /// of normally distributed numbers. Thus y is the real sampled point.
217  template<class RngType>
218  result_type operator()(RngType& rng) const {
219  result_type result;
220  RealVector& z = result.second;
221  RealVector& y = result.first;
222  generate(rng,y,z);
223  return result;
224  }
225 
226 private:
227  blas::matrix<double,blas::column_major> m_lowerCholesky; ///< The lower cholesky factor (actually any is okay as long as it is the left)
228 };
229 
230 
231 }
232 
233 #endif