SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2013 Heiko Strathmann 00008 */ 00009 00010 #ifdef HAVE_EIGEN3 00011 00012 #include <shogun/distributions/classical/GaussianDistribution.h> 00013 #include <shogun/base/Parameter.h> 00014 #include <shogun/mathematics/eigen3.h> 00015 00016 using namespace shogun; 00017 using namespace Eigen; 00018 00019 CGaussianDistribution::CGaussianDistribution() : CProbabilityDistribution() 00020 { 00021 init(); 00022 } 00023 00024 CGaussianDistribution::CGaussianDistribution(SGVector<float64_t> mean, 00025 SGMatrix<float64_t> cov, bool cov_is_factor) : 00026 CProbabilityDistribution(mean.vlen) 00027 { 00028 REQUIRE(cov.num_rows==cov.num_cols, "Covariance must be square but is " 00029 "%dx%d\n", cov.num_rows, cov.num_cols); 00030 REQUIRE(mean.vlen==cov.num_cols, "Mean must have same dimension as " 00031 "covariance, which is %dx%d, but is %d\n", 00032 cov.num_rows, cov.num_cols, mean.vlen); 00033 00034 init(); 00035 00036 m_mean=mean; 00037 00038 if (!cov_is_factor) 00039 { 00040 Map<MatrixXd> eigen_cov(cov.matrix, cov.num_rows, cov.num_cols); 00041 m_L=SGMatrix<float64_t>(cov.num_rows, cov.num_cols); 00042 Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols); 00043 00044 /* compute cholesky */ 00045 LLT<MatrixXd> llt(eigen_cov); 00046 if (llt.info()==NumericalIssue) 00047 { 00048 /* try to compute smalles eigenvalue for information */ 00049 SelfAdjointEigenSolver<MatrixXd> solver(eigen_cov); 00050 if (solver.info() == Success) 00051 { 00052 VectorXd ev=solver.eigenvalues(); 00053 SG_ERROR("Error computing Cholesky of Gaussian's covariance. " 00054 "Smallest Eigenvalue is %f.\n", ev[0]); 00055 } 00056 } 00057 00058 eigen_L=llt.matrixL(); 00059 } 00060 else 00061 m_L=cov; 00062 } 00063 00064 CGaussianDistribution::~CGaussianDistribution() 00065 { 00066 00067 } 00068 00069 SGMatrix<float64_t> CGaussianDistribution::sample(int32_t num_samples, 00070 SGMatrix<float64_t> pre_samples) const 00071 { 00072 REQUIRE(num_samples>0, "Number of samples (%d) must be positive\n", 00073 num_samples); 00074 00075 /* use pre-allocated samples? */ 00076 SGMatrix<float64_t> samples; 00077 if (pre_samples.matrix) 00078 { 00079 REQUIRE(pre_samples.num_rows==m_dimension, "Dimension of pre-samples" 00080 " (%d) does not match dimension of Gaussian (%d)\n", 00081 pre_samples.num_rows, m_dimension); 00082 00083 REQUIRE(pre_samples.num_cols==num_samples, "Number of pre-samples" 00084 " (%d) does not desired number of samples (%d)\n", 00085 pre_samples.num_cols, num_samples); 00086 00087 samples=pre_samples; 00088 } 00089 else 00090 { 00091 /* allocate memory and sample from std normal */ 00092 samples=SGMatrix<float64_t>(m_dimension, num_samples); 00093 for (index_t i=0; i<m_dimension*num_samples; ++i) 00094 samples.matrix[i]=sg_rand->std_normal_distrib(); 00095 } 00096 00097 /* map into desired Gaussian covariance */ 00098 Map<MatrixXd> eigen_samples(samples.matrix, samples.num_rows, 00099 samples.num_cols); 00100 Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols); 00101 eigen_samples=eigen_L*eigen_samples; 00102 00103 /* add mean */ 00104 Map<VectorXd> eigen_mean(m_mean.vector, m_mean.vlen); 00105 eigen_samples.colwise()+=eigen_mean; 00106 00107 return samples; 00108 } 00109 00110 SGVector<float64_t> CGaussianDistribution::log_pdf_multiple(SGMatrix<float64_t> samples) const 00111 { 00112 REQUIRE(samples.num_cols>0, "Number of samples must be positive, but is %d\n", 00113 samples.num_cols); 00114 REQUIRE(samples.num_rows==m_dimension, "Sample dimension (%d) does not match" 00115 "Gaussian dimension (%d)\n", samples.num_rows, m_dimension); 00116 00117 /* for easier to read code */ 00118 index_t num_samples=samples.num_cols; 00119 00120 float64_t const_part=-0.5 * m_dimension * CMath::log(2 * CMath::PI); 00121 00122 /* determinant is product of diagonal elements of triangular matrix */ 00123 float64_t log_det_part=0; 00124 Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols); 00125 VectorXd diag=eigen_L.diagonal(); 00126 log_det_part=-diag.array().log().sum(); 00127 00128 /* sample based part */ 00129 Map<MatrixXd> eigen_samples(samples.matrix, samples.num_rows, 00130 samples.num_cols); 00131 Map<VectorXd> eigen_mean(m_mean.vector, m_mean.vlen); 00132 00133 /* substract mean from samples (create copy) */ 00134 SGMatrix<float64_t> centred(m_dimension, num_samples); 00135 Map<MatrixXd> eigen_centred(centred.matrix, centred.num_rows, 00136 centred.num_cols); 00137 for (index_t dim=0; dim<m_dimension; ++dim) 00138 { 00139 for (index_t sample_idx=0; sample_idx<num_samples; ++sample_idx) 00140 centred(dim,sample_idx)=samples(dim,sample_idx)-m_mean[dim]; 00141 } 00142 00143 /* solve the linear system based on factorization */ 00144 MatrixXd solved=eigen_L.triangularView<Lower>().solve(eigen_centred); 00145 solved=eigen_L.transpose().triangularView<Upper>().solve(solved); 00146 00147 /* one quadratic part x^T C^-1 x for each sample x */ 00148 SGVector<float64_t> result(num_samples); 00149 Map<VectorXd> eigen_result(result.vector, result.vlen); 00150 for (index_t i=0; i<num_samples; ++i) 00151 { 00152 /* i-th centred sample */ 00153 VectorXd left=eigen_centred.block(0, i, m_dimension, 1); 00154 00155 /* inverted covariance times i-th centred sample */ 00156 VectorXd right=solved.block(0,i,m_dimension,1); 00157 result[i]=-0.5*left.dot(right); 00158 } 00159 00160 /* combine and return */ 00161 eigen_result=eigen_result.array()+(log_det_part+const_part); 00162 00163 /* contains everything */ 00164 return result; 00165 00166 } 00167 00168 void CGaussianDistribution::init() 00169 { 00170 SG_ADD(&m_mean, "mean", "Mean of the Gaussian.", MS_NOT_AVAILABLE); 00171 SG_ADD(&m_L, "L", "Lower factor of covariance matrix, " 00172 "depending on the factorization type.", MS_NOT_AVAILABLE); 00173 } 00174 00175 #endif // HAVE_EIGEN3