VarDTCInferenceMethod.cpp

/*
 * Copyright (c) The Shogun Machine Learning Toolbox
 * Written (W) 2015 Wu Lin
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Shogun Development Team.
 *
 * This code specifically adapted from function in varsgpLikelihood.m and varsgpPredict.m
 *
 * The reference paper is
 * Titsias, Michalis K.
 * "Variational learning of inducing variables in sparse Gaussian processes."
 * International Conference on Artificial Intelligence and Statistics. 2009.
 *
 */

#include <shogun/machine/gp/VarDTCInferenceMethod.h>
#include <shogun/machine/gp/GaussianLikelihood.h>
#include <shogun/mathematics/Math.h>
#include <shogun/labels/RegressionLabels.h>
#include <shogun/mathematics/eigen3.h>

using namespace shogun;
using namespace Eigen;

CVarDTCInferenceMethod::CVarDTCInferenceMethod() : CSingleSparseInference()
{
	init();
}

CVarDTCInferenceMethod::CVarDTCInferenceMethod(CKernel* kern, CFeatures* feat,
		CMeanFunction* m, CLabels* lab, CLikelihoodModel* mod, CFeatures* lat)
		: CSingleSparseInference(kern, feat, m, lab, mod, lat)
{
	init();
}

void CVarDTCInferenceMethod::init()
{
	m_yy=0.0;
	m_f3=0.0;
	m_sigma2=0.0;
	m_trk=0.0;
	m_Tmm=SGMatrix<float64_t>();
	m_Tnm=SGMatrix<float64_t>();
	m_inv_Lm=SGMatrix<float64_t>();
	m_inv_La=SGMatrix<float64_t>();
	m_Knm_inv_Lm=SGMatrix<float64_t>();

	SG_ADD(&m_yy, "yy", "yy", MS_NOT_AVAILABLE);
	SG_ADD(&m_f3, "f3", "f3", MS_NOT_AVAILABLE);
	SG_ADD(&m_sigma2, "sigma2", "sigma2", MS_NOT_AVAILABLE);
	SG_ADD(&m_trk, "trk", "trk", MS_NOT_AVAILABLE);
	SG_ADD(&m_Tmm, "Tmm", "Tmm", MS_NOT_AVAILABLE);
	SG_ADD(&m_Tnm, "Tnm", "Tnm", MS_NOT_AVAILABLE);
	SG_ADD(&m_inv_Lm, "inv_Lm", "inv_Lm", MS_NOT_AVAILABLE);
	SG_ADD(&m_inv_La, "inv_La", "inv_La", MS_NOT_AVAILABLE);
	SG_ADD(&m_Knm_inv_Lm, "Knm_Inv_Lm", "Knm_Inv_Lm", MS_NOT_AVAILABLE);
}

CVarDTCInferenceMethod::~CVarDTCInferenceMethod()
{
}

void CVarDTCInferenceMethod::compute_gradient()
{
	CInference::compute_gradient();

	if (!m_gradient_update)
	{
		update_deriv();
		m_gradient_update=true;
		update_parameter_hash();
	}
}

void CVarDTCInferenceMethod::update()
{
	SG_DEBUG("entering\n");

	CInference::update();
	update_chol();
	update_alpha();
	m_gradient_update=false;
	update_parameter_hash();

	SG_DEBUG("leaving\n");
}

CVarDTCInferenceMethod* CVarDTCInferenceMethod::obtain_from_generic(
		CInference* inference)
{
	if (inference==NULL)
		return NULL;

	if (inference->get_inference_type()!=INF_KL_SPARSE_REGRESSION)
		SG_SERROR("Provided inference is not of type CVarDTCInferenceMethod!\n")

	SG_REF(inference);
	return (CVarDTCInferenceMethod*)inference;
}

void CVarDTCInferenceMethod::check_members() const
{
	CSingleSparseInference::check_members();

	REQUIRE(m_model->get_model_type()==LT_GAUSSIAN,
			"VarDTC inference method can only use Gaussian likelihood function\n")
	REQUIRE(m_labels->get_label_type()==LT_REGRESSION, "Labels must be type "
			"of CRegressionLabels\n")
}

SGVector<float64_t> CVarDTCInferenceMethod::get_diagonal_vector()
{
	SG_NOTIMPLEMENTED
	//the inference method does not need to use this 
	return SGVector<float64_t>();
}

