diff --git a/src/shogun/mathematics/linalg/LinalgBackendEigen.cpp b/src/shogun/mathematics/linalg/LinalgBackendEigen.cpp new file mode 100644 index 00000000000..5af56548c58 --- /dev/null +++ b/src/shogun/mathematics/linalg/LinalgBackendEigen.cpp @@ -0,0 +1,683 @@ +/* + * Copyright (c) 2016, Shogun-Toolbox e.V. + * 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. + * + * 3. Neither the name of the copyright holder nor the names of its + * contributors may be used to endorse or promote products derived from + * this software without specific prior written permission. + * + * 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 HOLDER 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. + * + * Authors: 2016 Pan Deng, Soumyajit De, Heiko Strathmann, Viktor Gal + */ + +#include +#include + +using namespace shogun; + +#define BACKEND_GENERIC_IN_PLACE_ADD(Type, Container) \ +void LinalgBackendEigen::add(Container& a, Container& b, Type alpha, \ +Type beta, Container& result) const \ +{ \ +add_impl(a, b, alpha, beta, result); \ +} +DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_IN_PLACE_ADD, SGVector) +DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_IN_PLACE_ADD, SGMatrix) +#undef BACKEND_GENERIC_IN_PLACE_ADD + +#define BACKEND_GENERIC_ADD_COL_VEC(Type, Container) \ +void LinalgBackendEigen::add_col_vec(const SGMatrix& A, index_t i, \ + const SGVector& b, Container& result, Type alpha, Type beta) const \ +{ \ + add_col_vec_impl(A, i, b, result, alpha, beta); \ +} +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_CHOLESKY_FACTOR(Type, Container) \ +Container LinalgBackendEigen::cholesky_factor(const Container& A, \ + const bool lower) const \ +{ \ + return cholesky_factor_impl(A, lower); \ +} +DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_FACTOR, SGMatrix) +#undef BACKEND_GENERIC_CHOLESKY_FACTOR + +#define BACKEND_GENERIC_CHOLESKY_SOLVER(Type, Container) \ +SGVector LinalgBackendEigen::cholesky_solver(const Container& L, \ + const SGVector& b, const bool lower) const \ +{ \ + return cholesky_solver_impl(L, b, lower); \ +} +DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_SOLVER, SGMatrix) +#undef BACKEND_GENERIC_CHOLESKY_SOLVER + +#define BACKEND_GENERIC_DOT(Type, Container) \ +Type LinalgBackendEigen::dot(const Container& a, const Container& b) const \ +{ \ + return dot_impl(a, b); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_DOT, SGVector) +#undef BACKEND_GENERIC_DOT + +#define BACKEND_GENERIC_IN_PLACE_ELEMENT_PROD(Type, Container) \ +void LinalgBackendEigen::element_prod(Container& a, Container& b,\ + Container& result) const \ +{ \ + element_prod_impl(a, b, result); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_ELEMENT_PROD, SGMatrix) +#undef BACKEND_GENERIC_IN_PLACE_ELEMENT_PROD + +#define BACKEND_GENERIC_IN_PLACE_BLOCK_ELEMENT_PROD(Type, Container) \ +void LinalgBackendEigen::element_prod(linalg::Block>& a, \ + linalg::Block>& b, Container& result) const \ +{ \ + element_prod_impl(a, b, result); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_BLOCK_ELEMENT_PROD, SGMatrix) +#undef BACKEND_GENERIC_IN_PLACE_BLOCK_ELEMENT_PROD + +#define BACKEND_GENERIC_IDENTITY(Type, Container) \ +void LinalgBackendEigen::identity(Container& I) const \ +{ \ + identity_impl(I); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IDENTITY, SGMatrix) +#undef BACKEND_GENERIC_IDENTITY + +#define BACKEND_GENERIC_LOGISTIC(Type, Container) \ +void LinalgBackendEigen::logistic(Container& a, Container& result) const \ +{ \ + logistic_impl(a, result); \ +} +DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_LOGISTIC, SGMatrix) +#undef BACKEND_GENERIC_LOGISTIC + +#define BACKEND_GENERIC_IN_PLACE_MATRIX_PROD(Type, Container) \ +void LinalgBackendEigen::matrix_prod(SGMatrix& a, Container& b,\ + Container& result, bool transpose_A, bool transpose_B) const \ +{ \ + matrix_prod_impl(a, b, result, transpose_A, transpose_B); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_MATRIX_PROD, SGVector) +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_MATRIX_PROD, SGMatrix) +#undef BACKEND_GENERIC_IN_PLACE_MATRIX_PROD + +#define BACKEND_GENERIC_MAX(Type, Container) \ +Type LinalgBackendEigen::max(const Container& a) const \ +{ \ + return max_impl(a); \ +} +DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_MAX, SGVector) +DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_MAX, SGMatrix) +#undef BACKEND_GENERIC_MAX + +#define BACKEND_GENERIC_REAL_MEAN(Type, Container) \ +float64_t LinalgBackendEigen::mean(const Container& a) const \ +{ \ + return mean_impl(a); \ +} +DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGVector) +DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGMatrix) +#undef BACKEND_GENERIC_REAL_MEAN + +#define BACKEND_GENERIC_COMPLEX_MEAN(Container) \ +complex128_t LinalgBackendEigen::mean(const Container& a) const \ +{ \ + return mean_impl(a); \ +} + BACKEND_GENERIC_COMPLEX_MEAN(SGVector) + BACKEND_GENERIC_COMPLEX_MEAN(SGMatrix) +#undef BACKEND_GENERIC_COMPLEX_MEAN + +#define BACKEND_GENERIC_RANGE_FILL(Type, Container) \ +void LinalgBackendEigen::range_fill(Container& a, const Type start) const \ +{ \ + range_fill_impl(a, start); \ +} +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_IN_PLACE_SCALE(Type, Container) \ +void LinalgBackendEigen::scale(Container& a, Type alpha, Container& result) const \ +{ \ + scale_impl(a, alpha, result); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_SCALE, SGVector) +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_SCALE, SGMatrix) +#undef BACKEND_GENERIC_IN_PLACE_SCALE + +#define BACKEND_GENERIC_SET_CONST(Type, Container) \ +void LinalgBackendEigen::set_const(Container& a, const Type value) const \ +{ \ + set_const_impl(a, value); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_SET_CONST, SGVector) +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_SET_CONST, SGMatrix) +#undef BACKEND_GENERIC_SET_CONST + +#define BACKEND_GENERIC_SUM(Type, Container) \ +Type LinalgBackendEigen::sum(const Container& a, bool no_diag) const \ +{ \ + return sum_impl(a, no_diag); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_SUM, SGVector) +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_SUM, SGMatrix) +#undef BACKEND_GENERIC_SUM + +#define BACKEND_GENERIC_BLOCK_SUM(Type, Container) \ +Type LinalgBackendEigen::sum(const linalg::Block>& a, bool no_diag) const \ +{ \ + return sum_impl(a, no_diag); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_BLOCK_SUM, SGMatrix) +#undef BACKEND_GENERIC_BLOCK_SUM + +#define BACKEND_GENERIC_SYMMETRIC_SUM(Type, Container) \ +Type LinalgBackendEigen::sum_symmetric(const Container& a, bool no_diag) const \ +{ \ + return sum_symmetric_impl(a, no_diag); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_SYMMETRIC_SUM, SGMatrix) +#undef BACKEND_GENERIC_SYMMETRIC_SUM + +#define BACKEND_GENERIC_SYMMETRIC_BLOCK_SUM(Type, Container) \ +Type LinalgBackendEigen::sum_symmetric(const linalg::Block>& a, bool no_diag) const \ +{ \ + return sum_symmetric_impl(a, no_diag); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_SYMMETRIC_BLOCK_SUM, SGMatrix) +#undef BACKEND_GENERIC_SYMMETRIC_BLOCK_SUM + +#define BACKEND_GENERIC_COLWISE_SUM(Type, Container) \ +SGVector LinalgBackendEigen::colwise_sum(const Container& a, bool no_diag) const \ +{ \ + return colwise_sum_impl(a, no_diag); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_COLWISE_SUM, SGMatrix) +#undef BACKEND_GENERIC_COLWISE_SUM + +#define BACKEND_GENERIC_BLOCK_COLWISE_SUM(Type, Container) \ +SGVector LinalgBackendEigen::colwise_sum(const linalg::Block>& a, bool no_diag) const \ +{ \ + return colwise_sum_impl(a, no_diag); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_BLOCK_COLWISE_SUM, SGMatrix) +#undef BACKEND_GENERIC_BLOCK_COLWISE_SUM + +#define BACKEND_GENERIC_ROWWISE_SUM(Type, Container) \ +SGVector LinalgBackendEigen::rowwise_sum(const Container& a, bool no_diag) const \ +{ \ + return rowwise_sum_impl(a, no_diag); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ROWWISE_SUM, SGMatrix) +#undef BACKEND_GENERIC_ROWWISE_SUM + +#define BACKEND_GENERIC_BLOCK_ROWWISE_SUM(Type, Container) \ +SGVector LinalgBackendEigen::rowwise_sum(const linalg::Block>& a, bool no_diag) const \ +{ \ + return rowwise_sum_impl(a, no_diag); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_BLOCK_ROWWISE_SUM, SGMatrix) +#undef BACKEND_GENERIC_BLOCK_ROWWISE_SUM + +#define BACKEND_GENERIC_TRACE(Type, Container) \ +Type LinalgBackendEigen::trace(const Container& A) const \ +{ \ + return trace_impl(A); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_TRACE, SGMatrix) +#undef BACKEND_GENERIC_TRACE + +#define BACKEND_GENERIC_ZERO(Type, Container) \ +void LinalgBackendEigen::zero(Container& a) const \ +{ \ + zero_impl(a); \ +} +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ZERO, SGVector) +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ZERO, SGMatrix) +#undef BACKEND_GENERIC_ZERO + +#undef DEFINE_FOR_ALL_PTYPE +#undef DEFINE_FOR_REAL_PTYPE +#undef DEFINE_FOR_NON_INTEGER_PTYPE +#undef DEFINE_FOR_NUMERIC_PTYPE + +template +void LinalgBackendEigen::add_impl(SGVector& a, SGVector& b, T alpha, T beta, SGVector& result) const +{ + typename SGVector::EigenVectorXtMap a_eig = a; + typename SGVector::EigenVectorXtMap b_eig = b; + typename SGVector::EigenVectorXtMap result_eig = result; + + result_eig = alpha * a_eig + beta * b_eig; +} + +template +void LinalgBackendEigen::add_impl(SGMatrix& a, SGMatrix& b, T alpha, T beta, SGMatrix& result) const +{ + typename SGMatrix::EigenMatrixXtMap a_eig = a; + typename SGMatrix::EigenMatrixXtMap b_eig = b; + typename SGMatrix::EigenMatrixXtMap result_eig = result; + + result_eig = alpha * a_eig + beta * b_eig; +} + +template +void LinalgBackendEigen::add_col_vec_impl(const SGMatrix& A, index_t i, const SGVector& b, + SGMatrix& result, T alpha, T beta) const +{ + typename SGMatrix::EigenMatrixXtMap A_eig = A; + typename SGVector::EigenVectorXtMap b_eig = b; + typename SGMatrix::EigenMatrixXtMap result_eig = result; + + result_eig.col(i) = alpha * A_eig.col(i) + beta * b_eig; +} + +template +void LinalgBackendEigen::add_col_vec_impl(const SGMatrix& A, index_t i, const SGVector& b, + SGVector& result, T alpha, T beta) const +{ + typename SGMatrix::EigenMatrixXtMap A_eig = A; + typename SGVector::EigenVectorXtMap b_eig = b; + typename SGVector::EigenVectorXtMap result_eig = result; + + result_eig = alpha * A_eig.col(i) + beta * b_eig; +} + +template +SGMatrix LinalgBackendEigen::cholesky_factor_impl(const SGMatrix& A, const bool lower) const +{ + SGMatrix c(A.num_rows, A.num_cols); + set_const_impl(c, 0); + typename SGMatrix::EigenMatrixXtMap A_eig = A; + typename SGMatrix::EigenMatrixXtMap c_eig = c; + + Eigen::LLT > llt(A_eig); + + //compute matrix L or U + if(lower==false) + c_eig = llt.matrixU(); + else + c_eig = llt.matrixL(); + + /* + * checking for success + * + * 0: Eigen::Success. Decomposition was successful + * 1: Eigen::NumericalIssue. The provided data did not satisfy the prerequisites. + */ + REQUIRE(llt.info()!=Eigen::NumericalIssue, "Matrix is not Hermitian positive definite!\n"); + + return c; +} + +template +SGVector LinalgBackendEigen::cholesky_solver_impl(const SGMatrix& L, const SGVector& b, + const bool lower) const +{ + SGVector x(b.size()); + set_const_impl(x, 0); + typename SGMatrix::EigenMatrixXtMap L_eig = L; + typename SGVector::EigenVectorXtMap b_eig = b; + typename SGVector::EigenVectorXtMap x_eig = x; + + if (lower == false) + { + Eigen::TriangularView::EigenMatrixXt, + 0, Eigen::Stride<0,0> >, Eigen::Upper> tlv(L_eig); + + x_eig = (tlv.transpose()).solve(tlv.solve(b_eig)); + } + else + { + Eigen::TriangularView::EigenMatrixXt, + 0, Eigen::Stride<0,0> >, Eigen::Lower> tlv(L_eig); + x_eig = (tlv.transpose()).solve(tlv.solve(b_eig)); + } + + return x; +} + +template +T LinalgBackendEigen::dot_impl(const SGVector& a, const SGVector& b) const +{ + return (typename SGVector::EigenVectorXtMap(a)).dot(typename SGVector::EigenVectorXtMap(b)); +} + +template +void LinalgBackendEigen::element_prod_impl(SGMatrix& a, SGMatrix& b, SGMatrix& result) const +{ + typename SGMatrix::EigenMatrixXtMap a_eig = a; + typename SGMatrix::EigenMatrixXtMap b_eig = b; + typename SGMatrix::EigenMatrixXtMap result_eig = result; + + result_eig = a_eig.array() * b_eig.array(); +} + +template +void LinalgBackendEigen::identity_impl(SGMatrix& I) const +{ + typename SGMatrix::EigenMatrixXtMap I_eig = I; + I_eig.setIdentity(); +} + +template +void LinalgBackendEigen::logistic_impl(SGMatrix& a, SGMatrix& result) const +{ + typename SGMatrix::EigenMatrixXtMap a_eig = a; + typename SGMatrix::EigenMatrixXtMap result_eig = result; + + result_eig = (T)1 / (1 + ((-1 * a_eig).array()).exp()); +} + +template +void LinalgBackendEigen::element_prod_impl(linalg::Block>& a, + linalg::Block>& b, SGMatrix& result) const +{ + typename SGMatrix::EigenMatrixXtMap a_eig = a.m_matrix; + typename SGMatrix::EigenMatrixXtMap b_eig = b.m_matrix; + typename SGMatrix::EigenMatrixXtMap result_eig = result; + + Eigen::Block::EigenMatrixXtMap> a_block = + a_eig.block(a.m_row_begin, a.m_col_begin, a.m_row_size, a.m_col_size); + Eigen::Block::EigenMatrixXtMap> b_block = + b_eig.block(b.m_row_begin, b.m_col_begin, b.m_row_size, b.m_col_size); + + result_eig = a_block.array() * b_block.array(); +} + +template +void LinalgBackendEigen::matrix_prod_impl(SGMatrix& a, SGVector& b, SGVector& result, + bool transpose, bool transpose_B) const +{ + typename SGMatrix::EigenMatrixXtMap a_eig = a; + typename SGVector::EigenVectorXtMap b_eig = b; + typename SGVector::EigenVectorXtMap result_eig = result; + + if (transpose) + result_eig = a_eig.transpose() * b_eig; + else + result_eig = a_eig * b_eig; +} + +template +void LinalgBackendEigen::matrix_prod_impl(SGMatrix& a, SGMatrix& b, SGMatrix& result, + bool transpose_A, bool transpose_B) const +{ + typename SGMatrix::EigenMatrixXtMap a_eig = a; + typename SGMatrix::EigenMatrixXtMap b_eig = b; + typename SGMatrix::EigenMatrixXtMap result_eig = result; + + if (transpose_A && transpose_B) + result_eig = a_eig.transpose() * b_eig.transpose(); + + else if (transpose_A) + result_eig = a_eig.transpose() * b_eig; + + else if (transpose_B) + result_eig = a_eig * b_eig.transpose(); + + else + result_eig = a_eig * b_eig; +} + +template +T LinalgBackendEigen::max_impl(const SGVector& vec) const +{ + return (typename SGVector::EigenVectorXtMap(vec)).maxCoeff(); +} + +template +T LinalgBackendEigen::max_impl(const SGMatrix& mat) const +{ + return (typename SGMatrix::EigenMatrixXtMap(mat)).maxCoeff(); +} + +template class Container> +typename std::enable_if::value, float64_t>::type +LinalgBackendEigen::mean_impl(const Container& a) const +{ + return sum_impl(a)/(float64_t(a.size())); +} + +template