SHOGUN
v2.0.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 * Copyright (C) 2012 Jacob Walker 00008 * 00009 * Code adapted from Gaussian Process Machine Learning Toolbox 00010 * http://www.gaussianprocess.org/gpml/code/matlab/doc/ 00011 * 00012 */ 00013 #include <shogun/lib/config.h> 00014 00015 #ifdef HAVE_LAPACK 00016 #ifdef HAVE_EIGEN3 00017 00018 #include <shogun/regression/gp/FITCInferenceMethod.h> 00019 #include <shogun/regression/gp/GaussianLikelihood.h> 00020 #include <shogun/mathematics/lapack.h> 00021 #include <shogun/mathematics/Math.h> 00022 #include <shogun/labels/RegressionLabels.h> 00023 #include <shogun/kernel/GaussianKernel.h> 00024 #include <shogun/features/CombinedFeatures.h> 00025 #include <shogun/mathematics/eigen3.h> 00026 00027 00028 using namespace shogun; 00029 using namespace Eigen; 00030 00031 CFITCInferenceMethod::CFITCInferenceMethod() : CInferenceMethod() 00032 { 00033 init(); 00034 update_all(); 00035 update_parameter_hash(); 00036 } 00037 00038 CFITCInferenceMethod::CFITCInferenceMethod(CKernel* kern, CFeatures* feat, 00039 CMeanFunction* m, CLabels* lab, CLikelihoodModel* mod, CFeatures* lat) : 00040 CInferenceMethod(kern, feat, m, lab, mod) 00041 { 00042 init(); 00043 set_latent_features(lat); 00044 update_all(); 00045 } 00046 00047 void CFITCInferenceMethod::init() 00048 { 00049 m_latent_features = NULL; 00050 m_ind_noise = 1e-10; 00051 } 00052 00053 CFITCInferenceMethod::~CFITCInferenceMethod() 00054 { 00055 } 00056 00057 void CFITCInferenceMethod::update_all() 00058 { 00059 if (m_labels) 00060 m_label_vector = 00061 ((CRegressionLabels*) m_labels)->get_labels().clone(); 00062 00063 if (m_features && m_features->has_property(FP_DOT) 00064 && m_features->get_num_vectors()) 00065 { 00066 m_feature_matrix = 00067 ((CDotFeatures*)m_features)->get_computed_dot_feature_matrix(); 00068 } 00069 00070 else if (m_features && m_features->get_feature_class() == C_COMBINED) 00071 { 00072 CDotFeatures* feat = 00073 (CDotFeatures*)((CCombinedFeatures*)m_features)-> 00074 get_first_feature_obj(); 00075 00076 if (feat->get_num_vectors()) 00077 m_feature_matrix = feat->get_computed_dot_feature_matrix(); 00078 00079 SG_UNREF(feat); 00080 } 00081 00082 if (m_latent_features && m_latent_features->has_property(FP_DOT) && 00083 m_latent_features->get_num_vectors()) 00084 { 00085 m_latent_matrix = 00086 ((CDotFeatures*)m_latent_features)-> 00087 get_computed_dot_feature_matrix(); 00088 } 00089 00090 else if (m_latent_features && 00091 m_latent_features->get_feature_class() == C_COMBINED) 00092 { 00093 CDotFeatures* subfeat = 00094 (CDotFeatures*)((CCombinedFeatures*)m_latent_features)-> 00095 get_first_feature_obj(); 00096 00097 if (m_latent_features->get_num_vectors()) 00098 m_latent_matrix = subfeat->get_computed_dot_feature_matrix(); 00099 00100 SG_UNREF(subfeat); 00101 } 00102 00103 update_data_means(); 00104 00105 if (m_kernel) 00106 update_train_kernel(); 00107 00108 if (m_ktrtr.num_cols*m_ktrtr.num_rows && 00109 m_kuu.num_rows*m_kuu.num_cols && 00110 m_ktru.num_cols*m_ktru.num_rows) 00111 { 00112 update_chol(); 00113 update_alpha(); 00114 } 00115 } 00116 00117 void CFITCInferenceMethod::check_members() 00118 { 00119 if (!m_labels) 00120 SG_ERROR("No labels set\n"); 00121 00122 if (m_labels->get_label_type() != LT_REGRESSION) 00123 SG_ERROR("Expected RegressionLabels\n"); 00124 00125 if (!m_features) 00126 SG_ERROR("No features set!\n"); 00127 00128 if (!m_latent_features) 00129 SG_ERROR("No latent features set!\n"); 00130 00131 if (m_labels->get_num_labels() != m_features->get_num_vectors()) 00132 SG_ERROR("Number of training vectors does not match number of labels\n"); 00133 00134 if(m_features->get_feature_class() == C_COMBINED) 00135 { 00136 CDotFeatures* feat = 00137 (CDotFeatures*)((CCombinedFeatures*)m_features)-> 00138 get_first_feature_obj(); 00139 00140 if (!feat->has_property(FP_DOT)) 00141 SG_ERROR("Specified features are not of type CFeatures\n"); 00142 00143 if (feat->get_feature_class() != C_DENSE) 00144 SG_ERROR("Expected Simple Features\n"); 00145 00146 if (feat->get_feature_type() != F_DREAL) 00147 SG_ERROR("Expected Real Features\n"); 00148 00149 SG_UNREF(feat); 00150 } 00151 00152 else 00153 { 00154 if (!m_features->has_property(FP_DOT)) 00155 SG_ERROR("Specified features are not of type CFeatures\n"); 00156 00157 if (m_features->get_feature_class() != C_DENSE) 00158 SG_ERROR("Expected Simple Features\n"); 00159 00160 if (m_features->get_feature_type() != F_DREAL) 00161 SG_ERROR("Expected Real Features\n"); 00162 } 00163 00164 if(m_latent_features->get_feature_class() == C_COMBINED) 00165 { 00166 CDotFeatures* feat = 00167 (CDotFeatures*)((CCombinedFeatures*)m_latent_features)-> 00168 get_first_feature_obj(); 00169 00170 if (!feat->has_property(FP_DOT)) 00171 SG_ERROR("Specified features are not of type CFeatures\n"); 00172 00173 if (feat->get_feature_class() != C_DENSE) 00174 SG_ERROR("Expected Simple Features\n"); 00175 00176 if (feat->get_feature_type() != F_DREAL) 00177 SG_ERROR("Expected Real Features\n"); 00178 00179 SG_UNREF(feat); 00180 } 00181 00182 else 00183 { 00184 if (!m_latent_features->has_property(FP_DOT)) 00185 SG_ERROR("Specified features are not of type CFeatures\n"); 00186 00187 if (m_latent_features->get_feature_class() != C_DENSE) 00188 SG_ERROR("Expected Simple Features\n"); 00189 00190 if (m_latent_features->get_feature_type() != F_DREAL) 00191 SG_ERROR("Expected Real Features\n"); 00192 } 00193 00194 if (m_latent_matrix.num_rows != m_feature_matrix.num_rows) 00195 SG_ERROR("Regular and Latent Features do not match in dimensionality!\n"); 00196 00197 if (!m_kernel) 00198 SG_ERROR( "No kernel assigned!\n"); 00199 00200 if (!m_mean) 00201 SG_ERROR( "No mean function assigned!\n"); 00202 00203 if (m_model->get_model_type() != LT_GAUSSIAN) 00204 { 00205 SG_ERROR("FITC Inference Method can only use " \ 00206 "Gaussian Likelihood Function.\n"); 00207 } 00208 } 00209 00210 CMap<TParameter*, SGVector<float64_t> > CFITCInferenceMethod:: 00211 get_marginal_likelihood_derivatives(CMap<TParameter*, 00212 CSGObject*>& para_dict) 00213 { 00214 check_members(); 00215 00216 if(update_parameter_hash()) 00217 update_all(); 00218 00219 //Get the sigma variable from the likelihood model 00220 float64_t m_sigma = 00221 dynamic_cast<CGaussianLikelihood*>(m_model)->get_sigma(); 00222 00223 Map<MatrixXd> eigen_ktru(m_ktru.matrix, m_ktru.num_rows, m_ktru.num_cols); 00224 00225 MatrixXd W = eigen_ktru; 00226 00227 Map<VectorXd> eigen_dg(m_dg.vector, m_dg.vlen); 00228 00229 for (index_t j = 0; j < eigen_ktru.rows(); j++) 00230 { 00231 for (index_t i = 0; i < eigen_ktru.cols(); i++) 00232 W(i,j) = eigen_ktru(i,j) / sqrt(eigen_dg[j]); 00233 } 00234 00235 Map<MatrixXd> eigen_uu(m_kuu.matrix, m_kuu.num_rows, m_kuu.num_cols); 00236 LLT<MatrixXd> CholW(eigen_uu + W*W.transpose() + 00237 m_ind_noise*MatrixXd::Identity(eigen_uu.rows(), eigen_uu.cols())); 00238 W = CholW.matrixL(); 00239 00240 00241 W = W.colPivHouseholderQr().solve(eigen_ktru); 00242 00243 VectorXd true_lab(m_data_means.vlen); 00244 00245 for (index_t j = 0; j < m_data_means.vlen; j++) 00246 true_lab[j] = m_label_vector[j] - m_data_means[j]; 00247 00248 VectorXd al = W*true_lab.cwiseQuotient(eigen_dg); 00249 00250 al = W.transpose()*al; 00251 00252 al = true_lab - al; 00253 00254 al = al.cwiseQuotient(eigen_dg); 00255 00256 MatrixXd iKuu = eigen_uu.selfadjointView<Eigen::Upper>().llt() 00257 .solve(MatrixXd::Identity(eigen_uu.rows(), eigen_uu.cols())); 00258 00259 MatrixXd B = iKuu*eigen_ktru; 00260 00261 MatrixXd Wdg = W; 00262 00263 for (index_t j = 0; j < eigen_ktru.rows(); j++) 00264 { 00265 for (index_t i = 0; i < eigen_ktru.cols(); i++) 00266 Wdg(i,j) = Wdg(i,j) / eigen_dg[j]; 00267 } 00268 00269 VectorXd w = B*al; 00270 00271 VectorXd sum(1); 00272 sum[0] = 0; 00273 00274 m_kernel->build_parameter_dictionary(para_dict); 00275 m_mean->build_parameter_dictionary(para_dict); 00276 00277 //This will be the vector we return 00278 CMap<TParameter*, SGVector<float64_t> > gradient( 00279 3+para_dict.get_num_elements(), 00280 3+para_dict.get_num_elements()); 00281 00282 for (index_t i = 0; i < para_dict.get_num_elements(); i++) 00283 { 00284 shogun::CMapNode<TParameter*, CSGObject*>* node = 00285 para_dict.get_node_ptr(i); 00286 00287 TParameter* param = node->key; 00288 CSGObject* obj = node->data; 00289 00290 index_t length = 1; 00291 00292 if ((param->m_datatype.m_ctype== CT_VECTOR || 00293 param->m_datatype.m_ctype == CT_SGVECTOR) && 00294 param->m_datatype.m_length_y != NULL) 00295 length = *(param->m_datatype.m_length_y); 00296 00297 SGVector<float64_t> variables(length); 00298 00299 bool deriv_found = false; 00300 00301 for (index_t g = 0; g < length; g++) 00302 { 00303 00304 SGMatrix<float64_t> deriv; 00305 SGMatrix<float64_t> derivtru; 00306 SGMatrix<float64_t> derivuu; 00307 SGVector<float64_t> mean_derivatives; 00308 VectorXd mean_dev_temp; 00309 00310 if (param->m_datatype.m_ctype == CT_VECTOR || 00311 param->m_datatype.m_ctype == CT_SGVECTOR) 00312 { 00313 m_kernel->init(m_features, m_features); 00314 deriv = m_kernel->get_parameter_gradient(param, obj); 00315 00316 m_kernel->init(m_latent_features, m_features); 00317 derivtru = m_kernel->get_parameter_gradient(param, obj); 00318 00319 m_kernel->init(m_latent_features, m_latent_features); 00320 derivuu = m_kernel->get_parameter_gradient(param, obj); 00321 00322 m_kernel->remove_lhs_and_rhs(); 00323 00324 mean_derivatives = m_mean->get_parameter_derivative( 00325 param, obj, m_feature_matrix, g); 00326 00327 for (index_t d = 0; d < mean_derivatives.vlen; d++) 00328 mean_dev_temp[d] = mean_derivatives[d]; 00329 } 00330 00331 else 00332 { 00333 mean_derivatives = m_mean->get_parameter_derivative( 00334 param, obj, m_feature_matrix); 00335 00336 for (index_t d = 0; d < mean_derivatives.vlen; d++) 00337 mean_dev_temp[d] = mean_derivatives[d]; 00338 00339 m_kernel->init(m_features, m_features); 00340 deriv = m_kernel->get_parameter_gradient(param, obj); 00341 00342 m_kernel->init(m_latent_features, m_features); 00343 derivtru = m_kernel->get_parameter_gradient(param, obj); 00344 00345 m_kernel->init(m_latent_features, m_latent_features); 00346 derivuu = m_kernel->get_parameter_gradient(param, obj); 00347 00348 m_kernel->remove_lhs_and_rhs(); 00349 } 00350 00351 sum[0] = 0; 00352 00353 00354 if (deriv.num_cols*deriv.num_rows > 0) 00355 { 00356 MatrixXd ddiagKi(deriv.num_cols, deriv.num_rows); 00357 MatrixXd dKuui(derivuu.num_cols, derivuu.num_rows); 00358 MatrixXd dKui(derivtru.num_cols, derivtru.num_rows); 00359 00360 for (index_t d = 0; d < deriv.num_rows; d++) 00361 { 00362 for (index_t s = 0; s < deriv.num_cols; s++) 00363 ddiagKi(d,s) = deriv(d,s)*m_scale*m_scale; 00364 } 00365 00366 for (index_t d = 0; d < derivuu.num_rows; d++) 00367 { 00368 for (index_t s = 0; s < derivuu.num_cols; s++) 00369 dKuui(d,s) = derivuu(d,s)*m_scale*m_scale; 00370 } 00371 00372 for (index_t d = 0; d < derivtru.num_rows; d++) 00373 { 00374 for (index_t s = 0; s < derivtru.num_cols; s++) 00375 dKui(d,s) = derivtru(d,s)*m_scale*m_scale; 00376 } 00377 00378 MatrixXd R = 2*dKui-dKuui*B; 00379 MatrixXd v = ddiagKi; 00380 MatrixXd temp = R.cwiseProduct(B); 00381 00382 for (index_t d = 0; d < ddiagKi.rows(); d++) 00383 v(d,d) = v(d,d) - temp.col(d).sum(); 00384 00385 sum = sum + ddiagKi.diagonal().transpose()* 00386 VectorXd::Ones(eigen_dg.rows()).cwiseQuotient(eigen_dg); 00387 00388 sum = sum + w.transpose()*(dKuui*w-2*(dKui*al)); 00389 00390 sum = sum - al.transpose()*(v.diagonal().cwiseProduct(al)); 00391 00392 MatrixXd Wdg_temp = Wdg.cwiseProduct(Wdg); 00393 00394 VectorXd Wdg_sum(Wdg.rows()); 00395 00396 for (index_t d = 0; d < Wdg.rows(); d++) 00397 Wdg_sum[d] = Wdg_temp.col(d).sum(); 00398 00399 sum = sum - v.diagonal().transpose()*Wdg_sum; 00400 00401 Wdg_temp = (R*Wdg.transpose()).cwiseProduct(B*Wdg.transpose()); 00402 00403 sum[0] = sum[0] - Wdg_temp.sum(); 00404 00405 sum /= 2.0; 00406 00407 variables[g] = sum[0]; 00408 deriv_found = true; 00409 } 00410 00411 else if (mean_derivatives.vlen > 0) 00412 { 00413 sum = mean_dev_temp*al; 00414 variables[g] = sum[0]; 00415 deriv_found = true; 00416 } 00417 00418 00419 } 00420 00421 if (deriv_found) 00422 gradient.add(param, variables); 00423 00424 } 00425 00426 //Here we take the kernel scale derivative. 00427 { 00428 TParameter* param; 00429 index_t index = get_modsel_param_index("scale"); 00430 param = m_model_selection_parameters->get_parameter(index); 00431 00432 SGVector<float64_t> variables(1); 00433 00434 SGMatrix<float64_t> deriv; 00435 SGMatrix<float64_t> derivtru; 00436 SGMatrix<float64_t> derivuu; 00437 00438 m_kernel->init(m_features, m_features); 00439 deriv = m_kernel->get_kernel_matrix(); 00440 00441 m_kernel->init(m_latent_features, m_features); 00442 derivtru = m_kernel->get_kernel_matrix(); 00443 00444 m_kernel->init(m_latent_features, m_latent_features); 00445 derivuu = m_kernel->get_kernel_matrix(); 00446 00447 m_kernel->remove_lhs_and_rhs(); 00448 00449 MatrixXd ddiagKi(deriv.num_cols, deriv.num_rows); 00450 MatrixXd dKuui(derivuu.num_cols, derivuu.num_rows); 00451 MatrixXd dKui(derivtru.num_cols, derivtru.num_rows); 00452 00453 for (index_t d = 0; d < deriv.num_rows; d++) 00454 { 00455 for (index_t s = 0; s < deriv.num_cols; s++) 00456 ddiagKi(d,s) = deriv(d,s)*m_scale*2.0; 00457 } 00458 00459 for (index_t d = 0; d < derivuu.num_rows; d++) 00460 { 00461 for (index_t s = 0; s < derivuu.num_cols; s++) 00462 dKuui(d,s) = derivuu(d,s)*m_scale*2.0; 00463 } 00464 00465 for (index_t d = 0; d < derivtru.num_rows; d++) 00466 { 00467 for (index_t s = 0; s < derivtru.num_cols; s++) 00468 dKui(d,s) = derivtru(d,s)*m_scale*2.0; 00469 } 00470 00471 MatrixXd R = 2*dKui-dKuui*B; 00472 MatrixXd v = ddiagKi; 00473 MatrixXd temp = R.cwiseProduct(B); 00474 00475 for (index_t d = 0; d < ddiagKi.rows(); d++) 00476 v(d,d) = v(d,d) - temp.col(d).sum(); 00477 00478 sum = sum + ddiagKi.diagonal().transpose()* 00479 00480 VectorXd::Ones(eigen_dg.rows()).cwiseQuotient(eigen_dg); 00481 00482 sum = sum + w.transpose()*(dKuui*w-2*(dKui*al)); 00483 00484 sum = sum - al.transpose()*(v.diagonal().cwiseProduct(al)); 00485 00486 MatrixXd Wdg_temp = Wdg.cwiseProduct(Wdg); 00487 00488 VectorXd Wdg_sum(Wdg.rows()); 00489 00490 for (index_t d = 0; d < Wdg.rows(); d++) 00491 Wdg_sum[d] = Wdg_temp.col(d).sum(); 00492 00493 sum = sum - v.diagonal().transpose()*Wdg_sum; 00494 00495 Wdg_temp = (R*Wdg.transpose()).cwiseProduct(B*Wdg.transpose()); 00496 00497 sum[0] = sum[0] - Wdg_temp.sum(); 00498 00499 sum /= 2.0; 00500 00501 variables[0] = sum[0]; 00502 00503 gradient.add(param, variables); 00504 para_dict.add(param, this); 00505 00506 } 00507 00508 TParameter* param; 00509 index_t index; 00510 00511 index = m_model->get_modsel_param_index("sigma"); 00512 param = m_model->m_model_selection_parameters->get_parameter(index); 00513 00514 sum[0] = 0; 00515 00516 MatrixXd W_temp = W.cwiseProduct(W); 00517 VectorXd W_sum(W_temp.rows()); 00518 00519 for (index_t d = 0; d < W_sum.rows(); d++) 00520 W_sum[d] = W_temp.col(d).sum(); 00521 00522 W_sum = W_sum.cwiseQuotient(eigen_dg.cwiseProduct(eigen_dg)); 00523 00524 sum[0] = W_sum.sum(); 00525 00526 sum = sum + al.transpose()*al; 00527 00528 sum[0] = VectorXd::Ones(eigen_dg.rows()).cwiseQuotient(eigen_dg).sum() - sum[0]; 00529 00530 sum = sum*m_sigma*m_sigma; 00531 float64_t dKuui = 2.0*m_ind_noise; 00532 00533 MatrixXd R = -dKuui*B; 00534 00535 MatrixXd temp = R.cwiseProduct(B); 00536 VectorXd v(temp.rows()); 00537 00538 for (index_t d = 0; d < temp.rows(); d++) 00539 v[d] = temp.col(d).sum(); 00540 00541 sum = sum + (w.transpose()*dKuui*w)/2.0; 00542 00543 sum = sum - al.transpose()*(v.cwiseProduct(al))/2.0; 00544 00545 MatrixXd Wdg_temp = Wdg.cwiseProduct(Wdg); 00546 VectorXd Wdg_sum(Wdg.rows()); 00547 00548 for (index_t d = 0; d < Wdg.rows(); d++) 00549 Wdg_sum[d] = Wdg_temp.col(d).sum(); 00550 00551 sum = sum - v.transpose()*Wdg_sum/2.0; 00552 00553 00554 Wdg_temp = (R*Wdg.transpose()).cwiseProduct(B*Wdg.transpose()); 00555 00556 sum[0] = sum[0] - Wdg_temp.sum()/2.0; 00557 00558 SGVector<float64_t> sigma(1); 00559 00560 sigma[0] = sum[0]; 00561 gradient.add(param, sigma); 00562 para_dict.add(param, m_model); 00563 00564 return gradient; 00565 00566 } 00567 00568 SGVector<float64_t> CFITCInferenceMethod::get_diagonal_vector() 00569 { 00570 SGVector<float64_t> result; 00571 00572 return result; 00573 } 00574 00575 float64_t CFITCInferenceMethod::get_negative_marginal_likelihood() 00576 { 00577 if(update_parameter_hash()) 00578 update_all(); 00579 00580 Map<MatrixXd> eigen_chol_utr(m_chol_utr.matrix, m_chol_utr.num_rows, m_chol_utr.num_cols); 00581 00582 Map<VectorXd> eigen_dg(m_dg.vector, m_dg.vlen); 00583 00584 Map<VectorXd> eigen_r(m_r.vector, m_r.vlen); 00585 00586 Map<VectorXd> eigen_be(m_be.vector, m_be.vlen); 00587 00588 VectorXd temp = eigen_dg; 00589 VectorXd temp2(m_chol_utr.num_cols); 00590 00591 for (index_t i = 0; i < eigen_dg.rows(); i++) 00592 temp[i] = log(eigen_dg[i]); 00593 00594 for (index_t j = 0; j < eigen_chol_utr.rows(); j++) 00595 temp2[j] = log(eigen_chol_utr(j,j)); 00596 00597 VectorXd sum(1); 00598 00599 sum[0] = temp.sum(); 00600 sum = sum + eigen_r.transpose()*eigen_r; 00601 sum = sum - eigen_be.transpose()*eigen_be; 00602 sum[0] += m_label_vector.vlen*log(2*CMath::PI); 00603 sum /= 2.0; 00604 sum[0] += temp2.sum(); 00605 00606 return sum[0]; 00607 } 00608 00609 SGVector<float64_t> CFITCInferenceMethod::get_alpha() 00610 { 00611 if(update_parameter_hash()) 00612 update_all(); 00613 00614 SGVector<float64_t> result(m_alpha); 00615 return result; 00616 } 00617 00618 SGMatrix<float64_t> CFITCInferenceMethod::get_cholesky() 00619 { 00620 if(update_parameter_hash()) 00621 update_all(); 00622 00623 SGMatrix<float64_t> result(m_L); 00624 return result; 00625 } 00626 00627 void CFITCInferenceMethod::update_train_kernel() 00628 { 00629 m_kernel->init(m_features, m_features); 00630 00631 //K(X, X) 00632 SGMatrix<float64_t> kernel_matrix = m_kernel->get_kernel_matrix(); 00633 00634 m_ktrtr=kernel_matrix.clone(); 00635 00636 m_kernel->init(m_latent_features, m_latent_features); 00637 00638 //K(X, X) 00639 kernel_matrix = m_kernel->get_kernel_matrix(); 00640 00641 m_kuu = kernel_matrix.clone(); 00642 for (index_t i = 0; i < kernel_matrix.num_rows; i++) 00643 { 00644 for (index_t j = 0; j < kernel_matrix.num_cols; j++) 00645 m_kuu(i,j) = kernel_matrix(i,j)*m_scale*m_scale; 00646 } 00647 00648 m_kernel->init(m_latent_features, m_features); 00649 00650 kernel_matrix = m_kernel->get_kernel_matrix(); 00651 00652 m_kernel->remove_lhs_and_rhs(); 00653 00654 m_ktru = SGMatrix<float64_t>(kernel_matrix.num_rows, kernel_matrix.num_cols); 00655 00656 for (index_t i = 0; i < kernel_matrix.num_rows; i++) 00657 { 00658 for (index_t j = 0; j < kernel_matrix.num_cols; j++) 00659 m_ktru(i,j) = kernel_matrix(i,j)*m_scale*m_scale; 00660 } 00661 } 00662 00663 00664 void CFITCInferenceMethod::update_chol() 00665 { 00666 check_members(); 00667 00668 //Get the sigma variable from the likelihood model 00669 float64_t m_sigma = 00670 dynamic_cast<CGaussianLikelihood*>(m_model)->get_sigma(); 00671 00672 Map<MatrixXd> eigen_uu(m_kuu.matrix, m_kuu.num_rows, m_kuu.num_cols); 00673 00674 Map<MatrixXd> eigen_ktru(m_ktru.matrix, m_ktru.num_rows, 00675 m_ktru.num_cols); 00676 00677 LLT<MatrixXd> Luu(eigen_uu + 00678 m_ind_noise*MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)); 00679 00680 MatrixXd eigen_chol_uu = Luu.matrixL(); 00681 00682 MatrixXd V = eigen_chol_uu.colPivHouseholderQr().solve(eigen_ktru); 00683 00684 m_chol_uu = SGMatrix<float64_t>(eigen_chol_uu.rows(), 00685 eigen_chol_uu.cols()); 00686 00687 for (index_t f = 0; f < eigen_chol_uu.rows(); f++) 00688 { 00689 for (index_t s = 0; s < eigen_chol_uu.cols(); s++) 00690 m_chol_uu(f,s) = eigen_chol_uu(f,s); 00691 } 00692 00693 MatrixXd temp_V = V.cwiseProduct(V); 00694 00695 VectorXd eigen_dg; 00696 00697 eigen_dg.resize(m_ktrtr.num_cols); 00698 VectorXd sqrt_dg(m_ktrtr.num_cols); 00699 00700 for (index_t i = 0; i < m_ktrtr.num_cols; i++) 00701 { 00702 eigen_dg[i] = m_ktrtr(i,i)*m_scale*m_scale + m_sigma*m_sigma - temp_V.col(i).sum(); 00703 sqrt_dg[i] = sqrt(eigen_dg[i]); 00704 } 00705 00706 m_dg = SGVector<float64_t>(eigen_dg.rows()); 00707 00708 for (index_t i = 0; i < eigen_dg.rows(); i++) 00709 m_dg[i] = eigen_dg[i]; 00710 00711 for (index_t i = 0; i < V.rows(); i++) 00712 { 00713 for (index_t j = 0; j < V.cols(); j++) 00714 V(i,j) /= sqrt_dg[j]; 00715 } 00716 00717 LLT<MatrixXd> Lu(V*V.transpose() + 00718 MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)); 00719 00720 MatrixXd eigen_chol_utr = Lu.matrixL(); 00721 00722 m_chol_utr = SGMatrix<float64_t>(eigen_chol_utr.rows(), 00723 eigen_chol_utr.cols()); 00724 00725 for (index_t f = 0; f < eigen_chol_utr.rows(); f++) 00726 { 00727 for (index_t s = 0; s < eigen_chol_utr.cols(); s++) 00728 m_chol_utr(f,s) = eigen_chol_utr(f,s); 00729 } 00730 00731 VectorXd true_lab(m_data_means.vlen); 00732 00733 for (index_t j = 0; j < m_data_means.vlen; j++) 00734 true_lab[j] = m_label_vector[j] - m_data_means[j]; 00735 00736 VectorXd eigen_r; 00737 00738 eigen_r = true_lab.cwiseQuotient(sqrt_dg); 00739 00740 m_r = SGVector<float64_t>(eigen_r.rows()); 00741 00742 for (index_t j = 0; j < eigen_r.rows(); j++) 00743 m_r[j] = eigen_r[j]; 00744 00745 VectorXd eigen_be = eigen_chol_utr.colPivHouseholderQr().solve(V*eigen_r); 00746 00747 00748 m_be = SGVector<float64_t>(eigen_be.rows()); 00749 00750 for (index_t j = 0; j < eigen_be.rows(); j++) 00751 m_be[j] = eigen_be[j]; 00752 00753 MatrixXd iKuu = eigen_chol_uu.llt().solve( 00754 MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)); 00755 00756 MatrixXd chol = eigen_chol_utr*eigen_chol_uu; 00757 00758 chol *= chol.transpose(); 00759 00760 chol = chol.llt().solve(MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)); 00761 00762 chol = chol - iKuu; 00763 00764 m_L = SGMatrix<float64_t>(chol.rows(), chol.cols()); 00765 00766 for (index_t i = 0; i < chol.rows(); i++) 00767 { 00768 for (index_t j = 0; j < chol.cols(); j++) 00769 m_L(i,j) = chol(i,j); 00770 } 00771 } 00772 00773 void CFITCInferenceMethod::update_alpha() 00774 { 00775 MatrixXd alph; 00776 00777 Map<MatrixXd> eigen_chol_uu(m_chol_uu.matrix, m_chol_uu.num_rows, m_chol_uu.num_cols); 00778 00779 Map<MatrixXd> eigen_chol_utr(m_chol_utr.matrix, m_chol_utr.num_rows, m_chol_utr.num_cols); 00780 00781 Map<VectorXd> eigen_be(m_be.vector, m_be.vlen); 00782 00783 alph = eigen_chol_utr.colPivHouseholderQr().solve(eigen_be); 00784 00785 alph = eigen_chol_uu.colPivHouseholderQr().solve(alph); 00786 00787 m_alpha = SGVector<float64_t>(alph.rows()); 00788 00789 for (index_t i = 0; i < alph.rows(); i++) 00790 m_alpha[i] = alph(i,0); 00791 } 00792 00793 #endif // HAVE_EIGEN3 00794 #endif // HAVE_LAPACK