24 using namespace shogun;
25 using namespace Eigen;
40 void CFITCInferenceMethod::init()
42 SG_ADD((
CSGObject**)&m_latent_features,
"latent_features",
"Latent features",
45 m_latent_features=NULL;
60 SG_SERROR(
"Provided inference is not of type CFITCInferenceMethod!\n")
84 "FITC inference method can only use Gaussian likelihood function\n")
86 "of CRegressionLabels\n")
87 REQUIRE(m_latent_features,
"Latent features should not be NULL\n")
89 "Number of latent features must be greater than zero\n")
115 Map<MatrixXd> eigen_chol_utr(m_chol_utr.
matrix, m_chol_utr.
num_rows,
117 Map<VectorXd> eigen_dg(m_dg.
vector, m_dg.
vlen);
119 Map<VectorXd> eigen_be(m_be.
vector, m_be.
vlen);
123 float64_t result=eigen_chol_utr.diagonal().array().log().sum()+
124 (eigen_dg.array().log().sum()+eigen_r.dot(eigen_r)-eigen_be.dot(eigen_be)+
187 LLT<MatrixXd> Luu(eigen_kuu*
CMath::sq(
m_scale)+m_ind_noise*MatrixXd::Identity(
195 eigen_chol_uu=Luu.matrixU();
198 MatrixXd V=eigen_chol_uu.triangularView<Upper>().adjoint().solve(eigen_ktru*
204 Map<VectorXd> eigen_dg(m_dg.
vector, m_dg.
vlen);
207 VectorXd::Ones(m_dg.
vlen)-(V.cwiseProduct(V)).colwise().sum().adjoint();
210 LLT<MatrixXd> Lu(V*((VectorXd::Ones(m_dg.
vlen)).cwiseQuotient(eigen_dg)).asDiagonal()*
216 Map<MatrixXd> eigen_chol_utr(m_chol_utr.matrix, m_chol_utr.num_rows,
217 m_chol_utr.num_cols);
218 eigen_chol_utr=Lu.matrixU();
227 VectorXd sqrt_dg=eigen_dg.array().
sqrt();
233 eigen_r=(eigen_y-eigen_m).cwiseQuotient(sqrt_dg);
237 Map<VectorXd> eigen_be(m_be.
vector, m_be.
vlen);
238 eigen_be=eigen_chol_utr.triangularView<Upper>().adjoint().solve(
239 V*eigen_r.cwiseQuotient(sqrt_dg));
242 MatrixXd iKuu=Luu.solve(MatrixXd::Identity(m_kuu.
num_rows, m_kuu.
num_cols));
245 MatrixXd eigen_prod=eigen_chol_utr*eigen_chol_uu;
249 eigen_chol=eigen_prod.triangularView<Upper>().adjoint().solve(
251 eigen_chol=eigen_prod.triangularView<Upper>().solve(eigen_chol)-iKuu;
258 Map<MatrixXd> eigen_chol_utr(m_chol_utr.
matrix, m_chol_utr.
num_rows,
260 Map<VectorXd> eigen_be(m_be.
vector, m_be.
vlen);
267 eigen_alpha=eigen_chol_utr.triangularView<Upper>().solve(eigen_be);
268 eigen_alpha=eigen_chol_uu.triangularView<Upper>().solve(eigen_alpha);
279 Map<VectorXd> eigen_dg(m_dg.
vector, m_dg.
vlen);
280 Map<VectorXd> eigen_be(m_be.
vector, m_be.
vlen);
291 MatrixXd V=eigen_Luu.triangularView<Upper>().adjoint().solve(eigen_Ktru*
296 Map<VectorXd> eigen_al(m_al.
vector, m_al.
vlen);
299 eigen_al=((eigen_y-eigen_m)-(V.adjoint()*
300 eigen_Lu.triangularView<Upper>().solve(eigen_be))).cwiseQuotient(eigen_dg);
303 MatrixXd iKuu=eigen_Luu.triangularView<Upper>().adjoint().solve(
305 iKuu=eigen_Luu.triangularView<Upper>().solve(iKuu);
317 eigen_w=eigen_B*eigen_al;
324 eigen_W=eigen_Lu.triangularView<Upper>().adjoint().solve(V*VectorXd::Ones(
325 m_dg.
vlen).cwiseQuotient(eigen_dg).asDiagonal());
331 REQUIRE(!strcmp(param->
m_name,
"scale"),
"Can't compute derivative of "
332 "the nagative log marginal likelihood wrt %s.%s parameter\n",
336 Map<VectorXd> eigen_dg(m_dg.
vector, m_dg.
vlen);
337 Map<VectorXd> eigen_al(m_al.
vector, m_al.
vlen);
348 Map<VectorXd> ddiagKi(deriv_trtr.
vector, deriv_trtr.
vlen);
358 MatrixXd R=2*dKui-dKuui*eigen_B;
361 VectorXd v=ddiagKi-R.cwiseProduct(eigen_B).colwise().sum().adjoint();
367 result[0]=(ddiagKi.dot(VectorXd::Ones(m_dg.
vlen).cwiseQuotient(eigen_dg))+
368 eigen_w.dot(dKuui*eigen_w-2*(dKui*eigen_al))-
369 eigen_al.dot(v.cwiseProduct(eigen_al))-
370 eigen_W.cwiseProduct(eigen_W).colwise().sum().dot(v)-
371 (R*eigen_W.adjoint()).cwiseProduct(eigen_B*eigen_W.adjoint()).sum())/2.0;
379 REQUIRE(!strcmp(param->
m_name,
"sigma"),
"Can't compute derivative of "
380 "the nagative log marginal likelihood wrt %s.%s parameter\n",
384 Map<VectorXd> eigen_dg(m_dg.
vector, m_dg.
vlen);
385 Map<VectorXd> eigen_al(m_al.
vector, m_al.
vlen);
397 result[0]=
CMath::sq(sigma)*(VectorXd::Ones(m_dg.
vlen).cwiseQuotient(
398 eigen_dg).sum()-eigen_W.cwiseProduct(eigen_W).sum()-eigen_al.dot(eigen_al));
401 MatrixXd R=-dKuui*eigen_B;
402 VectorXd v=-R.cwiseProduct(eigen_B).colwise().sum().adjoint();
404 result[0]=result[0]+((eigen_w.dot(dKuui*eigen_w))-eigen_al.dot(
405 v.cwiseProduct(eigen_al))-eigen_W.cwiseProduct(eigen_W).colwise().sum().dot(v)-
406 (R*eigen_W.adjoint()).cwiseProduct(eigen_B*eigen_W.adjoint()).sum())/2.0;
415 Map<VectorXd> eigen_dg(m_dg.
vector, m_dg.
vlen);
416 Map<VectorXd> eigen_al(m_al.
vector, m_al.
vlen);
427 "Length of the parameter %s should not be NULL\n", param->
m_name)
465 Map<VectorXd> ddiagKi(deriv_trtr.
vector, deriv_trtr.
vlen);
476 MatrixXd R=2*dKui-dKuui*eigen_B;
479 VectorXd v=ddiagKi-R.cwiseProduct(eigen_B).colwise().sum().adjoint();
483 result[i]=(ddiagKi.dot(VectorXd::Ones(m_dg.
vlen).cwiseQuotient(eigen_dg))+
484 eigen_w.dot(dKuui*eigen_w-2*(dKui*eigen_al))-
485 eigen_al.dot(v.cwiseProduct(eigen_al))-
486 eigen_W.cwiseProduct(eigen_W).colwise().sum().dot(v)-
487 (R*eigen_W.adjoint()).cwiseProduct(eigen_B*eigen_W.adjoint()).sum())/2.0;
497 Map<VectorXd> eigen_al(m_al.
vector, m_al.
vlen);
505 "Length of the parameter %s should not be NULL\n", param->
m_name)
523 Map<VectorXd> eigen_dmu(dmu.
vector, dmu.
vlen);
526 result[i]=-eigen_dmu.
dot(eigen_al);