From f1e5729b30d0529aa1273fc5e2336de6537a16af Mon Sep 17 00:00:00 2001 From: Heiko Strathmann Date: Tue, 31 Jul 2018 12:18:27 +0100 Subject: [PATCH] add plain and cholesky rank updates --- src/shogun/clustering/GMM.cpp | 4 +- src/shogun/distributions/Gaussian.cpp | 2 +- .../mathematics/linalg/LinalgBackendBase.h | 54 +++++++- .../mathematics/linalg/LinalgBackendEigen.h | 48 ++++++- src/shogun/mathematics/linalg/LinalgMacros.h | 2 +- .../mathematics/linalg/LinalgNamespace.h | 117 ++++++++++++++---- .../linalg/backend/eigen/BasicOps.cpp | 46 ++++++- .../linalg/backend/eigen/Decompositions.cpp | 30 ++++- .../linalg/backend/eigen/EigenSolvers.cpp | 2 +- .../linalg/backend/eigen/Functions.cpp | 10 +- .../mathematics/linalg/backend/eigen/Misc.cpp | 2 +- .../linalg/backend/eigen/Solvers.cpp | 2 +- .../linalg/backend/eigen/SpecialPurposes.cpp | 4 +- src/shogun/metric/LMNNImpl.cpp | 10 +- .../operations/Eigen3_operations_unittest.cc | 116 +++++++++++++++++ 15 files changed, 394 insertions(+), 55 deletions(-) diff --git a/src/shogun/clustering/GMM.cpp b/src/shogun/clustering/GMM.cpp index 67186c05e0c..01fb85d791d 100644 --- a/src/shogun/clustering/GMM.cpp +++ b/src/shogun/clustering/GMM.cpp @@ -591,8 +591,8 @@ void CGMM::max_likelihood(SGMatrix alpha, float64_t min_cov) switch (cov_type) { case FULL: - linalg::dger( - alpha.matrix[j * alpha.num_cols + i], v, v, cov_sum); + linalg::rank_update( + cov_sum, v, alpha.matrix[j * alpha.num_cols + i]); break; case DIAG: { diff --git a/src/shogun/distributions/Gaussian.cpp b/src/shogun/distributions/Gaussian.cpp index 9ee3c0b9d8f..99b04b540d9 100644 --- a/src/shogun/distributions/Gaussian.cpp +++ b/src/shogun/distributions/Gaussian.cpp @@ -167,7 +167,7 @@ float64_t CGaussian::update_params_em(const SGVector alpha_k) CblasRowMajor, num_dim, num_dim, alpha_k[j], v.vector, 1, v.vector, 1, (double*)cov_sum.matrix, num_dim); #else - linalg::dger(alpha_k[j], v, v, cov_sum); + linalg::rank_update(cov_sum, v, alpha_k[j]); #endif break; case DIAG: diff --git a/src/shogun/mathematics/linalg/LinalgBackendBase.h b/src/shogun/mathematics/linalg/LinalgBackendBase.h index 9487f9deece..5964ce9a761 100644 --- a/src/shogun/mathematics/linalg/LinalgBackendBase.h +++ b/src/shogun/mathematics/linalg/LinalgBackendBase.h @@ -69,7 +69,7 @@ namespace shogun METHODNAME(floatmax_t, Container); \ METHODNAME(complex128_t, Container); -#define DEFINE_FOR_REAL_PTYPE(METHODNAME, Container) \ +#define DEFINE_FOR_NON_COMPLEX_PTYPE(METHODNAME, Container) \ METHODNAME(bool, Container); \ METHODNAME(char, Container); \ METHODNAME(int8_t, Container); \ @@ -84,6 +84,11 @@ namespace shogun METHODNAME(float64_t, Container); \ METHODNAME(floatmax_t, Container); +#define DEFINE_FOR_REAL_PTYPE(METHODNAME, Container) \ + METHODNAME(float32_t, Container); \ + METHODNAME(float64_t, Container); \ + METHODNAME(floatmax_t, Container); + #define DEFINE_FOR_NON_INTEGER_PTYPE(METHODNAME, Container) \ METHODNAME(float32_t, Container); \ METHODNAME(float64_t, Container); \ @@ -140,6 +145,20 @@ namespace shogun DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD_DIAG, SGMatrix); #undef BACKEND_GENERIC_ADD_DIAG +/** + * Wrapper method of add diagonal vector A.diagonal = A.diagonal + beta * b. + * + * @see linalg::add_ridge + */ +#define BACKEND_GENERIC_ADD_RIDGE(Type, Container) \ + virtual void add_ridge(SGMatrix& A, Type beta) const \ + { \ + SG_SNOTIMPLEMENTED; \ + return; \ + } + DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD_RIDGE, SGMatrix); +#undef BACKEND_GENERIC_ADD_RIDGE + /** * Wrapper method of add vector to each column of matrix. * @@ -214,6 +233,33 @@ namespace shogun DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_SOLVER, SGMatrix) #undef BACKEND_GENERIC_CHOLESKY_SOLVER +/** + * Wrapper for rank one update of Cholesky decomposition + * + * @see linalg::cholesky_factor + */ +#define BACKEND_GENERIC_CHOLESKY_RANK_UPDATE(Type, Container) \ + virtual void cholesky_rank_update( \ + Container& L, const SGVector& b, Type alpha, \ + const bool lower) const \ + { \ + SG_SNOTIMPLEMENTED; \ + } + DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_CHOLESKY_RANK_UPDATE, SGMatrix) +#undef BACKEND_GENERIC_CHOLESKY_RANK_UPDATE + +/** + * Wrapper for rank one update of a square matrix + */ +#define BACKEND_GENERIC_RANK_UPDATE(Type, Container) \ + virtual void rank_update( \ + Container& A, const SGVector& b, Type alpha) const \ + { \ + SG_SNOTIMPLEMENTED; \ + } + DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_RANK_UPDATE, SGMatrix) +#undef BACKEND_GENERIC_RANK_UPDATE + /** * Wrapper method of LDLT Cholesky decomposition * @@ -469,8 +515,8 @@ namespace shogun SG_SNOTIMPLEMENTED; \ return 0; \ } - DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGVector) - DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGMatrix) + DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGVector) + DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGMatrix) #undef BACKEND_GENERIC_REAL_MEAN /** @@ -849,7 +895,7 @@ namespace shogun #undef BACKEND_GENERIC_FROM_GPU #undef DEFINE_FOR_ALL_PTYPE -#undef DEFINE_FOR_REAL_PTYPE +#undef DEFINE_FOR_NON_COMPLEX_PTYPE #undef DEFINE_FOR_NON_INTEGER_PTYPE }; } diff --git a/src/shogun/mathematics/linalg/LinalgBackendEigen.h b/src/shogun/mathematics/linalg/LinalgBackendEigen.h index eca6a8afa6d..4ce8ee44d17 100644 --- a/src/shogun/mathematics/linalg/LinalgBackendEigen.h +++ b/src/shogun/mathematics/linalg/LinalgBackendEigen.h @@ -73,6 +73,12 @@ namespace shogun DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_DIAG, SGMatrix) #undef BACKEND_GENERIC_ADD_DIAG +/** Implementation of @see LinalgBackendBase::add_ridge */ +#define BACKEND_GENERIC_ADD_RIDGE(Type, Container) \ + virtual void add_ridge(SGMatrix& A, Type beta) const; + DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_RIDGE, SGMatrix) +#undef BACKEND_GENERIC_ADD_RIDGE + /** Implementation of @see LinalgBackendBase::add_vector */ #define BACKEND_GENERIC_ADD(Type, Container) \ virtual void add_vector( \ @@ -101,6 +107,21 @@ namespace shogun DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_FACTOR, SGMatrix) #undef BACKEND_GENERIC_CHOLESKY_FACTOR +/** Implementation of @see LinalgBackendBase::cholesky_rank_update */ +#define BACKEND_GENERIC_CHOLESKY_RANK_UPDATE(Type, Container) \ + virtual void cholesky_rank_update( \ + Container& L, const SGVector& b, Type alpha, \ + const bool lower) const; + DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_CHOLESKY_RANK_UPDATE, SGMatrix) +#undef BACKEND_GENERIC_CHOLESKY_RANK_UPDATE + +/** Implementation of @see LinalgBackendBase::rank_update */ +#define BACKEND_GENERIC_RANK_UPDATE(Type, Container) \ + virtual void rank_update( \ + Container& A, const SGVector& b, Type alpha) const; + DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_RANK_UPDATE, SGMatrix) +#undef BACKEND_GENERIC_RANK_UPDATE + /** Implementation of @see LinalgBackendBase::cholesky_solver */ #define BACKEND_GENERIC_CHOLESKY_SOLVER(Type, Container) \ virtual SGVector cholesky_solver( \ @@ -230,15 +251,15 @@ namespace shogun /** Implementation of @see LinalgBackendBase::max */ #define BACKEND_GENERIC_MAX(Type, Container) \ virtual Type max(const Container& a) const; - DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_MAX, SGVector) - DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_MAX, SGMatrix) + DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_MAX, SGVector) + DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_MAX, SGMatrix) #undef BACKEND_GENERIC_MAX /** Implementation of @see LinalgBackendBase::mean */ #define BACKEND_GENERIC_REAL_MEAN(Type, Container) \ virtual float64_t mean(const Container& a) const; - DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGVector) - DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGMatrix) + DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGVector) + DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGMatrix) #undef BACKEND_GENERIC_REAL_MEAN /** Implementation of @see LinalgBackendBase::mean */ @@ -283,7 +304,7 @@ namespace shogun #define BACKEND_GENERIC_RECTIFIED_LINEAR(Type, Container) \ virtual void rectified_linear( \ const Container& a, Container& result) const; - DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_RECTIFIED_LINEAR, SGMatrix) + DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_RECTIFIED_LINEAR, SGMatrix) #undef BACKEND_GENERIC_RECTIFIED_LINEAR /** Implementation of @see linalg::scale */ @@ -409,7 +430,7 @@ namespace shogun #undef BACKEND_GENERIC_ZERO #undef DEFINE_FOR_ALL_PTYPE -#undef DEFINE_FOR_REAL_PTYPE +#undef DEFINE_FOR_NON_COMPLEX_PTYPE #undef DEFINE_FOR_NON_INTEGER_PTYPE #undef DEFINE_FOR_NUMERIC_PTYPE @@ -443,6 +464,10 @@ namespace shogun void add_diag_impl( SGMatrix& A, const SGVector& b, T alpha, T beta) const; + /** Eigen3 add diagonal scalar method */ + template + void add_ridge_impl(SGMatrix& A, T beta) const; + /** Eigen3 add vector to each column of matrix method */ template void add_vector_impl( @@ -466,6 +491,17 @@ namespace shogun SGMatrix cholesky_factor_impl(const SGMatrix& A, const bool lower) const; + /** Eigen3 Cholesky rank one update */ + template + void cholesky_rank_update_impl( + SGMatrix& L, const SGVector& b, T alpha, + const bool lower) const; + + /** Eigen3 rank one update */ + template + void + rank_update_impl(SGMatrix& A, const SGVector& b, T alpha) const; + /** Eigen3 Cholesky solver */ template SGVector cholesky_solver_impl( diff --git a/src/shogun/mathematics/linalg/LinalgMacros.h b/src/shogun/mathematics/linalg/LinalgMacros.h index 0cdbcd90da4..4ba9033689a 100644 --- a/src/shogun/mathematics/linalg/LinalgMacros.h +++ b/src/shogun/mathematics/linalg/LinalgMacros.h @@ -46,7 +46,7 @@ METHODNAME(floatmax_t, Container); \ METHODNAME(complex128_t, Container); -#define DEFINE_FOR_REAL_PTYPE(METHODNAME, Container) \ +#define DEFINE_FOR_NON_COMPLEX_PTYPE(METHODNAME, Container) \ METHODNAME(bool, Container); \ METHODNAME(char, Container); \ METHODNAME(int8_t, Container); \ diff --git a/src/shogun/mathematics/linalg/LinalgNamespace.h b/src/shogun/mathematics/linalg/LinalgNamespace.h index 049b4d131d5..f14b8b11cea 100644 --- a/src/shogun/mathematics/linalg/LinalgNamespace.h +++ b/src/shogun/mathematics/linalg/LinalgNamespace.h @@ -583,6 +583,21 @@ namespace shogun infer_backend(A, SGMatrix(b))->add_diag(A, b, alpha, beta); } + /** + * Performs the operation A.diagonal = A.diagonal + beta. + * The matrix is not required to be square. + * + * @param A The matrix + * @param beta Constant to be multiplied by the vector + */ + template + void add_ridge(SGMatrix& A, T beta) + { + auto diagonal_len = CMath::min(A.num_cols, A.num_rows); + REQUIRE(diagonal_len > 0, "Matrix can't be empty.\n"); + infer_backend(A)->add_ridge(A, beta); + } + /** * Performs the operation result = alpha * A.col(i) + beta * b, * for each column of A. @@ -614,7 +629,7 @@ namespace shogun } /** - * Adds a scalar to each element of a vector or a matrix in-place. + * Adds a scalar to the matrix/vector in place. * * @param a Vector or matrix * @param b Scalar to be added @@ -644,34 +659,12 @@ namespace shogun infer_backend(A)->center_matrix(A); } - /** - * Performs the operation A = alpha * x * y' + A - * - * @param alpha scaling factor for vector x - * @param x vector - * @param y vector - * @param A m artix - */ - template - void - dger(T alpha, const SGVector x, const SGVector y, SGMatrix& A) - { - auto x_martix = - SGVector::convert_to_matrix(x, A.num_rows, 1, false); - auto y_martix = - SGVector::convert_to_matrix(y, A.num_cols, 1, false); - - auto temp_martix = SGMatrix::matrix_multiply( - x_martix, y_martix, false, true, alpha); - add(A, temp_martix, A); - } - /** * Compute the cholesky decomposition \f$A = L L^{*}\f$ or \f$A = U^{*} * U\f$ * of a Hermitian positive definite matrix * - * @param A The matrix whose cholesky decomposition is to be computed + * @param A The matrix whose Cholesky decomposition is to be computed * @param lower Whether to compute the upper or lower triangular * Cholesky factorization (default: lower) * @return The upper or lower triangular Cholesky factorization @@ -680,9 +673,63 @@ namespace shogun SGMatrix cholesky_factor(const SGMatrix& A, 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); return infer_backend(A)->cholesky_factor(A, lower); } + /** + * Updates the Cholesky factorization \f$A = L L^{*}\f$ with a rank one + * term in-place. + * If A = LL^T before the update, then after it + * LL^{*} = A + alpha * b b^{*} + * + * @param L Triangular matrix, Cholesky factorization of A + * @param b Vector whose outer product with itself is the update + * @param alpha Scaling factor + * @param lower Whether to use L as the upper or lower triangular + * Cholesky factorization (default:lower) + */ + template + void cholesky_rank_update( + SGMatrix& L, const SGVector& b, T alpha = 1, + bool lower = true) + { + REQUIRE( + L.num_rows == L.num_cols, "Matrix (%dx%d) is not square\n", + L.num_rows, L.num_cols); + REQUIRE( + L.num_rows == b.size(), + "Vector size (%d) must match matrix size (%dx%d)\n", b.size(), + L.num_rows, L.num_rows); + return infer_backend(L, SGMatrix(b)) + ->cholesky_rank_update(L, b, alpha, lower); + } + + /** + * Updates a matrix \f$A\f$ with a rank one term in-place,If A = LL^T + * before the update, then after it + * A = A + alpha * b b^{*} + * + * @param A square matrix + * @param b Vector whose outer product with itself is the update + * @param alpha Scaling factor + */ + template + void rank_update(SGMatrix& A, const SGVector& b, T alpha = 1) + { + REQUIRE( + A.num_rows == A.num_cols, "Matrix (%dx%d) is not square\n", + A.num_rows, A.num_cols); + REQUIRE( + A.num_rows == b.size(), + "Vector size (%d) must match matrix size (%dx%d)\n", b.size(), + A.num_rows, A.num_rows); + return infer_backend(A, SGMatrix(b))->rank_update(A, b, alpha); + } + /** * Solve the linear equations \f$Ax=b\f$, given the Cholesky * factorization of A, @@ -698,6 +745,14 @@ namespace shogun SGVector cholesky_solver( const SGMatrix& L, const SGVector& b, const bool lower = true) { + REQUIRE( + L.num_rows == L.num_cols, + "Matrix dimensions (%dx%d) are not square\n", L.num_rows, + L.num_cols); + REQUIRE( + L.num_rows == b.size(), + "Vector size (%d) must match matrix size (%dx%d)\n", b.size(), + L.num_rows); return infer_backend(L, SGMatrix(b)) ->cholesky_solver(L, b, lower); } @@ -763,6 +818,20 @@ namespace shogun const SGMatrix& L, const SGVector& d, SGVector& p, const SGVector& b, const bool lower = true) { + REQUIRE( + L.num_rows == L.num_cols, + "Matrix dimensions (%dx%d) are not square\n", L.num_rows, + L.num_cols); + REQUIRE( + d.vlen == L.num_rows, "Length of vector d (%d) doesn't match " + "length of diagonal of matrix L (%d)\n", + d.vlen, L.num_rows); + REQUIRE( + p.vlen = L.num_rows, "Length of transpositions vector p (%d) " + "doesn't match length of diagonal of " + "matrix L (%d)\n", + p.vlen, L.num_rows); + return infer_backend(L, SGMatrix(d), SGMatrix(b)) ->ldlt_solver(L, d, p, b, lower); } diff --git a/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp b/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp index b0488468839..bfd7e6afb6b 100644 --- a/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp +++ b/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp @@ -30,6 +30,7 @@ * Authors: 2016 Pan Deng, Soumyajit De, Heiko Strathmann, Viktor Gal */ +#include #include #include @@ -67,6 +68,14 @@ DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_COL_VEC, SGMatrix) DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_DIAG, SGMatrix) #undef BACKEND_GENERIC_ADD_DIAG +#define BACKEND_GENERIC_ADD_RIDGE(Type, Container) \ + void LinalgBackendEigen::add_ridge(SGMatrix& A, Type beta) const \ + { \ + add_ridge_impl(A, beta); \ + } +DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_RIDGE, SGMatrix) +#undef BACKEND_GENERIC_ADD_RIDGE + #define BACKEND_GENERIC_PINVH(Type, Container) \ void LinalgBackendEigen::pinvh( \ const SGMatrix& A, SGMatrix& result) const \ @@ -76,6 +85,15 @@ DEFINE_FOR_NUMERIC_PTYPE(BACKEND_GENERIC_ADD_DIAG, SGMatrix) DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_PINVH, SGMatrix) #undef BACKEND_GENERIC_PINVH +#define BACKEND_GENERIC_RANK_UPDATE(Type, Container) \ + void LinalgBackendEigen::rank_update( \ + SGMatrix& A, const SGVector& b, Type alpha) const \ + { \ + rank_update_impl(A, b, alpha); \ + } +DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_RANK_UPDATE, SGMatrix) +#undef BACKEND_GENERIC_RANK_UPDATE + #define BACKEND_GENERIC_PINV(Type, Container) \ void LinalgBackendEigen::pinv( \ const SGMatrix& A, SGMatrix& result) const \ @@ -176,7 +194,7 @@ DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_SCALE, SGMatrix) #undef BACKEND_GENERIC_IN_PLACE_SCALE #undef DEFINE_FOR_ALL_PTYPE -#undef DEFINE_FOR_REAL_PTYPE +#undef DEFINE_FOR_NON_COMPLEX_PTYPE #undef DEFINE_FOR_NON_INTEGER_PTYPE #undef DEFINE_FOR_NUMERIC_PTYPE @@ -226,6 +244,13 @@ void LinalgBackendEigen::add_diag_impl( A_eig.diagonal() = alpha * A_eig.diagonal() + beta * b_eig; } +template +void LinalgBackendEigen::add_ridge_impl(SGMatrix& A, T beta) const +{ + for (auto i : range(std::min(A.num_rows, A.num_cols))) + A(i, i) += beta; +} + template void LinalgBackendEigen::pinvh_impl( const SGMatrix& A, SGMatrix& result) const @@ -266,6 +291,25 @@ void LinalgBackendEigen::pinv_impl( result_eig = V_eig * s_eig.asDiagonal() * U_eig.transpose(); } +template +void LinalgBackendEigen::rank_update_impl( + SGMatrix& A, const SGVector& b, T alpha) const +{ + T x; + T update; + for (auto j : range(A.num_cols)) + { + x = b[j]; + for (auto i : range(j)) + { + update = alpha * x * b[i]; + A(i, j) += update; + A(j, i) += update; + } + A(j, j) += alpha * x * x; + } +} + template void LinalgBackendEigen::add_col_vec_impl( const SGMatrix& A, index_t i, const SGVector& b, SGVector& result, diff --git a/src/shogun/mathematics/linalg/backend/eigen/Decompositions.cpp b/src/shogun/mathematics/linalg/backend/eigen/Decompositions.cpp index 44f1ce9abc5..260acdfa3d6 100644 --- a/src/shogun/mathematics/linalg/backend/eigen/Decompositions.cpp +++ b/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_CHOLESKY_RANK_UPDATE(Type, Container) \ + void LinalgBackendEigen::cholesky_rank_update( \ + Container& L, const SGVector& b, Type alpha, \ + const bool lower) const \ + { \ + cholesky_rank_update_impl(L, b, alpha, lower); \ + } +DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_CHOLESKY_RANK_UPDATE, SGMatrix) +#undef BACKEND_GENERIC_CHOLESKY_RANK_UPDATE + #define BACKEND_GENERIC_LDLT_FACTOR(Type, Container) \ void LinalgBackendEigen::ldlt_factor( \ const Container& A, Container& L, SGVector& d, \ @@ -66,7 +76,7 @@ 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_COMPLEX_PTYPE #undef DEFINE_FOR_NON_INTEGER_PTYPE #undef DEFINE_FOR_NUMERIC_PTYPE @@ -101,6 +111,24 @@ SGMatrix LinalgBackendEigen::cholesky_factor_impl( return c; } +#include + +template +void LinalgBackendEigen::cholesky_rank_update_impl( + SGMatrix& L, const SGVector& b, T alpha, bool lower) const +{ + typename SGMatrix::EigenMatrixXtMap L_eig = L; + typename SGVector::EigenVectorXtMap b_eig = b; + + if (lower == false) + { + auto U = L_eig.transpose(); + Eigen::internal::llt_rank_update_lower(U, b_eig, alpha); + } + else + Eigen::internal::llt_rank_update_lower(L_eig, b_eig, alpha); +} + template void LinalgBackendEigen::ldlt_factor_impl( const SGMatrix& A, SGMatrix& L, SGVector& d, SGVector& p, diff --git a/src/shogun/mathematics/linalg/backend/eigen/EigenSolvers.cpp b/src/shogun/mathematics/linalg/backend/eigen/EigenSolvers.cpp index 6a98bf6a491..61691a70e58 100644 --- a/src/shogun/mathematics/linalg/backend/eigen/EigenSolvers.cpp +++ b/src/shogun/mathematics/linalg/backend/eigen/EigenSolvers.cpp @@ -57,7 +57,7 @@ DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC, SGMatrix) #undef BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC #undef DEFINE_FOR_ALL_PTYPE -#undef DEFINE_FOR_REAL_PTYPE +#undef DEFINE_FOR_NON_COMPLEX_PTYPE #undef DEFINE_FOR_NON_INTEGER_PTYPE #undef DEFINE_FOR_NUMERIC_PTYPE diff --git a/src/shogun/mathematics/linalg/backend/eigen/Functions.cpp b/src/shogun/mathematics/linalg/backend/eigen/Functions.cpp index 12f3de5d182..13a42928dd1 100644 --- a/src/shogun/mathematics/linalg/backend/eigen/Functions.cpp +++ b/src/shogun/mathematics/linalg/backend/eigen/Functions.cpp @@ -40,8 +40,8 @@ using namespace shogun; { \ return max_impl(a); \ } -DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_MAX, SGVector) -DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_MAX, SGMatrix) +DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_MAX, SGVector) +DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_MAX, SGMatrix) #undef BACKEND_GENERIC_MAX #define BACKEND_GENERIC_REAL_MEAN(Type, Container) \ @@ -49,8 +49,8 @@ DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_MAX, SGMatrix) { \ return mean_impl(a); \ } -DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGVector) -DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGMatrix) +DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGVector) +DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_REAL_MEAN, SGMatrix) #undef BACKEND_GENERIC_REAL_MEAN #define BACKEND_GENERIC_COMPLEX_MEAN(Container) \ @@ -144,7 +144,7 @@ DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_TRACE, SGMatrix) #undef BACKEND_GENERIC_TRACE #undef DEFINE_FOR_ALL_PTYPE -#undef DEFINE_FOR_REAL_PTYPE +#undef DEFINE_FOR_NON_COMPLEX_PTYPE #undef DEFINE_FOR_NON_INTEGER_PTYPE #undef DEFINE_FOR_NUMERIC_PTYPE diff --git a/src/shogun/mathematics/linalg/backend/eigen/Misc.cpp b/src/shogun/mathematics/linalg/backend/eigen/Misc.cpp index a7ba5fa3528..182280f6e90 100644 --- a/src/shogun/mathematics/linalg/backend/eigen/Misc.cpp +++ b/src/shogun/mathematics/linalg/backend/eigen/Misc.cpp @@ -90,7 +90,7 @@ 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_COMPLEX_PTYPE #undef DEFINE_FOR_NON_INTEGER_PTYPE #undef DEFINE_FOR_NUMERIC_PTYPE diff --git a/src/shogun/mathematics/linalg/backend/eigen/Solvers.cpp b/src/shogun/mathematics/linalg/backend/eigen/Solvers.cpp index fca9065cad5..5c34accede7 100644 --- a/src/shogun/mathematics/linalg/backend/eigen/Solvers.cpp +++ b/src/shogun/mathematics/linalg/backend/eigen/Solvers.cpp @@ -80,7 +80,7 @@ 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_COMPLEX_PTYPE #undef DEFINE_FOR_NON_INTEGER_PTYPE #undef DEFINE_FOR_NUMERIC_PTYPE diff --git a/src/shogun/mathematics/linalg/backend/eigen/SpecialPurposes.cpp b/src/shogun/mathematics/linalg/backend/eigen/SpecialPurposes.cpp index 425406b901e..316722382cd 100644 --- a/src/shogun/mathematics/linalg/backend/eigen/SpecialPurposes.cpp +++ b/src/shogun/mathematics/linalg/backend/eigen/SpecialPurposes.cpp @@ -78,7 +78,7 @@ DEFINE_FOR_NON_INTEGER_REAL_PTYPE( { \ rectified_linear_impl(a, result); \ } -DEFINE_FOR_REAL_PTYPE(BACKEND_GENERIC_RECTIFIED_LINEAR, SGMatrix) +DEFINE_FOR_NON_COMPLEX_PTYPE(BACKEND_GENERIC_RECTIFIED_LINEAR, SGMatrix) #undef BACKEND_GENERIC_RECTIFIED_LINEAR #define BACKEND_GENERIC_SOFTMAX(Type, Container) \ @@ -99,7 +99,7 @@ DEFINE_FOR_NON_INTEGER_REAL_PTYPE(BACKEND_GENERIC_SQUARED_ERROR, SGMatrix) #undef BACKEND_GENERIC_SQUARED_ERROR #undef DEFINE_FOR_ALL_PTYPE -#undef DEFINE_FOR_REAL_PTYPE +#undef DEFINE_FOR_NON_COMPLEX_PTYPE #undef DEFINE_FOR_NON_INTEGER_PTYPE #undef DEFINE_FOR_NUMERIC_PTYPE diff --git a/src/shogun/metric/LMNNImpl.cpp b/src/shogun/metric/LMNNImpl.cpp index abf0e02c212..23d2c5b9cd9 100644 --- a/src/shogun/metric/LMNNImpl.cpp +++ b/src/shogun/metric/LMNNImpl.cpp @@ -186,7 +186,7 @@ SGMatrix CLMNNImpl::sum_outer_products( { auto dx = linalg::add( X.get_column(i), X.get_column(target_nn(j, i)), 1.0, -1.0); - linalg::dger(1.0, dx, dx, sop); + linalg::rank_update(sop, dx); } } @@ -256,8 +256,8 @@ void CLMNNImpl::update_gradient( X.get_column(it->example), X.get_column(it->target), 1.0, -1.0); auto dx2 = linalg::add( X.get_column(it->example), X.get_column(it->impostor), 1.0, -1.0); - linalg::dger(-regularization, dx1, dx1, G); - linalg::dger(regularization, dx2, dx2, G); + linalg::rank_update(G, dx1, -regularization); + linalg::rank_update(G, dx2, regularization); } // add the gradient contributions of the new impostors @@ -268,8 +268,8 @@ void CLMNNImpl::update_gradient( X.get_column(it->example), X.get_column(it->target), 1.0, -1.0); auto dx2 = linalg::add( X.get_column(it->example), X.get_column(it->impostor), 1.0, -1.0); - linalg::dger(regularization, dx1, dx1, G); - linalg::dger(-regularization, dx2, dx2, G); + linalg::rank_update(G, dx1, regularization); + linalg::rank_update(G, dx2, -regularization); } } diff --git a/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc b/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc index 309117fcd1a..4e01a70c147 100644 --- a/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc +++ b/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc @@ -225,6 +225,31 @@ TEST(LinalgBackendEigen, add_diag) EXPECT_THROW(add_diag(A3, b3), ShogunException); } +TEST(LinalgBackendEigen, add_ridge) +{ + SGMatrix A1(2, 3); + + A1(0, 0) = 1; + A1(0, 1) = 2; + A1(0, 2) = 3; + A1(1, 0) = 4; + A1(1, 1) = 5; + A1(1, 2) = 6; + + add_ridge(A1, 1.0); + + EXPECT_NEAR(A1(0, 0), 2, 1e-15); + EXPECT_NEAR(A1(0, 1), 2, 1e-15); + EXPECT_NEAR(A1(0, 2), 3, 1e-15); + EXPECT_NEAR(A1(1, 0), 4, 1e-15); + EXPECT_NEAR(A1(1, 1), 6, 1e-15); + EXPECT_NEAR(A1(1, 2), 6, 1e-15); + + // test error cases + SGMatrix A2; + EXPECT_THROW(add_ridge(A2, 1.0), ShogunException); +} + TEST(LinalgBackendEigen, add_vector) { const float64_t alpha = 0.7; @@ -342,6 +367,68 @@ TEST(LinalgBackendEigen, SGMatrix_cholesky_llt_upper) EXPECT_EQ(m.num_cols, U.num_cols); } +TEST(LinalgBackendEigen, SGMatrix_cholesky_rank_update_upper) +{ + const index_t size = 2; + SGMatrix A(size, size); + SGMatrix U(size, size); + SGVector b(size); + Map A_eig(A.matrix, size, size); + Map U_eig(U.matrix, size, size); + Map b_eig(b.vector, size); + + U(0, 0) = 2.0; + U(0, 1) = 1.0; + U(1, 1) = 2.5; + b[0] = 2; + b[1] = 3; + A_eig = U_eig.transpose() * U_eig; + + auto A2 = A.clone(); + Map A2_eig(A2.matrix, A2.num_rows, A2.num_cols); + A2(0, 0) += b[0] * b[0]; + A2(0, 1) += b[0] * b[1]; + A2(1, 0) += b[1] * b[0]; + A2(1, 1) += b[1] * b[1]; + + cholesky_rank_update(U, b, 1.0, false); + EXPECT_NEAR((A2_eig - U_eig.transpose() * U_eig).norm(), 0.0, 1e-14); + + cholesky_rank_update(U, b, -1., false); + EXPECT_NEAR((A_eig - U_eig.transpose() * U_eig).norm(), 0.0, 1e-14); +} + +TEST(LinalgBackendEigen, SGMatrix_cholesky_rank_update_lower) +{ + const index_t size = 2; + SGMatrix A(size, size); + SGMatrix L(size, size); + SGVector b(size); + Map A_eig(A.matrix, size, size); + Map L_eig(L.matrix, size, size); + Map b_eig(b.vector, size); + + L(0, 0) = 2.0; + L(1, 0) = 1.0; + L(1, 1) = 2.5; + b[0] = 2; + b[1] = 3; + A_eig = L_eig * L_eig.transpose(); + + auto A2 = A.clone(); + Map A2_eig(A2.matrix, A2.num_rows, A2.num_cols); + A2(0, 0) += b[0] * b[0]; + A2(0, 1) += b[0] * b[1]; + A2(1, 0) += b[1] * b[0]; + A2(1, 1) += b[1] * b[1]; + + cholesky_rank_update(L, b); + EXPECT_NEAR((A2_eig - L_eig * L_eig.transpose()).norm(), 0.0, 1e-14); + + cholesky_rank_update(L, b, -1.); + EXPECT_NEAR((A_eig - L_eig * L_eig.transpose()).norm(), 0.0, 1e-14); +} + TEST(LinalgBackendEigen, SGMatrix_cholesky_ldlt_lower) { const index_t size = 3; @@ -1973,3 +2060,32 @@ TEST(LinalgBackendEigen, SGMatrix_zero) for (index_t i = 0; i < nrows*ncols; ++i) EXPECT_EQ(A[i], 0); } + +TEST(LinalgBackendEigen, SGMatrix_rank_update) +{ + const index_t size = 2; + SGMatrix A(size, size); + SGVector b(size); + Map A_eig(A.matrix, size, size); + Map b_eig(b.vector, size); + + A(0, 0) = 2.0; + A(1, 0) = 1.0; + A(0, 1) = 1.0; + A(1, 1) = 2.5; + b[0] = 2; + b[1] = 3; + + auto A2 = A.clone(); + Map A2_eig(A2.matrix, size, size); + A2(0, 0) += b[0] * b[0]; + A2(0, 1) += b[0] * b[1]; + A2(1, 0) += b[1] * b[0]; + A2(1, 1) += b[1] * b[1]; + + rank_update(A, b); + EXPECT_NEAR((A2_eig - A_eig).norm(), 0.0, 1e-14); + + rank_update(A, b, -1.); + EXPECT_NEAR((A_eig - A_eig).norm(), 0.0, 1e-14); +}