Source

Bayesian-Optimization / src / hierarchical_gaussian_process.cpp

/*
-------------------------------------------------------------------------
   This file is part of BayesOpt, an efficient C++ library for 
   Bayesian optimization.

   Copyright (C) 2011-2013 Ruben Martinez-Cantin <rmcantin@unizar.es>
 
   BayesOpt is free software: you can redistribute it and/or modify it 
   under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   BayesOpt is distributed in the hope that it will be useful, but 
   WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with BayesOpt.  If not, see <http://www.gnu.org/licenses/>.
------------------------------------------------------------------------
*/
#include "hierarchical_gaussian_process.hpp"
#include "cholesky.hpp"
#include "trace_ublas.hpp"

namespace bayesopt
{
  namespace ublas = boost::numeric::ublas;

  HierarchicalGaussianProcess::HierarchicalGaussianProcess(size_t dim, 
							   bopt_params params):
    NonParametricProcess(dim, params) {};

  double HierarchicalGaussianProcess::negativeTotalLogLikelihood()
  {
    /*This is the restricted version. For the unrestricted, make p=0
      and remove the last term of loglik*/

    const matrixd K = computeCorrMatrix();
    const size_t n = K.size1();
    const size_t p = mFeatM.size1(); 
    matrixd L(n,n);
    utils::cholesky_decompose(K,L);

    matrixd KF(ublas::trans(mFeatM));
    inplace_solve(L,KF,ublas::lower_tag());
    
    matrixd FKF = prod(ublas::trans(KF),KF);
    matrixd L2(p,p);
    utils::cholesky_decompose(FKF,L2);

    vectord Ky(mGPY);
    inplace_solve(L,Ky,ublas::lower_tag());

    vectord wML = prod(Ky,KF);
    utils::cholesky_solve(L2,wML,ublas::lower());

    vectord alpha = mGPY - prod(wML,mFeatM);
    inplace_solve(L,alpha,ublas::lower_tag());
    double sigma = ublas::inner_prod(alpha,alpha)/(n-p);

    double loglik = .5*(n-p)*log(ublas::inner_prod(alpha,alpha)) 
      + utils::log_trace(L) + utils::log_trace(L2);
    return loglik;
  }
}
  
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.