28 using namespace shogun;
29 using namespace Eigen;
47 void CFITCInferenceMethod::init()
74 get_first_feature_obj();
87 get_computed_dot_feature_matrix();
95 get_first_feature_obj();
117 void CFITCInferenceMethod::check_members()
123 SG_ERROR(
"Expected RegressionLabels\n");
129 SG_ERROR(
"No latent features set!\n");
132 SG_ERROR(
"Number of training vectors does not match number of labels\n");
138 get_first_feature_obj();
141 SG_ERROR(
"Specified features are not of type CFeatures\n");
144 SG_ERROR(
"Expected Simple Features\n");
147 SG_ERROR(
"Expected Real Features\n");
155 SG_ERROR(
"Specified features are not of type CFeatures\n");
158 SG_ERROR(
"Expected Simple Features\n");
161 SG_ERROR(
"Expected Real Features\n");
168 get_first_feature_obj();
171 SG_ERROR(
"Specified features are not of type CFeatures\n");
174 SG_ERROR(
"Expected Simple Features\n");
177 SG_ERROR(
"Expected Real Features\n");
185 SG_ERROR(
"Specified features are not of type CFeatures\n");
188 SG_ERROR(
"Expected Simple Features\n");
191 SG_ERROR(
"Expected Real Features\n");
195 SG_ERROR(
"Regular and Latent Features do not match in dimensionality!\n");
201 SG_ERROR(
"No mean function assigned!\n");
205 SG_ERROR(
"FITC Inference Method can only use " \
206 "Gaussian Likelihood Function.\n");
225 MatrixXd W = eigen_ktru;
227 Map<VectorXd> eigen_dg(m_dg.
vector, m_dg.
vlen);
229 for (
index_t j = 0; j < eigen_ktru.rows(); j++)
231 for (
index_t i = 0; i < eigen_ktru.cols(); i++)
232 W(i,j) = eigen_ktru(i,j) / sqrt(eigen_dg[j]);
236 LLT<MatrixXd> CholW(eigen_uu + W*W.transpose() +
237 m_ind_noise*MatrixXd::Identity(eigen_uu.rows(), eigen_uu.cols()));
241 W = W.colPivHouseholderQr().solve(eigen_ktru);
248 VectorXd al = W*true_lab.cwiseQuotient(eigen_dg);
250 al = W.transpose()*al;
254 al = al.cwiseQuotient(eigen_dg);
256 MatrixXd iKuu = eigen_uu.selfadjointView<Eigen::Upper>().llt()
257 .solve(MatrixXd::Identity(eigen_uu.rows(), eigen_uu.cols()));
259 MatrixXd B = iKuu*eigen_ktru;
263 for (
index_t j = 0; j < eigen_ktru.rows(); j++)
265 for (
index_t i = 0; i < eigen_ktru.cols(); i++)
266 Wdg(i,j) = Wdg(i,j) / eigen_dg[j];
279 3+para_dict.get_num_elements(),
280 3+para_dict.get_num_elements());
282 for (
index_t i = 0; i < para_dict.get_num_elements(); i++)
284 shogun::CMapNode<TParameter*, CSGObject*>*
node =
285 para_dict.get_node_ptr(i);
299 bool deriv_found =
false;
301 for (
index_t g = 0; g < length; g++)
308 VectorXd mean_dev_temp;
327 for (
index_t d = 0; d < mean_derivatives.
vlen; d++)
328 mean_dev_temp[d] = mean_derivatives[d];
336 for (
index_t d = 0; d < mean_derivatives.
vlen; d++)
337 mean_dev_temp[d] = mean_derivatives[d];
378 MatrixXd R = 2*dKui-dKuui*B;
379 MatrixXd v = ddiagKi;
380 MatrixXd temp = R.cwiseProduct(B);
382 for (
index_t d = 0; d < ddiagKi.rows(); d++)
383 v(d,d) = v(d,d) - temp.col(d).sum();
385 sum = sum + ddiagKi.diagonal().transpose()*
386 VectorXd::Ones(eigen_dg.rows()).cwiseQuotient(eigen_dg);
388 sum = sum + w.transpose()*(dKuui*w-2*(dKui*al));
390 sum = sum - al.transpose()*(v.diagonal().cwiseProduct(al));
392 MatrixXd Wdg_temp = Wdg.cwiseProduct(Wdg);
394 VectorXd Wdg_sum(Wdg.rows());
396 for (
index_t d = 0; d < Wdg.rows(); d++)
397 Wdg_sum[d] = Wdg_temp.col(d).sum();
399 sum = sum - v.diagonal().transpose()*Wdg_sum;
401 Wdg_temp = (R*Wdg.transpose()).cwiseProduct(B*Wdg.transpose());
403 sum[0] = sum[0] - Wdg_temp.sum();
407 variables[g] = sum[0];
411 else if (mean_derivatives.
vlen > 0)
413 sum = mean_dev_temp*al;
414 variables[g] = sum[0];
422 gradient.
add(param, variables);
456 ddiagKi(d,s) = deriv(d,s)*
m_scale*2.0;
462 dKuui(d,s) = derivuu(d,s)*
m_scale*2.0;
468 dKui(d,s) = derivtru(d,s)*
m_scale*2.0;
471 MatrixXd R = 2*dKui-dKuui*B;
472 MatrixXd v = ddiagKi;
473 MatrixXd temp = R.cwiseProduct(B);
475 for (
index_t d = 0; d < ddiagKi.rows(); d++)
476 v(d,d) = v(d,d) - temp.col(d).sum();
478 sum = sum + ddiagKi.diagonal().transpose()*
480 VectorXd::Ones(eigen_dg.rows()).cwiseQuotient(eigen_dg);
482 sum = sum + w.transpose()*(dKuui*w-2*(dKui*al));
484 sum = sum - al.transpose()*(v.diagonal().cwiseProduct(al));
486 MatrixXd Wdg_temp = Wdg.cwiseProduct(Wdg);
488 VectorXd Wdg_sum(Wdg.rows());
490 for (
index_t d = 0; d < Wdg.rows(); d++)
491 Wdg_sum[d] = Wdg_temp.col(d).sum();
493 sum = sum - v.diagonal().transpose()*Wdg_sum;
495 Wdg_temp = (R*Wdg.transpose()).cwiseProduct(B*Wdg.transpose());
497 sum[0] = sum[0] - Wdg_temp.sum();
501 variables[0] = sum[0];
503 gradient.
add(param, variables);
504 para_dict.add(param,
this);
516 MatrixXd W_temp = W.cwiseProduct(W);
517 VectorXd W_sum(W_temp.rows());
519 for (
index_t d = 0; d < W_sum.rows(); d++)
520 W_sum[d] = W_temp.col(d).sum();
522 W_sum = W_sum.cwiseQuotient(eigen_dg.cwiseProduct(eigen_dg));
524 sum[0] = W_sum.sum();
526 sum = sum + al.transpose()*al;
528 sum[0] = VectorXd::Ones(eigen_dg.rows()).cwiseQuotient(eigen_dg).sum() - sum[0];
530 sum = sum*m_sigma*m_sigma;
533 MatrixXd R = -dKuui*B;
535 MatrixXd temp = R.cwiseProduct(B);
536 VectorXd v(temp.rows());
538 for (
index_t d = 0; d < temp.rows(); d++)
539 v[d] = temp.col(d).sum();
541 sum = sum + (w.transpose()*dKuui*w)/2.0;
543 sum = sum - al.transpose()*(v.cwiseProduct(al))/2.0;
545 MatrixXd Wdg_temp = Wdg.cwiseProduct(Wdg);
546 VectorXd Wdg_sum(Wdg.rows());
548 for (
index_t d = 0; d < Wdg.rows(); d++)
549 Wdg_sum[d] = Wdg_temp.col(d).sum();
551 sum = sum - v.transpose()*Wdg_sum/2.0;
554 Wdg_temp = (R*Wdg.transpose()).cwiseProduct(B*Wdg.transpose());
556 sum[0] = sum[0] - Wdg_temp.sum()/2.0;
561 gradient.
add(param, sigma);
582 Map<VectorXd> eigen_dg(m_dg.
vector, m_dg.
vlen);
586 Map<VectorXd> eigen_be(m_be.
vector, m_be.
vlen);
588 VectorXd temp = eigen_dg;
589 VectorXd temp2(m_chol_utr.
num_cols);
591 for (
index_t i = 0; i < eigen_dg.rows(); i++)
592 temp[i] = log(eigen_dg[i]);
594 for (
index_t j = 0; j < eigen_chol_utr.rows(); j++)
595 temp2[j] = log(eigen_chol_utr(j,j));
600 sum = sum + eigen_r.transpose()*eigen_r;
601 sum = sum - eigen_be.transpose()*eigen_be;
604 sum[0] += temp2.sum();
641 m_kuu = kernel_matrix.
clone();
677 LLT<MatrixXd> Luu(eigen_uu +
680 MatrixXd eigen_chol_uu = Luu.matrixL();
682 MatrixXd V = eigen_chol_uu.colPivHouseholderQr().solve(eigen_ktru);
685 eigen_chol_uu.cols());
687 for (
index_t f = 0; f < eigen_chol_uu.rows(); f++)
689 for (
index_t s = 0; s < eigen_chol_uu.cols(); s++)
690 m_chol_uu(f,s) = eigen_chol_uu(f,s);
693 MatrixXd temp_V = V.cwiseProduct(V);
703 sqrt_dg[i] = sqrt(eigen_dg[i]);
708 for (
index_t i = 0; i < eigen_dg.rows(); i++)
709 m_dg[i] = eigen_dg[i];
711 for (
index_t i = 0; i < V.rows(); i++)
713 for (
index_t j = 0; j < V.cols(); j++)
714 V(i,j) /= sqrt_dg[j];
717 LLT<MatrixXd> Lu(V*V.transpose() +
720 MatrixXd eigen_chol_utr = Lu.matrixL();
723 eigen_chol_utr.cols());
725 for (
index_t f = 0; f < eigen_chol_utr.rows(); f++)
727 for (
index_t s = 0; s < eigen_chol_utr.cols(); s++)
728 m_chol_utr(f,s) = eigen_chol_utr(f,s);
738 eigen_r = true_lab.cwiseQuotient(sqrt_dg);
742 for (
index_t j = 0; j < eigen_r.rows(); j++)
745 VectorXd eigen_be = eigen_chol_utr.colPivHouseholderQr().solve(V*eigen_r);
750 for (
index_t j = 0; j < eigen_be.rows(); j++)
751 m_be[j] = eigen_be[j];
753 MatrixXd iKuu = eigen_chol_uu.llt().solve(
756 MatrixXd chol = eigen_chol_utr*eigen_chol_uu;
758 chol *= chol.transpose();
760 chol = chol.llt().solve(MatrixXd::Identity(m_kuu.
num_rows, m_kuu.
num_cols));
766 for (
index_t i = 0; i < chol.rows(); i++)
768 for (
index_t j = 0; j < chol.cols(); j++)
769 m_L(i,j) = chol(i,j);
781 Map<VectorXd> eigen_be(m_be.
vector, m_be.
vlen);
783 alph = eigen_chol_utr.colPivHouseholderQr().solve(eigen_be);
785 alph = eigen_chol_uu.colPivHouseholderQr().solve(alph);
789 for (
index_t i = 0; i < alph.rows(); i++)
793 #endif // HAVE_EIGEN3
794 #endif // HAVE_LAPACK