Add LDLT Cholesky decomposition (#4130)

* Add LDLT decomposition

* Add LDLT Cholesky solver

* Polish error message

* Add ref to ldlt_factor

src/shogun/mathematics/linalg/LinalgBackendBase.h

@@ -214,6 +214,38 @@ namespace shogun
 		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_SOLVER, SGMatrix)
 #undef BACKEND_GENERIC_CHOLESKY_SOLVER
 
+/**
+ * Wrapper method of LDLT Cholesky decomposition
+ *
+ * @see linalg::ldlt_factor
+ */
+#define BACKEND_GENERIC_LDLT_FACTOR(Type, Container)                           \
+	virtual void ldlt_factor(                                                  \
+	    const Container<Type>& A, Container<Type>& L, SGVector<Type>& d,       \
+	    SGVector<index_t>& p, const bool lower) const                          \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+	}
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_LDLT_FACTOR, SGMatrix)
+#undef BACKEND_GENERIC_LDLT_FACTOR
+
+/**
+ * Wrapper method of LDLT Cholesky solver
+ *
+ * @see linalg::ldlt_solver
+ */
+#define BACKEND_GENERIC_LDLT_SOLVER(Type, Container)                           \
+	virtual SGVector<Type> ldlt_solver(                                        \
+	    const Container<Type>& L, const SGVector<Type>& d,                     \
+	    const SGVector<index_t>& p, const SGVector<Type>& b, const bool lower) \
+	    const                                                                  \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+		return 0;                                                              \
+	}
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_LDLT_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_LDLT_SOLVER
+
 /**
  * Wrapper method of cross entropy.
  *

src/shogun/mathematics/linalg/LinalgBackendEigen.h

@@ -109,6 +109,23 @@ namespace shogun
 		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_SOLVER, SGMatrix)
 #undef BACKEND_GENERIC_CHOLESKY_SOLVER
 
+/** Implementation of @see LinalgBackendBase::ldlt_factor */
+#define BACKEND_GENERIC_LDLT_FACTOR(Type, Container)                           \
+	virtual void ldlt_factor(                                                  \
+	    const Container<Type>& A, Container<Type>& L, SGVector<Type>& d,       \
+	    SGVector<index_t>& p, const bool lower) const;
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_LDLT_FACTOR, SGMatrix)
+#undef BACKEND_GENERIC_LDLT_FACTOR
+
+/** Implementation of @see LinalgBackendBase::ldlt_solver */
+#define BACKEND_GENERIC_LDLT_SOLVER(Type, Container)                           \
+	virtual SGVector<Type> ldlt_solver(                                        \
+	    const Container<Type>& L, const SGVector<Type>& d,                     \
+	    const SGVector<index_t>& p, const SGVector<Type>& b, const bool lower) \
+	    const;
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_LDLT_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_LDLT_SOLVER
+
 /** Implementation of @see linalg::cross_entropy */
 #define BACKEND_GENERIC_CROSS_ENTROPY(Type, Container)                         \
 	virtual Type cross_entropy(                                                \
@@ -432,6 +449,19 @@ namespace shogun
 		SGVector<T> cholesky_solver_impl(
 		    const SGMatrix<T>& L, const SGVector<T>& b, const bool lower) const;
 
+		/** Eigen3 LDLT Cholesky decomposition */
+		template <typename T>
+		void ldlt_factor_impl(
+		    const SGMatrix<T>& A, SGMatrix<T>& L, SGVector<T>& d,
+		    SGVector<index_t>& p, const bool lower) const;
+
+		/** Eigen3 LDLT Cholesky solver */
+		template <typename T>
+		SGVector<T> ldlt_solver_impl(
+		    const SGMatrix<T>& L, const SGVector<T>& d,
+		    const SGVector<index_t>& p, const SGVector<T>& b,
+		    const bool lower) const;
+
 		/** Eigen3 cross_entropy method
 		 * The cross entropy is defined as \f$ H(P,Q) = - \sum_{ij}
 		 * P[i,j]log(Q[i,j]) \f$

src/shogun/mathematics/linalg/LinalgNamespace.h

@@ -106,6 +106,45 @@ namespace shogun
 				return sg_linalg->get_cpu_backend();
 		}
 
+		/** Infer the appropriate backend for linalg operations
+		 * from the input SGVector or SGMatrix (Container).
+		 * Raise error if the backends of the Containers conflict.
+		 *
+		 * @param a The first SGVector/SGMatrix
+		 * @param b The second SGVector/SGMatrix
+		 * @param c The third SGVector/SGMatrix
+		 * @return @see LinalgBackendBase pointer
+		 */
+		template <typename T, template <typename> class Container>
+		LinalgBackendBase* infer_backend(
+		    const Container<T>& a, const Container<T>& b, const Container<T>& c)
+		{
+			if (a.on_gpu() && b.on_gpu() && c.on_gpu())
+			{
+				if (sg_linalg->get_gpu_backend())
+					return sg_linalg->get_gpu_backend();
+				else
+				{
+					SG_SERROR(
+					    "Vector or matrix is on GPU but no GPU backend registered. \
+					  This can happen if the GPU backend was de-activated \
+					  after memory has been transferred to GPU.\n");
+					return NULL;
+				}
+			}
+			else if (a.on_gpu() || b.on_gpu() || c.on_gpu())
+			{
+				SG_SERROR(
+				    "Cannot operate with first vector/matrix on_gpu flag(%d),\
+					second vector/matrix on_gpu flag (%d) and third vector/matrix \
+					on_gpu flag (%d).\n",
+				    a.on_gpu(), b.on_gpu(), c.on_gpu());
+				return NULL;
+			}
+			else
+				return sg_linalg->get_cpu_backend();
+		}
+
 		/**
 		 * Transfers data to GPU memory.
 		 * Shallow-copy of SGVector with vector on CPU if GPU backend not
@@ -590,6 +629,71 @@ namespace shogun
 			    ->cholesky_solver(L, b, lower);
 		}
 
+		/**
+		 * Compute the LDLT cholesky decomposition \f$A = P^{T} L D L^{*} P\f$
+		 * or \f$A = P^{T} U^{*} D U P\f$
+		 * of a positive semidefinite or negative semidefinite Hermitan matrix
+		 *
+		 * @param A The matrix whose LDLT cholesky decomposition is to be
+		 *  computed
+		 * @param L The matrix that saves the triangular LDLT
+		 *  Cholesky factorization (default: lower)
+		 * @param d The vector that saves the diagonal of the diagonal matrix D
+		 * @param p The vector that saves the permutation matrix P as a
+		 * transposition sequence
+		 * @param lower Whether to use L as the upper or lower triangular
+		 *  Cholesky factorization (default:lower)
+		 */
+		template <typename T>
+		void ldlt_factor(
+		    const SGMatrix<T>& A, SGMatrix<T>& L, SGVector<T>& d,
+		    SGVector<index_t>& p, const bool lower = true)
+		{
+			REQUIRE(
+			    A.num_rows == A.num_cols,
+			    "Matrix dimensions (%dx%d) are not square\n", A.num_rows,
+			    A.num_cols);
+			REQUIRE(
+			    A.num_rows == L.num_rows && A.num_cols == L.num_cols,
+			    "Shape of matrix A (%d, %d) doesn't match matrix L (%d, %d)\n",
+			    A.num_rows, A.num_cols, L.num_rows, L.num_rows);
+			REQUIRE(
+			    d.vlen == A.num_rows, "Length of vector d (%d) doesn't match "
+			                          "length of diagonal of matrix L (%d)\n",
+			    d.vlen, A.num_rows);
+			REQUIRE(
+			    p.vlen = A.num_rows, "Length of transpositions vector p (%d) "
+			                         "doesn't match length of diagonal of "
+			                         "matrix L (%d)\n",
+			    p.vlen, A.num_rows);
+
+			infer_backend(A)->ldlt_factor(A, L, d, p, lower);
+		}
+
+		/**
+		 * Solve the linear equations \f$Ax=b\f$, given the LDLT Cholesky
+		 * factorization of A,
+		 * where \f$A\f$ is a positive semidefinite or negative semidefinite
+		 * Hermitan matrix @see ldlt_factor
+		 *
+		 * @param L Triangular matrix, LDLT Cholesky factorization of A
+		 * @param d The diagonal of the diagonal matrix D
+		 * @param p The permuattion matrix P as a
+		 * transposition sequence
+		 * @param b Right-hand side array
+		 * @param lower Whether to use L as the upper or lower triangular
+		 *  Cholesky factorization (default:lower)
+		 * @return \f$\x\f$
+		 */
+		template <typename T>
+		SGVector<T> ldlt_solver(
+		    const SGMatrix<T>& L, const SGVector<T>& d, SGVector<index_t>& p,
+		    const SGVector<T>& b, const bool lower = true)
+		{
+			return infer_backend(L, SGMatrix<T>(d), SGMatrix<T>(b))
+			    ->ldlt_solver(L, d, p, b, lower);
+		}
+
 		/**
 		 * Vector dot-product that works with generic vectors.
 		 *

src/shogun/mathematics/linalg/backend/eigen/Decompositions.cpp

@@ -45,6 +45,16 @@ using namespace shogun;
 DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_FACTOR, SGMatrix)
 #undef BACKEND_GENERIC_CHOLESKY_FACTOR
 
+#define BACKEND_GENERIC_LDLT_FACTOR(Type, Container)                           \
+	void LinalgBackendEigen::ldlt_factor(                                      \
+	    const Container<Type>& A, Container<Type>& L, SGVector<Type>& d,       \
+	    SGVector<index_t>& p, const bool lower) const                          \
+	{                                                                          \
+		return ldlt_factor_impl(A, L, d, p, lower);                            \
+	}
+DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_LDLT_FACTOR, SGMatrix)
+#undef BACKEND_GENERIC_LDLT_FACTOR
+
 #define BACKEND_GENERIC_SVD(Type, Container)                                   \
 	void LinalgBackendEigen::svd(                                              \
 	    const Container<Type>& A, SGVector<Type> s, Container<Type> U,         \
@@ -91,6 +101,33 @@ SGMatrix<T> LinalgBackendEigen::cholesky_factor_impl(
 	return c;
 }
 
+template <typename T>
+void LinalgBackendEigen::ldlt_factor_impl(
+    const SGMatrix<T>& A, SGMatrix<T>& L, SGVector<T>& d, SGVector<index_t>& p,
+    const bool lower) const
+{
+	set_const(L, 0);
+	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
+	typename SGMatrix<T>::EigenMatrixXtMap L_eig = L;
+	typename SGVector<T>::EigenVectorXtMap d_eig = d;
+	typename SGVector<index_t>::EigenVectorXtMap p_eig = p;
+
+	Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>> ldlt(A_eig);
+
+	d_eig = ldlt.vectorD().template cast<T>();
+	if (lower)
+		L_eig = ldlt.matrixL();
+	else
+		L_eig = ldlt.matrixU();
+
+	// flatten N*1 matrix into vector
+	p_eig = ldlt.transpositionsP().indices().template cast<index_t>();
+
+	REQUIRE(
+	    ldlt.info() != Eigen::NumericalIssue,
+	    "The factorization failed because of a zero pivot.\n");
+}
+
 template <typename T>
 void LinalgBackendEigen::svd_impl(
     const SGMatrix<T>& A, SGVector<T>& s, SGMatrix<T>& U, bool thin_U,

src/shogun/mathematics/linalg/backend/eigen/Solvers.cpp

@@ -30,6 +30,8 @@
  * Authors: 2016 Pan Deng, Soumyajit De, Heiko Strathmann, Viktor Gal
  */
 
+#include <cmath>
+#include <shogun/base/range.h>
 #include <shogun/mathematics/linalg/LinalgBackendEigen.h>
 #include <shogun/mathematics/linalg/LinalgMacros.h>
 
@@ -45,6 +47,17 @@ using namespace shogun;
 DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_SOLVER, SGMatrix)
 #undef BACKEND_GENERIC_CHOLESKY_SOLVER
 
