From 333cf7142fa00d7c97a91baad0294b952b798b8b Mon Sep 17 00:00:00 2001
From: Michele Mazzoni <mic.mazzoni@gmail.com>
Date: Mon, 26 Jun 2017 16:23:43 +0200
Subject: [PATCH] Add to linalg (CPU-only): add_scalar, center_matrix,
 eigen_solver, eigen_solver_symmetric, qr_solver, svd, transpose_matrix,
 triangular_solver.

---
 .../mathematics/linalg/LinalgBackendBase.h    | 125 +++++++
 .../mathematics/linalg/LinalgBackendEigen.h   | 126 +++++++
 .../mathematics/linalg/LinalgNamespace.h      | 206 ++++++++++++
 .../linalg/backend/eigen/BasicOps.cpp         |  23 ++
 .../linalg/backend/eigen/Decompositions.cpp   |  89 +++++
 .../mathematics/linalg/backend/eigen/Misc.cpp |  45 +++
 .../linalg/backend/eigen/Solvers.cpp          | 107 ++++++
 .../operations/Eigen3_operations_unittest.cc  | 313 ++++++++++++++++++
 8 files changed, 1034 insertions(+)

diff --git a/src/shogun/mathematics/linalg/LinalgBackendBase.h b/src/shogun/mathematics/linalg/LinalgBackendBase.h
index 8dd02e1ccf5..8f4967b70f9 100644
--- a/src/shogun/mathematics/linalg/LinalgBackendBase.h
+++ b/src/shogun/mathematics/linalg/LinalgBackendBase.h
@@ -39,6 +39,7 @@
 #include <shogun/lib/SGVector.h>
 #include <shogun/lib/common.h>
 #include <shogun/lib/config.h>
+#include <shogun/mathematics/Math.h>
 #include <shogun/mathematics/linalg/GPUMemoryBase.h>
 #include <shogun/mathematics/linalg/internal/Block.h>
 
@@ -121,6 +122,33 @@ namespace shogun
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD_COL_VEC, SGMatrix)
 #undef BACKEND_GENERIC_ADD_COL_VEC
 
+/**
+ * Wrapper method of add scalar operation.
+ *
+ * @see linalg::add_scalar
+ */
+#define BACKEND_GENERIC_ADD_SCALAR(Type, Container)                            \
+	virtual void add_scalar(Container<Type>& a, Type b) const                  \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+	}
+		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD_SCALAR, SGVector)
+		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD_SCALAR, SGMatrix)
+#undef BACKEND_GENERIC_ADD_SCALAR
+
+/**
+ * Wrapper method of center matrix operation.
+ *
+ * @see linalg::center_matrix
+ */
+#define BACKEND_GENERIC_CENTER_MATRIX(Type, Container)                         \
+	virtual void center_matrix(Container<Type>& A) const                       \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+	}
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CENTER_MATRIX, SGMatrix)
+#undef BACKEND_GENERIC_CENTER_MATRIX
+
 /**
  * Wrapper method of Cholesky decomposition.
  *
@@ -166,6 +194,38 @@ namespace shogun
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_DOT, SGVector)
 #undef BACKEND_GENERIC_DOT
 
+/**
+ * Wrapper method of eigenvalues and eigenvectors computation.
+ *
+ * @see linalg::eigen_solver
+ */
+#define BACKEND_GENERIC_EIGEN_SOLVER(Type, Container)                          \
+	virtual void eigen_solver(                                                 \
+	    const Container<Type>& A, SGVector<Type>& eigenvalues,                 \
+	    SGMatrix<Type>& eigenvectors) const                                    \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+	}
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_EIGEN_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_EIGEN_SOLVER
+
+/**
+ * Wrapper method of eigenvalues and eigenvectors computation
+ * for symmetric matrices.
+ *
+ * @see linalg::eigen_solver_symmetric
+ */
+#define BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC(Type, Container)                \
+	virtual void eigen_solver_symmetric(                                       \
+	    const Container<Type>& A, SGVector<Type>& eigenvalues,                 \
+	    SGMatrix<Type>& eigenvectors) const                                    \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+	}
+		DEFINE_FOR_NON_INTEGER_PTYPE(
+		    BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC, SGMatrix)
+#undef BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC
+
 /**
  * Wrapper method of in-place matrix elementwise product.
  *
@@ -301,6 +361,23 @@ namespace shogun
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_RANGE_FILL, SGMatrix)
 #undef BACKEND_GENERIC_RANGE_FILL
 
+/**
+ * Wrapper method that solves a system of linear equations
+ * using QR decomposition.
+ *
+ * @see linalg::qr_solver
+ */
+#define BACKEND_GENERIC_QR_SOLVER(Type, Container)                             \
+	virtual Container<Type> qr_solver(                                         \
+	    const SGMatrix<Type>& A, const Container<Type>& b) const               \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+		return 0;                                                              \
+	}
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_QR_SOLVER, SGVector)
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_QR_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_QR_SOLVER
+
 /**
  * Wrapper method of scale operation the operation result = alpha*A.
  *
@@ -449,6 +526,21 @@ namespace shogun
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_BLOCK_ROWWISE_SUM, SGMatrix)
 #undef BACKEND_GENERIC_BLOCK_ROWWISE_SUM
 
+/**
+ * Wrapper method of svd computation.
+ *
+ * @see linalg::svd
+ */
+#define BACKEND_GENERIC_SVD(Type, Container)                                   \
+	virtual void svd(                                                          \
+	    const Container<Type>& A, SGVector<Type> s, SGMatrix<Type> U,          \
+	    bool thin_U) const                                                     \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+	}
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_SVD, SGMatrix)
+#undef BACKEND_GENERIC_SVD
+
 /**
  * Wrapper method of trace computation.
  *
@@ -463,6 +555,39 @@ namespace shogun
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_TRACE, SGMatrix)
 #undef BACKEND_GENERIC_TRACE
 
+/**
+ * Wrapper method of trace computation.
+ *
+ * @see linalg::transpose_matrix
+ */
+#define BACKEND_GENERIC_TRANSPOSE_MATRIX(Type, Container)                      \
+	virtual Container<Type> transpose_matrix(const Container<Type>& A) const   \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+		return 0;                                                              \
+	}
+		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_TRANSPOSE_MATRIX, SGMatrix)
+#undef BACKEND_GENERIC_TRANSPOSE_MATRIX
+
+/**
+ * Wrapper method of triangular solver.
+ *
+ * @see linalg::triangular_solver
+ */
+#define BACKEND_GENERIC_TRIANGULAR_SOLVER(Type, Container)                     \
+	virtual Container<Type> triangular_solver(                                 \
+	    const SGMatrix<Type>& L, const Container<Type>& b,                     \
+	    const bool lower = true) const                                         \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+		return 0;                                                              \
+	}
+		DEFINE_FOR_NON_INTEGER_PTYPE(
+		    BACKEND_GENERIC_TRIANGULAR_SOLVER, SGVector)
+		DEFINE_FOR_NON_INTEGER_PTYPE(
+		    BACKEND_GENERIC_TRIANGULAR_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_TRIANGULAR_SOLVER
+
 /**
  * Wrapper method of set vector or matrix to zero.
  *
diff --git a/src/shogun/mathematics/linalg/LinalgBackendEigen.h b/src/shogun/mathematics/linalg/LinalgBackendEigen.h
index 519ec9fe463..875d7d13910 100644
--- a/src/shogun/mathematics/linalg/LinalgBackendEigen.h
+++ b/src/shogun/mathematics/linalg/LinalgBackendEigen.h
@@ -63,6 +63,19 @@ namespace shogun
 		DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_COL_VEC, SGMatrix)
 #undef BACKEND_GENERIC_ADD_COL_VEC
 
+/** Implementation of @see LinalgBackendBase::add_scalar */
+#define BACKEND_GENERIC_ADD_SCALAR(Type, Container)                            \
+	virtual void add_scalar(Container<Type>& a, Type b) const;
+		DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_SCALAR, SGVector)
+		DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_SCALAR, SGMatrix)
+#undef BACKEND_GENERIC_ADD_SCALAR
+
+/** Implementation of @see LinalgBackendBase::center_matrix */
+#define BACKEND_GENERIC_CENTER_MATRIX(Type, Container)                         \
+	virtual void center_matrix(Container<Type>& A) const;
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CENTER_MATRIX, SGMatrix)
+#undef BACKEND_GENERIC_CENTER_MATRIX
+
 /** Implementation of @see LinalgBackendBase::cholesky_factor */
 #define BACKEND_GENERIC_CHOLESKY_FACTOR(Type, Container)                       \
 	virtual Container<Type> cholesky_factor(                                   \
@@ -84,6 +97,23 @@ namespace shogun
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_DOT, SGVector)
 #undef BACKEND_GENERIC_DOT
 
