LLT linear solver using eigen3

src/shogun/mathematics/linalg/linsolver/LLTLinearSolver.cpp

@@ -0,0 +1,99 @@
+/*
+ * 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) 2016 Kunal Arora
+ */
+
+#include <shogun/lib/config.h>
+
+#include <shogun/lib/SGVector.h>
+#include <shogun/lib/SGMatrix.h>
+#include <shogun/mathematics/eigen3.h>
+#include <shogun/mathematics/linalg/linop/DenseMatrixOperator.h>
+#include <shogun/mathematics/linalg/linop/SparseMatrixOperator.h>
+#include <shogun/mathematics/linalg/linsolver/LLTLinearSolver.h>
+
+using namespace Eigen;
+
+namespace shogun
+{
+
+CLLTLinearSolver::CLLTLinearSolver()
+	: CLinearSolver<float64_t, float64_t>()
+{
+}
+
+CLLTLinearSolver::~CLLTLinearSolver()
+{
+}
+
+SGMatrix<float64_t> CLLTLinearSolver::compute_cholesky(
+	CLinearOperator<float64_t>* A)
+{
+	//sanity check
+	REQUIRE(A, "Operator is NULL!\n");
+	CDenseMatrixOperator<float64_t>* op
+		=dynamic_cast<CDenseMatrixOperator<float64_t>*>(A);
+	REQUIRE(op, "Operator is not SparseMatrixOperator type!\n");
+
+	//creating eigen3 dense matrices
+	SGMatrix<float64_t> mat_A=op->get_matrix_operator();
+	Map<MatrixXd> map_A(mat_A.matrix, mat_A.num_rows, mat_A.num_cols);
+
+	SGMatrix<float64_t> L(mat_A.num_rows,mat_A.num_cols);
+	L.set_const(0.0);
+	Map<MatrixXd> map_L(L.matrix, L.num_rows, L.num_cols);
+
+	LLT<MatrixXd> llt(map_A);
+
+	//compute matrix L
+	map_L= llt.matrixL();
+
+	// checking for success
+	if (llt.info()==NumericalIssue)
+		SG_WARNING("Matrix is not Hermitian positive definite!\n");
+
+	//return matrix L
+	return L;
+
+}
+
+SGVector<float64_t> CLLTLinearSolver::triangular_solve(
+	SGMatrix<float64_t> L, SGVector<float64_t> b)
+{
+	//sanity check
+	REQUIRE(L.num_cols==b.vlen, "Dimension mismatch!\n");
+
+	// creating eigen3 maps for vectors
+	SGVector<float64_t> x(L.num_cols);
+	x.set_const(0.0);
+
+	//creating eigen3 matrix
+	Map<MatrixXd> map_L(L.matrix, L.num_rows, L.num_cols);
+
+	//creating eigen3 vectors
+	Map<VectorXd> map_x(x.vector, x.vlen);
+	Map<VectorXd> map_b(b.vector, b.vlen);
+
+	// x=L*^{-1}(L^{1}(b))
+	map_x=map_L.triangularView<Eigen::Lower>().adjoint().solve(map_L.triangularView<Eigen::Lower>().solve(map_b));
+
+	return x;
+
+}
+SGVector<float64_t> CLLTLinearSolver::solve(
+		CLinearOperator<float64_t>* A, SGVector<float64_t> b)
+{
+	//compute the matrix L
+	SGMatrix<float64_t> L=compute_cholesky(A);
+
+	//solve the linear equation using triangular solve
+	SGVector<float64_t>x =triangular_solve(L, b);
+
+	return x;
+}
+
+}

src/shogun/mathematics/linalg/linsolver/LLTLinearSolver.h

@@ -0,0 +1,70 @@
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General turalPublic License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Written (W) 2016 Kunal Arora
+ */
+
+#ifndef LLT_LINEAR_SOLVER_H_
+#define LLT_LINEAR_SOLVER_H_
+
+#include <shogun/lib/config.h>
+
+#include <shogun/mathematics/linalg/linsolver/LinearSolver.h>
+
+namespace shogun
+{
+
+/** @brief Class that provides a solve method for real matrix
+ * linear systems using LLT
+ */
+class CLLTLinearSolver : public CLinearSolver<float64_t, float64_t>
+{
+public:
+	/** default constructor */
+	CLLTLinearSolver();
+
+	/** destructor */
+	virtual ~CLLTLinearSolver();
+
+	/**
+	 * compute_cholesky method for computing the cholesky factor of linear systems
+	 *
+	 * @param A the linear operator of the system
+	 * @return the cholesky factor
+	 */
+	virtual SGMatrix<float64_t> compute_cholesky(
+		CLinearOperator<float64_t>* A);
+
+	/**
+	 * triangular_solve method for solving linear systems given their cholesky factor
+	 *
+	 * @param L the lower triangular matrix L i.e. the cholesky factor
+	 * @param b the vector of the system
+	 * @return the solution vector
+	 */
+	virtual SGVector<float64_t> triangular_solve(
+		SGMatrix<float64_t> L, SGVector<float64_t> b);
+	/**
+	 * solve method for solving inear systems using LLT
+	 *
+	 * @param A the linear operator of the system
+	 * @param b the vector of the system
+	 * @return the solution vector
+	 */
+	virtual SGVector<float64_t> solve(CLinearOperator<float64_t>* A,
+		SGVector<float64_t> b);
+
+	/** @return object name */
+	virtual const char* get_name() const
+	{
+		return "CLLTLinearSolver";
+	}
+
+};
+
+}
+
+#endif // DIRECT_SPARSE_LINEAR_SOLVER_H_