+#define BACKEND_GENERIC_LDLT_SOLVER(Type, Container)                           \
+	SGVector<Type> LinalgBackendEigen::ldlt_solver(                            \
+	    const Container<Type>& L, const SGVector<Type>& d,                     \
+	    const SGVector<index_t>& p, const SGVector<Type>& b, const bool lower) \
+	    const                                                                  \
+	{                                                                          \
+		return ldlt_solver_impl(L, d, p, b, lower);                            \
+	}
+DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_LDLT_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_LDLT_SOLVER
+
 #define BACKEND_GENERIC_QR_SOLVER(Type, Container)                             \
 	Container<Type> LinalgBackendEigen::qr_solver(                             \
 	    const SGMatrix<Type>& A, const Container<Type>& b) const               \
@@ -102,6 +115,67 @@ SGVector<T> LinalgBackendEigen::cholesky_solver_impl(
 	return x;
 }
 
+template <typename T>
+SGVector<T> LinalgBackendEigen::ldlt_solver_impl(
+    const SGMatrix<T>& L, const SGVector<T>& d, const SGVector<index_t>& p,
+    const SGVector<T>& b, const bool lower) const
+{
+	SGVector<T> result(b.vlen);
+	set_const(result, 0);
+
+	typename SGMatrix<T>::EigenMatrixXtMap L_eig = L;
+	typename SGVector<T>::EigenVectorXtMap b_eig = b;
+	typename SGVector<T>::EigenVectorXtMap result_eig = result;
+	typename SGVector<index_t>::EigenVectorXtMap p_eig = p;
+	Eigen::Transpositions<Eigen::Dynamic> transpositions(p_eig);
+
+	// result = P b
+	result_eig = transpositions * b_eig;
+
+	// result = L^-1 (P b)
+	if (lower)
+		Eigen::TriangularView<Eigen::Map<typename SGMatrix<T>::EigenMatrixXt, 0,
+		                                 Eigen::Stride<0, 0>>,
+		                      Eigen::Lower>(L_eig)
+		    .solveInPlace(result_eig);
+	else
+		Eigen::TriangularView<Eigen::Map<typename SGMatrix<T>::EigenMatrixXt, 0,
+		                                 Eigen::Stride<0, 0>>,
+		                      Eigen::Upper>(L_eig)
+		    .transpose()
+		    .solveInPlace(result_eig);
+
+	auto tolerance =
+	    1.0 / Eigen::NumTraits<typename Eigen::NumTraits<T>::Real>::highest();
+
+	// result = D^-1 L^-1 P b
+	for (auto i : range(d.vlen))
+	{
+		if (std::abs(d[i]) > tolerance)
+			result_eig.row(i) /= d[i];
+		else
+			result_eig.row(i).setZero();
+	}
+
+	// result = U^-1 (D^-1 L^-1 P b)
+	if (lower)
+		Eigen::TriangularView<Eigen::Map<typename SGMatrix<T>::EigenMatrixXt, 0,
+		                                 Eigen::Stride<0, 0>>,
+		                      Eigen::Lower>(L_eig)
+		    .transpose()
+		    .solveInPlace(result_eig);
+	else
+		Eigen::TriangularView<Eigen::Map<typename SGMatrix<T>::EigenMatrixXt, 0,
+		                                 Eigen::Stride<0, 0>>,
+		                      Eigen::Upper>(L_eig)
+		    .solveInPlace(result_eig);
+
+	// result = P^-1 (U^-1 D^-1 L^-1 P b) = A^-1 b
+	result_eig = transpositions.transpose() * result_eig;
+
+	return result;
+}
+
 template <typename T>
 SGVector<T> LinalgBackendEigen::qr_solver_impl(
     const SGMatrix<T>& A, const SGVector<T>& b) const