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 Roman Votyakov 00008 */ 00009 00010 #include <shogun/lib/config.h> 00011 00012 #ifdef HAVE_EIGEN3 00013 00014 #include <shogun/machine/GaussianProcessMachine.h> 00015 #include <shogun/mathematics/Math.h> 00016 #include <shogun/kernel/Kernel.h> 00017 #include <shogun/machine/gp/FITCInferenceMethod.h> 00018 00019 #include <shogun/mathematics/eigen3.h> 00020 00021 using namespace shogun; 00022 using namespace Eigen; 00023 00024 CGaussianProcessMachine::CGaussianProcessMachine() 00025 { 00026 init(); 00027 } 00028 00029 CGaussianProcessMachine::CGaussianProcessMachine(CInferenceMethod* method) 00030 { 00031 init(); 00032 set_inference_method(method); 00033 } 00034 00035 void CGaussianProcessMachine::init() 00036 { 00037 m_method=NULL; 00038 00039 SG_ADD((CSGObject**) &m_method, "inference_method", "Inference method", 00040 MS_AVAILABLE); 00041 } 00042 00043 CGaussianProcessMachine::~CGaussianProcessMachine() 00044 { 00045 SG_UNREF(m_method); 00046 } 00047 00048 SGVector<float64_t> CGaussianProcessMachine::get_posterior_means(CFeatures* data) 00049 { 00050 REQUIRE(m_method, "Inference method should not be NULL\n") 00051 00052 CFeatures* feat; 00053 00054 // use latent features for FITC inference method 00055 if (m_method->get_inference_type()==INF_FITC) 00056 { 00057 CFITCInferenceMethod* fitc_method= 00058 CFITCInferenceMethod::obtain_from_generic(m_method); 00059 feat=fitc_method->get_latent_features(); 00060 SG_UNREF(fitc_method); 00061 } 00062 else 00063 feat=m_method->get_features(); 00064 00065 // get kernel and compute kernel matrix: K(feat, data)*scale^2 00066 CKernel* kernel=m_method->get_kernel(); 00067 kernel->init(feat, data); 00068 00069 // get kernel matrix and create eigen representation of it 00070 SGMatrix<float64_t> k_trts=kernel->get_kernel_matrix(); 00071 Map<MatrixXd> eigen_Ks(k_trts.matrix, k_trts.num_rows, k_trts.num_cols); 00072 00073 // compute Ks=Ks*scale^2 00074 eigen_Ks*=CMath::sq(m_method->get_scale()); 00075 00076 // cleanup 00077 SG_UNREF(feat); 00078 SG_UNREF(kernel); 00079 00080 // get alpha and create eigen representation of it 00081 SGVector<float64_t> alpha=m_method->get_alpha(); 00082 Map<VectorXd> eigen_alpha(alpha.vector, alpha.vlen); 00083 00084 // get mean and create eigen representation of it 00085 CMeanFunction* mean_function=m_method->get_mean(); 00086 SGVector<float64_t> m=mean_function->get_mean_vector(data); 00087 Map<VectorXd> eigen_m(m.vector, m.vlen); 00088 SG_UNREF(mean_function); 00089 00090 // compute mean: mu=Ks'*alpha+m 00091 SGVector<float64_t> mu(m.vlen); 00092 Map<VectorXd> eigen_mu(mu.vector, mu.vlen); 00093 eigen_mu=eigen_Ks.adjoint()*eigen_alpha+eigen_m; 00094 00095 return mu; 00096 } 00097 00098 SGVector<float64_t> CGaussianProcessMachine::get_posterior_variances( 00099 CFeatures* data) 00100 { 00101 REQUIRE(m_method, "Inference method should not be NULL\n") 00102 00103 CFeatures* feat; 00104 00105 // use latent features for FITC inference method 00106 if (m_method->get_inference_type()==INF_FITC) 00107 { 00108 CFITCInferenceMethod* fitc_method= 00109 CFITCInferenceMethod::obtain_from_generic(m_method); 00110 feat=fitc_method->get_latent_features(); 00111 SG_UNREF(fitc_method); 00112 } 00113 else 00114 feat=m_method->get_features(); 00115 00116 SG_REF(data); 00117 00118 // get kernel and compute kernel matrix: K(data, data)*scale^2 00119 CKernel* kernel=m_method->get_kernel(); 00120 kernel->init(data, data); 00121 00122 // get kernel matrix and create eigen representation of it 00123 SGMatrix<float64_t> k_tsts=kernel->get_kernel_matrix(); 00124 Map<MatrixXd> eigen_Kss(k_tsts.matrix, k_tsts.num_rows, k_tsts.num_cols); 00125 00126 // compute Kss=Kss*scale^2 00127 eigen_Kss*=CMath::sq(m_method->get_scale()); 00128 00129 // compute kernel matrix: K(feat, data)*scale^2 00130 kernel->init(feat, data); 00131 00132 // get kernel matrix and create eigen representation of it 00133 SGMatrix<float64_t> k_trts=kernel->get_kernel_matrix(); 00134 Map<MatrixXd> eigen_Ks(k_trts.matrix, k_trts.num_rows, k_trts.num_cols); 00135 00136 // compute Ks=Ks*scale^2 00137 eigen_Ks*=CMath::sq(m_method->get_scale()); 00138 00139 // cleanup 00140 SG_UNREF(kernel); 00141 SG_UNREF(feat); 00142 SG_UNREF(data); 00143 00144 // get shogun representation of cholesky and create eigen representation 00145 SGMatrix<float64_t> L=m_method->get_cholesky(); 00146 Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols); 00147 00148 // result variance vector 00149 SGVector<float64_t> s2(k_tsts.num_cols); 00150 Map<VectorXd> eigen_s2(s2.vector, s2.vlen); 00151 00152 if (eigen_L.isUpperTriangular()) 00153 { 00154 // get shogun of diagonal sigma vector and create eigen representation 00155 SGVector<float64_t> sW=m_method->get_diagonal_vector(); 00156 Map<VectorXd> eigen_sW(sW.vector, sW.vlen); 00157 00158 // solve L' * V = sW * Ks and compute V.^2 00159 MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve( 00160 eigen_sW.asDiagonal()*eigen_Ks); 00161 MatrixXd eigen_sV=eigen_V.cwiseProduct(eigen_V); 00162 00163 eigen_s2=eigen_Kss.diagonal()-eigen_sV.colwise().sum().adjoint(); 00164 } 00165 else 00166 { 00167 // M = Ks .* (L * Ks) 00168 MatrixXd eigen_M=eigen_Ks.cwiseProduct(eigen_L*eigen_Ks); 00169 eigen_s2=eigen_Kss.diagonal()+eigen_M.colwise().sum().adjoint(); 00170 } 00171 00172 return s2; 00173 } 00174 00175 #endif /* HAVE_EIGEN3 */