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 * 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 */