+/** Implementation of @see LinalgBackendBase::eigen_solver */
+#define BACKEND_GENERIC_EIGEN_SOLVER(Type, Container)                          \
+	virtual void eigen_solver(                                                 \
+	    const Container<Type>& A, SGVector<Type>& eigenvalues,                 \
+	    SGMatrix<Type>& eigenvectors) const;
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_EIGEN_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_EIGEN_SOLVER
+
+/** Implementation of @see LinalgBackendBase::eigen_solver_symmetric */
+#define BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC(Type, Container)                \
+	virtual void eigen_solver_symmetric(                                       \
+	    const Container<Type>& A, SGVector<Type>& eigenvalues,                 \
+	    SGMatrix<Type>& eigenvectors) const;
+		DEFINE_FOR_NON_INTEGER_PTYPE(
+		    BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC, SGMatrix)
+#undef BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC
+
 /** Implementation of @see LinalgBackendBase::element_prod */
 #define BACKEND_GENERIC_IN_PLACE_ELEMENT_PROD(Type, Container)                 \
 	virtual void element_prod(                                                 \
@@ -143,6 +173,14 @@ namespace shogun
 		BACKEND_GENERIC_COMPLEX_MEAN(SGMatrix)
 #undef BACKEND_GENERIC_COMPLEX_MEAN
 
+/** Implementation of @see LinalgBackendBase::qr_solver */
+#define BACKEND_GENERIC_QR_SOLVER(Type, Container)                             \
+	virtual Container<Type> qr_solver(                                         \
+	    const SGMatrix<Type>& A, const Container<Type>& b) const;
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_QR_SOLVER, SGVector)
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_QR_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_QR_SOLVER
+
 /** Implementation of @see LinalgBackendBase::range_fill */
 #define BACKEND_GENERIC_RANGE_FILL(Type, Container)                            \
 	virtual void range_fill(Container<Type>& a, const Type start) const;
@@ -220,12 +258,37 @@ namespace shogun
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_BLOCK_ROWWISE_SUM, SGMatrix)
 #undef BACKEND_GENERIC_BLOCK_ROWWISE_SUM
 
+/** Implementation of @see LinalgBackendBase::svd */
+#define BACKEND_GENERIC_SVD(Type, Container)                                   \
+	virtual void svd(                                                          \
+	    const Container<Type>& A, SGVector<Type> s, Container<Type> U,         \
+	    bool thin_U) const;
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_SVD, SGMatrix)
+#undef BACKEND_GENERIC_SVD
+
 /** Implementation of @see LinalgBackendBase::trace */
 #define BACKEND_GENERIC_TRACE(Type, Container)                                 \
 	virtual Type trace(const Container<Type>& A) const;
 		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_TRACE, SGMatrix)
 #undef BACKEND_GENERIC_TRACE
 
+/** Implementation of @see LinalgBackendBase::transpose_matrix */
+#define BACKEND_GENERIC_TRANSPOSE_MATRIX(Type, Container)                      \
+	virtual Container<Type> transpose_matrix(const Container<Type>& A) const;
+		DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_TRANSPOSE_MATRIX, SGMatrix)
+#undef BACKEND_GENERIC_TRANSPOSE_MATRIX
+
+/** Implementation of @see LinalgBackendBase::triangular_solver */
+#define BACKEND_GENERIC_TRIANGULAR_SOLVER(Type, Container)                     \
+	virtual Container<Type> triangular_solver(                                 \
+	    const SGMatrix<Type>& L, const Container<Type>& b, const bool lower)   \
+	    const;
+		DEFINE_FOR_NON_INTEGER_PTYPE(
+		    BACKEND_GENERIC_TRIANGULAR_SOLVER, SGVector)
+		DEFINE_FOR_NON_INTEGER_PTYPE(
+		    BACKEND_GENERIC_TRIANGULAR_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_TRIANGULAR_SOLVER
+
 /** Implementation of @see LinalgBackendBase::zero */
 #define BACKEND_GENERIC_ZERO(Type, Container)                                  \
 	virtual void zero(Container<Type>& a) const;
@@ -263,6 +326,18 @@ namespace shogun
 		    const SGMatrix<T>& A, index_t i, const SGVector<T>& b,
 		    SGVector<T>& result, T alpha, T beta) const;
 
+		/** Eigen3 vector add scalar method */
+		template <typename T>
+		void add_scalar_impl(SGVector<T>& a, T b) const;
+
+		/** Eigen3 matrix add scalar method */
+		template <typename T>
+		void add_scalar_impl(SGMatrix<T>& a, T b) const;
+
+		/** Eigen3 center matrix method */
+		template <typename T>
+		void center_matrix_impl(SGMatrix<T>& A) const;
+
 		/** Eigen3 Cholesky decomposition */
 		template <typename T>
 		SGMatrix<T>
@@ -277,6 +352,27 @@ namespace shogun
 		template <typename T>
 		T dot_impl(const SGVector<T>& a, const SGVector<T>& b) const;
 