float64_t CVarDTCInferenceMethod::get_negative_log_marginal_likelihood()
{
	if (parameter_hash_changed())
		update();

	Map<MatrixXd> eigen_inv_La(m_inv_La.matrix, m_inv_La.num_rows, m_inv_La.num_cols);
	Map<VectorXd> eigen_ktrtr_diag(m_ktrtr_diag.vector, m_ktrtr_diag.vlen);

	//F012 =-(model.n-model.m)*model.Likelihood.logtheta-0.5*model.n*log(2*pi)-(0.5/sigma2)*(model.yy)-sum(log(diag(La)));
	float64_t neg_f012 =
	    (m_ktru.num_cols - m_ktru.num_rows) * std::log(m_sigma2) / 2.0 +
	    0.5 * m_ktru.num_cols * std::log(2 * CMath::PI) +
	    0.5 * m_yy / (m_sigma2)-eigen_inv_La.diagonal().array().log().sum();

	//F3 = (0.5/sigma2)*(yKnmInvLmInvLa*yKnmInvLmInvLa');
	float64_t neg_f3=-m_f3;

	//model.TrKnn = sum(model.diagKnn);
	//TrK = - (0.5/sigma2)*(model.TrKnn  - sum(diag(C)) );
	float64_t neg_trk=-m_trk;

	//F = F012 + F3 + TrK;
	//F = - F;
	return neg_f012+neg_f3+neg_trk;
}

void CVarDTCInferenceMethod::update_chol()
{
	// get the sigma variable from the Gaussian likelihood model
	CGaussianLikelihood* lik = m_model->as<CGaussianLikelihood>();
	float64_t sigma=lik->get_sigma();
	m_sigma2=sigma*sigma;

	//m-by-m matrix
	Map<MatrixXd> eigen_kuu(m_kuu.matrix, m_kuu.num_rows, m_kuu.num_cols);

	//m-by-n matrix
	Map<MatrixXd> eigen_ktru(m_ktru.matrix, m_ktru.num_rows, m_ktru.num_cols);

	Map<VectorXd> eigen_ktrtr_diag(m_ktrtr_diag.vector, m_ktrtr_diag.vlen);

	//Lm = chol(model.Kmm + model.jitter*eye(model.m))
	LLT<MatrixXd> Luu(
	    eigen_kuu * std::exp(m_log_scale * 2.0) +
	    std::exp(m_log_ind_noise) *
	        MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols));
	m_inv_Lm=SGMatrix<float64_t>(Luu.rows(), Luu.cols());
	Map<MatrixXd> eigen_inv_Lm(m_inv_Lm.matrix, m_inv_Lm.num_rows, m_inv_Lm.num_cols);
	//invLm = Lm\eye(model.m); 
	eigen_inv_Lm=Luu.matrixU().solve(MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols));

	m_Knm_inv_Lm=SGMatrix<float64_t>(m_ktru.num_cols, m_ktru.num_rows);
	Map<MatrixXd> eigen_Knm_inv_Lm(m_Knm_inv_Lm.matrix, m_Knm_inv_Lm.num_rows, m_Knm_inv_Lm.num_cols);
	// KnmInvLm = model.Knm*invLm;
	eigen_Knm_inv_Lm =
	    (eigen_ktru.transpose() * std::exp(m_log_scale * 2.0)) * eigen_inv_Lm;

	m_Tmm=SGMatrix<float64_t>(m_kuu.num_rows, m_kuu.num_cols);
	Map<MatrixXd> eigen_C(m_Tmm.matrix, m_Tmm.num_rows, m_Tmm.num_cols);
	//C = KnmInvLm'*KnmInvLm; 
	eigen_C=eigen_Knm_inv_Lm.transpose()*eigen_Knm_inv_Lm;

	m_inv_La=SGMatrix<float64_t>(m_kuu.num_rows, m_kuu.num_cols);
	Map<MatrixXd> eigen_inv_La(m_inv_La.matrix, m_inv_La.num_rows, m_inv_La.num_cols);
	//A = sigma2*eye(model.m) + C;    
	LLT<MatrixXd> chol_A(m_sigma2*MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)+eigen_C);

	//La = chol(A);
	//invLa =  La\eye(model.m);  
	eigen_inv_La=chol_A.matrixU().solve(MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols));

	//L=-invLm*invLm' +  sigma2*(invLm*invLa*invLa'*invLm');
	m_L=SGMatrix<float64_t>(m_kuu.num_rows, m_kuu.num_cols);
	Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
	eigen_L=eigen_inv_Lm*(
		m_sigma2*eigen_inv_La*eigen_inv_La.transpose()-MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)
		)*eigen_inv_Lm.transpose();

	//TrK = - (0.5/sigma2)*(model.TrKnn  - sum(diag(C)) );
	m_trk = -0.5 / (m_sigma2) *
	        (eigen_ktrtr_diag.array().sum() * std::exp(m_log_scale * 2.0) -
	         eigen_C.diagonal().array().sum());
}

