diff --git a/src/shogun/classifier/svm/NewtonSVM.cpp b/src/shogun/classifier/svm/NewtonSVM.cpp index ba633f3250e..f95e450e3b5 100644 --- a/src/shogun/classifier/svm/NewtonSVM.cpp +++ b/src/shogun/classifier/svm/NewtonSVM.cpp @@ -126,12 +126,10 @@ void CNewtonSVM::iteration() Xsv2sum(x_d, x_d) = size_sv; linalg::add(Xsv2sum, lcrossdiag, Xsv2sum); - SGMatrix inverse((x_d + 1), (x_d + 1)); - SGMatrix::pinv(Xsv2sum.data(), x_d + 1, x_d + 1, inverse.data()); SGVector step(x_d + 1); SGVector s2(x_d + 1); - s2 = linalg::matrix_prod(inverse, grad); + s2 = linalg::matrix_prod(linalg::pinvh(Xsv2sum), grad); for (int32_t i = 0; i < x_d + 1; i++) step[i] = -s2[i]; diff --git a/src/shogun/lib/SGMatrix.cpp b/src/shogun/lib/SGMatrix.cpp index 6b29892917e..c8fc0874834 100644 --- a/src/shogun/lib/SGMatrix.cpp +++ b/src/shogun/lib/SGMatrix.cpp @@ -948,54 +948,7 @@ SGMatrix SGMatrix::create_identity_matrix(index_t si return identity_matrix; } -//Howto construct the pseudo inverse (from "The Matrix Cookbook") -// -//Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m). -// -//The matrix A+ known as the pseudo inverse is unique and does always exist. -// -//The pseudo inverse A+ can be constructed from the singular value -//decomposition A = UDV^T , by A^+ = V(D+)U^T. - #ifdef HAVE_LAPACK -template -float64_t* SGMatrix::pinv( - float64_t* matrix, int32_t rows, int32_t cols, float64_t* target) -{ - if (!target) - target=SG_MALLOC(float64_t, rows*cols); - - char jobu='A'; - char jobvt='A'; - int m=rows; /* for calling external lib */ - int n=cols; /* for calling external lib */ - int lda=m; /* for calling external lib */ - int ldu=m; /* for calling external lib */ - int ldvt=n; /* for calling external lib */ - int info=-1; /* for calling external lib */ - int32_t lsize=CMath::min((int32_t) m, (int32_t) n); - double* s=SG_MALLOC(double, lsize); - double* u=SG_MALLOC(double, m*m); - double* vt=SG_MALLOC(double, n*n); - - wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info); - ASSERT(info==0) - - for (int64_t i=0; i void SGMatrix::inverse(SGMatrix matrix) diff --git a/src/shogun/lib/SGMatrix.h b/src/shogun/lib/SGMatrix.h index e00d551fadf..0b036f3c762 100644 --- a/src/shogun/lib/SGMatrix.h +++ b/src/shogun/lib/SGMatrix.h @@ -393,13 +393,6 @@ template class SGMatrix : public SGReferencedData /** Inverses square matrix in-place */ static void inverse(SGMatrix matrix); - /** Return the pseudo inverse for matrix - * when matrix has shape (rows, cols) the pseudo inverse has (cols, rows) - */ - static float64_t* pinv( - float64_t* matrix, int32_t rows, int32_t cols, - float64_t* target=NULL); - #endif /** Compute trace */ diff --git a/src/shogun/mathematics/linalg/LinalgBackendBase.h b/src/shogun/mathematics/linalg/LinalgBackendBase.h index 4660a188cc0..9487f9deece 100644 --- a/src/shogun/mathematics/linalg/LinalgBackendBase.h +++ b/src/shogun/mathematics/linalg/LinalgBackendBase.h @@ -229,6 +229,32 @@ namespace shogun DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_LDLT_FACTOR, SGMatrix) #undef BACKEND_GENERIC_LDLT_FACTOR +/** + * Wrapper method of pseudo inverse for self adjoint matrices + * + * @see linalg::pinvh + */ +#define BACKEND_GENERIC_PINVH(Type, Container) \ + virtual void pinvh(const SGMatrix& A, SGMatrix& result) const \ + { \ + SG_SNOTIMPLEMENTED; \ + } + DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_PINVH, SGMatrix) +#undef BACKEND_GENERIC_PINVH + +/** + * Wrapper method of pseudo inverse + * + * @see linalg::pinvh + */ +#define BACKEND_GENERIC_PINV(Type, Container) \ + virtual void pinv(const SGMatrix& A, SGMatrix& result) const \ + { \ + SG_SNOTIMPLEMENTED; \ + } + DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_PINV, SGMatrix) +#undef BACKEND_GENERIC_PINV + /** * Wrapper method of LDLT Cholesky solver * diff --git a/src/shogun/mathematics/linalg/LinalgBackendEigen.h b/src/shogun/mathematics/linalg/LinalgBackendEigen.h index 45d926acf56..eca6a8afa6d 100644 --- a/src/shogun/mathematics/linalg/LinalgBackendEigen.h +++ b/src/shogun/mathematics/linalg/LinalgBackendEigen.h @@ -117,6 +117,18 @@ namespace shogun DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_LDLT_FACTOR, SGMatrix) #undef BACKEND_GENERIC_LDLT_FACTOR +/** Implementation of @see LinalgBackendBase::pinvh */ +#define BACKEND_GENERIC_PINVH(Type, Container) \ + virtual void pinvh(const SGMatrix& A, SGMatrix& result) const; + DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_PINVH, SGMatrix) +#undef BACKEND_GENERIC_PINVH + +/** Implementation of @see LinalgBackendBase::pinv */ +#define BACKEND_GENERIC_PINV(Type, Container) \ + virtual void pinv(const SGMatrix& A, SGMatrix& result) const; + DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_PINV, SGMatrix) +#undef BACKEND_GENERIC_PINV + /** Implementation of @see LinalgBackendBase::ldlt_solver */ #define BACKEND_GENERIC_LDLT_SOLVER(Type, Container) \ virtual SGVector ldlt_solver( \ @@ -465,6 +477,14 @@ namespace shogun const SGMatrix& A, SGMatrix& L, SGVector& d, SGVector& p, const bool lower) const; + /** Eigen3 Pseudo inverse of self adjoint matrix */ + template + void pinvh_impl(const SGMatrix& A, SGMatrix& result) const; + + /** Eigen3 Pseudo inverse of self adjoint matrix */ + template + void pinv_impl(const SGMatrix& A, SGMatrix& result) const; + /** Eigen3 LDLT Cholesky solver */ template SGVector ldlt_solver_impl( diff --git a/src/shogun/mathematics/linalg/LinalgNamespace.h b/src/shogun/mathematics/linalg/LinalgNamespace.h index c04389f9c24..049b4d131d5 100644 --- a/src/shogun/mathematics/linalg/LinalgNamespace.h +++ b/src/shogun/mathematics/linalg/LinalgNamespace.h @@ -484,6 +484,79 @@ namespace shogun ->add_col_vec(A, i, b, result, alpha, beta); } + /** Calculates pseudo inverse A+ from eigen values solved by a self + * adjoint eigen solver. + * User should pass an appropriately pre-allocated memory matrix + * or pass the operand argument A as a result. + * + * @param A The matrix + * @param result The matrix that saves the result + */ + template + void pinvh(const SGMatrix& A, SGMatrix& result) + { + + REQUIRE( + result.num_rows == A.num_rows && result.num_cols == A.num_cols, + "Dimension mismatch! A (%d x %d) vs result (%d x %d).\n", + A.num_rows, A.num_cols, result.num_rows, result.num_cols); + + REQUIRE( + A.num_rows == A.num_cols, + "Given matrix is not square (%d x %d)", A.num_rows, A.num_cols); + + infer_backend(A)->pinvh(A, result); + } + + /** + * Calculates pseudo inverse A+ from eigen values solved by a self + * adjoint eigen solver. + * This version returns the result in a newly created matrix. + * + * @param A The matrix + * @return The result matrix + */ + template + SGMatrix pinvh(const SGMatrix& A) + { + SGMatrix result(A.num_rows, A.num_rows); + pinvh(A, result); + return result; + } + + /** Calculates pseudo inverse A+ using singular value decomposition. + * User should pass an appropriately pre-allocated memory matrix + * + * @param A The matrix + * @param result The matrix that saves the result + */ + template + void pinv(const SGMatrix& A, SGMatrix& result) + { + REQUIRE( + result.num_rows == A.num_cols && result.num_cols == A.num_rows, + "Dimension mismatch! Result must be of (%d x %d) provided is " + "(%d x %d).\n", + A.num_cols, A.num_rows, result.num_rows, result.num_cols); + + infer_backend(A)->pinv(A, result); + } + + /** + * Calculates pseudo inverse A+ using singular value decomposition. + * This version returns the result in a newly created matrix. + * + * @param A The matrix + * @return The result matrix + */ + template + SGMatrix pinv(const SGMatrix& A) + { + SGMatrix result(A.num_cols, A.num_rows); + pinv(A, result); + return result; + } + /** * Performs the operation A.diagonal = alpha * A.diagonal + beta * b. * The matrix is not required to be square. diff --git a/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp b/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp index ec01019ff6e..b0488468839 100644 --- a/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp +++ b/src/shogun/mathematics/linalg/backend/eigen/BasicOps.cpp @@ -67,6 +67,24 @@ 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_PINVH(Type, Container) \ + void LinalgBackendEigen::pinvh( \ + const SGMatrix& A, SGMatrix& result) const \ + { \ + pinvh_impl(A, result); \ + } +DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_PINVH, SGMatrix) +#undef BACKEND_GENERIC_PINVH + +#define BACKEND_GENERIC_PINV(Type, Container) \ + void LinalgBackendEigen::pinv( \ + const SGMatrix& A, SGMatrix& result) const \ + { \ + pinv_impl(A, result); \ + } +DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_PINV, SGMatrix) +#undef BACKEND_GENERIC_PINV + #define BACKEND_GENERIC_ADD_VECTOR(Type, Container) \ void LinalgBackendEigen::add_vector( \ const SGMatrix& A, const SGVector& b, \ @@ -208,6 +226,46 @@ void LinalgBackendEigen::add_diag_impl( A_eig.diagonal() = alpha * A_eig.diagonal() + beta * b_eig; } +template +void LinalgBackendEigen::pinvh_impl( + const SGMatrix& A, SGMatrix& result) const +{ + auto m = A.num_rows; + SGVector s(m); + SGMatrix V(m, m); + typename SGMatrix::EigenMatrixXtMap A_eig = A; + typename SGMatrix::EigenMatrixXtMap result_eig = result; + + eigen_solver_symmetric_impl(A, s, V, m); + + typename SGVector::EigenVectorXtMap s_eig = s; + typename SGMatrix::EigenMatrixXtMap V_eig = V; + + float64_t pinv_tol = CMath::MACHINE_EPSILON * m * s_eig.real().maxCoeff(); + s_eig.array() = (s_eig.real().array() > pinv_tol) + .select(s_eig.real().array().inverse(), 0); + result_eig = V_eig * s_eig.asDiagonal() * V_eig.transpose(); +} + +template +void LinalgBackendEigen::pinv_impl( + const SGMatrix& A, SGMatrix& result) const +{ + auto m = std::min(A.num_cols, A.num_rows); + typename SGMatrix::EigenMatrixXtMap A_eig = A; + typename SGMatrix::EigenMatrixXtMap result_eig = result; + + auto svd_eig = A_eig.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV); + auto s_eig = svd_eig.singularValues().real(); + auto U_eig = svd_eig.matrixU().real(); + auto V_eig = svd_eig.matrixV().real(); + + float64_t pinv_tol = CMath::MACHINE_EPSILON * m * s_eig.maxCoeff(); + s_eig.array() = + (s_eig.array() > pinv_tol).select(s_eig.array().inverse(), 0); + result_eig = V_eig * s_eig.asDiagonal() * U_eig.transpose(); +} + template void LinalgBackendEigen::add_col_vec_impl( const SGMatrix& A, index_t i, const SGVector& b, SGVector& result, diff --git a/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc b/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc index 2597f97a652..309117fcd1a 100644 --- a/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc +++ b/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc @@ -465,6 +465,116 @@ TEST(LinalgBackendEigen, SGMatrix_cross_entropy) EXPECT_NEAR(ref, result, 1e-15); } +TEST(LinalgBackendEigen, SGMatrix_pinv_psd) +{ + float64_t A_data[] = {2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0}; + // inverse generated by scipy pinv + float64_t scipy_result_data[] = {0.75, 0.5, 0.25, 0.5, 1.0, + 0.5, 0.25, 0.5, 0.75}; + + SGMatrix A(A_data, 3, 3, false); + SGMatrix result(scipy_result_data, 3, 3, false); + + SGMatrix identity_matrix(3, 3); + linalg::identity(identity_matrix); + // using symmetric eigen solver + SGMatrix A_pinvh(3, 3); + linalg::pinvh(A, A_pinvh); + // using singular value decomposition + SGMatrix A_pinv(3, 3); + linalg::pinv(A, A_pinv); + SGMatrix I_check = linalg::matrix_prod(A, A_pinvh); + for (auto i : range(3)) + { + for (auto j : range(3)) + { + EXPECT_NEAR(identity_matrix(i, j), I_check(i, j), 1e-14); + EXPECT_NEAR(result(i, j), A_pinvh(i, j), 1e-14); + EXPECT_NEAR(result(i, j), A_pinv(i, j), 1e-14); + } + } + // no memory errors + EXPECT_NO_THROW(linalg::pinvh(A, A)); +} + +TEST(LinalgBackendEigen, SGMatrix_pinv_2x4) +{ + SGMatrix A(2, 4); + A(0, 0) = 1; + A(0, 1) = 1; + A(0, 2) = 1; + A(0, 3) = 1; + A(1, 0) = 5; + A(1, 1) = 7; + A(1, 2) = 7; + A(1, 3) = 9; + + SGMatrix identity_matrix(2, 2); + linalg::identity(identity_matrix); + SGMatrix A_pinverse(4, 2); + linalg::pinv(A, A_pinverse); + SGMatrix I_check = linalg::matrix_prod(A, A_pinverse); + for (auto i : range(2)) + { + for (auto j : range(2)) + { + EXPECT_NEAR(identity_matrix(i, j), I_check(i, j), 1e-14); + } + } + + // compare result with scipy pinv + EXPECT_NEAR(A_pinverse(0, 0), 2.0, 1e-14); + EXPECT_NEAR(A_pinverse(0, 1), -0.25, 1e-14); + + EXPECT_NEAR(A_pinverse(1, 0), 0.25, 1e-14); + EXPECT_NEAR(A_pinverse(1, 1), 0.0, 1e-14); + + EXPECT_NEAR(A_pinverse(2, 0), 0.25, 1e-14); + EXPECT_NEAR(A_pinverse(2, 1), 0.0, 1e-14); + + EXPECT_NEAR(A_pinverse(3, 0), -1.5, 1e-14); + EXPECT_NEAR(A_pinverse(3, 1), 0.25, 1e-14); + + // incorrect dimension + EXPECT_THROW(linalg::pinv(A, A), ShogunException); +} + +TEST(LinalgBackendEigen, SGMatrix_pinv_4x2) +{ + SGMatrix A(4, 2); + A(0, 0) = 2.0; + A(0, 1) = -0.25; + A(1, 0) = 0.25; + A(1, 1) = 0.0; + A(2, 0) = 0.25; + A(2, 1) = 0.0; + A(3, 0) = -1.5; + A(3, 1) = 0.25; + + SGMatrix identity_matrix(2, 2); + linalg::identity(identity_matrix); + SGMatrix A_pinverse(2, 4); + linalg::pinv(A, A_pinverse); + SGMatrix I_check = linalg::matrix_prod(A_pinverse, A); + for (auto i : range(2)) + { + for (auto j : range(2)) + { + EXPECT_NEAR(identity_matrix(i, j), I_check(i, j), 1e-14); + } + } + // compare with results from scipy + EXPECT_NEAR(A_pinverse(0, 0), 1.0, 1e-14); + EXPECT_NEAR(A_pinverse(0, 1), 1.0, 1e-14); + EXPECT_NEAR(A_pinverse(0, 2), 1.0, 1e-14); + EXPECT_NEAR(A_pinverse(0, 3), 1.0, 1e-14); + + EXPECT_NEAR(A_pinverse(1, 0), 5.0, 1e-14); + EXPECT_NEAR(A_pinverse(1, 1), 7.0, 1e-14); + EXPECT_NEAR(A_pinverse(1, 2), 7.0, 1e-14); + EXPECT_NEAR(A_pinverse(1, 3), 9.0, 1e-14); +} + TEST(LinalgBackendEigen, SGVector_dot) { const index_t size = 3;