+		/** Eigen3 eigenvalues and eigenvectors computation for real matrices.
+		 */
+		template <typename T>
+		void eigen_solver_impl(
+		    const SGMatrix<T>& A, SGVector<T>& eigenvalues,
+		    SGMatrix<T>& eigenvectors) const;
+
+		/** Eigen3 eigenvalues and eigenvectors computation for complex
+		 * matrices. */
+		void eigen_solver_impl(
+		    const SGMatrix<complex128_t>& A,
+		    SGVector<complex128_t>& eigenvalues,
+		    SGMatrix<complex128_t>& eigenvectors) const;
+
+		/** Eigen3 eigenvalues and eigenvectors computation of symmetric
+		 * matrices */
+		template <typename T>
+		void eigen_solver_symmetric_impl(
+		    const SGMatrix<T>& A, SGVector<T>& eigenvalues,
+		    SGMatrix<T>& eigenvectors) const;
+
 		/** Eigen3 matrix in-place elementwise product method */
 		template <typename T>
 		void element_prod_impl(
@@ -326,6 +422,16 @@ namespace shogun
 		template <template <typename> class Container>
 		complex128_t mean_impl(const Container<complex128_t>& a) const;
 
+		/** Eigen3 vector QR solver. */
+		template <typename T>
+		SGVector<T>
+		qr_solver_impl(const SGMatrix<T>& A, const SGVector<T>& b) const;
+
+		/** Eigen3 matrix QR solver. */
+		template <typename T>
+		SGMatrix<T>
+		qr_solver_impl(const SGMatrix<T>& A, const SGMatrix<T> b) const;
+
 		/** Range fill a vector or matrix with start...start+len-1. */
 		template <typename T, template <typename> class Container>
 		void range_fill_impl(Container<T>& a, const T start) const;
@@ -385,10 +491,30 @@ namespace shogun
 		SGVector<T> rowwise_sum_impl(
 		    const linalg::Block<SGMatrix<T>>& mat, bool no_diag) const;
 
+		/** Eigen3 compute svd method */
+		template <typename T>
+		void svd_impl(
+		    const SGMatrix<T>& A, SGVector<T>& s, SGMatrix<T>& U,
+		    bool thin_U) const;
+
 		/** Eigen3 compute trace method */
 		template <typename T>
 		T trace_impl(const SGMatrix<T>& A) const;
 
+		/** Eigen3 transpose matrix method */
+		template <typename T>
+		SGMatrix<T> transpose_matrix_impl(const SGMatrix<T>& A) const;
+
+		/** Eigen3 triangular solver method */
+		template <typename T>
+		SGMatrix<T> triangular_solver_impl(
+		    const SGMatrix<T>& L, const SGMatrix<T>& b, const bool lower) const;
+
+		/** Eigen3 triangular solver method */
+		template <typename T>
+		SGVector<T> triangular_solver_impl(
+		    const SGMatrix<T>& L, const SGVector<T>& b, const bool lower) const;
+
 		/** Eigen3 set vector to zero method */
 		template <typename T>
 		void zero_impl(SGVector<T>& a) const;
diff --git a/src/shogun/mathematics/linalg/LinalgNamespace.h b/src/shogun/mathematics/linalg/LinalgNamespace.h
index a55738b6f65..97c8625f2c7 100644
--- a/src/shogun/mathematics/linalg/LinalgNamespace.h
+++ b/src/shogun/mathematics/linalg/LinalgNamespace.h
@@ -443,6 +443,37 @@ namespace shogun
 			    ->add_col_vec(A, i, b, result, alpha, beta);
 		}
 