void CVarDTCInferenceMethod::update_alpha()
{
	Map<MatrixXd> eigen_Knm_inv_Lm(m_Knm_inv_Lm.matrix, m_Knm_inv_Lm.num_rows, m_Knm_inv_Lm.num_cols);
	Map<MatrixXd> eigen_inv_La(m_inv_La.matrix, m_inv_La.num_rows, m_inv_La.num_cols);
	Map<MatrixXd> eigen_inv_Lm(m_inv_Lm.matrix, m_inv_Lm.num_rows, m_inv_Lm.num_cols);

	SGVector<float64_t> y=((CRegressionLabels*) m_labels)->get_labels();
	Map<VectorXd> eigen_y(y.vector, y.vlen);
	SGVector<float64_t> m=m_mean->get_mean_vector(m_features);
	Map<VectorXd> eigen_m(m.vector, m.vlen);

	//yKnmInvLm = (model.y'*KnmInvLm);  
	//yKnmInvLmInvLa = yKnmInvLm*invLa;     
	VectorXd y_cor=eigen_y-eigen_m;
	VectorXd eigen_y_Knm_inv_Lm_inv_La_transpose=eigen_inv_La.transpose()*(
		eigen_Knm_inv_Lm.transpose()*y_cor);
	//alpha = invLm*invLa*yKnmInvLmInvLa'; 
	m_alpha=SGVector<float64_t>(m_kuu.num_rows);
	Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
	eigen_alpha=eigen_inv_Lm*eigen_inv_La*eigen_y_Knm_inv_Lm_inv_La_transpose;

	m_yy=y_cor.dot(y_cor);
	//F3 = (0.5/sigma2)*(yKnmInvLmInvLa*yKnmInvLmInvLa');
	m_f3=0.5*eigen_y_Knm_inv_Lm_inv_La_transpose.dot(eigen_y_Knm_inv_Lm_inv_La_transpose)/m_sigma2;
}

void CVarDTCInferenceMethod::update_deriv()
{
	Map<MatrixXd> eigen_inv_La(m_inv_La.matrix, m_inv_La.num_rows, m_inv_La.num_cols);
	Map<MatrixXd> eigen_inv_Lm(m_inv_Lm.matrix, m_inv_Lm.num_rows, m_inv_Lm.num_cols);
	Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
	Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
	//m-by-n matrix
	Map<MatrixXd> eigen_ktru(m_ktru.matrix, m_ktru.num_rows, m_ktru.num_cols);
	Map<MatrixXd> eigen_C(m_Tmm.matrix, m_Tmm.num_rows, m_Tmm.num_cols);

	//m_Tmm=SGMatrix<float64_t>(m_kuu.num_rows, m_kuu.num_cols);
	m_Tnm=SGMatrix<float64_t>(m_ktru.num_cols, m_ktru.num_rows);

	Map<MatrixXd> eigen_Tmm(m_Tmm.matrix, m_Tmm.num_rows, m_Tmm.num_cols);
	Map<MatrixXd> eigen_Tnm(m_Tnm.matrix, m_Tnm.num_rows, m_Tnm.num_cols);

	CGaussianLikelihood* lik = m_model->as<CGaussianLikelihood>();
	float64_t sigma=lik->get_sigma();
	m_sigma2=sigma*sigma;

	//invLmInvLa = invLm*invLa;  
	//invA = invLmInvLa*invLmInvLa';    
	//yKnmInvA = yKnmInvLmInvLa*invLmInvLa'; 
	//invKmm = invLm*invLm'; 
	//Tmm = sigma2*invA + yKnmInvA'*yKnmInvA;  
	//Tmm = invKmm - Tmm;
	MatrixXd Tmm=-eigen_L-eigen_alpha*eigen_alpha.transpose();

	// Tnm = model.Knm*Tmm;
	eigen_Tnm = (eigen_ktru.transpose() * std::exp(m_log_scale * 2.0)) * Tmm;

	//Tmm = Tmm - (invLm*(C*invLm'))/sigma2; 
	eigen_Tmm = Tmm - (eigen_inv_Lm*eigen_C*eigen_inv_Lm.transpose()/m_sigma2);

	SGVector<float64_t> y=((CRegressionLabels*) m_labels)->get_labels();
	Map<VectorXd> eigen_y(y.vector, y.vlen);
	SGVector<float64_t> m=m_mean->get_mean_vector(m_features);
	Map<VectorXd> eigen_m(m.vector, m.vlen);

	//Tnm = Tnm + (model.y*yKnmInvA); 
	eigen_Tnm += (eigen_y-eigen_m)*eigen_alpha.transpose();
}

