SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
LaplacianInferenceMethod.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  * Copyright (C) 2012 Jacob Walker
00009  * Copyright (C) 2013 Roman Votyakov
00010  *
00011  * Code adapted from Gaussian Process Machine Learning Toolbox
00012  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
00013  * This code specifically adapted from infLaplace.m
00014  */
00015 
00016 #include <shogun/machine/gp/LaplacianInferenceMethod.h>
00017 
00018 #ifdef HAVE_EIGEN3
00019 
00020 #include <shogun/machine/gp/StudentsTLikelihood.h>
00021 #include <shogun/mathematics/Math.h>
00022 #include <shogun/lib/external/brent.h>
00023 #include <shogun/mathematics/eigen3.h>
00024 
00025 using namespace shogun;
00026 using namespace Eigen;
00027 
00028 namespace shogun
00029 {
00030 
00031 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00032 
00034 class CPsiLine : public func_base
00035 {
00036 public:
00037     float64_t scale;
00038     MatrixXd K;
00039     VectorXd dalpha;
00040     VectorXd start_alpha;
00041     Map<VectorXd>* alpha;
00042     SGVector<float64_t>* dlp;
00043     SGVector<float64_t>* W;
00044     SGVector<float64_t>* f;
00045     SGVector<float64_t>* m;
00046     CLikelihoodModel* lik;
00047     CLabels* lab;
00048 
00049     virtual double operator() (double x)
00050     {
00051         Map<VectorXd> eigen_f(f->vector, f->vlen);
00052         Map<VectorXd> eigen_m(m->vector, m->vlen);
00053 
00054         // compute alpha=alpha+x*dalpha and f=K*alpha+m
00055         (*alpha)=start_alpha+x*dalpha;
00056         eigen_f=K*(*alpha)*CMath::sq(scale)+eigen_m;
00057 
00058         // get first and second derivatives of log likelihood
00059         (*dlp)=lik->get_log_probability_derivative_f(lab, (*f), 1);
00060 
00061         (*W)=lik->get_log_probability_derivative_f(lab, (*f), 2);
00062         W->scale(-1.0);
00063 
00064         // compute psi=alpha'*(f-m)/2-lp
00065         float64_t result = (*alpha).dot(eigen_f-eigen_m)/2.0-
00066             SGVector<float64_t>::sum(lik->get_log_probability_f(lab, *f));
00067 
00068         return result;
00069     }
00070 };
00071 
00072 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
00073 
00074 CLaplacianInferenceMethod::CLaplacianInferenceMethod() : CInferenceMethod()
00075 {
00076     init();
00077 }
00078 
00079 CLaplacianInferenceMethod::CLaplacianInferenceMethod(CKernel* kern,
00080         CFeatures* feat, CMeanFunction* m, CLabels* lab, CLikelihoodModel* mod)
00081         : CInferenceMethod(kern, feat, m, lab, mod)
00082 {
00083     init();
00084 }
00085 
00086 void CLaplacianInferenceMethod::init()
00087 {
00088     m_iter=20;
00089     m_tolerance=1e-6;
00090     m_opt_tolerance=1e-10;
00091     m_opt_max=10;
00092 }
00093 
00094 CLaplacianInferenceMethod::~CLaplacianInferenceMethod()
00095 {
00096 }
00097 
00098 void CLaplacianInferenceMethod::update()
00099 {
00100     CInferenceMethod::update();
00101     update_alpha();
00102     update_chol();
00103     update_approx_cov();
00104     update_deriv();
00105 }
00106 
00107 SGVector<float64_t> CLaplacianInferenceMethod::get_diagonal_vector()
00108 {
00109     if (update_parameter_hash())
00110         update();
00111 
00112     return SGVector<float64_t>(sW);
00113 }
00114 
00115 float64_t CLaplacianInferenceMethod::get_negative_log_marginal_likelihood()
00116 {
00117     if (update_parameter_hash())
00118         update();
00119 
00120     // create eigen representations alpha, f, W, L
00121     Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
00122     Map<VectorXd> eigen_mu(m_mu.vector, m_mu.vlen);
00123     Map<VectorXd> eigen_W(W.vector, W.vlen);
00124     Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
00125 
00126     // get mean vector and create eigen representation of it
00127     SGVector<float64_t> mean=m_mean->get_mean_vector(m_features);
00128     Map<VectorXd> eigen_mean(mean.vector, mean.vlen);
00129 
00130     // get log likelihood
00131     float64_t lp=SGVector<float64_t>::sum(m_model->get_log_probability_f(m_labels,
00132         m_mu));
00133 
00134     float64_t result;
00135 
00136     if (eigen_W.minCoeff()<0)
00137     {
00138         Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
00139         Map<MatrixXd> eigen_ktrtr(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
00140 
00141         FullPivLU<MatrixXd> lu(MatrixXd::Identity(m_ktrtr.num_rows, m_ktrtr.num_cols)+
00142             eigen_ktrtr*CMath::sq(m_scale)*eigen_sW.asDiagonal());
00143 
00144         result=(eigen_alpha.dot(eigen_mu-eigen_mean))/2.0-
00145             lp+log(lu.determinant())/2.0;
00146     }
00147     else
00148     {
00149         result=eigen_alpha.dot(eigen_mu-eigen_mean)/2.0-lp+
00150             eigen_L.diagonal().array().log().sum();
00151     }
00152 
00153     return result;
00154 }
00155 
00156 SGVector<float64_t> CLaplacianInferenceMethod::get_alpha()
00157 {
00158     if (update_parameter_hash())
00159         update();
00160 
00161     SGVector<float64_t> result(m_alpha);
00162     return result;
00163 }
00164 
00165 SGMatrix<float64_t> CLaplacianInferenceMethod::get_cholesky()
00166 {
00167     if (update_parameter_hash())
00168         update();
00169 
00170     return SGMatrix<float64_t>(m_L);
00171 }
00172 
00173 SGVector<float64_t> CLaplacianInferenceMethod::get_posterior_mean()
00174 {
00175     if (update_parameter_hash())
00176         update();
00177 
00178     return SGVector<float64_t>(m_mu);
00179 }
00180 
00181 SGMatrix<float64_t> CLaplacianInferenceMethod::get_posterior_covariance()
00182 {
00183     if (update_parameter_hash())
00184         update();
00185 
00186     return SGMatrix<float64_t>(m_Sigma);
00187 }
00188 
00189 void CLaplacianInferenceMethod::update_approx_cov()
00190 {
00191     Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
00192     Map<MatrixXd> eigen_K(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
00193     Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
00194 
00195     m_Sigma=SGMatrix<float64_t>(m_ktrtr.num_rows, m_ktrtr.num_cols);
00196     Map<MatrixXd> eigen_Sigma(m_Sigma.matrix, m_Sigma.num_rows, m_Sigma.num_cols);
00197 
00198     // compute V = L^(-1) * W^(1/2) * K, using upper triangular factor L^T
00199     MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
00200             eigen_sW.asDiagonal()*eigen_K*CMath::sq(m_scale));
00201 
00202     // compute covariance matrix of the posterior:
00203     // Sigma = K - K * W^(1/2) * (L * L^T)^(-1) * W^(1/2) * K =
00204     // K - (K * W^(1/2)) * (L^T)^(-1) * L^(-1) * W^(1/2) * K =
00205     // K - (W^(1/2) * K)^T * (L^(-1))^T * L^(-1) * W^(1/2) * K = K - V^T * V
00206     eigen_Sigma=eigen_K*CMath::sq(m_scale)-eigen_V.adjoint()*eigen_V;
00207 }
00208 
00209 void CLaplacianInferenceMethod::update_chol()
00210 {
00211     Map<VectorXd> eigen_W(W.vector, W.vlen);
00212 
00213     // create eigen representation of kernel matrix
00214     Map<MatrixXd> eigen_ktrtr(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
00215 
00216     // create shogun and eigen representation of posterior cholesky
00217     m_L=SGMatrix<float64_t>(m_ktrtr.num_rows, m_ktrtr.num_cols);
00218     Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
00219 
00220     if (eigen_W.minCoeff() < 0)
00221     {
00222         // compute inverse of diagonal noise: iW = 1/W
00223         VectorXd eigen_iW = (VectorXd::Ones(W.vlen)).cwiseQuotient(eigen_W);
00224 
00225         FullPivLU<MatrixXd> lu(
00226             eigen_ktrtr*CMath::sq(m_scale)+MatrixXd(eigen_iW.asDiagonal()));
00227 
00228         // compute cholesky: L = -(K + iW)^-1
00229         eigen_L = -lu.inverse();
00230     }
00231     else
00232     {
00233         Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
00234 
00235         // compute cholesky: L = chol(sW * sW' .* K + I)
00236         LLT<MatrixXd> L(
00237             (eigen_sW*eigen_sW.transpose()).cwiseProduct(eigen_ktrtr*CMath::sq(m_scale))+
00238             MatrixXd::Identity(m_ktrtr.num_rows, m_ktrtr.num_cols));
00239 
00240         eigen_L = L.matrixU();
00241     }
00242 }
00243 
00244 void CLaplacianInferenceMethod::update_alpha()
00245 {
00246     float64_t Psi_Old = CMath::INFTY;
00247     float64_t Psi_New;
00248     float64_t Psi_Def;
00249 
00250     // get mean vector and create eigen representation of it
00251     SGVector<float64_t> mean=m_mean->get_mean_vector(m_features);
00252     Map<VectorXd> eigen_mean(mean.vector, mean.vlen);
00253 
00254     // create eigen representation of kernel matrix
00255     Map<MatrixXd> eigen_ktrtr(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
00256 
00257     // create shogun and eigen representation of function vector
00258     m_mu=SGVector<float64_t>(mean.vlen);
00259     Map<VectorXd> eigen_mu(m_mu, m_mu.vlen);
00260 
00261     if (m_alpha.vlen!=m_labels->get_num_labels())
00262     {
00263         // set alpha a zero vector
00264         m_alpha=SGVector<float64_t>(m_labels->get_num_labels());
00265         m_alpha.zero();
00266 
00267         // f = mean, if length of alpha and length of y doesn't match
00268         eigen_mu=eigen_mean;
00269 
00270         // compute W = -d2lp
00271         W=m_model->get_log_probability_derivative_f(m_labels, m_mu, 2);
00272         W.scale(-1.0);
00273 
00274         Psi_New=-SGVector<float64_t>::sum(m_model->get_log_probability_f(
00275             m_labels, m_mu));
00276     }
00277     else
00278     {
00279         Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
00280 
00281         // compute f = K * alpha + m
00282         eigen_mu=eigen_ktrtr*CMath::sq(m_scale)*eigen_alpha+eigen_mean;
00283 
00284         // compute W = -d2lp
00285         W=m_model->get_log_probability_derivative_f(m_labels, m_mu, 2);
00286         W.scale(-1.0);
00287 
00288         Psi_New=eigen_alpha.dot(eigen_mu-eigen_mean)/2.0-
00289             SGVector<float64_t>::sum(m_model->get_log_probability_f(m_labels, m_mu));
00290 
00291         Psi_Def=-SGVector<float64_t>::sum(m_model->get_log_probability_f(m_labels, mean));
00292 
00293         // if default is better, then use it
00294         if (Psi_Def < Psi_New)
00295         {
00296             m_alpha.zero();
00297             eigen_mu=eigen_mean;
00298             Psi_New=-SGVector<float64_t>::sum(m_model->get_log_probability_f(
00299                 m_labels, m_mu));
00300         }
00301     }
00302 
00303     Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
00304 
00305     // get first derivative of log probability function
00306     dlp=m_model->get_log_probability_derivative_f(m_labels, m_mu, 1);
00307 
00308     // create shogun and eigen representation of sW
00309     sW=SGVector<float64_t>(W.vlen);
00310     Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
00311 
00312     index_t iter=0;
00313 
00314     while (Psi_Old-Psi_New>m_tolerance && iter<m_iter)
00315     {
00316         Map<VectorXd> eigen_W(W.vector, W.vlen);
00317         Map<VectorXd> eigen_dlp(dlp.vector, dlp.vlen);
00318 
00319         Psi_Old = Psi_New;
00320         iter++;
00321 
00322         if (eigen_W.minCoeff() < 0)
00323         {
00324             // Suggested by Vanhatalo et. al.,
00325             // Gaussian Process Regression with Student's t likelihood, NIPS 2009
00326             // Quoted from infLaplace.m
00327             float64_t df;
00328 
00329             if (m_model->get_model_type()==LT_STUDENTST)
00330             {
00331                 CStudentsTLikelihood* lik=CStudentsTLikelihood::obtain_from_generic(m_model);
00332                 df=lik->get_degrees_freedom();
00333                 SG_UNREF(lik);
00334             }
00335             else
00336                 df=1;
00337 
00338             eigen_W+=(2.0/df)*eigen_dlp.cwiseProduct(eigen_dlp);
00339         }
00340 
00341         // compute sW = sqrt(W)
00342         eigen_sW=eigen_W.cwiseSqrt();
00343 
00344         LLT<MatrixXd> L((eigen_sW*eigen_sW.transpose()).cwiseProduct(eigen_ktrtr*CMath::sq(m_scale))+
00345             MatrixXd::Identity(m_ktrtr.num_rows, m_ktrtr.num_cols));
00346 
00347         VectorXd b=eigen_W.cwiseProduct(eigen_mu - eigen_mean)+eigen_dlp;
00348 
00349         VectorXd dalpha=b-eigen_sW.cwiseProduct(
00350             L.solve(eigen_sW.cwiseProduct(eigen_ktrtr*b*CMath::sq(m_scale))))-eigen_alpha;
00351 
00352         // perform Brent's optimization
00353         CPsiLine func;
00354 
00355         func.scale=m_scale;
00356         func.K=eigen_ktrtr;
00357         func.dalpha=dalpha;
00358         func.start_alpha=eigen_alpha;
00359         func.alpha=&eigen_alpha;
00360         func.dlp=&dlp;
00361         func.f=&m_mu;
00362         func.m=&mean;
00363         func.W=&W;
00364         func.lik=m_model;
00365         func.lab=m_labels;
00366 
00367         float64_t x;
00368         Psi_New=local_min(0, m_opt_max, m_opt_tolerance, func, x);
00369     }
00370 
00371     // compute f = K * alpha + m
00372     eigen_mu=eigen_ktrtr*CMath::sq(m_scale)*eigen_alpha+eigen_mean;
00373 
00374     // get log probability derivatives
00375     dlp=m_model->get_log_probability_derivative_f(m_labels, m_mu, 1);
00376     d2lp=m_model->get_log_probability_derivative_f(m_labels, m_mu, 2);
00377     d3lp=m_model->get_log_probability_derivative_f(m_labels, m_mu, 3);
00378 
00379     // W = -d2lp
00380     W=d2lp.clone();
00381     W.scale(-1.0);
00382 
00383     // compute sW
00384     Map<VectorXd> eigen_W(W.vector, W.vlen);
00385 
00386     if (eigen_W.minCoeff()>0)
00387         eigen_sW=eigen_W.cwiseSqrt();
00388     else
00389         eigen_sW.setZero();
00390 }
00391 
00392 void CLaplacianInferenceMethod::update_deriv()
00393 {
00394     // create eigen representation of W, sW, dlp, d3lp, K, alpha and L
00395     Map<VectorXd> eigen_W(W.vector, W.vlen);
00396     Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
00397     Map<VectorXd> eigen_dlp(dlp.vector, dlp.vlen);
00398     Map<VectorXd> eigen_d3lp(d3lp.vector, d3lp.vlen);
00399     Map<MatrixXd> eigen_K(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
00400     Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
00401     Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
00402 
00403     // create shogun and eigen representation of matrix Z
00404     m_Z=SGMatrix<float64_t>(m_L.num_rows, m_L.num_cols);
00405     Map<MatrixXd> eigen_Z(m_Z.matrix, m_Z.num_rows, m_Z.num_cols);
00406 
00407     // create shogun and eigen representation of the vector g
00408     m_g=SGVector<float64_t>(m_Z.num_rows);
00409     Map<VectorXd> eigen_g(m_g.vector, m_g.vlen);
00410 
00411     if (eigen_W.minCoeff()<0)
00412     {
00413         eigen_Z=-eigen_L;
00414 
00415         // compute iA = (I + K * diag(W))^-1
00416         FullPivLU<MatrixXd> lu(MatrixXd::Identity(m_ktrtr.num_rows, m_ktrtr.num_cols)+
00417                 eigen_K*CMath::sq(m_scale)*eigen_W.asDiagonal());
00418         MatrixXd iA=lu.inverse();
00419 
00420         // compute derivative ln|L'*L| wrt W: g=sum(iA.*K,2)/2
00421         eigen_g=(iA.cwiseProduct(eigen_K*CMath::sq(m_scale))).rowwise().sum()/2.0;
00422     }
00423     else
00424     {
00425         // solve L'*L*Z=diag(sW) and compute Z=diag(sW)*Z
00426         eigen_Z=eigen_L.triangularView<Upper>().adjoint().solve(
00427                 MatrixXd(eigen_sW.asDiagonal()));
00428         eigen_Z=eigen_L.triangularView<Upper>().solve(eigen_Z);
00429         eigen_Z=eigen_sW.asDiagonal()*eigen_Z;
00430 
00431         // solve L'*C=diag(sW)*K
00432         MatrixXd C=eigen_L.triangularView<Upper>().adjoint().solve(
00433                 eigen_sW.asDiagonal()*eigen_K*CMath::sq(m_scale));
00434 
00435         // compute derivative ln|L'*L| wrt W: g=(diag(K)-sum(C.^2,1)')/2
00436         eigen_g=(eigen_K.diagonal()*CMath::sq(m_scale)-
00437                 (C.cwiseProduct(C)).colwise().sum().adjoint())/2.0;
00438     }
00439 
00440     // create shogun and eigen representation of the vector dfhat
00441     m_dfhat=SGVector<float64_t>(m_g.vlen);
00442     Map<VectorXd> eigen_dfhat(m_dfhat.vector, m_dfhat.vlen);
00443 
00444     // compute derivative of nlZ wrt fhat
00445     eigen_dfhat=eigen_g.cwiseProduct(eigen_d3lp);
00446 }
00447 
00448 SGVector<float64_t> CLaplacianInferenceMethod::get_derivative_wrt_inference_method(
00449         const TParameter* param)
00450 {
00451     REQUIRE(!strcmp(param->m_name, "scale"), "Can't compute derivative of "
00452             "the nagative log marginal likelihood wrt %s.%s parameter\n",
00453             get_name(), param->m_name)
00454 
00455     // create eigen representation of K, Z, dfhat, dlp and alpha
00456     Map<MatrixXd> eigen_K(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
00457     Map<MatrixXd> eigen_Z(m_Z.matrix, m_Z.num_rows, m_Z.num_cols);
00458     Map<VectorXd> eigen_dfhat(m_dfhat.vector, m_dfhat.vlen);
00459     Map<VectorXd> eigen_dlp(dlp.vector, dlp.vlen);
00460     Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
00461 
00462     SGVector<float64_t> result(1);
00463 
00464     // compute derivative K wrt scale
00465     MatrixXd dK=eigen_K*m_scale*2.0;
00466 
00467     // compute dnlZ=sum(sum(Z.*dK))/2-alpha'*dK*alpha/2
00468     result[0]=(eigen_Z.cwiseProduct(dK)).sum()/2.0-
00469         (eigen_alpha.adjoint()*dK).dot(eigen_alpha)/2.0;
00470 
00471     // compute b=dK*dlp
00472     VectorXd b=dK*eigen_dlp;
00473 
00474     // compute dnlZ=dnlZ-dfhat'*(b-K*(Z*b))
00475     result[0]=result[0]-eigen_dfhat.dot(b-eigen_K*CMath::sq(m_scale)*(eigen_Z*b));
00476 
00477     return result;
00478 }
00479 
00480 SGVector<float64_t> CLaplacianInferenceMethod::get_derivative_wrt_likelihood_model(
00481         const TParameter* param)
00482 {
00483     // create eigen representation of K, Z, g and dfhat
00484     Map<MatrixXd> eigen_K(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
00485     Map<MatrixXd> eigen_Z(m_Z.matrix, m_Z.num_rows, m_Z.num_cols);
00486     Map<VectorXd> eigen_g(m_g.vector, m_g.vlen);
00487     Map<VectorXd> eigen_dfhat(m_dfhat.vector, m_dfhat.vlen);
00488 
00489     // get derivatives wrt likelihood model parameters
00490     SGVector<float64_t> lp_dhyp=m_model->get_first_derivative(m_labels,
00491             m_mu, param);
00492     SGVector<float64_t> dlp_dhyp=m_model->get_second_derivative(m_labels,
00493             m_mu, param);
00494     SGVector<float64_t> d2lp_dhyp=m_model->get_third_derivative(m_labels,
00495             m_mu, param);
00496 
00497     // create eigen representation of the derivatives
00498     Map<VectorXd> eigen_lp_dhyp(lp_dhyp.vector, lp_dhyp.vlen);
00499     Map<VectorXd> eigen_dlp_dhyp(dlp_dhyp.vector, dlp_dhyp.vlen);
00500     Map<VectorXd> eigen_d2lp_dhyp(d2lp_dhyp.vector, d2lp_dhyp.vlen);
00501 
00502     SGVector<float64_t> result(1);
00503 
00504     // compute b vector
00505     VectorXd b=eigen_K*eigen_dlp_dhyp;
00506 
00507     // compute dnlZ=-g'*d2lp_dhyp-sum(lp_dhyp)-dfhat'*(b-K*(Z*b))
00508     result[0]=-eigen_g.dot(eigen_d2lp_dhyp)-eigen_lp_dhyp.sum()-
00509         eigen_dfhat.dot(b-eigen_K*CMath::sq(m_scale)*(eigen_Z*b));
00510 
00511     return result;
00512 }
00513 
00514 SGVector<float64_t> CLaplacianInferenceMethod::get_derivative_wrt_kernel(
00515         const TParameter* param)
00516 {
00517     // create eigen representation of K, Z, dfhat, dlp and alpha
00518     Map<MatrixXd> eigen_K(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
00519     Map<MatrixXd> eigen_Z(m_Z.matrix, m_Z.num_rows, m_Z.num_cols);
00520     Map<VectorXd> eigen_dfhat(m_dfhat.vector, m_dfhat.vlen);
00521     Map<VectorXd> eigen_dlp(dlp.vector, dlp.vlen);
00522     Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
00523 
00524     SGVector<float64_t> result;
00525 
00526     if (param->m_datatype.m_ctype==CT_VECTOR ||
00527             param->m_datatype.m_ctype==CT_SGVECTOR)
00528     {
00529         REQUIRE(param->m_datatype.m_length_y,
00530                 "Length of the parameter %s should not be NULL\n", param->m_name)
00531         result=SGVector<float64_t>(*(param->m_datatype.m_length_y));
00532     }
00533     else
00534     {
00535         result=SGVector<float64_t>(1);
00536     }
00537 
00538     for (index_t i=0; i<result.vlen; i++)
00539     {
00540         SGMatrix<float64_t> dK;
00541 
00542         if (result.vlen==1)
00543             dK=m_kernel->get_parameter_gradient(param);
00544         else
00545             dK=m_kernel->get_parameter_gradient(param, i);
00546 
00547         Map<MatrixXd> eigen_dK(dK.matrix, dK.num_cols, dK.num_rows);
00548 
00549         // compute dnlZ=sum(sum(Z.*dK))/2-alpha'*dK*alpha/2
00550         result[i]=(eigen_Z.cwiseProduct(eigen_dK)).sum()/2.0-
00551             (eigen_alpha.adjoint()*eigen_dK).dot(eigen_alpha)/2.0;
00552 
00553         // compute b=dK*dlp
00554         VectorXd b=eigen_dK*eigen_dlp;
00555 
00556         // compute dnlZ=dnlZ-dfhat'*(b-K*(Z*b))
00557         result[i]=result[i]-eigen_dfhat.dot(b-eigen_K*CMath::sq(m_scale)*
00558                 (eigen_Z*b));
00559     }
00560 
00561     return result;
00562 }
00563 
00564 SGVector<float64_t> CLaplacianInferenceMethod::get_derivative_wrt_mean(
00565         const TParameter* param)
00566 {
00567     // create eigen representation of K, Z, dfhat and alpha
00568     Map<MatrixXd> eigen_K(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
00569     Map<MatrixXd> eigen_Z(m_Z.matrix, m_Z.num_rows, m_Z.num_cols);
00570     Map<VectorXd> eigen_dfhat(m_dfhat.vector, m_dfhat.vlen);
00571     Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
00572 
00573     SGVector<float64_t> result;
00574 
00575     if (param->m_datatype.m_ctype==CT_VECTOR ||
00576             param->m_datatype.m_ctype==CT_SGVECTOR)
00577     {
00578         REQUIRE(param->m_datatype.m_length_y,
00579                 "Length of the parameter %s should not be NULL\n", param->m_name)
00580 
00581         result=SGVector<float64_t>(*(param->m_datatype.m_length_y));
00582     }
00583     else
00584     {
00585         result=SGVector<float64_t>(1);
00586     }
00587 
00588     for (index_t i=0; i<result.vlen; i++)
00589     {
00590         SGVector<float64_t> dmu;
00591 
00592         if (result.vlen==1)
00593             dmu=m_mean->get_parameter_derivative(m_features, param);
00594         else
00595             dmu=m_mean->get_parameter_derivative(m_features, param, i);
00596 
00597         Map<VectorXd> eigen_dmu(dmu.vector, dmu.vlen);
00598 
00599         // compute dnlZ=-alpha'*dm-dfhat'*(dm-K*(Z*dm))
00600         result[i]=-eigen_alpha.dot(eigen_dmu)-eigen_dfhat.dot(eigen_dmu-
00601                 eigen_K*CMath::sq(m_scale)*(eigen_Z*eigen_dmu));
00602     }
00603 
00604     return result;
00605 }
00606 }
00607 
00608 #endif /* HAVE_EIGEN3 */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation