Simplify LRR via using CDenseFeatures interface for cov, gram, sum

make bias computation optional

doc/cookbook/source/examples/regression/linear_ridge_regression.rst

@@ -10,6 +10,10 @@ We aim to find the linear function that best explains the data, i.e. minimizes t
 
 where :math:`X=\left[{\bf x}_{1},\dots{\bf x}_{N}\right]\in\mathbb{R}^{D\times N}` is the training data matrix, containing :math:`N` training samples of dimension :math:`D`, :math:`y=[y_{1},\dots,y_{N}]^{\top}\in\mathbb{R}^{N}` are the labels, and :math:`\tau>0` scales the regularization term.
 
+Alternatively if :math:`D>N`, the solution can be written as
+.. math::
+    {\bf w}=X\left(\tau I_{N}+X^{\top}X\right)^{-1}y
+
 The bias term is computed as :math:`b=\frac{1}{N}\sum_{i=1}^{N}y_{i}-{\bf w}\cdot\bar{\mathbf{x}}`, where :math:`\bar{\mathbf{x}}=\frac{1}{N}\sum_{i=1}^{N}{\bf x}_{i}`.
 
 For the special case when :math:`\tau = 0`, a wrapper class :sgclass:`CLeastSquaresRegression` is available.
@@ -34,6 +38,10 @@ After training, we can extract :math:`{\bf w}` and the bias.
 
 .. sgexample:: linear_ridge_regression.sg:extract_w
 
+We could also have trained without bias and set it manually.
+
+.. sgexample:: linear_ridge_regression.sg:manual_bias
+
 Finally, we can evaluate the :sgclass:`CMeanSquaredError`.
 
 .. sgexample:: linear_ridge_regression.sg:evaluate_error

examples/meta/src/regression/linear_ridge_regression.sg

@@ -24,10 +24,19 @@ real b = lrr.get_real("bias")
 RealVector w = lrr.get_real_vector("w")
 #[!extract_w]
 
+#[!manual_bias]
+Machine lrr2 = machine("LinearRidgeRegression", tau=0.001, labels=labels_train, use_bias=False)
+lrr2.train(features_train)
+real my_bias = 0.1
+lrr2.put("bias", my_bias)
+Labels labels_predict2 = lrr2.apply(features_test)
+#[!manual_bias]
+
 #![evaluate_error]
 Evaluation eval = evaluation("MeanSquaredError")
 real mse = eval.evaluate(labels_predict, labels_test)
 #![evaluate_error]
 
 # integration testing variables
-RealVector output = labels_test.get_real_vector("labels")
+RealVector output = labels_predict.get_real_vector("labels")
+RealVector output2 = labels_predict2.get_real_vector("labels")

src/interfaces/swig/factory.i

@@ -4,6 +4,8 @@
 %include <shogun/util/factory.h>
 
 %template(features) shogun::features<float64_t>;
+%template(labels) shogun::labels<float64_t>;
+
 
 %newobject shogun::string_features(CFile*, EAlphabet alpha = DNA, EPrimitiveType primitive_type = PT_CHAR);
 %newobject shogun::transformer(const std::string&);

src/shogun/features/DenseFeatures.cpp

@@ -19,6 +19,19 @@
 #include <algorithm>
 #include <string.h>
 
+#define ASSERT_FLOATING_POINT                                                  \
+	switch (get_feature_type())                                                \
+	{                                                                          \
+	case F_SHORTREAL:                                                          \
+	case F_DREAL:                                                              \
+	case F_LONGREAL:                                                           \
+		break;                                                                 \
+	default:                                                                   \
+		REQUIRE(                                                               \
+		    false, "Only defined for %s with real type, not for %s.\n",        \
+		    get_name(), demangled_type<ST>().c_str());                         \
+	}
+
 namespace shogun {
 
 template<class ST> CDenseFeatures<ST>::CDenseFeatures(int32_t size) : CDotFeatures(size)
@@ -1001,6 +1014,53 @@ template< class ST > CDenseFeatures< ST >* CDenseFeatures< ST >::obtain_from_gen
 	return (CDenseFeatures< ST >*) base_features;
 }
 
+template <typename ST>
+SGVector<ST> CDenseFeatures<ST>::sum() const
+{
+	// TODO optimize non batch mode, but get_feature_vector is non const :(
+	SGVector<ST> result = linalg::rowwise_sum(get_feature_matrix());
+	return result;
+}
+
+template <typename ST>
+SGVector<ST> CDenseFeatures<ST>::mean() const
+{
+	ASSERT_FLOATING_POINT
+
+	auto result = sum();
+	ST scale = ((ST)1.0) / get_num_vectors();
+	linalg::scale(result, result, scale);
+	return result;
+}
+
+template <typename ST>
+SGMatrix<ST> CDenseFeatures<ST>::cov() const
+{
+	// TODO optimize non batch mode, but get_feature_vector is non const :(
+	auto mat = get_feature_matrix();
+	return linalg::matrix_prod(mat, mat, false, true);
+}
+
+template <typename ST>
+SGMatrix<ST> CDenseFeatures<ST>::gram() const
+{
+	// TODO optimize non batch mode, but get_feature_vector is non const :(
+	auto mat = get_feature_matrix();
+	return linalg::matrix_prod(mat, mat, true, false);
+}
+
+template <typename ST>
+SGVector<ST> CDenseFeatures<ST>::dot(const SGVector<ST>& other) const
+{
+	REQUIRE(
+		get_num_vectors() == other.size(), "Number of feature vectors (%d) "
+		                                   "must match provided vector's size "
+		                                   "(%d).\n",
+		get_num_features(), other.size());
+	// TODO optimize non batch mode, but get_feature_vector is non const :(
+	return linalg::matrix_prod(get_feature_matrix(), other, false);
+}
+
 template class CDenseFeatures<bool>;
 template class CDenseFeatures<char>;
 template class CDenseFeatures<int8_t>;

