Add Pseudo Inverse to linalg (#4356)

* pinv impl for hermitian and general matrix
* pinv in newton svm
* remove pinv from SGMatrix

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<float64_t> inverse((x_d + 1), (x_d + 1));
-	SGMatrix<float64_t>::pinv(Xsv2sum.data(), x_d + 1, x_d + 1, inverse.data());
 
 	SGVector<float64_t> step(x_d + 1);
 	SGVector<float64_t> 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];

src/shogun/lib/SGMatrix.cpp

@@ -948,54 +948,7 @@ SGMatrix<complex128_t> SGMatrix<complex128_t>::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 <class T>
-float64_t* SGMatrix<T>::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<n; i++)
-	{
-		for (int64_t j=0; j<lsize; j++)
-			vt[i*n+j]=vt[i*n+j]/s[j];
-	}
-
-	cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);
-
-	SG_FREE(u);
-	SG_FREE(vt);
-	SG_FREE(s);
-
-	return target;
-}
-
 /// inverses square matrix in-place
 template <class T>
 void SGMatrix<T>::inverse(SGMatrix<float64_t> matrix)

src/shogun/lib/SGMatrix.h

@@ -393,13 +393,6 @@ template<class T> class SGMatrix : public SGReferencedData
 		/** Inverses square matrix in-place */
 		static void inverse(SGMatrix<float64_t> 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 */

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<Type>& A, SGMatrix<Type>& 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<Type>& A, SGMatrix<Type>& result) const   \
+	{                                                                          \
+		SG_SNOTIMPLEMENTED;                                                    \
+	}
+		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_PINV, SGMatrix)
+#undef BACKEND_GENERIC_PINV
+
 /**
  * Wrapper method of LDLT Cholesky solver
  *

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<Type>& A, SGMatrix<Type>& 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<Type>& A, SGMatrix<Type>& 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<Type> ldlt_solver(                                        \
@@ -465,6 +477,14 @@ namespace shogun
 		    const SGMatrix<T>& A, SGMatrix<T>& L, SGVector<T>& d,
 		    SGVector<index_t>& p, const bool lower) const;
 
+		/** Eigen3 Pseudo inverse of self adjoint matrix */
+		template <typename T>
+		void pinvh_impl(const SGMatrix<T>& A, SGMatrix<T>& result) const;
+
+		/** Eigen3 Pseudo inverse of self adjoint matrix */
+		template <typename T>
+		void pinv_impl(const SGMatrix<T>& A, SGMatrix<T>& result) const;
+
 		/** Eigen3 LDLT Cholesky solver */
 		template <typename T>
 		SGVector<T> ldlt_solver_impl(

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 <typename T>
+		void pinvh(const SGMatrix<T>& A, SGMatrix<T>& 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 <typename T>
+		SGMatrix<T> pinvh(const SGMatrix<T>& A)
+		{
+			SGMatrix<T> 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 <typename T>
+		void pinv(const SGMatrix<T>& A, SGMatrix<T>& 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 <typename T>
+		SGMatrix<T> pinv(const SGMatrix<T>& A)
+		{
+			SGMatrix<T> 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.

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<Type>& A, SGMatrix<Type>& 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<Type>& A, SGMatrix<Type>& 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<Type>& A, const SGVector<Type>& b,                      \
@@ -208,6 +226,46 @@ void LinalgBackendEigen::add_diag_impl(
 	A_eig.diagonal() = alpha * A_eig.diagonal() + beta * b_eig;
 }
 
+template <typename T>
+void LinalgBackendEigen::pinvh_impl(
+    const SGMatrix<T>& A, SGMatrix<T>& result) const
+{
+	auto m = A.num_rows;
+	SGVector<T> s(m);
+	SGMatrix<T> V(m, m);
+	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
+	typename SGMatrix<T>::EigenMatrixXtMap result_eig = result;
+
+	eigen_solver_symmetric_impl<T>(A, s, V, m);
+
+	typename SGVector<T>::EigenVectorXtMap s_eig = s;
+	typename SGMatrix<T>::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 <typename T>
+void LinalgBackendEigen::pinv_impl(
+    const SGMatrix<T>& A, SGMatrix<T>& result) const
+{
+	auto m = std::min(A.num_cols, A.num_rows);
+	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
+	typename SGMatrix<T>::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 <typename T>
 void LinalgBackendEigen::add_col_vec_impl(
     const SGMatrix<T>& A, index_t i, const SGVector<T>& b, SGVector<T>& result,