inverse_benchmark.cpp

/*
 * This program 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.
 *
 * Written (W) 2013 Soumyajit De
 */

#include <shogun/base/init.h>
#include <shogun/lib/config.h>
#include <shogun/io/SGIO.h>
#include <shogun/lib/SGMatrix.h>
#include <shogun/lib/SGVector.h>
#include <shogun/mathematics/eigen3.h>
#include <shogun/mathematics/Math.h>
#include <shogun/lib/Time.h>
#include <shogun/mathematics/Statistics.h>
#include <iostream>
#include <map>

using namespace shogun;
using namespace Eigen;
using namespace std;

SGMatrix<float64_t> gen(int size)
{
	const index_t n=size;
	SGMatrix<float64_t> m(n, n);

	// creating a pd matrix
	for (index_t i=0; i<n*n; ++i)
		m.matrix[i]=CMath::randn_double();
	Map<MatrixXd> eig_m(m.matrix, m.num_rows, m.num_cols);
	eig_m=eig_m*eig_m.transpose()+MatrixXd::Identity(n, n);
	m.center();
	return m;
}

std::pair<float64_t,float64_t> test1(int size)
{
	SGMatrix<float64_t> m=gen(size);
	const index_t n=size;
	float64_t eps=1e-2;

	CTime* time=new CTime();

	Map<MatrixXd> mat(m.matrix, n, n);
	SGMatrix<float64_t> res(n, n);
	res.zero();
	Map<MatrixXd> eig_res(res.matrix, n, n);
	MatrixXd eig_m=mat+MatrixXd::Identity(n,n)*n*eps;

	time->start();
	eig_res=mat*eig_m.inverse();
	float64_t elapsed_b=time->cur_time_diff();
//	cout << "naive" << endl << eig_res << endl;

	// compute cholesky decomposition
	time->start();
	LLT<MatrixXd> chol(eig_m);

	// compute the inverse by solving systems
	VectorXd e=VectorXd::Zero(n);
	for (index_t i=0; i<n; ++i)
	{
		e(i)=1;
		const VectorXd& x=chol.solve(e);
#pragma omp parallel for shared (eig_res, mat, x, i)
		for (index_t j=0; j<n; ++j)
		{
			// since mat is symmetric we can use mat.col instead of mat.row here
			// for faster execution since matrices are col-major
			eig_res(j,i)=x.dot(mat.col(j));
		}
		e(i)=0;
	}
	float64_t elapsed_a=time->cur_time_diff();
//	cout << "cholesky" << endl << eig_res << endl;

	//SG_SPRINT("Time taken: %f (cholesky) %f (naive)\n", elapsed_a, elapsed_b);
	return make_pair(elapsed_a, elapsed_b);
}

std::pair<float64_t,float64_t> test2(int size)
{
	SGMatrix<float64_t> m=gen(size);
	const index_t n=size;
	float64_t eps=1e-2;

	CTime* time=new CTime();

	Map<MatrixXd> mat(m.matrix, n, n);
	SGMatrix<float64_t> res(n, n);
	res.zero();
	Map<MatrixXd> eig_res(res.matrix, n, n);
	MatrixXd eig_m=mat+MatrixXd::Identity(n,n)*n*eps;

	time->start();
	eig_res=mat*eig_m.inverse();
	float64_t elapsed_b=time->cur_time_diff();
//	cout << "naive" << endl << eig_res << endl;

	// compute cholesky decomposition
	time->start();
	LLT<MatrixXd> chol(eig_m);

	// compute the inverse by solving systems
	VectorXd e=VectorXd::Zero(n);
	for (index_t i=0; i<n; ++i)
	{
		e(i)=1;
		eig_res.col(i)=chol.solve(e);
		e(i)=0;
	}
	eig_res=mat*eig_res;
	float64_t elapsed_a=time->cur_time_diff();
//	cout << "cholesky" << endl << eig_res << endl;

	//SG_SPRINT("Time taken: %f (cholesky) %f (naive)\n", elapsed_a, elapsed_b);
	return make_pair(elapsed_a, elapsed_b);
}

int main(int argc, char** argv)
{
	init_shogun_with_defaults();
//	sg_io->set_loglevel(MSG_DEBUG);
//	sg_io->set_location_info(MSG_FUNCTION);
	sg_rand->set_seed(123);
	const index_t iter=100;
	const index_t size=1000;
	SGVector<float64_t> time_a(iter), time_b(iter);
	for (index_t i=0; i<iter; ++i)
	{
		std::pair<float64_t,float64_t> time=test1(size);
		time_a[i]=time.first;
		time_b[i]=time.second;
	}
	SG_SPRINT("Without storing the inverse\n");
	SG_SPRINT("Avg time: %f (cholesky) %f (naive)\n", CStatistics::mean(time_a), CStatistics::mean(time_b));
	SG_SPRINT("Variance: %f (cholesky) %f (naive)\n", CStatistics::variance(time_a), CStatistics::variance(time_b));

	for (index_t i=0; i<iter; ++i)
	{
		std::pair<float64_t,float64_t> time=test2(size);
		time_a[i]=time.first;
		time_b[i]=time.second;
	}
	SG_SPRINT("Storing the inverse\n");
	SG_SPRINT("Avg time: %f (cholesky) %f (naive)\n", CStatistics::mean(time_a), CStatistics::mean(time_b));
	SG_SPRINT("Variance: %f (cholesky) %f (naive)\n", CStatistics::variance(time_a), CStatistics::variance(time_b));
	exit_shogun();
	return 0;
}

inverse_benchmark.out

Without storing the inverse
Avg time: 1.323329 (cholesky) 0.380563 (naive)
Variance: 0.036054 (cholesky) 0.000549 (naive)
Storing the inverse
Avg time: 0.662864 (cholesky) 0.366209 (naive)
Variance: 0.000435 (cholesky) 0.000333 (naive)