30 using namespace shogun;
31 using namespace Eigen;
41 Eigen::Map<Eigen::VectorXd>*
alpha;
43 Eigen::Map<Eigen::MatrixXd>*
K;
46 Eigen::Map<Eigen::VectorXd>*
dl2;
56 virtual double operator() (
double x)
58 Eigen::Map<Eigen::VectorXd> eigen_f(f->vector, f->vlen);
59 Eigen::Map<Eigen::VectorXd> eigen_m(m->vector, m->vlen);
61 *alpha = start_alpha + x*(*dalpha);
62 (eigen_f) = (*K)*(*alpha)*scale*scale+(eigen_m);
65 for (
index_t i = 0; i < eigen_f.rows(); i++)
68 (*dl1) = lik->get_log_probability_derivative_f(lab, (*f), 1);
69 (*mW) = lik->get_log_probability_derivative_f(lab, (*f), 2);
70 float64_t result = ((*alpha).dot(((eigen_f)-(eigen_m))))/2.0;
72 for (
index_t i = 0; i < (*mW).vlen; i++)
77 result -= lik->get_log_probability_f(lab, *f);
100 void CLaplacianInferenceMethod::init()
104 m_opt_tolerance = 1e-6;
131 get_first_feature_obj();
151 void CLaplacianInferenceMethod::check_members()
157 SG_ERROR(
"Expected RegressionLabels\n");
163 SG_ERROR(
"Number of training vectors does not match number of labels\n");
169 get_first_feature_obj();
172 SG_ERROR(
"Specified features are not of type CFeatures\n");
175 SG_ERROR(
"Expected Simple Features\n");
178 SG_ERROR(
"Expected Real Features\n");
186 SG_ERROR(
"Specified features are not of type CFeatures\n");
189 SG_ERROR(
"Expected Simple Features\n");
192 SG_ERROR(
"Expected Real Features\n");
199 SG_ERROR(
"No mean function assigned!\n");
214 Map<VectorXd> eigen_dlp(dlp.
vector, dlp.
vlen);
230 sW_temp(i,j) = sW[i];
237 Map<MatrixXd> eigen_temp_kernel(temp_kernel.
matrix,
240 if (eigen_W.minCoeff() < 0)
246 MatrixXd temp_diagonal(sW.
vlen, sW.
vlen);
247 temp_diagonal.setZero(sW.
vlen, sW.
vlen);
249 for (
index_t s = 0; s < temp_diagonal.rows(); s++)
251 for (
index_t r = 0; r < temp_diagonal.cols(); r++)
252 temp_diagonal(r,s) = W[s];
257 FullPivLU<MatrixXd> lu(A);
259 MatrixXd temp_matrix =
262 VectorXd temp_sum(temp_matrix.rows());
264 for (
index_t i = 0; i < temp_matrix.rows(); i++)
266 for (
index_t j = 0; j < temp_matrix.cols(); j++)
267 temp_sum[i] += temp_matrix(j,i);
275 MatrixXd C = Z.transpose().colPivHouseholderQr().solve(
278 MatrixXd temp_diagonal(sW.
vlen, sW.
vlen);
279 temp_diagonal.setZero(sW.
vlen, sW.
vlen);
282 temp_diagonal(s,s) = sW[s];
284 MatrixXd temp = Z.transpose();
286 Z = Z.transpose().colPivHouseholderQr().solve(temp_diagonal);
288 Z = temp.transpose().colPivHouseholderQr().solve(Z);
290 for (
index_t s = 0; s < Z.rows(); s++)
292 for (
index_t r = 0; r < Z.cols(); r++)
296 VectorXd temp_sum(C.rows());
298 temp_sum.setZero(C.rows());
300 for (
index_t i = 0; i < C.rows(); i++)
302 for (
index_t j = 0; j < C.cols(); j++)
303 temp_sum[i] += C(j,i)*C(j,i);
306 g = (eigen_temp_kernel.diagonal()*m_scale*m_scale-temp_sum)/2.0;
309 Map<VectorXd> eigen_d3lp(d3lp.
vector, d3lp.
vlen);
311 VectorXd dfhat = g.cwiseProduct(eigen_d3lp);
318 3+para_dict.get_num_elements(),
319 3+para_dict.get_num_elements());
321 for (
index_t i = 0; i < para_dict.get_num_elements(); i++)
323 shogun::CMapNode<TParameter*, CSGObject*>*
node =
324 para_dict.get_node_ptr(i);
338 bool deriv_found =
false;
340 Map<VectorXd> eigen_temp_alpha(temp_alpha.
vector,
343 for (
index_t h = 0; h < length; h++)
348 VectorXd mean_dev_temp;
366 for (
index_t d = 0; d < mean_derivatives.
vlen; d++)
367 mean_dev_temp[d] = mean_derivatives[d];
375 for (
index_t d = 0; d < mean_derivatives.
vlen; d++)
376 mean_dev_temp[d] = mean_derivatives[d];
394 dK(d,s) = deriv(d,s);
398 sum[0] = (Z.cwiseProduct(dK)).sum()/2.0;
401 sum = sum - eigen_temp_alpha.transpose()*dK*eigen_temp_alpha/2.0;
403 VectorXd b = dK*eigen_dlp;
408 variables[h] = sum[0];
413 else if (mean_derivatives.
vlen > 0)
415 sum = -eigen_temp_alpha.transpose()*mean_dev_temp;
416 sum = sum - dfhat.transpose()*(mean_dev_temp-eigen_temp_kernel*
418 variables[h] = sum[0];
422 else if (lik_first_deriv[0]+lik_second_deriv[0] !=
CMath::INFTY)
424 Map<VectorXd> eigen_fd(lik_first_deriv.
vector, lik_first_deriv.
vlen);
426 Map<VectorXd> eigen_sd(lik_second_deriv.
vector, lik_second_deriv.
vlen);
428 sum[0] = -g.dot(eigen_sd);
429 sum[0] = sum[0] - eigen_fd.sum();
430 variables[h] = sum[0];
437 gradient.
add(param, variables);
453 Map<VectorXd> eigen_temp_alpha(temp_alpha.
vector,
456 sum[0] = (Z.cwiseProduct(dK)).sum()/2.0;
458 sum = sum - eigen_temp_alpha.transpose()*dK*eigen_temp_alpha/2.0;
460 VectorXd b = dK*eigen_dlp;
462 sum = sum - dfhat.transpose()*(b-eigen_temp_kernel*(Z*b)*
m_scale*
m_scale);
468 gradient.
add(param, scale);
469 para_dict.add(param,
this);
489 Map<VectorXd> eigen_temp_alpha(temp_alpha.
vector, temp_alpha.
vlen);
491 Map<MatrixXd> eigen_temp_kernel(temp_kernel.
matrix,
494 Map<VectorXd> eigen_function(
function.vector,
499 Map<VectorXd> eigen_m_means(m_means.
vector, m_means.
vlen);
501 if (eigen_W.minCoeff() < 0)
505 MatrixXd temp_diagonal(sW.
vlen, sW.
vlen);
506 temp_diagonal.setZero(sW.
vlen, sW.
vlen);
509 temp_diagonal(s,s) = sW[s];
513 FullPivLU<MatrixXd> lu(A);
515 float64_t result = (eigen_temp_alpha.transpose()*(eigen_function-eigen_m_means))[0]/2.0 -
516 lp + log(lu.determinant())/2.0;
527 (eigen_sW*eigen_sW.transpose()).cwiseProduct(eigen_temp_kernel*
m_scale*
m_scale) +
530 MatrixXd chol = L.matrixL();
535 sum += log(
m_L(i,i));
537 float64_t result = eigen_temp_alpha.dot(eigen_function-eigen_m_means)/2.0 -
579 temp_kernel(i,j) = kernel_matrix(i,j);
590 Map<MatrixXd> eigen_temp_kernel(temp_kernel.
matrix,
594 if (eigen_W.minCoeff() < 0)
596 MatrixXd temp_diagonal(sW.
vlen, sW.
vlen);
597 temp_diagonal.setZero(sW.
vlen, sW.
vlen);
599 for (
index_t s = 0; s < temp_diagonal.rows(); s++)
601 for (
index_t r = 0; r < temp_diagonal.cols(); r++)
602 temp_diagonal(s,s) = 1.0/W[s];
608 FullPivLU<MatrixXd> lu(A);
610 MatrixXd chol = -lu.inverse();
614 for (
index_t i = 0; i < chol.rows(); i++)
616 for (
index_t j = 0; j < chol.cols(); j++)
617 m_L(i,j) = chol(i,j);
628 (eigen_sW*eigen_sW.transpose()).cwiseProduct((eigen_temp_kernel*
m_scale*
m_scale)) +
631 MatrixXd chol = L.matrixU();
635 for (
index_t i = 0; i < chol.rows(); i++)
637 for (
index_t j = 0; j < chol.cols(); j++)
638 m_L(i,j) = chol(i,j);
654 VectorXd first_derivative;
657 m_means[i] = temp_mean[i];
665 temp_kernel(i,j) =
m_ktrtr(i,j);
668 Map<MatrixXd> eigen_temp_kernel(temp_kernel.matrix, temp_kernel.num_rows, temp_kernel.num_cols);
670 VectorXd eigen_function;
672 Map<VectorXd> eigen_m_means(m_means.
vector, m_means.
vlen);
680 for (
index_t i = 0; i < temp_alpha.vlen; i++)
683 Map<VectorXd> eigen_temp_alpha2(temp_alpha.vector, temp_alpha.vlen);
685 eigen_function = eigen_temp_kernel*eigen_temp_alpha2*
m_scale*
m_scale+eigen_m_means;
689 for (
index_t i = 0; i < eigen_function.rows(); i++)
690 function[i] = eigen_function[i];
705 Map<VectorXd> eigen_temp_alpha2(temp_alpha.vector, temp_alpha.vlen);
707 eigen_function = eigen_temp_kernel*
m_scale*
m_scale*eigen_temp_alpha2+eigen_m_means;
711 for (
index_t i = 0; i < eigen_function.rows(); i++)
712 function[i] = eigen_function[i];
723 Psi_New = (eigen_temp_alpha2.transpose()*(eigen_function-eigen_m_means))[0]/2.0;
733 if (Psi_Def < Psi_New)
738 temp_alpha[i] = eigen_temp_alpha2[i];
751 Map<VectorXd> eigen_temp_alpha(temp_alpha.vector, temp_alpha.vlen);
764 Map<VectorXd> eigen_d2lp(d2lp.
vector, d2lp.
vlen);
768 while (Psi_Old - Psi_New > m_tolerance && itr < m_max_itr)
776 if (eigen_W.minCoeff() < 0)
783 for (
index_t i = 0; i < eigen_W.rows(); i++)
784 eigen_W[i] += 2.0/(df)*dlp[i]*dlp[i];
788 for (
index_t i = 0; i < eigen_W.rows(); i++)
796 LLT<MatrixXd> L((eigen_sW*eigen_sW.transpose()).cwiseProduct(
800 MatrixXd chol = L.matrixL();
802 MatrixXd temp = L.matrixL();
804 Map<VectorXd> temp_eigen_dlp(dlp.
vector, dlp.
vlen);
806 VectorXd b = eigen_W.cwiseProduct((eigen_function - eigen_m_means)) + temp_eigen_dlp;
808 chol = chol.colPivHouseholderQr().solve(eigen_sW.cwiseProduct(
809 (eigen_temp_kernel*b*m_scale*m_scale)));
811 chol = temp.transpose().colPivHouseholderQr().solve(chol);
813 VectorXd dalpha = b - eigen_sW.cwiseProduct(chol) - eigen_temp_alpha;
818 func.
K = &eigen_temp_kernel;
820 func.
alpha = &eigen_temp_alpha;
824 func.
dl2 = &eigen_d2lp;
830 local_min(0, m_max, m_opt_tolerance, func, Psi_New);
834 temp_alpha[i] = eigen_temp_alpha[i];
836 Map<VectorXd> eigen_dlp(dlp.
vector, dlp.
vlen);
838 eigen_function = eigen_temp_kernel*
m_scale*
m_scale*eigen_temp_alpha+eigen_m_means;
842 for (
index_t i = 0; i < eigen_function.rows(); i++)
843 function[i] = eigen_function[i];
867 if (eigen_W2.minCoeff() > 0)
876 m_alpha[i] = eigen_temp_alpha[i];
880 #endif // HAVE_EIGEN3
881 #endif // HAVE_LAPACK