SGVector<float64_t> CVarDTCInferenceMethod::get_posterior_mean()
{
	SG_NOTIMPLEMENTED
	//TODO: implement this method once I get time
	return SGVector<float64_t>();
}

SGMatrix<float64_t> CVarDTCInferenceMethod::get_posterior_covariance()
{
	SG_NOTIMPLEMENTED
	//TODO: implement this method once I get time
	return SGMatrix<float64_t>();
}

SGVector<float64_t> CVarDTCInferenceMethod::get_derivative_wrt_likelihood_model(
		const TParameter* param)
{
	REQUIRE(!strcmp(param->m_name, "log_sigma"), "Can't compute derivative of "
			"the nagative log marginal likelihood wrt %s.%s parameter\n",
			m_model->get_name(), param->m_name)

	SGVector<float64_t> dlik(1);

	Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
	Map<MatrixXd> eigen_inv_La(m_inv_La.matrix, m_inv_La.num_rows, m_inv_La.num_cols);
	Map<MatrixXd> eigen_ktru(m_ktru.matrix, m_ktru.num_rows, m_ktru.num_cols);
	Map<MatrixXd> eigen_kuu(m_kuu.matrix, m_kuu.num_rows, m_kuu.num_cols);
	//yKnmInvLmInvLainvLa = yKnmInvLmInvLa*invLa';
	//sigma2aux = sigma2*sum(sum(invLa.*invLa))  + yKnmInvLmInvLainvLa*yKnmInvLmInvLainvLa';
	float64_t sigma2aux =
	    m_sigma2 * eigen_inv_La.cwiseProduct(eigen_inv_La).array().sum() +
	    eigen_alpha.transpose() *
	        (eigen_kuu * std::exp(m_log_scale * 2.0) +
	         std::exp(m_log_ind_noise) *
	             MatrixXd::Identity(m_kuu.num_rows, m_kuu.num_cols)) *
	        eigen_alpha;
	//Dlik_neg = - (model.n-model.m) + model.yy/sigma2 - 2*F3 - sigma2aux - 2*TrK;
	
	dlik[0]=(m_ktru.num_cols-m_ktru.num_rows)-m_yy/m_sigma2+2.0*m_f3+sigma2aux+2.0*m_trk;
	return dlik;
}

SGVector<float64_t> CVarDTCInferenceMethod::get_derivative_wrt_inducing_features(
	const TParameter* param)
{
	//[DXu DXunm] = kernelSparseGradInd(model, Tmm, Tnm);
    //DXu_neg = DXu + DXunm/model.sigma2;
	
	Map<MatrixXd> eigen_Tmm(m_Tmm.matrix, m_Tmm.num_rows, m_Tmm.num_cols);
	Map<MatrixXd> eigen_Tnm(m_Tnm.matrix, m_Tnm.num_rows, m_Tnm.num_cols);

	int32_t dim=m_inducing_features.num_rows;
	int32_t num_samples=m_inducing_features.num_cols;
	SGVector<float64_t>deriv_lat(dim*num_samples);
	deriv_lat.zero();

	m_lock->lock();
	CFeatures *inducing_features=get_inducing_features();
	//asymtric part (related to xu and x)
	m_kernel->init(inducing_features, m_features);
	for(int32_t lat_idx=0; lat_idx<eigen_Tnm.cols(); lat_idx++)
	{
		Map<VectorXd> deriv_lat_col_vec(deriv_lat.vector+lat_idx*dim,dim);
		//p by n
		SGMatrix<float64_t> deriv_mat=m_kernel->get_parameter_gradient(param, lat_idx);
		Map<MatrixXd> eigen_deriv_mat(deriv_mat.matrix, deriv_mat.num_rows, deriv_mat.num_cols);
		//DXunm/model.sigma2;
		deriv_lat_col_vec +=
		    eigen_deriv_mat *
		    (-std::exp(m_log_scale * 2.0) / m_sigma2 * eigen_Tnm.col(lat_idx));
	}

	//symtric part (related to xu and xu)
	m_kernel->init(inducing_features, inducing_features);
	for(int32_t lat_lidx=0; lat_lidx<eigen_Tmm.cols(); lat_lidx++)
	{
		Map<VectorXd> deriv_lat_col_vec(deriv_lat.vector+lat_lidx*dim,dim);
		//p by n
		SGMatrix<float64_t> deriv_mat=m_kernel->get_parameter_gradient(param, lat_lidx);
		Map<MatrixXd> eigen_deriv_mat(deriv_mat.matrix, deriv_mat.num_rows, deriv_mat.num_cols);
		//DXu
		deriv_lat_col_vec += eigen_deriv_mat * (-std::exp(m_log_scale * 2.0) *
		                                        eigen_Tmm.col(lat_lidx));
	}
	SG_UNREF(inducing_features);
	m_lock->unlock();
	return deriv_lat;
}

