SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GaussianDistribution.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Heiko Strathmann
8  */
9 
10 #ifdef HAVE_EIGEN3
11 
13 #include <shogun/base/Parameter.h>
15 
16 using namespace shogun;
17 using namespace Eigen;
18 
20 {
21  init();
22 }
23 
25  SGMatrix<float64_t> cov, bool cov_is_factor) :
26  CProbabilityDistribution(mean.vlen)
27 {
28  REQUIRE(cov.num_rows==cov.num_cols, "Covariance must be square but is "
29  "%dx%d\n", cov.num_rows, cov.num_cols);
30  REQUIRE(mean.vlen==cov.num_cols, "Mean must have same dimension as "
31  "covariance, which is %dx%d, but is %d\n",
32  cov.num_rows, cov.num_cols, mean.vlen);
33 
34  init();
35 
36  m_mean=mean;
37 
38  if (!cov_is_factor)
39  {
40  Map<MatrixXd> eigen_cov(cov.matrix, cov.num_rows, cov.num_cols);
42  Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
43 
44  /* compute cholesky */
45  LLT<MatrixXd> llt(eigen_cov);
46  if (llt.info()==NumericalIssue)
47  {
48  /* try to compute smalles eigenvalue for information */
49  SelfAdjointEigenSolver<MatrixXd> solver(eigen_cov);
50  if (solver.info() == Success)
51  {
52  VectorXd ev=solver.eigenvalues();
53  SG_ERROR("Error computing Cholesky of Gaussian's covariance. "
54  "Smallest Eigenvalue is %f.\n", ev[0]);
55  }
56  }
57 
58  eigen_L=llt.matrixL();
59  }
60  else
61  m_L=cov;
62 }
63 
65 {
66 
67 }
68 
70  SGMatrix<float64_t> pre_samples) const
71 {
72  REQUIRE(num_samples>0, "Number of samples (%d) must be positive\n",
73  num_samples);
74 
75  /* use pre-allocated samples? */
76  SGMatrix<float64_t> samples;
77  if (pre_samples.matrix)
78  {
79  REQUIRE(pre_samples.num_rows==m_dimension, "Dimension of pre-samples"
80  " (%d) does not match dimension of Gaussian (%d)\n",
81  pre_samples.num_rows, m_dimension);
82 
83  REQUIRE(pre_samples.num_cols==num_samples, "Number of pre-samples"
84  " (%d) does not desired number of samples (%d)\n",
85  pre_samples.num_cols, num_samples);
86 
87  samples=pre_samples;
88  }
89  else
90  {
91  /* allocate memory and sample from std normal */
92  samples=SGMatrix<float64_t>(m_dimension, num_samples);
93  for (index_t i=0; i<m_dimension*num_samples; ++i)
94  samples.matrix[i]=sg_rand->std_normal_distrib();
95  }
96 
97  /* map into desired Gaussian covariance */
98  Map<MatrixXd> eigen_samples(samples.matrix, samples.num_rows,
99  samples.num_cols);
100  Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
101  eigen_samples=eigen_L*eigen_samples;
102 
103  /* add mean */
104  Map<VectorXd> eigen_mean(m_mean.vector, m_mean.vlen);
105  eigen_samples.colwise()+=eigen_mean;
106 
107  return samples;
108 }
109 
111 {
112  REQUIRE(samples.num_cols>0, "Number of samples must be positive, but is %d\n",
113  samples.num_cols);
114  REQUIRE(samples.num_rows==m_dimension, "Sample dimension (%d) does not match"
115  "Gaussian dimension (%d)\n", samples.num_rows, m_dimension);
116 
117  /* for easier to read code */
118  index_t num_samples=samples.num_cols;
119 
120  float64_t const_part=-0.5 * m_dimension * CMath::log(2 * CMath::PI);
121 
122  /* determinant is product of diagonal elements of triangular matrix */
123  float64_t log_det_part=0;
124  Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
125  VectorXd diag=eigen_L.diagonal();
126  log_det_part=-diag.array().log().sum();
127 
128  /* sample based part */
129  Map<MatrixXd> eigen_samples(samples.matrix, samples.num_rows,
130  samples.num_cols);
131  Map<VectorXd> eigen_mean(m_mean.vector, m_mean.vlen);
132 
133  /* substract mean from samples (create copy) */
134  SGMatrix<float64_t> centred(m_dimension, num_samples);
135  Map<MatrixXd> eigen_centred(centred.matrix, centred.num_rows,
136  centred.num_cols);
137  for (index_t dim=0; dim<m_dimension; ++dim)
138  {
139  for (index_t sample_idx=0; sample_idx<num_samples; ++sample_idx)
140  centred(dim,sample_idx)=samples(dim,sample_idx)-m_mean[dim];
141  }
142 
143  /* solve the linear system based on factorization */
144  MatrixXd solved=eigen_L.triangularView<Lower>().solve(eigen_centred);
145  solved=eigen_L.transpose().triangularView<Upper>().solve(solved);
146 
147  /* one quadratic part x^T C^-1 x for each sample x */
148  SGVector<float64_t> result(num_samples);
149  Map<VectorXd> eigen_result(result.vector, result.vlen);
150  for (index_t i=0; i<num_samples; ++i)
151  {
152  /* i-th centred sample */
153  VectorXd left=eigen_centred.block(0, i, m_dimension, 1);
154 
155  /* inverted covariance times i-th centred sample */
156  VectorXd right=solved.block(0,i,m_dimension,1);
157  result[i]=-0.5*left.dot(right);
158  }
159 
160  /* combine and return */
161  eigen_result=eigen_result.array()+(log_det_part+const_part);
162 
163  /* contains everything */
164  return result;
165 
166 }
167 
168 void CGaussianDistribution::init()
169 {
170  SG_ADD(&m_mean, "mean", "Mean of the Gaussian.", MS_NOT_AVAILABLE);
171  SG_ADD(&m_L, "L", "Lower factor of covariance matrix, "
172  "depending on the factorization type.", MS_NOT_AVAILABLE);
173 }
174 
175 #endif // HAVE_EIGEN3

SHOGUN Machine Learning Toolbox - Documentation