tests/unit/mathematics/linalg/LLTLinearSolver_unittest.cc

@@ -0,0 +1,118 @@
+/*
+ * 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) 2016 Kunal Arora
+ */
+
+#include <shogun/lib/common.h>
+
+#include <shogun/lib/SGVector.h>
+#include <shogun/lib/SGMatrix.h>
+#include <shogun/mathematics/eigen3.h>
+#include <shogun/mathematics/Math.h>
+#include <shogun/mathematics/linalg/linop/DenseMatrixOperator.h>
+#include <shogun/mathematics/linalg/linop/MatrixOperator.h>
+#include <shogun/mathematics/linalg/linsolver/LLTLinearSolver.h>
+#include <gtest/gtest.h>
+
+using namespace shogun;
+using namespace Eigen;
+
+TEST(LLTLinearSolver, solve)
+{
+	const index_t size=2;
+	SGMatrix<float64_t> m(size, size);
+
+	// LLT doesn't work on non-symmetric matrices
+	m(0,0)=2.0;
+	m(0,1)=1.0;
+	m(1,0)=1.0;
+	m(1,1)=2.5;
+
+	CDenseMatrixOperator<float64_t>* A
+		=new CDenseMatrixOperator<float64_t>(m);
+
+	SGVector<float64_t> b(size);
+	b.set_const(1.0);
+
+	CLLTLinearSolver* linear_solver=new CLLTLinearSolver();
+	SGVector<float64_t> x
+		=linear_solver->solve((CLinearOperator<float64_t>*)A, b);
+
+	SGVector<float64_t> Ax=A->apply(x);
+	Map<VectorXd> map_b(b.vector, b.vlen);
+	Map<VectorXd> map_Ax(Ax.vector, Ax.vlen);
+
+	EXPECT_NEAR((map_Ax-map_b).norm(),
+		0.0, 1E-15);
+
+	SG_UNREF(linear_solver);
+	SG_UNREF(A);
+}
+
+TEST(LLTLinearSolver, compute_cholesky)
+{
+	const index_t size=2;
+	SGMatrix<float64_t> m(size, size);
+
+	// LLT doesn't work on non-symmetric matrices
+	m(0,0)=2.0;
+	m(0,1)=1.0;
+	m(1,0)=1.0;
+	m(1,1)=2.5;
+
+	CDenseMatrixOperator<float64_t>* A
+		=new CDenseMatrixOperator<float64_t>(m);
+
+	CLLTLinearSolver* linear_solver=new CLLTLinearSolver();
+	SGMatrix<float64_t> L
+		=linear_solver->compute_cholesky((CLinearOperator<float64_t>*)A);
+
+
+	Map<MatrixXd> map_A(m.matrix,m.num_rows,m.num_cols);
+	Map<MatrixXd> map_L(L.matrix,L.num_rows,L.num_cols);
+	EXPECT_NEAR((map_A-map_L*map_L.transpose()).norm(),
+		0.0, 1E-15);
+
+	SG_UNREF(linear_solver);
+	SG_UNREF(A);
+}
+
+TEST(LLTLinearSolver, triangular_solve)
+{
+	const index_t size=2;
+	SGMatrix<float64_t> m(size, size);
+
+	// LLT doesn't work on non-symmetric matrices
+	m(0,0)=2.0;
+	m(0,1)=1.0;
+	m(1,0)=1.0;
+	m(1,1)=2.5;
+
+	CDenseMatrixOperator<float64_t>* A
+		=new CDenseMatrixOperator<float64_t>(m);
+
+	CLLTLinearSolver* linear_solver=new CLLTLinearSolver();
+	SGMatrix<float64_t> L
+		=linear_solver->compute_cholesky((CLinearOperator<float64_t>*)A);
+
+	SGVector<float64_t> b(size);
+	b.set_const(1.0);
+
+	SGVector<float64_t> x
+		=linear_solver->triangular_solve(L, b);
+
+	SGVector<float64_t> Ax=A->apply(x);
+	Map<VectorXd> map_b(b.vector, b.vlen);
+	Map<VectorXd> map_Ax(Ax.vector, Ax.vlen);
+
+	EXPECT_NEAR((map_Ax-map_b).norm(),
+		0.0, 1E-15);
+
+	SG_UNREF(linear_solver);
+	SG_UNREF(A);
+}
+