SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
GaussianProcessMachine.cpp
Go to the documentation of this file.
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 */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation