Merge pull request #4387 from karlnapf/feature/rank_updates

add plain and cholesky rank updates

src/shogun/clustering/GMM.cpp

@@ -591,8 +591,8 @@ void CGMM::max_likelihood(SGMatrix<float64_t> 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:
 			    {

src/shogun/distributions/Gaussian.cpp

@@ -167,7 +167,7 @@ float64_t CGaussian::update_params_em(const SGVector<float64_t> 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<float64_t>(alpha_k[j], v, v, cov_sum);
+			linalg::rank_update(cov_sum, v, alpha_k[j]);
 #endif
 			break;
 		case DIAG:

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<Type>& 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<Type>& L, const SGVector<Type>& 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<Type>& A, const SGVector<Type>& 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
 	};
 }

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<Type>& 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<Type>& L, const SGVector<Type>& 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<Type>& A, const SGVector<Type>& 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<Type> cholesky_solver(                                    \
@@ -230,15 +251,15 @@ namespace shogun
 /** Implementation of @see LinalgBackendBase::max */
 #define BACKEND_GENERIC_MAX(Type, Container)                                   \
 	virtual Type max(const Container<Type>& 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<Type>& 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<Type>& a, Container<Type>& 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<T>& A, const SGVector<T>& b, T alpha, T beta) const;
 
+		/** Eigen3 add diagonal scalar method */
+		template <typename T>
+		void add_ridge_impl(SGMatrix<T>& A, T beta) const;
+
 		/** Eigen3 add vector to each column of matrix method */
 		template <typename T>
 		void add_vector_impl(
@@ -466,6 +491,17 @@ namespace shogun
 		SGMatrix<T>
 		cholesky_factor_impl(const SGMatrix<T>& A, const bool lower) const;
 
+		/** Eigen3 Cholesky rank one update */
+		template <typename T>
+		void cholesky_rank_update_impl(
+		    SGMatrix<T>& L, const SGVector<T>& b, T alpha,
+		    const bool lower) const;
+
+		/** Eigen3 rank one update */
+		template <typename T>
+		void
+		rank_update_impl(SGMatrix<T>& A, const SGVector<T>& b, T alpha) const;
+
 		/** Eigen3 Cholesky solver */
 		template <typename T>
 		SGVector<T> cholesky_solver_impl(

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);                                             \

src/shogun/mathematics/linalg/LinalgNamespace.h

@@ -583,6 +583,21 @@ namespace shogun
 			infer_backend(A, SGMatrix<T>(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 <typename T>
+		void add_ridge(SGMatrix<T>& 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 <typename T>
-		void
-		dger(T alpha, const SGVector<T> x, const SGVector<T> y, SGMatrix<T>& A)
-		{
-			auto x_martix =
-			    SGVector<T>::convert_to_matrix(x, A.num_rows, 1, false);
-			auto y_martix =
-			    SGVector<T>::convert_to_matrix(y, A.num_cols, 1, false);
-
-			auto temp_martix = SGMatrix<T>::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<T>
 		cholesky_factor(const SGMatrix<T>& 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 <typename T>
+		void cholesky_rank_update(
+		    SGMatrix<T>& L, const SGVector<T>& 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<T>(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 <typename T>
+		void rank_update(SGMatrix<T>& A, const SGVector<T>& 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<T>(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<T> cholesky_solver(
 		    const SGMatrix<T>& L, const SGVector<T>& 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<T>(b))
 			    ->cholesky_solver(L, b, lower);
 		}
@@ -763,6 +818,20 @@ namespace shogun
 		    const SGMatrix<T>& L, const SGVector<T>& d, SGVector<index_t>& p,
 		    const SGVector<T>& 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<T>(d), SGMatrix<T>(b))
 			    ->ldlt_solver(L, d, p, b, lower);
 		}