27 using namespace shogun;
28 using namespace Eigen;
31 #define CREATE_SGVECTOR(vec, len, sg_type) \
33 if (!vec.vector || vec.vlen!=len) \
34 vec=SGVector<sg_type>(len); \
38 #define CREATE_SGMATRIX(mat, rows, cols, sg_type) \
40 if (!mat.matrix || mat.num_rows!=rows || mat.num_cols!=cols) \
41 mat=SGMatrix<sg_type>(rows, cols); \
60 void CEPInferenceMethod::init()
146 if (m_ttau.
vlen!=n || m_nlZ>nlZ0)
185 while ((
CMath::abs(m_nlZ-nlZ_old)>m_tol && sweep<m_max_sweep) ||
199 tau_n[i]=1.0/m_Sigma(i,i)-m_ttau[i];
200 nu_n[i]=m_mu[i]/m_Sigma(i,i)+mean[i]*tau_n[i]-m_tnu[i];
203 mu_n[i]=nu_n[i]/tau_n[i];
204 s2_n[i]=1.0/tau_n[i];
218 m_tnu[i]=mu/s2-nu_n[i];
226 Map<VectorXd> eigen_tnu(m_tnu.vector, m_tnu.vlen);
227 Map<VectorXd> eigen_mu(m_mu.
vector, m_mu.
vlen);
229 VectorXd eigen_si=eigen_Sigma.col(i);
232 eigen_Sigma=eigen_Sigma-ds2/(1.0+ds2*eigen_si(i))*eigen_si*
236 eigen_mu=eigen_Sigma*eigen_tnu;
248 if (sweep==m_max_sweep &&
CMath::abs(m_nlZ-nlZ_old)>m_tol)
250 SG_ERROR(
"Maximum number (%d) of sweeps reached, but tolerance (%f) was "
251 "not yet reached. You can manually set maximum number of sweeps "
252 "or tolerance to fix this problem.\n", m_max_sweep, m_tol);
271 Map<VectorXd> eigen_tnu(m_tnu.
vector, m_tnu.
vlen);
272 Map<VectorXd> eigen_sttau(m_sttau.
vector, m_sttau.
vlen);
280 VectorXd eigen_v=eigen_L.triangularView<Upper>().adjoint().solve(
282 eigen_v=eigen_L.triangularView<Upper>().solve(eigen_v);
287 eigen_alpha=eigen_tnu-eigen_sttau.cwiseProduct(eigen_v);
294 Map<VectorXd> eigen_sttau(m_sttau.
vector, m_sttau.
vlen);
303 LLT<MatrixXd> eigen_chol((eigen_sttau*eigen_sttau.adjoint()).cwiseProduct(
307 eigen_L=eigen_chol.matrixU();
315 Map<VectorXd> eigen_sttau(m_sttau.
vector, m_sttau.
vlen);
323 MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
337 Map<VectorXd> eigen_tnu(m_tnu.
vector, m_tnu.
vlen);
341 Map<VectorXd> eigen_mu(m_mu.
vector, m_mu.
vlen);
344 eigen_mu=eigen_Sigma*eigen_tnu;
352 Map<VectorXd> eigen_mu(m_mu.
vector, m_mu.
vlen);
353 Map<VectorXd> eigen_tnu(m_tnu.
vector, m_tnu.
vlen);
354 Map<VectorXd> eigen_ttau(m_ttau.
vector, m_ttau.
vlen);
361 VectorXd eigen_tau_n=(VectorXd::Ones(m_ttau.
vlen)).cwiseQuotient(
362 eigen_Sigma.diagonal())-eigen_ttau;
365 VectorXd eigen_nu_n=eigen_mu.cwiseQuotient(eigen_Sigma.diagonal())-
366 eigen_tnu+eigen_m.cwiseProduct(eigen_tau_n);
370 Map<VectorXd> eigen_mu_n(mu_n.
vector, mu_n.
vlen);
372 eigen_mu_n=eigen_nu_n.cwiseQuotient(eigen_tau_n);
376 Map<VectorXd> eigen_s2_n(s2_n.
vector, s2_n.
vlen);
378 eigen_s2_n=(VectorXd::Ones(m_ttau.
vlen)).cwiseQuotient(eigen_tau_n);
384 float64_t nlZ_part1=eigen_L.diagonal().array().log().sum()-lZ-
385 (eigen_tnu.adjoint()*eigen_Sigma).dot(eigen_tnu)/2.0;
388 float64_t nlZ_part2=(eigen_tnu.array().square()/
389 (eigen_tau_n+eigen_ttau).array()).sum()/2.0-(1.0+eigen_ttau.array()/
390 eigen_tau_n.array()).log().sum()/2.0;
394 float64_t nlZ_part3=-(eigen_nu_n-eigen_m.cwiseProduct(eigen_tau_n)).dot(
395 ((eigen_ttau.array()/eigen_tau_n.array()*(eigen_nu_n.array()-
396 eigen_m.array()*eigen_tau_n.array())-2*eigen_tnu.array())/
397 (eigen_ttau.array()+eigen_tau_n.array())).matrix())/2.0;
400 m_nlZ=nlZ_part1+nlZ_part2+nlZ_part3;
407 Map<VectorXd> eigen_sttau(m_sttau.
vector, m_sttau.
vlen);
415 MatrixXd V=eigen_L.triangularView<Upper>().adjoint().solve(
416 MatrixXd(eigen_sttau.asDiagonal()));
417 V=eigen_L.triangularView<Upper>().solve(V);
420 eigen_F=eigen_alpha*eigen_alpha.adjoint()-eigen_sttau.asDiagonal()*V;
426 REQUIRE(!strcmp(param->
m_name,
"scale"),
"Can't compute derivative of "
427 "the nagative log marginal likelihood wrt %s.%s parameter\n",
436 result[0]=-(eigen_F.cwiseProduct(eigen_K)*
m_scale*2.0).
sum()/2.0;
460 "Length of the parameter %s should not be NULL\n", param->
m_name)