+		/**
+		 * Adds a scalar to each element of a vector or a matrix in-place.
+		 *
+		 * @param a Vector or matrix
+		 * @param b Scalar to be added
+		 */
+		template <typename T, template <typename> class Container>
+		void add_scalar(Container<T>& a, T b)
+		{
+			infer_backend(a)->add_scalar(a, b);
+		}
+
+		/**
+		 * Centers a square matrix in-place, i.e. removes column/row mean from
+		 * columns/rows.
+		 * In particular it computes A = A - 1N*A - A*1N + 1N*A*1N
+		 * where 1N denotes the matrix of the same size as A for which each
+		 * element
+		 * takes value 1/n, where n is the number of columns and rows of A.
+		 *
+		 * @param A The matrix to be centered
+		 */
+		template <typename T>
+		void center_matrix(SGMatrix<T>& A)
+		{
+			REQUIRE(
+			    A.num_rows == A.num_cols, "Matrix A (%d x% d) is not square!\n",
+			    A.num_rows, A.num_cols);
+			infer_backend(A)->center_matrix(A);
+		}
+
 		/**
 		 * Performs the operation A = alpha * x * y' + A
 		 *
@@ -520,6 +551,84 @@ namespace shogun
 			return infer_backend(a, b)->dot(a, b);
 		}
 
+		/**
+		 * Compute the eigenvalues and eigenvectors of a matrix.
+		 * Note that the type of the computed values is the same
+		 * as the matrix's type, i.e. for real matrices it returns
+		 * only the real part of the eigenvalues/vectors.
+		 *
+		 * User should pass an appropriately pre-allocated memory vector
+		 * to store the eigenvalues and an appropriately pre-allocated memory
+		 * matrix to store the eigenvectors.
+		 *
+		 * @param A The matrix whose eigenvalues and eigenvectors are to be
+		 * computed
+		 * @param eigenvalues Eigenvalues result vector
+		 * @param eigenvectors Eigenvectors result matrix
+		 */
+		template <typename T>
+		void eigen_solver(
+		    const SGMatrix<T>& A, SGVector<T>& eigenvalues,
+		    SGMatrix<T>& eigenvectors)
+		{
+			REQUIRE(
+			    A.num_rows == A.num_cols, "Matrix A (%d x% d) is not square!\n",
+			    A.num_rows, A.num_cols);
+			REQUIRE(
+			    A.num_rows == eigenvectors.num_rows,
+			    "Number of rows of A (%d) doesn't match eigenvectors' matrix "
+			    "(%d).\n",
+			    A.num_rows, eigenvectors.num_rows);
+			REQUIRE(
+			    A.num_cols == eigenvectors.num_cols,
+			    "Number of columns of A (%d) doesn't match eigenvectors' "
+			    "matrix (%d).\n",
+			    A.num_cols, eigenvectors.num_cols);
+			REQUIRE(
+			    A.num_cols == eigenvalues.vlen,
+			    "Length of eigenvalues' vector doesn't match matrix A");
+
+			infer_backend(A)->eigen_solver(A, eigenvalues, eigenvectors);
+		}
+
+		/**
+		 * Compute the eigenvalues and eigenvectors of a symmetric matrix.
+		 *
+		 * User should pass an appropriately pre-allocated memory vector
+		 * to store the eigenvalues and an appropriately pre-allocated memory
+		 * matrix to store the eigenvectors.
+		 *
+		 * @param A The matrix whose eigenvalues and eigenvectors are to be
+		 * computed
+		 * @param eigenvalues Eigenvalues result vector
+		 * @param eigenvectors Eigenvectors result matrix
+		 */
+		template <typename T>
+		void eigen_solver_symmetric(
+		    const SGMatrix<T>& A, SGVector<T>& eigenvalues,
+		    SGMatrix<T>& eigenvectors)
+		{
+			REQUIRE(
+			    A.num_rows == A.num_cols, "Matrix A (%d x% d) is not square!\n",
+			    A.num_rows, A.num_cols);
+			REQUIRE(
+			    A.num_rows == eigenvectors.num_rows,
+			    "Number of rows of A (%d) doesn't match eigenvectors' matrix "
+			    "(%d).\n",
+			    A.num_rows, eigenvectors.num_rows);
+			REQUIRE(
+			    A.num_cols == eigenvectors.num_cols,
+			    "Number of columns of A (%d) doesn't match eigenvectors' "
+			    "matrix (%d).\n",
+			    A.num_cols, eigenvectors.num_cols);
+			REQUIRE(
+			    A.num_cols == eigenvalues.vlen,
+			    "Length of eigenvalues' vector doesn't match matrix A");
+
+			infer_backend(A)->eigen_solver_symmetric(
+			    A, eigenvalues, eigenvectors);
+		}
+
 		/** Performs the operation C = A .* B where ".*" denotes elementwise
 		 * multiplication
 		 * on matrix blocks.
@@ -1006,6 +1115,24 @@ namespace shogun
 			return infer_backend(a)->mean(a);
 		}
 
+		/**
+		 * Solve the linear equations \f$Ax=b\f$ through the
+		 * QR decomposition of A.
+		 *
+		 * @param A The matrix
+		 * @param b Right-hand side vector or matrix
+		 * @return \f$\x\f$
+		 */
+		template <typename T, template <typename> class Container>
+		Container<T> qr_solver(const SGMatrix<T>& A, const Container<T>& b)
+		{
+			REQUIRE(
+			    A.num_rows == A.num_cols, "Matrix A (%d x% d) is not square!\n",
+			    A.num_rows, A.num_cols);
+
+			return infer_backend(A, SGMatrix<T>(b))->qr_solver(A, b);
+		}
+
 		/**
 		 * Range fill a vector or matrix with start...start+len-1
 		 *
@@ -1222,6 +1349,55 @@ namespace shogun
 			return sg_linalg->get_cpu_backend()->rowwise_sum(a, no_diag);
 		}
 
+		/**
+		 * Compute the singular value decomposition \f$A = U S V^{*}\f$ of a
+		 * matrix.
+		 * Given the \f$m \times n\f$ matrix A with \f$r = min(m,n)\f$, user
+		 * should pass
+		 * an appropriately pre-allocated vector of length r and a pre-allocated
+		 * matrix of size \f$m \times r\f$ for thin U or \f$m \times m\f$
+		 * otherwise.
+		 *
+		 * @param A The matrix whose svd is to be computed
+		 * @param s The vector that stores the resulting singular values
+		 * @param U The matrix that stores the resulting unitary matrix U
+		 * @param thin_U Whether to compute the full or thin matrix U
+		 * (default:thin)
+		 */
+		template <typename T>
+		void
+		svd(const SGMatrix<T>& A, SGVector<T>& s, SGMatrix<T>& U,
+		    bool thin_U = true)
+		{
+			auto r = CMath::min(A.num_cols, A.num_rows);
+			REQUIRE(
+			    (A.num_rows == U.num_rows),
+			    "Number of rows of matrix A (%d) must match matrix U (%d).\n",
+			    A.num_rows, U.num_rows);
+			if (thin_U)
+			{
+				REQUIRE(
+				    (U.num_cols == r), "Number of columns of matrix A (%d) "
+				                       "must match A's smaller dimension "
+				                       "(%d).\n",
+				    A.num_cols, r);
+			}
+			else
+			{
+				REQUIRE(
+				    (A.num_rows == U.num_cols), "Number of rows of matrix A "
+				                                "(%d) must match number of "
+				                                "columns of U (%d).\n",
+				    A.num_cols, r);
+			}
+			REQUIRE(
+			    (s.vlen == r), "Length of vector s (%d) doesn't match A's "
+			                   "smaller dimension (%d).\n",
+			    s.vlen, r);
+
+			infer_backend(A)->svd(A, s, U, thin_U);
+		}
+
 		/**
 		 * Method that computes the trace of square matrix.
 		 *
@@ -1235,6 +1411,36 @@ namespace shogun
 			return infer_backend(A)->trace(A);
 		}
 
+		/**
+		 * Method that computes the transpose of a matrix.
+		 *
+		 * @param A The matrix
+		 * @return The transposed matrix
+		 */
+		template <typename T>
+		SGMatrix<T> transpose_matrix(const SGMatrix<T>& A)
+		{
+			return infer_backend(A)->transpose_matrix(A);
+		}
+
+		/**
+		 * Solve the linear equations \f$Lx=b\f$,
+		 * where \f$L\f$ is a triangular matrix.
+		 *
+		 * @param L Triangular matrix
+		 * @param b Right-hand side array
+		 * @param lower Whether L is upper or lower triangular (default:lower)
+		 * @return \f$\x\f$
+		 */
+		template <typename T, template <typename> class Container>
+		Container<T> triangular_solver(
+		    const SGMatrix<T>& L, const Container<T>& b,
+		    const bool lower = true)
+		{
+			return infer_backend(L, SGMatrix<T>(b))
+			    ->triangular_solver(L, b, lower);
+		}
+
 		/**
 		 * Method that fills with zero a vector or a matrix.
 		 *
diff --git a/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp b/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp
index dc4559ac533..1d8ae3fb4d1 100644
--- a/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp
+++ b/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp
@@ -57,6 +57,15 @@ DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_COL_VEC, SGVector)
 DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_COL_VEC, SGMatrix)
 #undef BACKEND_GENERIC_ADD_COL_VEC
 
+#define BACKEND_GENERIC_ADD_SCALAR(Type, Container)                            \
+	void LinalgBackendEigen::add_scalar(Container<Type>& a, Type b) const      \
+	{                                                                          \
+		add_scalar_impl(a, b);                                                 \
+	}
+DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_SCALAR, SGVector)
+DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_SCALAR, SGMatrix)
+#undef BACKEND_GENERIC_ADD_SCALAR
+
 #define BACKEND_GENERIC_DOT(Type, Container)                                   \
 	Type LinalgBackendEigen::dot(                                              \
 	    const Container<Type>& a, const Container<Type>& b) const              \
@@ -157,6 +166,20 @@ void LinalgBackendEigen::add_col_vec_impl(
 	result_eig = alpha * A_eig.col(i) + beta * b_eig;
 }
 
+template <typename T>
+void LinalgBackendEigen::add_scalar_impl(SGVector<T>& a, T b) const
+{
+	typename SGVector<T>::EigenVectorXtMap a_eig = a;
+	a_eig = a_eig.array() + b;
+}
+
+template <typename T>
+void LinalgBackendEigen::add_scalar_impl(SGMatrix<T>& a, T b) const
+{
+	typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
+	a_eig = a_eig.array() + b;
+}
+
 template <typename T>
 T LinalgBackendEigen::dot_impl(const SGVector<T>& a, const SGVector<T>& b) const
 {
diff --git a/src/shogun/mathematics/linalg/backend/eigen/Decompositions.cpp b/src/shogun/mathematics/linalg/backend/eigen/Decompositions.cpp
index 6edaa88464b..e57634ee4cd 100644
--- a/src/shogun/mathematics/linalg/backend/eigen/Decompositions.cpp
+++ b/src/shogun/mathematics/linalg/backend/eigen/Decompositions.cpp
@@ -44,6 +44,36 @@ using namespace shogun;
 DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_FACTOR, SGMatrix)
 #undef BACKEND_GENERIC_CHOLESKY_FACTOR
 
+#define BACKEND_GENERIC_EIGEN_SOLVER(Type, Container)                          \
+	void LinalgBackendEigen::eigen_solver(                                     \
+	    const Container<Type>& A, SGVector<Type>& eigenvalues,                 \
+	    SGMatrix<Type>& eigenvectors) const                                    \
+	{                                                                          \
+		eigen_solver_impl(A, eigenvalues, eigenvectors);                       \
+	}
+DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_EIGEN_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_EIGEN_SOLVER
+
+#define BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC(Type, Container)                \
+	void LinalgBackendEigen::eigen_solver_symmetric(                           \
+	    const Container<Type>& A, SGVector<Type>& eigenvalues,                 \
+	    SGMatrix<Type>& eigenvectors) const                                    \
+	{                                                                          \
+		eigen_solver_symmetric_impl(A, eigenvalues, eigenvectors);             \
+	}
+DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC, SGMatrix)
+#undef BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC
+
+#define BACKEND_GENERIC_SVD(Type, Container)                                   \
+	void LinalgBackendEigen::svd(                                              \
+	    const Container<Type>& A, SGVector<Type> s, Container<Type> U,         \
+	    bool thin_U) const                                                     \
+	{                                                                          \
+		return svd_impl(A, s, U, thin_U);                                      \
+	}
+DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_SVD, SGMatrix)
+#undef BACKEND_GENERIC_SVD
+
 #undef DEFINE_FOR_ALL_PTYPE
 #undef DEFINE_FOR_REAL_PTYPE
 #undef DEFINE_FOR_NON_INTEGER_PTYPE
@@ -79,3 +109,62 @@ SGMatrix<T> LinalgBackendEigen::cholesky_factor_impl(
 
 	return c;
 }
+
+template <typename T>
+void LinalgBackendEigen::eigen_solver_impl(
+    const SGMatrix<T>& A, SGVector<T>& eigenvalues,
+    SGMatrix<T>& eigenvectors) const
+{
+	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
+	typename SGMatrix<T>::EigenMatrixXtMap eigenvectors_eig = eigenvectors;
+	typename SGVector<T>::EigenVectorXtMap eigenvalues_eig = eigenvalues;
+
+	Eigen::EigenSolver<typename SGMatrix<T>::EigenMatrixXt> solver(A_eig);
+	eigenvalues_eig = solver.eigenvalues().real();
+	eigenvectors_eig = solver.eigenvectors().real();
+}
+
+void LinalgBackendEigen::eigen_solver_impl(
+    const SGMatrix<complex128_t>& A, SGVector<complex128_t>& eigenvalues,
+    SGMatrix<complex128_t>& eigenvectors) const
+{
+	typename SGMatrix<complex128_t>::EigenMatrixXtMap A_eig = A;
+	typename SGMatrix<complex128_t>::EigenMatrixXtMap eigenvectors_eig =
+	    eigenvectors;
+	typename SGVector<complex128_t>::EigenVectorXtMap eigenvalues_eig =
+	    eigenvalues;
+
+	Eigen::ComplexEigenSolver<typename SGMatrix<complex128_t>::EigenMatrixXt>
+	    solver(A_eig);
+	eigenvalues_eig = solver.eigenvalues();
+	eigenvectors_eig = solver.eigenvectors();
+}
+
+template <typename T>
+void LinalgBackendEigen::eigen_solver_symmetric_impl(
+    const SGMatrix<T>& A, SGVector<T>& eigenvalues,
+    SGMatrix<T>& eigenvectors) const
+{
+	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
+	typename SGMatrix<T>::EigenMatrixXtMap eigenvectors_eig = eigenvectors;
+	typename SGVector<T>::EigenVectorXtMap eigenvalues_eig = eigenvalues;
+
+	Eigen::SelfAdjointEigenSolver<typename SGMatrix<T>::EigenMatrixXt> solver(
+	    A_eig);
+	eigenvalues_eig = solver.eigenvalues().template cast<T>();
+	eigenvectors_eig = solver.eigenvectors().template cast<T>();
+}
+
+template <typename T>
+void LinalgBackendEigen::svd_impl(
+    const SGMatrix<T>& A, SGVector<T>& s, SGMatrix<T>& U, bool thin_U) const
+{
+	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
+	typename SGVector<T>::EigenVectorXtMap s_eig = s;
+	typename SGMatrix<T>::EigenMatrixXtMap U_eig = U;
+
+	auto svd_eig =
+	    A_eig.jacobiSvd(thin_U ? Eigen::ComputeThinU : Eigen::ComputeFullU);
+	s_eig = svd_eig.singularValues().template cast<T>();
+	U_eig = svd_eig.matrixU().template cast<T>();
+}
\ No newline at end of file
diff --git a/src/shogun/mathematics/linalg/backend/eigen/Misc.cpp b/src/shogun/mathematics/linalg/backend/eigen/Misc.cpp
index eb64e3f978f..f392d9616fa 100644
--- a/src/shogun/mathematics/linalg/backend/eigen/Misc.cpp
+++ b/src/shogun/mathematics/linalg/backend/eigen/Misc.cpp
@@ -35,6 +35,14 @@
 
 using namespace shogun;
 
+#define BACKEND_GENERIC_CENTER_MATRIX(Type, Container)                         \
+	void LinalgBackendEigen::center_matrix(Container<Type>& A) const           \
+	{                                                                          \
+		center_matrix_impl(A);                                                 \
+	}
+DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CENTER_MATRIX, SGMatrix)
+#undef BACKEND_GENERIC_CENTER_MATRIX
+
 #define BACKEND_GENERIC_IDENTITY(Type, Container)                              \
 	void LinalgBackendEigen::identity(Container<Type>& I) const                \
 	{                                                                          \
@@ -63,6 +71,15 @@ DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_RANGE_FILL, SGVector)
 DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_RANGE_FILL, SGMatrix)
 #undef BACKEND_GENERIC_RANGE_FILL
 
+#define BACKEND_GENERIC_TRANSPOSE_MATRIX(Type, Container)                      \
+	Container<Type> LinalgBackendEigen::transpose_matrix(                      \
+	    const Container<Type>& A) const                                        \
+	{                                                                          \
+		return transpose_matrix_impl(A);                                       \
+	}
+DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_TRANSPOSE_MATRIX, SGMatrix)
+#undef BACKEND_GENERIC_TRANSPOSE_MATRIX
+
 #define BACKEND_GENERIC_ZERO(Type, Container)                                  \
 	void LinalgBackendEigen::zero(Container<Type>& a) const                    \
 	{                                                                          \
@@ -77,6 +94,22 @@ DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ZERO, SGMatrix)
 #undef DEFINE_FOR_NON_INTEGER_PTYPE
 #undef DEFINE_FOR_NUMERIC_PTYPE
 
+template <typename T>
+void LinalgBackendEigen::center_matrix_impl(SGMatrix<T>& A) const
+{
+	index_t n = A.num_cols;
+	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
+
+	typename SGVector<T>::EigenVectorXt rows_sum =
+	    (A_eig.rowwise().sum()).array() / (T)n;
+	typename SGVector<T>::EigenRowVectorXt cols_sum =
+	    (A_eig.colwise().sum()).array() / (T)n;
+	T m = rows_sum.sum() / (T)n;
+
+	A_eig = ((A_eig.array() + m).matrix().rowwise() - cols_sum).colwise() -
+	        rows_sum;
+}
+
 template <typename T>
 void LinalgBackendEigen::identity_impl(SGMatrix<T>& I) const
 {
@@ -98,6 +131,18 @@ void LinalgBackendEigen::set_const_impl(Container<T>& a, T value) const
 		a[i] = value;
 }
 
+template <typename T>
+SGMatrix<T>
+LinalgBackendEigen::transpose_matrix_impl(const SGMatrix<T>& A) const
+{
+	SGMatrix<T> tr(A.num_cols, A.num_rows);
+	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
+	typename SGMatrix<T>::EigenMatrixXtMap tr_eig = tr;
+
+	tr_eig = A_eig.transpose();
+	return tr;
+}
+
 template <typename T>
 void LinalgBackendEigen::zero_impl(SGVector<T>& a) const
 {
diff --git a/src/shogun/mathematics/linalg/backend/eigen/Solvers.cpp b/src/shogun/mathematics/linalg/backend/eigen/Solvers.cpp
index 1ba41cceec3..47d8b0f3de7 100644
--- a/src/shogun/mathematics/linalg/backend/eigen/Solvers.cpp
+++ b/src/shogun/mathematics/linalg/backend/eigen/Solvers.cpp
@@ -45,6 +45,27 @@ using namespace shogun;
 DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_SOLVER, SGMatrix)
 #undef BACKEND_GENERIC_CHOLESKY_SOLVER
 
+#define BACKEND_GENERIC_QR_SOLVER(Type, Container)                             \
+	Container<Type> LinalgBackendEigen::qr_solver(                             \
+	    const SGMatrix<Type>& A, const Container<Type>& b) const               \
+	{                                                                          \
+		return qr_solver_impl(A, b);                                           \
+	}
+DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_QR_SOLVER, SGVector)
+DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_QR_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_QR_SOLVER
+
+#define BACKEND_GENERIC_TRIANGULAR_SOLVER(Type, Container)                     \
+	Container<Type> LinalgBackendEigen::triangular_solver(                     \
+	    const SGMatrix<Type>& L, const Container<Type>& b, const bool lower)   \
+	    const                                                                  \
+	{                                                                          \
+		return triangular_solver_impl(L, b, lower);                            \
+	}
+DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_TRIANGULAR_SOLVER, SGVector)
+DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_TRIANGULAR_SOLVER, SGMatrix)
+#undef BACKEND_GENERIC_TRIANGULAR_SOLVER
+
 #undef DEFINE_FOR_ALL_PTYPE
 #undef DEFINE_FOR_REAL_PTYPE
 #undef DEFINE_FOR_NON_INTEGER_PTYPE
@@ -80,3 +101,89 @@ SGVector<T> LinalgBackendEigen::cholesky_solver_impl(
 
 	return x;
 }
+
+template <typename T>
+SGVector<T> LinalgBackendEigen::qr_solver_impl(
+    const SGMatrix<T>& A, const SGVector<T>& b) const
+{
+	SGVector<T> result(b.vlen);
+	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
+	typename SGVector<T>::EigenVectorXtMap b_eig = b;
+	typename SGVector<T>::EigenVectorXtMap result_eig = result;
+
+	result_eig = (A_eig.householderQr().solve(b_eig));
+
+	return result;
+}
+
+template <typename T>
+SGMatrix<T> LinalgBackendEigen::qr_solver_impl(
+    const SGMatrix<T>& A, const SGMatrix<T> b) const
+{
+	SGMatrix<T> result(b.num_rows, b.num_cols);
+	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
+	typename SGMatrix<T>::EigenMatrixXtMap b_eig = b;
+	typename SGMatrix<T>::EigenMatrixXtMap result_eig = result;
+
+	result_eig = (A_eig.householderQr().solve(b_eig));
+
+	return result;
+}
+
+template <typename T>
+SGMatrix<T> LinalgBackendEigen::triangular_solver_impl(
+    const SGMatrix<T>& L, const SGMatrix<T>& b, const bool lower) const
+{
+	SGMatrix<T> x(b.num_rows, b.num_cols);
+	typename SGMatrix<T>::EigenMatrixXtMap L_eig = L;
+	typename SGMatrix<T>::EigenMatrixXtMap b_eig = b;
+	typename SGMatrix<T>::EigenMatrixXtMap x_eig = x;
+
+	if (lower == false)
+	{
+		Eigen::TriangularView<Eigen::Map<typename SGMatrix<T>::EigenMatrixXt, 0,
+		                                 Eigen::Stride<0, 0>>,
+		                      Eigen::Upper>
+		    tlv(L_eig);
+		x_eig = tlv.solve(b_eig);
+	}
+	else
+	{
+		Eigen::TriangularView<Eigen::Map<typename SGMatrix<T>::EigenMatrixXt, 0,
+		                                 Eigen::Stride<0, 0>>,
+		                      Eigen::Lower>
+		    tlv(L_eig);
+		x_eig = tlv.solve(b_eig);
+	}
+
+	return x;
+}
+
+template <typename T>
+SGVector<T> LinalgBackendEigen::triangular_solver_impl(
+    const SGMatrix<T>& L, const SGVector<T>& b, const bool lower) const
+{
+	SGVector<T> x(b.size());
+	typename SGMatrix<T>::EigenMatrixXtMap L_eig = L;
+	typename SGVector<T>::EigenVectorXtMap b_eig = b;
+	typename SGVector<T>::EigenVectorXtMap x_eig = x;
+
+	if (lower == false)
+	{
+		Eigen::TriangularView<Eigen::Map<typename SGMatrix<T>::EigenMatrixXt, 0,
+		                                 Eigen::Stride<0, 0>>,
+		                      Eigen::Upper>
+		    tlv(L_eig);
+		x_eig = tlv.solve(b_eig);
+	}
+	else
+	{
+		Eigen::TriangularView<Eigen::Map<typename SGMatrix<T>::EigenMatrixXt, 0,
+		                                 Eigen::Stride<0, 0>>,
+		                      Eigen::Lower>
+		    tlv(L_eig);
+		x_eig = tlv.solve(b_eig);
+	}
+
+	return x;
+}
diff --git a/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc b/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc
index 4ea962f32e2..5ed92dcfee1 100644
--- a/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc
+++ b/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc
@@ -187,6 +187,56 @@ TEST(LinalgBackendEigen, SGMatrix_add_col_vec_in_place)
 		}
 }
 
+TEST(LinalgBackendEigen, SGVector_add_scalar)
+{
+	const index_t n = 4;
+	float64_t s = -0.3;
+
+	SGVector<float64_t> a(n);
+	for (index_t i = 0; i < (index_t)a.size(); ++i)
+		a[i] = i;
+	SGVector<float64_t> orig = a.clone();
+
+	add_scalar(a, s);
+
+	for (index_t i = 0; i < (index_t)a.size(); ++i)
+		EXPECT_NEAR(a[i], orig[i] + s, 1e-15);
+}
+
+TEST(LinalgBackendEigen, SGMatrix_add_scalar)
+{
+	const index_t r = 4, c = 3;
+	float64_t s = 0.4;
+
+	SGMatrix<float64_t> a(r, c);
+	for (index_t i = 0; i < (index_t)a.size(); ++i)
+		a[i] = i;
+	SGMatrix<float64_t> orig = a.clone();
+
+	add_scalar(a, s);
+
+	for (index_t i = 0; i < (index_t)a.size(); ++i)
+		EXPECT_NEAR(a[i], orig[i] + s, 1e-15);
+}
+
+TEST(LinalgBackendEigen, center_matrix)
+{
+	const index_t n = 3;
+	float64_t data[] = {0.8192343,  0.13191962, 0.50888604,
+	                    0.16857468, 0.24107738, 0.89455301,
+	                    0.40657379, 0.07902286, 0.24319651};
+	float64_t result[] = {0.25587541,  -0.1173183, -0.13855711,
+	                      -0.34283925, 0.04378442, 0.29905482,
+	                      0.08696383,  0.07353387, -0.16049771};
+
+	SGMatrix<float64_t> m(data, n, n, false);
+
+	center_matrix(m);
+
+	for (index_t i = 0; i < (index_t)m.size(); ++i)
+		EXPECT_NEAR(m[i], result[i], 1e-8);
+}
+
 TEST(LinalgBackendEigen, SGMatrix_cholesky_llt_lower)
 {
 	const index_t size=2;
@@ -266,6 +316,78 @@ TEST(LinalgBackendEigen, SGVector_dot)
 	EXPECT_NEAR(result, 5, 1E-15);
 }
 
+TEST(LinalgBackendEigen, eigensolver)
+{
+	const index_t n = 4;
+	float64_t data[] = {0.09987322, 0.80575314, 0.79068641, 0.69989667,
+	                    0.62323516, 0.16837367, 0.85027625, 0.60165948,
+	                    0.04898732, 0.96701123, 0.51683275, 0.51116495,
+	                    0.18277926, 0.6179262,  0.43745891, 0.63685464};
+	float64_t result_eigenvectors[] = {
+	    -0.63494074, 0.75831593,  -0.14014031, 0.04656076,
+	    0.82257205,  -0.28671857, -0.44196422, -0.21409185,
+	    -0.005932,   -0.20233724, -0.52285555, 0.82803776,
+	    -0.23930111, -0.56199714, -0.57298901, -0.54642272};
+	float64_t result_eigenvalues[] = {-0.6470538, -0.19125664, 0.16205101,
+	                                  2.0981937};
+
+	SGMatrix<float64_t> m(data, n, n, false);
+	SGMatrix<float64_t> eigenvectors(n, n);
+	SGVector<float64_t> eigenvalues(n);
+
+	eigen_solver(m, eigenvalues, eigenvectors);
+
+	auto args = CMath::argsort(eigenvalues);
+	for (index_t i = 0; i < n; ++i)
+	{
+		index_t idx = args[i];
+		EXPECT_NEAR(eigenvalues[idx], result_eigenvalues[i], 1e-7);
+
+		auto s =
+		    CMath::sign(eigenvectors[idx * n] * result_eigenvectors[i * n]);
+		for (index_t j = 0; j < n; ++j)
+			EXPECT_NEAR(
+			    eigenvectors[idx * n + j], s * result_eigenvectors[i * n + j],
+			    1e-7);
+	}
+}
+
+TEST(LinalgBackendEigen, eigensolver_symmetric)
+{
+	const index_t n = 4;
+	float64_t data[] = {0.09987322, 0.80575314, 0.04898732, 0.69989667,
+	                    0.80575314, 0.16837367, 0.96701123, 0.6179262,
+	                    0.04898732, 0.96701123, 0.51683275, 0.43745891,
+	                    0.69989667, 0.6179262,  0.43745891, 0.63685464};
+	float64_t result_eigenvectors[] = {
+	    -0.54618542, 0.69935447,  -0.45219663, 0.09001671,
+	    -0.56171388, -0.41397154, 0.17642953,  0.69424612,
+	    -0.46818396, 0.16780603,  0.73247599,  -0.46489119,
+	    0.40861077,  0.55800718,  0.47735703,  0.542029037};
+	float64_t result_eigenvalues[] = {-1.00663298, -0.18672196, 0.42940933,
+	                                  2.18587989};
+
+	SGMatrix<float64_t> m(data, n, n, false);
+	SGMatrix<float64_t> eigenvectors(n, n);
+	SGVector<float64_t> eigenvalues(n);
+
+	eigen_solver(m, eigenvalues, eigenvectors);
+
+	auto args = CMath::argsort(eigenvalues);
+	for (index_t i = 0; i < n; ++i)
+	{
+		index_t idx = args[i];
+		EXPECT_NEAR(eigenvalues[idx], result_eigenvalues[i], 1e-7);
+
+		auto s =
+		    CMath::sign(eigenvectors[idx * n] * result_eigenvectors[i * n]);
+		for (index_t j = 0; j < n; ++j)
+			EXPECT_NEAR(
+			    eigenvectors[idx * n + j], s * result_eigenvectors[i * n + j],
+			    1e-7);
+	}
+}
+
 TEST(LinalgBackendEigen, SGMatrix_elementwise_product)
 {
 	const index_t nrows = 3;
@@ -678,6 +800,44 @@ TEST(LinalgBackendEigen, SGMatrix_mean)
 	EXPECT_NEAR(result, 2.5, 1E-15);
 }
 
+TEST(LinalgBackendEigen, SGVector_qr_solver)
+{
+	const index_t n = 3;
+	float64_t data_A[] = {0.02800922, 0.99326012, 0.15204902,
+	                      0.30492837, 0.39708534, 0.40466969,
+	                      0.36415317, 0.04407589, 0.9095746};
+	float64_t data_b[] = {0.39461571, 0.6816856, 0.43323709};
+	float64_t result[] = {0.07135206, 1.56393127, -0.23141312};
+
+	SGMatrix<float64_t> A(data_A, n, n, false);
+	SGVector<float64_t> b(data_b, n, false);
+
+	auto x = qr_solver(A, b);
+
+	for (index_t i = 0; i < x.size(); ++i)
+		EXPECT_NEAR(x[i], result[i], 1E-7);
+}
+
+TEST(LinalgBackendEigen, SGMatrix_qr_solver)
+{
+	const index_t n = 3, m = 2;
+	float64_t data_A[] = {0.02800922, 0.99326012, 0.15204902,
+	                      0.30492837, 0.39708534, 0.40466969,
+	                      0.36415317, 0.04407589, 0.9095746};
+	float64_t data_B[] = {0.76775073, 0.88471312, 0.34795225,
+	                      0.94311546, 0.59630347, 0.65820143};
+	float64_t result[] = {-0.73834587, 4.22750496, -1.37484721,
+	                      -1.14718091, 4.49142548, -1.08282992};
+
+	SGMatrix<float64_t> A(data_A, n, n, false);
+	SGMatrix<float64_t> B(data_B, n, m, false);
+
+	auto X = qr_solver(A, B);
+
+	for (index_t i = 0; i < (index_t)X.size(); ++i)
+		EXPECT_NEAR(X[i], result[i], 1E-7);
+}
+
 TEST(LinalgBackendEigen, SGVector_range_fill)
 {
 	const index_t size = 5;
@@ -1035,6 +1195,66 @@ TEST(LinalgBackendEigen, SGMatrix_block_rowwise_sum)
 	}
 }
 
+TEST(LinalgBackendEigen, SGMatrix_svd_thinU)
+{
+	const index_t m = 5, n = 3;
+	float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194,
+	                    0.30786772, 0.25503552, 0.34367041, 0.66491478,
+	                    0.20488809, 0.5734351,  0.87179189, 0.07139643,
+	                    0.28540373, 0.06264684, 0.56204061};
+	float64_t result_s[] = {1.75382524, 0.56351367, 0.41124883};
+	float64_t result_U[] = {-0.60700926, -0.16647013, -0.56501385, -0.26696629,
+	                        -0.46186125, -0.69145782, 0.29548428,  0.5718984,
+	                        0.31771648,  -0.08101592, -0.27461424, 0.37170223,
+	                        -0.12681555, -0.53830325, 0.69323293};
+
+	SGMatrix<float64_t> A(data, m, n, false);
+	SGMatrix<float64_t> U(m, n);
+	SGVector<float64_t> s(n);
+
+	svd(A, s, U, true);
+
+	for (index_t i = 0; i < n; ++i)
+	{
+		auto c = CMath::sign(U[i * m] * result_U[i * m]);
+		for (index_t j = 0; j < m; ++j)
+			EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-7);
+	}
+	for (index_t i = 0; i < (index_t)s.size(); ++i)
+		EXPECT_NEAR(s[i], result_s[i], 1e-7);
+}
+
+TEST(LinalgBackendEigen, SGMatrix_svd_fullU)
+{
+	const index_t m = 5, n = 3;
+	float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194,
+	                    0.30786772, 0.25503552, 0.34367041, 0.66491478,
+	                    0.20488809, 0.5734351,  0.87179189, 0.07139643,
+	                    0.28540373, 0.06264684, 0.56204061};
+	float64_t result_s[] = {1.75382524, 0.56351367, 0.41124883};
+	float64_t result_U[] = {
+	    -0.60700926, -0.16647013, -0.56501385, -0.26696629, -0.46186125,
+	    -0.69145782, 0.29548428,  0.5718984,   0.31771648,  -0.08101592,
+	    -0.27461424, 0.37170223,  -0.12681555, -0.53830325, 0.69323293,
+	    -0.27809756, -0.68975171, -0.11662812, 0.38274703,  0.53554354,
+	    0.025973184, 0.520631112, -0.56921636, 0.62571522,  0.11287970};
+
+	SGMatrix<float64_t> A(data, m, n, false);
+	SGMatrix<float64_t> U(m, m);
+	SGVector<float64_t> s(n);
+
+	svd(A, s, U, false);
+
+	for (index_t i = 0; i < n; ++i)
+	{
+		auto c = CMath::sign(U[i * m] * result_U[i * m]);
+		for (index_t j = 0; j < m; ++j)
+			EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-7);
+	}
+	for (index_t i = 0; i < (index_t)s.size(); ++i)
+		EXPECT_NEAR(s[i], result_s[i], 1e-7);
+}
+
 TEST(LinalgBackendEigen, SGMatrix_trace)
 {
 	const index_t n = 4;
@@ -1050,6 +1270,99 @@ TEST(LinalgBackendEigen, SGMatrix_trace)
 	EXPECT_NEAR(trace(A), tr, 1e-15);
 }
 
+TEST(LinalgBackendEigen, SGMatrix_transpose_matrix)
+{
+	const index_t m = 5, n = 3;
+	float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194,
+	                    0.30786772, 0.25503552, 0.34367041, 0.66491478,
+	                    0.20488809, 0.5734351,  0.87179189, 0.07139643,
+	                    0.28540373, 0.06264684, 0.56204061};
+
+	SGMatrix<float64_t> A(data, m, n, false);
+
+	auto T = transpose_matrix(A);
+
+	for (index_t i = 0; i < m; ++i)
+		for (index_t j = 0; j < n; ++j)
+			EXPECT_NEAR(A.get_element(i, j), T.get_element(j, i), 1e-15);
+}
+
+TEST(LinalgBackendEigen, SGVector_triangular_solver_lower)
+{
+	const index_t n = 3;
+	float64_t data_L[] = {-0.92947874, -1.1432887,  -0.87119086,
+	                      0.,          -0.27048649, -0.05919915,
+	                      0.,          0.,          0.11869106};
+	float64_t data_b[] = {0.39461571, 0.6816856, 0.43323709};
+	float64_t result[] = {-0.42455592, -0.72571316, 0.17192745};
+
+	SGMatrix<float64_t> L(data_L, n, n, false);
+	SGVector<float64_t> b(data_b, n, false);
+
+	auto x = triangular_solver(L, b, true);
+
+	for (index_t i = 0; i < (index_t)x.size(); ++i)
+		EXPECT_NEAR(x[i], result[i], 1E-6);
+}
+
+TEST(LinalgBackendEigen, SGVector_triangular_solver_upper)
+{
+	const index_t n = 3;
+	float64_t data_U[] = {-0.92947874, 0.,          0.,
+	                      -1.1432887,  -0.27048649, 0.,
+	                      -0.87119086, -0.05919915, 0.11869106};
+	float64_t data_b[] = {0.39461571, 0.6816856, 0.43323709};
+	float64_t result[] = {0.23681135, -3.31909306, 3.65012412};
+
+	SGMatrix<float64_t> U(data_U, n, n, false);
+	SGVector<float64_t> b(data_b, n, false);
+
+	auto x = triangular_solver(U, b, false);
+
+	for (index_t i = 0; i < (index_t)x.size(); ++i)
+		EXPECT_NEAR(x[i], result[i], 1E-6);
+}
+
+TEST(LinalgBackendEigen, SGMatrix_triangular_solver_lower)
+{
+	const index_t n = 3, m = 2;
+	float64_t data_L[] = {-0.92947874, -1.1432887,  -0.87119086,
+	                      0.,          -0.27048649, -0.05919915,
+	                      0.,          0.,          0.11869106};
+	float64_t data_B[] = {0.76775073, 0.88471312, 0.34795225,
+	                      0.94311546, 0.59630347, 0.65820143};
+	float64_t result[] = {-0.82600139, 0.22050986, -3.02127745,
+	                      -1.01467136, 2.08424024, -0.86262387};
+
+	SGMatrix<float64_t> L(data_L, n, n, false);
+	SGMatrix<float64_t> B(data_B, n, m, false);
+
+	auto X = triangular_solver(L, B, true);
+
+	for (index_t i = 0; i < (index_t)X.size(); ++i)
+		EXPECT_NEAR(X[i], result[i], 1E-6);
+}
+
+TEST(LinalgBackendEigen, SGMatrix_triangular_solver_upper)
+{
+	const index_t n = 3, m = 2;
+	float64_t data_U[] = {-0.92947874, 0.,          0.,
+	                      -1.1432887,  -0.27048649, 0.,
+	                      -0.87119086, -0.05919915, 0.11869106};
+	float64_t data_B[] = {0.76775073, 0.88471312, 0.34795225,
+	                      0.94311546, 0.59630347, 0.65820143};
+	float64_t result[] = {1.238677,    -3.91243241, 2.9315793,
+	                      -2.00784647, -3.41825732, 5.54550138};
+
+	SGMatrix<float64_t> L(data_U, n, n, false);
+	SGMatrix<float64_t> B(data_B, n, m, false);
+
+	auto X = triangular_solver(L, B, false);
+
+	for (index_t i = 0; i < (index_t)X.size(); ++i)
+		EXPECT_NEAR(X[i], result[i], 1E-6);
+}
+
 TEST(LinalgBackendEigen, SGVector_zero)
 {
 	const index_t n = 16;