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 */ 00014 00015 #include <shogun/machine/gp/FITCInferenceMethod.h> 00016 00017 #ifdef HAVE_EIGEN3 00018 00019 #include <shogun/machine/gp/GaussianLikelihood.h> 00020 #include <shogun/mathematics/Math.h> 00021 #include <shogun/labels/RegressionLabels.h> 00022 #include <shogun/mathematics/eigen3.h> 00023 00024 using namespace shogun; 00025 using namespace Eigen; 00026 00027 CFITCInferenceMethod::CFITCInferenceMethod() : CInferenceMethod() 00028 { 00029 init(); 00030 } 00031 00032 CFITCInferenceMethod::CFITCInferenceMethod(CKernel* kern, CFeatures* feat, 00033 CMeanFunction* m, CLabels* lab, CLikelihoodModel* mod, CFeatures* lat) 00034 : CInferenceMethod(kern, feat, m, lab, mod) 00035 { 00036 init(); 00037 set_latent_features(lat); 00038 } 00039 00040 void CFITCInferenceMethod::init() 00041 { 00042 SG_ADD((CSGObject**)&m_latent_features, "latent_features", "Latent features", 00043 MS_NOT_AVAILABLE); 00044 00045 m_latent_features=NULL; 00046 m_ind_noise=1e-10; 00047 } 00048 00049 CFITCInferenceMethod::~CFITCInferenceMethod() 00050 { 00051 SG_UNREF(m_latent_features); 00052 } 00053 00054 CFITCInferenceMethod* CFITCInferenceMethod::obtain_from_generic( 00055 CInferenceMethod* inference) 00056 { 00057 ASSERT(inference!=NULL); 00058 00059 if (inference->get_inference_type()!=INF_FITC) 00060 SG_SERROR("Provided inference is not of type CFITCInferenceMethod!\n") 00061 00062 SG_REF(inference); 00063 return (CFITCInferenceMethod*)inference; 00064 } 00065 00066 void CFITCInferenceMethod::update() 00067 { 00068 CInferenceMethod::update(); 00069 update_chol(); 00070 update_alpha(); 00071 update_deriv(); 00072 } 00073 00074 void CFITCInferenceMethod::check_members() const 00075 { 00076 CInferenceMethod::check_members(); 00077 00078 REQUIRE(m_model->get_model_type()==LT_GAUSSIAN, 00079 "FITC inference method can only use Gaussian likelihood function\n") 00080 REQUIRE(m_labels->get_label_type()==LT_REGRESSION, "Labels must be type " 00081 "of CRegressionLabels\n") 00082 REQUIRE(m_latent_features, "Latent features should not be NULL\n") 00083 REQUIRE(m_latent_features->get_num_vectors(), 00084 "Number of latent features must be greater than zero\n") 00085 } 00086 00087 SGVector<float64_t> CFITCInferenceMethod::get_diagonal_vector() 00088 { 00089 if (update_parameter_hash()) 00090 update(); 00091 00092 // get the sigma variable from the Gaussian likelihood model 00093 CGaussianLikelihood* lik=CGaussianLikelihood::obtain_from_generic(m_model); 00094 float64_t sigma=lik->get_sigma(); 00095 SG_UNREF(lik); 00096 00097 // compute diagonal vector: sW=1/sigma 00098 SGVector<float64_t> result(m_features->get_num_vectors()); 00099 result.fill_vector(result.vector, m_features->get_num_vectors(), 1.0/sigma); 00100 00101 return result; 00102 } 00103 00104 float64_t CFITCInferenceMethod::get_negative_log_marginal_likelihood() 00105 { 00106 if (update_parameter_hash()) 00107 update(); 00108 00109 // create eigen representations of chol_utr, dg, r, be 00110 Map<MatrixXd> eigen_chol_utr(m_chol_utr.matrix, m_chol_utr.num_rows, 00111 m_chol_utr.num_cols); 00112 Map<VectorXd> eigen_dg(m_dg.vector, m_dg.vlen); 00113 Map<VectorXd> eigen_r(m_r.vector, m_r.vlen); 00114 Map<VectorXd> eigen_be(m_be.vector, m_be.vlen); 00115 00116 // compute negative log marginal likelihood: 00117 // nlZ=sum(log(diag(utr)))+(sum(log(dg))+r'*r-be'*be+n*log(2*pi))/2 00118 float64_t result=eigen_chol_utr.diagonal().array().log().sum()+ 00119 (eigen_dg.array().log().sum()+eigen_r.dot(eigen_r)-eigen_be.dot(eigen_be)+ 00120 m_ktrtr.num_rows*CMath::log(2*CMath::PI))/2.0; 00121 00122 return result; 00123 } 00124 00125 SGVector<float64_t> CFITCInferenceMethod::get_alpha() 00126 { 00127 if (update_parameter_hash()) 00128 update(); 00129 00130 SGVector<float64_t> result(m_alpha); 00131 return result; 00132 } 00133 00134 SGMatrix<float64_t> CFITCInferenceMethod::get_cholesky() 00135 { 00136 if (update_parameter_hash()) 00137 update(); 00138 00139 SGMatrix<float64_t> result(m_L); 00140 return result; 00141 } 00142 00143 SGVector<float64_t> CFITCInferenceMethod::get_posterior_mean() 00144 { 00145 SG_NOTIMPLEMENTED 00146 return SGVector<float64_t>(); 00147 } 00148 00149 SGMatrix<float64_t> CFITCInferenceMethod::get_posterior_covariance() 00150 { 00151 SG_NOTIMPLEMENTED 00152 return SGMatrix<float64_t>(); 00153 } 00154 00155 void CFITCInferenceMethod::update_train_kernel() 00156 { 00157 CInferenceMethod::update_train_kernel(); 00158 00159 // create kernel matrix for latent features 00160 m_kernel->init(m_latent_features, m_latent_features); 00161 m_kuu=m_kernel->get_kernel_matrix(); 00162 00163 // create kernel matrix for latent and training features 00164 m_kernel->init(m_latent_features, m_features); 00165 m_ktru=m_kernel->get_kernel_matrix(); 00166 } 00167 00168 void CFITCInferenceMethod::update_chol() 00169 { 00170 // get the sigma variable from the Gaussian likelihood model 00171 CGaussianLikelihood* lik=CGaussianLikelihood::obtain_from_generic(m_model); 00172 float64_t sigma=lik->get_sigma(); 00173 SG_UNREF(lik); 00174 00175 // eigen3 representation of covariance matrix of latent features (m_kuu) 00176 // and training features (m_ktru) 00177 Map<MatrixXd> eigen_kuu(m_kuu.matrix, m_kuu.num_rows, m_kuu.num_cols); 00178 Map<MatrixXd> eigen_ktru(m_ktru.matrix, m_ktru.num_rows, m_ktru.num_cols); 00179 Map<MatrixXd> eigen_ktrtr(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols); 00180 00181 // solve Luu' * Luu = Kuu + m_ind_noise * I 00182 LLT<MatrixXd> Luu(eigen_kuu*CMath::sq(m_scale)+m_ind_noise*MatrixXd::Identity( 00183 m_kuu.num_rows, m_kuu.num_cols)); 00184 00185 // create shogun and eigen3 representation of cholesky of covariance of 00186 // latent features Luu (m_chol_uu and eigen_chol_uu) 00187 m_chol_uu=SGMatrix<float64_t>(Luu.rows(), Luu.cols()); 00188 Map<MatrixXd> eigen_chol_uu(m_chol_uu.matrix, m_chol_uu.num_rows, 00189 m_chol_uu.num_cols); 00190 eigen_chol_uu=Luu.matrixU(); 00191 00192 // solve Luu' * V = Ktru 00193 MatrixXd V=eigen_chol_uu.triangularView<Upper>().adjoint().solve(eigen_ktru* 00194 CMath::sq(m_scale)); 00195 00196 // create shogun and eigen3 representation of 00197 // dg = diag(K) + sn2 - diag(Q) 00198 m_dg=SGVector<float64_t>(m_ktrtr.num_cols); 00199 Map<VectorXd> eigen_dg(m_dg.vector, m_dg.vlen); 00200 00201 eigen_dg=eigen_ktrtr.diagonal()*CMath::sq(m_scale)+CMath::sq(sigma)* 00202 VectorXd::Ones(m_dg.vlen)-(V.cwiseProduct(V)).colwise().sum().adjoint(); 00203 00204 // solve Lu' * Lu = V * diag(1/dg) * V' + I 00205 LLT<MatrixXd> Lu(V*((VectorXd::Ones(m_dg.vlen)).cwiseQuotient(eigen_dg)).asDiagonal()* 00206 V.adjoint()+MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)); 00207 00208 // create shogun and eigen3 representation of cholesky of covariance of 00209 // training features Luu (m_chol_utr and eigen_chol_utr) 00210 m_chol_utr=SGMatrix<float64_t>(Lu.rows(), Lu.cols()); 00211 Map<MatrixXd> eigen_chol_utr(m_chol_utr.matrix, m_chol_utr.num_rows, 00212 m_chol_utr.num_cols); 00213 eigen_chol_utr=Lu.matrixU(); 00214 00215 // create eigen representation of labels and mean vectors 00216 SGVector<float64_t> y=((CRegressionLabels*) m_labels)->get_labels(); 00217 Map<VectorXd> eigen_y(y.vector, y.vlen); 00218 SGVector<float64_t> m=m_mean->get_mean_vector(m_features); 00219 Map<VectorXd> eigen_m(m.vector, m.vlen); 00220 00221 // compute sgrt_dg = sqrt(dg) 00222 VectorXd sqrt_dg=eigen_dg.array().sqrt(); 00223 00224 // create shogun and eigen3 representation of labels adjusted for 00225 // noise and means (m_r) 00226 m_r=SGVector<float64_t>(y.vlen); 00227 Map<VectorXd> eigen_r(m_r.vector, m_r.vlen); 00228 eigen_r=(eigen_y-eigen_m).cwiseQuotient(sqrt_dg); 00229 00230 // compute be 00231 m_be=SGVector<float64_t>(m_chol_utr.num_cols); 00232 Map<VectorXd> eigen_be(m_be.vector, m_be.vlen); 00233 eigen_be=eigen_chol_utr.triangularView<Upper>().adjoint().solve( 00234 V*eigen_r.cwiseQuotient(sqrt_dg)); 00235 00236 // compute iKuu 00237 MatrixXd iKuu=Luu.solve(MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)); 00238 00239 // create shogun and eigen3 representation of posterior cholesky 00240 MatrixXd eigen_prod=eigen_chol_utr*eigen_chol_uu; 00241 m_L=SGMatrix<float64_t>(m_kuu.num_rows, m_kuu.num_cols); 00242 Map<MatrixXd> eigen_chol(m_L.matrix, m_L.num_rows, m_L.num_cols); 00243 00244 eigen_chol=eigen_prod.triangularView<Upper>().adjoint().solve( 00245 MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)); 00246 eigen_chol=eigen_prod.triangularView<Upper>().solve(eigen_chol)-iKuu; 00247 } 00248 00249 void CFITCInferenceMethod::update_alpha() 00250 { 00251 Map<MatrixXd> eigen_chol_uu(m_chol_uu.matrix, m_chol_uu.num_rows, 00252 m_chol_uu.num_cols); 00253 Map<MatrixXd> eigen_chol_utr(m_chol_utr.matrix, m_chol_utr.num_rows, 00254 m_chol_utr.num_cols); 00255 Map<VectorXd> eigen_be(m_be.vector, m_be.vlen); 00256 00257 // create shogun and eigen representations of alpha 00258 // and solve Luu * Lu * alpha = be 00259 m_alpha=SGVector<float64_t>(m_chol_uu.num_rows); 00260 Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen); 00261 00262 eigen_alpha=eigen_chol_utr.triangularView<Upper>().solve(eigen_be); 00263 eigen_alpha=eigen_chol_uu.triangularView<Upper>().solve(eigen_alpha); 00264 } 00265 00266 void CFITCInferenceMethod::update_deriv() 00267 { 00268 // create eigen representation of Ktru, Lu, Luu, dg, be 00269 Map<MatrixXd> eigen_Ktru(m_ktru.matrix, m_ktru.num_rows, m_ktru.num_cols); 00270 Map<MatrixXd> eigen_Lu(m_chol_utr.matrix, m_chol_utr.num_rows, 00271 m_chol_utr.num_cols); 00272 Map<MatrixXd> eigen_Luu(m_chol_uu.matrix, m_chol_uu.num_rows, 00273 m_chol_uu.num_cols); 00274 Map<VectorXd> eigen_dg(m_dg.vector, m_dg.vlen); 00275 Map<VectorXd> eigen_be(m_be.vector, m_be.vlen); 00276 00277 // get and create eigen representation of labels 00278 SGVector<float64_t> y=((CRegressionLabels*) m_labels)->get_labels(); 00279 Map<VectorXd> eigen_y(y.vector, y.vlen); 00280 00281 // get and create eigen representation of mean vector 00282 SGVector<float64_t> m=m_mean->get_mean_vector(m_features); 00283 Map<VectorXd> eigen_m(m.vector, m.vlen); 00284 00285 // compute V=inv(Luu')*Ku 00286 MatrixXd V=eigen_Luu.triangularView<Upper>().adjoint().solve(eigen_Ktru* 00287 CMath::sq(m_scale)); 00288 00289 // create shogun and eigen representation of al 00290 m_al=SGVector<float64_t>(m.vlen); 00291 Map<VectorXd> eigen_al(m_al.vector, m_al.vlen); 00292 00293 // compute al=(Kt+sn2*eye(n))\y 00294 eigen_al=((eigen_y-eigen_m)-(V.adjoint()* 00295 eigen_Lu.triangularView<Upper>().solve(eigen_be))).cwiseQuotient(eigen_dg); 00296 00297 // compute inv(Kuu+snu2*I)=iKuu 00298 MatrixXd iKuu=eigen_Luu.triangularView<Upper>().adjoint().solve( 00299 MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)); 00300 iKuu=eigen_Luu.triangularView<Upper>().solve(iKuu); 00301 00302 // create shogun and eigen representation of B 00303 m_B=SGMatrix<float64_t>(iKuu.rows(), eigen_Ktru.cols()); 00304 Map<MatrixXd> eigen_B(m_B.matrix, m_B.num_rows, m_B.num_cols); 00305 00306 eigen_B=iKuu*eigen_Ktru*CMath::sq(m_scale); 00307 00308 // create shogun and eigen representation of w 00309 m_w=SGVector<float64_t>(m_B.num_rows); 00310 Map<VectorXd> eigen_w(m_w.vector, m_w.vlen); 00311 00312 eigen_w=eigen_B*eigen_al; 00313 00314 // create shogun and eigen representation of W 00315 m_W=SGMatrix<float64_t>(m_chol_utr.num_cols, m_dg.vlen); 00316 Map<MatrixXd> eigen_W(m_W.matrix, m_W.num_rows, m_W.num_cols); 00317 00318 // compute W=Lu'\(V./repmat(g_sn2',nu,1)) 00319 eigen_W=eigen_Lu.triangularView<Upper>().adjoint().solve(V*VectorXd::Ones( 00320 m_dg.vlen).cwiseQuotient(eigen_dg).asDiagonal()); 00321 } 00322 00323 SGVector<float64_t> CFITCInferenceMethod::get_derivative_wrt_inference_method( 00324 const TParameter* param) 00325 { 00326 REQUIRE(!strcmp(param->m_name, "scale"), "Can't compute derivative of " 00327 "the nagative log marginal likelihood wrt %s.%s parameter\n", 00328 get_name(), param->m_name) 00329 00330 // create eigen representation of dg, al, B, W, w 00331 Map<VectorXd> eigen_dg(m_dg.vector, m_dg.vlen); 00332 Map<VectorXd> eigen_al(m_al.vector, m_al.vlen); 00333 Map<VectorXd> eigen_w(m_w.vector, m_w.vlen); 00334 Map<MatrixXd> eigen_B(m_B.matrix, m_B.num_rows, m_B.num_cols); 00335 Map<MatrixXd> eigen_W(m_W.matrix, m_W.num_rows, m_W.num_cols); 00336 00337 // clone kernel matrices 00338 SGVector<float64_t> deriv_trtr=m_ktrtr.get_diagonal_vector(); 00339 SGMatrix<float64_t> deriv_uu=m_kuu.clone(); 00340 SGMatrix<float64_t> deriv_tru=m_ktru.clone(); 00341 00342 // create eigen representation of kernel matrices 00343 Map<VectorXd> ddiagKi(deriv_trtr.vector, deriv_trtr.vlen); 00344 Map<MatrixXd> dKuui(deriv_uu.matrix, deriv_uu.num_rows, deriv_uu.num_cols); 00345 Map<MatrixXd> dKui(deriv_tru.matrix, deriv_tru.num_rows, deriv_tru.num_cols); 00346 00347 // compute derivatives wrt scale for each kernel matrix 00348 ddiagKi*=m_scale*2.0; 00349 dKuui*=m_scale*2.0; 00350 dKui*=m_scale*2.0; 00351 00352 // compute R=2*dKui-dKuui*B 00353 MatrixXd R=2*dKui-dKuui*eigen_B; 00354 00355 // compute v=ddiagKi-sum(R.*B, 1)' 00356 VectorXd v=ddiagKi-R.cwiseProduct(eigen_B).colwise().sum().adjoint(); 00357 00358 SGVector<float64_t> result(1); 00359 00360 // compute dnlZ=(ddiagKi'*(1./g_sn2)+w'*(dKuui*w-2*(dKui*al))-al'*(v.*al)- 00361 // sum(W.*W,1)*v- sum(sum((R*W').*(B*W'))))/2; 00362 result[0]=(ddiagKi.dot(VectorXd::Ones(m_dg.vlen).cwiseQuotient(eigen_dg))+ 00363 eigen_w.dot(dKuui*eigen_w-2*(dKui*eigen_al))- 00364 eigen_al.dot(v.cwiseProduct(eigen_al))- 00365 eigen_W.cwiseProduct(eigen_W).colwise().sum().dot(v)- 00366 (R*eigen_W.adjoint()).cwiseProduct(eigen_B*eigen_W.adjoint()).sum())/2.0; 00367 00368 return result; 00369 } 00370 00371 SGVector<float64_t> CFITCInferenceMethod::get_derivative_wrt_likelihood_model( 00372 const TParameter* param) 00373 { 00374 REQUIRE(!strcmp(param->m_name, "sigma"), "Can't compute derivative of " 00375 "the nagative log marginal likelihood wrt %s.%s parameter\n", 00376 m_model->get_name(), param->m_name) 00377 00378 // create eigen representation of dg, al, w, W and B 00379 Map<VectorXd> eigen_dg(m_dg.vector, m_dg.vlen); 00380 Map<VectorXd> eigen_al(m_al.vector, m_al.vlen); 00381 Map<VectorXd> eigen_w(m_w.vector, m_w.vlen); 00382 Map<MatrixXd> eigen_W(m_W.matrix, m_W.num_rows, m_W.num_cols); 00383 Map<MatrixXd> eigen_B(m_B.matrix, m_B.num_rows, m_B.num_cols); 00384 00385 // get the sigma variable from the Gaussian likelihood model 00386 CGaussianLikelihood* lik=CGaussianLikelihood::obtain_from_generic(m_model); 00387 float64_t sigma=lik->get_sigma(); 00388 SG_UNREF(lik); 00389 00390 SGVector<float64_t> result(1); 00391 00392 result[0]=CMath::sq(sigma)*(VectorXd::Ones(m_dg.vlen).cwiseQuotient( 00393 eigen_dg).sum()-eigen_W.cwiseProduct(eigen_W).sum()-eigen_al.dot(eigen_al)); 00394 00395 float64_t dKuui=2.0*m_ind_noise; 00396 MatrixXd R=-dKuui*eigen_B; 00397 VectorXd v=-R.cwiseProduct(eigen_B).colwise().sum().adjoint(); 00398 00399 result[0]=result[0]+((eigen_w.dot(dKuui*eigen_w))-eigen_al.dot( 00400 v.cwiseProduct(eigen_al))-eigen_W.cwiseProduct(eigen_W).colwise().sum().dot(v)- 00401 (R*eigen_W.adjoint()).cwiseProduct(eigen_B*eigen_W.adjoint()).sum())/2.0; 00402 00403 return result; 00404 } 00405 00406 SGVector<float64_t> CFITCInferenceMethod::get_derivative_wrt_kernel( 00407 const TParameter* param) 00408 { 00409 // create eigen representation of dg, al, w, W, B 00410 Map<VectorXd> eigen_dg(m_dg.vector, m_dg.vlen); 00411 Map<VectorXd> eigen_al(m_al.vector, m_al.vlen); 00412 Map<VectorXd> eigen_w(m_w.vector, m_w.vlen); 00413 Map<MatrixXd> eigen_W(m_W.matrix, m_W.num_rows, m_W.num_cols); 00414 Map<MatrixXd> eigen_B(m_B.matrix, m_B.num_rows, m_B.num_cols); 00415 00416 SGVector<float64_t> result; 00417 00418 if (param->m_datatype.m_ctype==CT_VECTOR || 00419 param->m_datatype.m_ctype==CT_SGVECTOR) 00420 { 00421 REQUIRE(param->m_datatype.m_length_y, 00422 "Length of the parameter %s should not be NULL\n", param->m_name) 00423 result=SGVector<float64_t>(*(param->m_datatype.m_length_y)); 00424 } 00425 else 00426 { 00427 result=SGVector<float64_t>(1); 00428 } 00429 00430 for (index_t i=0; i<result.vlen; i++) 00431 { 00432 SGVector<float64_t> deriv_trtr; 00433 SGMatrix<float64_t> deriv_uu; 00434 SGMatrix<float64_t> deriv_tru; 00435 00436 if (result.vlen==1) 00437 { 00438 m_kernel->init(m_features, m_features); 00439 deriv_trtr=m_kernel->get_parameter_gradient(param).get_diagonal_vector(); 00440 00441 m_kernel->init(m_latent_features, m_latent_features); 00442 deriv_uu=m_kernel->get_parameter_gradient(param); 00443 00444 m_kernel->init(m_latent_features, m_features); 00445 deriv_tru=m_kernel->get_parameter_gradient(param); 00446 } 00447 else 00448 { 00449 m_kernel->init(m_features, m_features); 00450 deriv_trtr=m_kernel->get_parameter_gradient(param, i).get_diagonal_vector(); 00451 00452 m_kernel->init(m_latent_features, m_latent_features); 00453 deriv_uu=m_kernel->get_parameter_gradient(param, i); 00454 00455 m_kernel->init(m_latent_features, m_features); 00456 deriv_tru=m_kernel->get_parameter_gradient(param, i); 00457 } 00458 00459 // create eigen representation of derivatives 00460 Map<VectorXd> ddiagKi(deriv_trtr.vector, deriv_trtr.vlen); 00461 Map<MatrixXd> dKuui(deriv_uu.matrix, deriv_uu.num_rows, 00462 deriv_uu.num_cols); 00463 Map<MatrixXd> dKui(deriv_tru.matrix, deriv_tru.num_rows, 00464 deriv_tru.num_cols); 00465 00466 ddiagKi*=CMath::sq(m_scale); 00467 dKuui*=CMath::sq(m_scale); 00468 dKui*=CMath::sq(m_scale); 00469 00470 // compute R=2*dKui-dKuui*B 00471 MatrixXd R=2*dKui-dKuui*eigen_B; 00472 00473 // compute v=ddiagKi-sum(R.*B, 1)' 00474 VectorXd v=ddiagKi-R.cwiseProduct(eigen_B).colwise().sum().adjoint(); 00475 00476 // compute dnlZ=(ddiagKi'*(1./g_sn2)+w'*(dKuui*w-2*(dKui*al))-al'* 00477 // (v.*al)-sum(W.*W,1)*v- sum(sum((R*W').*(B*W'))))/2; 00478 result[i]=(ddiagKi.dot(VectorXd::Ones(m_dg.vlen).cwiseQuotient(eigen_dg))+ 00479 eigen_w.dot(dKuui*eigen_w-2*(dKui*eigen_al))- 00480 eigen_al.dot(v.cwiseProduct(eigen_al))- 00481 eigen_W.cwiseProduct(eigen_W).colwise().sum().dot(v)- 00482 (R*eigen_W.adjoint()).cwiseProduct(eigen_B*eigen_W.adjoint()).sum())/2.0; 00483 } 00484 00485 return result; 00486 } 00487 00488 SGVector<float64_t> CFITCInferenceMethod::get_derivative_wrt_mean( 00489 const TParameter* param) 00490 { 00491 // create eigen representation of al vector 00492 Map<VectorXd> eigen_al(m_al.vector, m_al.vlen); 00493 00494 SGVector<float64_t> result; 00495 00496 if (param->m_datatype.m_ctype==CT_VECTOR || 00497 param->m_datatype.m_ctype==CT_SGVECTOR) 00498 { 00499 REQUIRE(param->m_datatype.m_length_y, 00500 "Length of the parameter %s should not be NULL\n", param->m_name) 00501 00502 result=SGVector<float64_t>(*(param->m_datatype.m_length_y)); 00503 } 00504 else 00505 { 00506 result=SGVector<float64_t>(1); 00507 } 00508 00509 for (index_t i=0; i<result.vlen; i++) 00510 { 00511 SGVector<float64_t> dmu; 00512 00513 if (result.vlen==1) 00514 dmu=m_mean->get_parameter_derivative(m_features, param); 00515 else 00516 dmu=m_mean->get_parameter_derivative(m_features, param, i); 00517 00518 Map<VectorXd> eigen_dmu(dmu.vector, dmu.vlen); 00519 00520 // compute dnlZ=-dm'*al 00521 result[i]=-eigen_dmu.dot(eigen_al); 00522 } 00523 00524 return result; 00525 } 00526 00527 #endif /* HAVE_EIGEN3 */