src/shogun/features/DenseFeatures.h

@@ -11,14 +11,14 @@
 
 #include <shogun/lib/config.h>
 
-#include <shogun/lib/common.h>
-#include <shogun/lib/Cache.h>
-#include <shogun/io/File.h>
 #include <shogun/features/DotFeatures.h>
 #include <shogun/features/StringFeatures.h>
+#include <shogun/io/File.h>
+#include <shogun/lib/Cache.h>
 #include <shogun/lib/DataType.h>
-
 #include <shogun/lib/SGMatrix.h>
+#include <shogun/lib/common.h>
+#include <shogun/mathematics/linalg/LinalgNamespace.h>
 
 namespace shogun {
 template<class ST> class CStringFeatures;
@@ -303,6 +303,46 @@ template<class ST> class CDenseFeatures: public CDotFeatures
 	virtual float64_t dot(int32_t vec_idx1, CDotFeatures* df,
 			int32_t vec_idx2);
 
+	/** Computes the sum of all feature vectors
+	 * @return Sum of all feature vectors
+	 */
+	SGVector<ST> sum() const;
+
+	/** Computes the empirical mean of all feature vectors
+	 * @return Mean of all feature vectors
+	 */
+	SGVector<ST> mean() const;
+
+	/** Computes the \f$DxD\f$ (uncentered, un-normalized) covariance matrix
+	 *
+	 *\f[
+	 * X X^\top
+	 * \f]
+	 *
+	 * where \f$X\f$ is the \f$DxN\f$ dimensional feature matrix with \f$N\f$
+	 * feature vectors of dimension \f$D\f$.
+	 */
+	SGMatrix<ST> cov() const;
+	/** Computes the \f$fNxN\f$ (uncentered, un-normalized) gram matrix of
+	 * pairwise dot products, that is
+	 *
+	 *\f[
+	 * X^\top X
+	 * \f]
+	 *
+	 * where \f$X\f$ is the \f$DxN\f$ dimensional feature matrix with \f$N\f$
+	 * feature vectors of dimension \f$D\f$.
+	 */
+	SGMatrix<ST> gram() const;
+
+	/** Computes the dot product of the feature matrix with a given vector.
+	 *
+	 * @param other Vector to compute dot products with, size must match number
+	 * of feature vectors
+	 * @return Vector as many entries as feature dimensions
+	 */
+	SGVector<ST> dot(const SGVector<ST>& other) const;
+
 	/** compute dot product between vector1 and a dense vector
 	 *
 	 * possible with subset

src/shogun/machine/LinearMachine.cpp

@@ -37,7 +37,9 @@ void CLinearMachine::init()
 
 	SG_ADD(&m_w, "w", "Parameter vector w.", MS_NOT_AVAILABLE);
 	SG_ADD(&bias, "bias", "Bias b.", MS_NOT_AVAILABLE);
-	SG_ADD(&features, "features", "Feature object.", MS_NOT_AVAILABLE);
+	SG_ADD(
+	    (CFeatures**)&features, "features", "Feature object.",
+	    MS_NOT_AVAILABLE);
 }
 
 

src/shogun/mathematics/linalg/LinalgBackendBase.h

@@ -140,6 +140,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 +228,22 @@ namespace shogun
 		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_CHOLESKY_SOLVER, SGMatrix)
 #undef BACKEND_GENERIC_CHOLESKY_SOLVER
 
+/**
+ * Wrapper 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_NON_INTEGER_PTYPE(
+		    BACKEND_GENERIC_CHOLESKY_RANK_UPDATE, SGMatrix)
+#undef BACKEND_GENERIC_CHOLESKY_RANK_UPDATE
+
 /**
  * Wrapper method of LDLT Cholesky decomposition
  *

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,15 @@ 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_NON_INTEGER_PTYPE(
+		    BACKEND_GENERIC_CHOLESKY_RANK_UPDATE, SGMatrix)
+#undef BACKEND_GENERIC_CHOLESKY_RANK_UPDATE
+
 /** Implementation of @see LinalgBackendBase::cholesky_solver */
 #define BACKEND_GENERIC_CHOLESKY_SOLVER(Type, Container)                       \
 	virtual SGVector<Type> cholesky_solver(                                    \
@@ -443,6 +458,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 +485,12 @@ namespace shogun
 		SGMatrix<T>
 		cholesky_factor_impl(const SGMatrix<T>& A, const bool lower) const;
 
+		/** Eigen3 Cholesky rank one update */
+		template <typename T>
+		SGMatrix<T> cholesky_rank_update_impl(
+		    SGMatrix<T>& L, const SGVector<T>& b, T alpha,
+		    const bool lower) const;
+
 		/** Eigen3 Cholesky solver */
 		template <typename T>
 		SGVector<T> cholesky_solver_impl(