SGVector<float64_t> CVarDTCInferenceMethod::get_derivative_wrt_inducing_noise(
	const TParameter* param)
{
	REQUIRE(param, "Param not set\n");
	REQUIRE(!strcmp(param->m_name, "log_inducing_noise"), "Can't compute derivative of "
			"the nagative log marginal likelihood wrt %s.%s parameter\n",
			get_name(), param->m_name)

	Map<MatrixXd> eigen_Tmm(m_Tmm.matrix, m_Tmm.num_rows, m_Tmm.num_cols);
	SGVector<float64_t> result(1);
	result[0] =
	    -0.5 * std::exp(m_log_ind_noise) * eigen_Tmm.diagonal().array().sum();
	return result;
}

float64_t CVarDTCInferenceMethod::get_derivative_related_cov(SGVector<float64_t> ddiagKi,
	SGMatrix<float64_t> dKuui, SGMatrix<float64_t> dKui)
{
	Map<VectorXd> eigen_ddiagKi(ddiagKi.vector, ddiagKi.vlen);
	Map<MatrixXd> eigen_dKuui(dKuui.matrix, dKuui.num_rows, dKuui.num_cols);
	Map<MatrixXd> eigen_dKui(dKui.matrix, dKui.num_rows, dKui.num_cols);

	Map<MatrixXd> eigen_Tmm(m_Tmm.matrix, m_Tmm.num_rows, m_Tmm.num_cols);
	Map<MatrixXd> eigen_Tnm(m_Tnm.matrix, m_Tnm.num_rows, m_Tnm.num_cols);

	//[Dkern Dkernnm DTrKnn] = kernelSparseGradHyp(model, Tmm, Tnm); 
	//Dkern_neg = 0.5*Dkern + Dkernnm/model.sigma2 - (0.5/model.sigma2)*DTrKnn;
	float64_t dkern= -0.5*eigen_dKuui.cwiseProduct(eigen_Tmm).sum()
		-eigen_dKui.cwiseProduct(eigen_Tnm.transpose()).sum()/m_sigma2
		+0.5*eigen_ddiagKi.array().sum()/m_sigma2;
	return dkern;
}

SGVector<float64_t> CVarDTCInferenceMethod::get_derivative_wrt_mean(
	const TParameter* param)
{
	REQUIRE(param, "Param not set\n");
	SGVector<float64_t> result;
	int64_t len=const_cast<TParameter *>(param)->m_datatype.get_num_elements();
	result=SGVector<float64_t>(len);

	SGVector<float64_t> y=((CRegressionLabels*) m_labels)->get_labels();
	Map<VectorXd> eigen_y(y.vector, y.vlen);
	SGVector<float64_t> m=m_mean->get_mean_vector(m_features);
	Map<VectorXd> eigen_m(m.vector, m.vlen);
	Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);

	//m-by-n matrix
	Map<MatrixXd> eigen_ktru(m_ktru.matrix, m_ktru.num_rows, m_ktru.num_cols);

	for (index_t i=0; i<result.vlen; i++)
	{
		SGVector<float64_t> dmu=m_mean->get_parameter_derivative(m_features, param, i);
		Map<VectorXd> eigen_dmu(dmu.vector, dmu.vlen);

		result[i] = eigen_dmu.dot(
		                eigen_ktru.transpose() * std::exp(m_log_scale * 2.0) *
		                    eigen_alpha +
		                (eigen_m - eigen_y)) /
		            m_sigma2;
	}
	return result;
}

void CVarDTCInferenceMethod::register_minimizer(Minimizer* minimizer)
{
	SG_WARNING("The method does not require a minimizer. The provided minimizer will not be used.\n");
}